/[escript]/trunk/downunder/py_src/regularizations.py
ViewVC logotype

Contents of /trunk/downunder/py_src/regularizations.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4097 - (show annotations)
Fri Dec 7 01:18:35 2012 UTC (6 years, 8 months ago) by caltinay
File MIME type: text/x-python
File size: 16353 byte(s)
Minor edits.

1
2 ##############################################################################
3 #
4 # Copyright (c) 2003-2012 by University of Queensland
5 # http://www.uq.edu.au
6 #
7 # Primary Business: Queensland, Australia
8 # Licensed under the Open Software License version 3.0
9 # http://www.opensource.org/licenses/osl-3.0.php
10 #
11 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 # Development since 2012 by School of Earth Sciences
13 #
14 ##############################################################################
15
16 __copyright__="""Copyright (c) 2003-2012 by University of Queensland
17 http://www.uq.edu.au
18 Primary Business: Queensland, Australia"""
19 __license__="""Licensed under the Open Software License version 3.0
20 http://www.opensource.org/licenses/osl-3.0.php"""
21 __url__="https://launchpad.net/escript-finley"
22
23 __all__ = ['Regularization']
24
25 from costfunctions import CostFunction
26
27 import numpy as np
28 from esys.escript import Data, grad, inner, integrate, interpolate, kronecker, boundingBoxEdgeLengths, vol
29 from esys.escript.linearPDEs import LinearPDE, IllegalCoefficientValue
30 from esys.escript.pdetools import ArithmeticTuple
31
32 class Regularization(CostFunction):
33 """
34 The regularization term for the level set function ``m`` within the cost
35 function J for an inversion:
36
37 *J(m)=1/2 * sum_k imtegrate( mu_0[k] * s0[k] * m_k**2 + mu_1[k] * s1[k,i] * m_{k,i}**2) + sum_l<k mu_c[l,k] sc[l,l] * | curl(m_k) x curl(m_l) |^2*
38
39 where s0[k], s1[k,i] and sc[k,l] are non-negative scaling factors and
40 mu_0[k], mu_1[k], mu_c[l,k] are weighting factors which may be altered
41 during the inversion. The scaling factors are normalized such that their
42 integrals over the domain are constant:
43
44 *integrate(s0[k])=1*
45 *integrate(inner(s1[k,:],L[:])=1*
46 *integrate(inner(sc[l,k]*L**4)=1*
47
48 """
49 def __init__(self, domain, numLevelSets=1,
50 s0=None, s1=None, sc=None,
51 location_of_set_m=Data(),
52 useDiagonalHessianApproximation=True, tol=1e-8):
53 """
54 initialization.
55
56 :param domain: domain
57 :type domain: `Domain`
58 :param numLevelSets: number of level sets
59 :type numLevelSets: ``int``
60 :param s0: scaling factor for the m**2 term. If not set zero is assumed.
61 :type s0: ``Scalar`` if ``numLevelSets`` == 1 or `Data` object of shape
62 (``numLevelSets`` ,) if ``numLevelSets`` > 1
63 :param s1: scaling factor for the grad(m_i) terms. If not set zero is assumed
64 :type s1: ``Vector`` if ``numLevelSets`` == 1 or `Data` object of shape
65 (``numLevelSets`` , DIM) if ``numLevelSets`` > 1
66 :param sc: scaling factor for the cross correlation terms. If not set
67 zero is assumed. Used for the case if ``numLevelSets`` > 1
68 only. Only values ``sc[l,k]`` in the lower triangle (l<k)
69 are used.
70 :type sc: `Data` object of shape (``numLevelSets`` , ``numLevelSets``)
71 :param location_of_set_m: marks location of zero values of the level
72 set function ``m`` by a positive entry.
73 :type location_of_set_m: ``Scalar`` if ``numLevelSets`` == 1 or `Data`
74 object of shape (``numLevelSets`` ,) if ``numLevelSets`` > 1
75 :param useDiagonalHessianApproximation: if True cross correlation terms
76 between level set components are ignored when calculating
77 approximations of the inverse of the Hessian Operator.
78 This can speed-up the calculation of the inverse but may
79 lead to an increase of the number of iteration steps in the
80 inversion.
81 :type useDiagonalHessianApproximation: ``bool``
82 :param tol: tolerance when solving the PDE for the inverse of the
83 Hessian Operator
84 :type tol: positive ``float``
85
86 """
87 if numLevelSets>1:
88 raise ValueError("Currently only numLevelSets<=1 is supported.")
89 if s0 == None and s1==None:
90 raise ValueError("Values for s0 or s1 must be given.")
91
92 self.__domain=domain
93 DIM=self.__domain.getDim()
94 self.__L2=np.asarray(boundingBoxEdgeLengths(domain))**2
95 self.__L4=np.sum(self.__L2)**2
96 self.__numLevelSets=numLevelSets
97
98 if self.__numLevelSets > 1:
99 self.__useDiagonalHessianApproximation=useDiagonalHessianApproximation
100 else:
101 self.__useDiagonalHessianApproximation=True
102 self._update_Hessian=True
103
104 self.__pde=LinearPDE(self.__domain, numEquations=self.__numLevelSets)
105 self.__pde.getSolverOptions().setTolerance(tol)
106 self.__pde.setSymmetryOn()
107 self.__pde_is_set=False
108 try:
109 self.__pde.setValue(q=location_of_set_m)
110 except IllegalCoefficientValue:
111 raise ValueError("Unable to set location of fixed level set function.")
112
113 self.__total_num_weights=2*numLevelSets+((numLevelSets-1)*numLevelSets)/2
114 self.__weight_index=[] # this is a mapping from the relevant mu-coefficients to the set of all mu-coefficients
115 # we count s0, then s1, then sc (k<l).
116 # THE S0 weighting factor
117 n=0
118 VV=vol(domain)
119 if s0 is not None:
120 s0 = interpolate(s0,self.__pde.getFunctionSpaceForCoefficient('D'))
121 s=s0.getShape()
122 if numLevelSets == 1 :
123 if s == () :
124 V=integrate(s0)
125 if V > 0:
126 self.__weight_index.append(n)
127 s0*=VV/V
128 else:
129 s0=None
130 else:
131 raise ValueError("Unexpected shape %s for weight s0."%s)
132 else:
133 if s == (numLevelSets,):
134 for k in xrange(numLevelSets):
135 V=integrate(s0[k])
136 if V > 0:
137 self.__weight_index.append(n+k)
138 s0[k]*=VV/V
139 else:
140 raise ValueError("Unexpected shape %s for weight s0."%s)
141 self.__s0=s0
142 n+=numLevelSets
143
144 # The S1 weighting factor
145 if not s1 is None:
146 s1 = interpolate(s1,self.__pde.getFunctionSpaceForCoefficient('A'))
147 s=s1.getShape()
148
149 if numLevelSets == 1 :
150 if s == (DIM,) :
151 V=integrate(inner(s1, 1/self.__L2))
152 if V > 0:
153 self.__weight_index.append(n)
154 s1*=VV/V
155 print "REG SCALE = ",s1
156 else:
157 raise ValueError("Unexpected shape %s for weight s1."%s)
158 else:
159 if s == (numLevelSets,DIM):
160 for k in xrange(numLevelSets):
161 for i in xrange(DIM):
162 ww=s1[k,:]
163 V=integrate(inner(ww,1/self.__L2))
164 if V > 0:
165 self.__weight_index.append(n+i)
166 s1[k,:]=ww*(VV/V)
167 else:
168 raise ValueError("Unexpected shape %s for weight s1."%s)
169
170 self.__s1=s1
171 n+=numLevelSets
172
173 # The SC weighting factor
174 if not sc is None:
175 if numLevelSets == 1 :
176 sc=None
177 else:
178 sc = interpolate(sc,self.__pde.getFunctionSpaceForCoefficient('A'))
179 s=sc.getShape()
180 if s == (numLevelSets,numLevelSets):
181 for k in xrange(numLevelSets):
182 sc[k,k]=0.
183 for l in xrange(k):
184 ww=sc[l,k]
185 V=integrate(ww)
186 if V > 0:
187 self.__weight_index.append(n+k*numLevelSets+l)
188 sc[l,k]=ww*VV/V*self.__L4
189 sc[k,l]=0
190 else:
191 raise ValueError("Unexpected shape %s for weight s0."%s)
192
193 self.__sc=sc
194 self.setWeights()
195
196 def getDomain(self):
197 """
198 returns the domain of the regularization term
199
200 :rtype: ``Domain``
201 """
202 return self.__domain
203
204 def getNumLevelSets(self):
205 """
206 returns the number of level set functions
207
208 :rtype: ``int``
209 """
210 return self.__numLevelSets
211
212 def getPDE(self):
213 """
214 returns the linear PDE to be solved for the Hessian Operator inverse
215
216 :rtype: `LinearPDE`
217 """
218 return self.__pde
219
220 def getDualProduct(self, m, r):
221 """
222 returns the dual product of a gradient represented by X=r[1] and Y=r[0]
223 with a level set function m:
224
225 *Y_i*m_i + X_ij*m_{i,j}*
226
227 :type m: `Data`
228 :type r: `ArithmeticTuple`
229 :rtype: ``float``
230 """
231 A=0
232 if not r[0].isEmpty(): A+=integrate(inner(r[0], m))
233 if not r[1].isEmpty(): A+=integrate(inner(r[1], grad(m)))
234 return A
235
236 def getValue(self, m, grad_m):
237 """
238 returns the value of the cost function J with respect to m.
239
240 :rtype: ``float``
241 """
242 mu=self.getWeights( uncompress=True)
243 DIM=self.getDomain().getDim()
244 numLS=self.getNumLevelSets()
245
246 A=0
247 n=0
248
249 if self.__s0 is not None:
250 A+=inner(integrate(m**2*self.__s0), mu[:numLS])
251 n+=numLS
252
253 if self.__s1 is not None:
254 if numLS == 1:
255 A+=integrate(inner(grad_m**2, self.__s1))*mu[n]
256 else:
257 for k in xrange(numLS):
258 A+=mu[n+k]*integrate(inner(grad_m[k,:]**2,self.__s1[k,:]))
259 n+=numLS
260
261 if self.__sc is not None:
262 for k in xrange(numLS):
263 gk=grad_m[k,:]
264 len_gk=length(gk)
265 for l in xrange(k):
266 gl=grad_m[l,:]
267 A+= mu[n+k*numLS+l] * integrate( self.__sc[l,k] * ( len_gk * length(gl) )**2 - inner(gk, gl)**2 )
268 return A/2
269
270 def getGradient(self, m, grad_m):
271 """
272 returns the gradient of the cost function J with respect to m.
273 The function returns Y_k=dPsi/dm_k and X_kj=dPsi/dm_kj
274 """
275
276 mu=self.getWeights( uncompress=True)
277 DIM=self.getDomain().getDim()
278 numLS=self.getNumLevelSets()
279
280 n=0
281
282 if self.__s0 is not None:
283 Y = m * self.__s0 * mu[:numLS]
284 else:
285 Y = Data()
286 n+=numLS
287
288 if self.__s1 is not None:
289 if numLS == 1:
290 X=grad_m* self.__s1*mu[n]
291 else:
292 X=self.getPDE().createCoefficient("X")
293 for k in xrange(numLS):
294 X[k,:]=mu[n+k]*grad_m[k,:]*self.__s1[k,:]
295 else:
296 X = Data()
297 n+=numLS
298 if self.__sc is not None:
299 raise NotImplementedError
300 return ArithmeticTuple(Y, X)
301
302 def getArguments(self, m):
303 """
304 """
305 return ( grad(m),)
306
307 def getInverseHessianApproximation(self, m, r, grad_m):
308 """
309 """
310 if self._new_mu or self._update_Hessian:
311
312 self._new_mu=False
313 self._update_Hessian=False
314
315 mu=self.getWeights( uncompress=True)
316 DIM=self.getDomain().getDim()
317 numLS=self.getNumLevelSets()
318 n=0
319 if self.__s0 is not None:
320 if numLS == 1:
321 D=self.__s0 * mu[n]
322 else:
323 D=self.getPDE().createCoefficient("D")
324 for k in xrange(numLS): D[k,k]=self.__s0[k] * mu[n+k]
325 self.getPDE().setValue(D=D)
326 n+=numLS
327
328 A=self.getPDE().createCoefficient("A")
329 if self.__s1 is not None:
330 if numLS == 1:
331 for i in xrange(DIM): A[i,i]=self.__s1[i] * mu[n]
332 else:
333 for k in xrange(numLS):
334 for i in xrange(DIM): A[k,i,k,i]=self.__s1[k,i] * mu[n+k]
335 n+=numLS
336 if self.__sc is not None:
337 raise NotImplementedError
338 print A
339 self.getPDE().setValue(A=A)
340 self.getPDE().resetRightHandSideCoefficients()
341 self.getPDE().setValue(X=r[1])
342 print "X only: ",self.getPDE().getSolution()
343 self.getPDE().resetRightHandSideCoefficients()
344 self.getPDE().setValue(Y=r[0])
345 print "Y only: ",self.getPDE().getSolution()
346
347 self.getPDE().resetRightHandSideCoefficients()
348 self.getPDE().setValue(X=r[1], Y=r[0])
349 return self.getPDE().getSolution()
350
351 def updateHessian(self):
352 """
353 notify the class to recalculate the Hessian operator
354 """
355 if not self.__useDiagonalHessianApproximation:
356 self._update_Hessian=True
357
358 # ================ we should factor these out ==============================
359
360 def getNumRelevantTerms(self):
361 """
362 returns the number of terms in the cost function of regularization
363 with non-zero scaling factors.
364
365 :rtype: ``int``
366 """
367 return len(self.__weight_index)
368
369 def getNumTerms(self):
370 """
371 returns the number of terms in the cost function of regularization
372 with non-zero scaling factors.
373
374 :rtype: ``int``
375 """
376 return len(self.__weight_index)
377
378 def setWeights(self, mu=None):
379 """
380 sets the values for the weighting terms in the cost function.
381 Note that only values corresponding to entries with non-negative
382 scaling factors are used.
383
384 :param mu: weights
385 :type mu: ``list`` of ``float`` or ``np,array``
386 """
387 if mu == None:
388 mu = np.ones(self.getNumRelevantTerms())
389 mu=np.asarray(mu)
390
391 if len(mu) == self.getNumRelevantTerms():
392 if not mu.shape == (self.getNumRelevantTerms(),):
393 raise ValueError("%s values are expected."%self.getNumRelevantTerms())
394 self.__mu=mu
395 else:
396 if not mu.shape == (self.__total_num_weights,):
397 raise ValueError("%s values are expected."%self.__total_num_weights)
398 self.__mu = np.zeros(self.getNumRelevantTerms())
399 for i in xrange(len(self.__weight_index)): self.__mu[i]=mu[self.__weight_index[i]]
400 self._new_mu=True
401
402 def setWeightsForS0(self, mu=None):
403 """
404 sets the weights for s0-terms
405 """
406 numLS=self.getNumLevelSets()
407 my_mu=self.getWeights(uncompress=True)
408 if mu is None:
409 my_mu[:numLS]=1
410 else:
411 my_mu[:numLS]=mu
412 self.setWeights(my_mu)
413
414 def setWeightsForS1(self, mu=None):
415 """
416 sets the weights for s1-terms
417 """
418 numLS=self.getNumLevelSets()
419 my_mu=self.getWeights(uncompress=True)
420 if mu is None:
421 my_mu[numLS:2*numLS]=1
422 else:
423 my_mu[numLS:2*numLS]=mu
424 self.setWeights(my_mu)
425
426 def setWeightsForSc(self, mu):
427 """
428 sets the weights for sc-terms
429 """
430 numLS=self.getNumLevelSets()
431 my_mu=self.getWeights(uncompress=True)
432 if mu is None:
433 my_mu[2*numLS:]=1
434 else:
435 my_mu[2*numLS:]=mu
436 self.setWeights(my_mu)
437
438 def getWeights(self, uncompress=False):
439 """
440 Returns the weights for the terms in the cost function.
441 The first ``numLevelSets`` values are used for the regularization terms
442 and the remaining values for the cross correlation terms.
443 """
444 if uncompress:
445 mu = np.zeros(self.__total_num_weights)
446 for i in xrange(len(self.__weight_index)):
447 mu[self.__weight_index[i]] = self.__mu[i]
448 return mu
449 else:
450 return self.__mu
451
452 def getWeightIndex(self):
453 """
454 returns an index to the contributions of terms with non-zero scaling
455 factor.
456 """
457 return self.__weight_index
458

  ViewVC Help
Powered by ViewVC 1.1.26