/[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 4122 - (show annotations)
Thu Dec 20 05:42:35 2012 UTC (7 years, 11 months ago) by gross
File MIME type: text/x-python
File size: 19148 byte(s)
More work towards joint inversion. There is now one class for inversion cost function which can handle all relevant cases 
for a single model inversion, cross gradient case and functional dependence of physical parameters.


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, Scalar, grad, inner, integrate, interpolate, kronecker, boundingBoxEdgeLengths, vol, sqrt, length
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[k] * ( w0[k] * m_k**2 * w1[k,i] * m_{k,i}**2) + sum_l<k mu_c[l,k] wc[l,k] * | curl(m_k) x curl(m_l) |^2*
38
39 where w0[k], w1[k,i] and wc[k,l] are non-negative weighting factors and
40 mu[k] and mu_c[l,k] are trade-off factors which may be altered
41 during the inversion. The weighting factors are normalized such that their
42 integrals over the domain are constant:
43
44 *integrate(w0[k] + inner(w1[k,:],1/L[:]**2))=scale[k]* volume(domain)*
45 *integrate(wc[l,k]*1/L**4)=scale_c[k]* volume(domain) *
46
47 """
48 def __init__(self, domain, numLevelSets=1,
49 w0=None, w1=None, wc=None,
50 location_of_set_m=Data(),
51 useDiagonalHessianApproximation=False, tol=1e-8,
52 scale=None, scale_c=None):
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 w0: weighting factor for the m**2 term. If not set zero is assumed.
61 :type w0: ``Scalar`` if ``numLevelSets`` == 1 or `Data` object of shape
62 (``numLevelSets`` ,) if ``numLevelSets`` > 1
63 :param w1: weighting factor for the grad(m_i) terms. If not set zero is assumed
64 :type w1: ``Vector`` if ``numLevelSets`` == 1 or `Data` object of shape
65 (``numLevelSets`` , DIM) if ``numLevelSets`` > 1
66 :param wc: weighting factor for the cross gradient terms. If not set
67 zero is assumed. Used for the case if ``numLevelSets`` > 1
68 only. Only values ``wc[l,k]`` in the lower triangle (l<k)
69 are used.
70 :type wc: `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 gradient 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 :param scale: weighting factor for level set function variation terms. If not set one is used.
87 :type scale: ``Scalar`` if ``numLevelSets`` == 1 or `Data` object of shape
88 (``numLevelSets`` ,) if ``numLevelSets`` > 1
89 :param scale_c: scale for the cross gradient terms. If not set
90 one is assumed. Used for the case if ``numLevelSets`` > 1
91 only. Only values ``scale_c[l,k]`` in the lower triangle (l<k)
92 are used.
93 :type scale_c: `Data` object of shape (``numLevelSets`` , ``numLevelSets``)
94
95
96 """
97 if w0 == None and w1==None:
98 raise ValueError("Values for w0 or for w1 must be given.")
99 if wc == None and numLevelSets>1:
100 raise ValueError("Values for wc must be given.")
101
102 self.__domain=domain
103 DIM=self.__domain.getDim()
104 self.__numLevelSets=numLevelSets
105
106 self.__pde=LinearPDE(self.__domain, numEquations=self.__numLevelSets)
107 self.__pde.getSolverOptions().setTolerance(tol)
108 self.__pde.setSymmetryOn()
109 try:
110 self.__pde.setValue(q=location_of_set_m)
111 except IllegalCoefficientValue:
112 raise ValueError("Unable to set location of fixed level set function.")
113
114 # =========== check the shape of the scales: =================================
115 if scale is None:
116 if numLevelSets == 1 :
117 scale = 1.
118 else:
119 scale = np.ones((numLevelSets,))
120 else:
121 scale=np.asarray(scale)
122 if numLevelSets == 1 :
123 if scale.shape == ():
124 if not scale > 0 :
125 raise ValueError("Value for scale must be positive.")
126 else:
127 raise ValueError("Unexpected shape %s for scale."%scale.shape)
128 else:
129 if scale.shape is (numLevelSets,):
130 if not min(scale) > 0:
131 raise ValueError("All value for scale must be positive.")
132 else:
133 raise ValueError("Unexpected shape %s for scale."%scale.shape)
134
135 if scale_c is None or numLevelSets < 2:
136 scale_c = np.ones((numLevelSets,numLevelSets))
137 else:
138 scale_c=np.asarray(scale_c)
139 if scale_c.shape == (numLevelSets,numLevelSets):
140 if not all( [ [ scale_c[l,k] > 0. for l in xrange(k) ] for k in xrange(1,numLevelSets) ]):
141 raise ValueError("All values in the lower triangle of scale_c must be positive.")
142 else:
143 raise ValueError("Unexpected shape %s for scale."%scale_c.shape)
144 # ===== check the shape of the weights: ============================================
145 if w0 is not None:
146 w0 = interpolate(w0,self.__pde.getFunctionSpaceForCoefficient('D'))
147 s0=w0.getShape()
148 if numLevelSets == 1 :
149 if not s0 == () :
150 raise ValueError("Unexpected shape %s for weight w0."%s0)
151 else:
152 if not s0 == (numLevelSets,):
153 raise ValueError("Unexpected shape %s for weight w0."%s0)
154 if not w1 is None:
155 w1 = interpolate(w1,self.__pde.getFunctionSpaceForCoefficient('A'))
156 s1=w1.getShape()
157 if numLevelSets is 1 :
158 if not s1 == (DIM,) :
159 raise ValueError("Unexpected shape %s for weight w1."%s1)
160 else:
161 if not s1 == (numLevelSets,DIM):
162 raise ValueError("Unexpected shape %s for weight w1."%s1)
163 if numLevelSets == 1 :
164 wc=None
165 else:
166 wc = interpolate(wc,self.__pde.getFunctionSpaceForCoefficient('A'))
167 sc=wc.getShape()
168 raise ValueError("Unexpected shape %s for weight wc."%sc)
169 # ============= now we rescale weights: ======================================
170 L2s=np.asarray(boundingBoxEdgeLengths(domain))**2
171 L4=1/np.sum(1/L2s)**2
172 if numLevelSets is 1 :
173 A=0
174 if w0 is not None:
175 A = integrate(w0)
176 if w1 is not None:
177 A += integrate(inner(w1, 1/L2s))
178 if A > 0:
179 f = scale/A
180 if w0 is not None:
181 w0*=f
182 if w1 is not None:
183 w1*=f
184 else:
185 raise ValueError("Non-positive weighting factor detected.")
186 else:
187 for k in xrange(numLevelSets):
188 A=0
189 if w0 is not None:
190 A = integrate(w0[k])
191 if w1 is not None:
192 A += integrate(inner(w1[k,:], 1/L2s))
193 if A > 0:
194 f = scale[k]/A
195 if w0 is not None:
196 w0[k]*=f
197 if w1 is not None:
198 w1[k,:]*=f
199 else:
200 raise ValueError("Non-positive weighting factor for level set component %d detected."%k)
201
202 # and now the cross-gradient:
203 for l in xrange(k):
204 A = integrate(wc[l,k])/L4
205 if A > 0:
206 f = scale_c[l,k]/A
207 wc[l,k]*=f
208 else:
209 raise ValueError("Non-positive weighting factor for cross-gradient level set components %d and %d detected."%(l,k))
210
211 self.__w0=w0
212 self.__w1=w1
213 self.__wc=wc
214
215
216
217 self.__pde_is_set=False
218 if self.__numLevelSets > 1:
219 self.__useDiagonalHessianApproximation=useDiagonalHessianApproximation
220 else:
221 self.__useDiagonalHessianApproximation=True
222 self._update_Hessian=True
223
224 self.__num_tradeoff_factors=numLevelSets+((numLevelSets-1)*numLevelSets)/2
225 self.setTradeOffFactors()
226 self.__vol_d=vol(self.__domain)
227
228 def getDomain(self):
229 """
230 returns the domain of the regularization term
231
232 :rtype: ``Domain``
233 """
234 return self.__domain
235
236 def getNumLevelSets(self):
237 """
238 returns the number of level set functions
239
240 :rtype: ``int``
241 """
242 return self.__numLevelSets
243
244 def getPDE(self):
245 """
246 returns the linear PDE to be solved for the Hessian Operator inverse
247
248 :rtype: `LinearPDE`
249 """
250 return self.__pde
251
252 def getDualProduct(self, m, r):
253 """
254 returns the dual product of a gradient represented by X=r[1] and Y=r[0]
255 with a level set function m:
256
257 *Y_i*m_i + X_ij*m_{i,j}*
258
259 :type m: `Data`
260 :type r: `ArithmeticTuple`
261 :rtype: ``float``
262 """
263 A=0
264 if not r[0].isEmpty(): A+=integrate(inner(r[0], m))
265 if not r[1].isEmpty(): A+=integrate(inner(r[1], grad(m)))
266 return A
267 def getNumTradeOffFactors(self):
268 """
269 returns the number of trade-off factors being used.
270
271 :rtype: ``int``
272 """
273 return self.__num_tradeoff_factors
274
275 def setTradeOffFactors(self, mu=None):
276 """
277 sets the trade-off factors for the level-set variation and the cross-gradient
278
279 :param mu: new values for the trade-off factors where values mu[:numLevelSets] are the
280 trade-off factors for the level-set variation and the remaining values for
281 the cross-gradient part with mu_c[l,k]=mu[numLevelSets+l+((k-1)*k)/2] (l<k).
282 If no values for mu is given ones are used. Values must be positive.
283 :type mu: ``list`` of ``float`` or ```numpy.array```
284 """
285 numLS=self.getNumLevelSets()
286 numTF=self.getNumTradeOffFactors()
287 if mu is None:
288 mu = np.ones((numTF,))
289 else:
290 mu = np.asarray(mu)
291
292 if mu.shape == (numTF,):
293 self.setTradeOffFactorsForVariation(mu[:numLS])
294 mu_c2=np.zeros((numLS,numLS))
295 for k in xrange(numLS):
296 for l in xrange(k):
297 mu_c2[l,k] = mu[numLS+l+((k-1)*k)/2]
298 self.setTradeOffFactorsForCrossGradient(mu_c2)
299 elif mu.shape == () and numLS ==1:
300 self.setTradeOffFactorsForVariation(mu)
301 else:
302 raise ValueError("Unexpected shape %s for mu."%(mu.shape,))
303
304 def setTradeOffFactorsForVariation(self, mu=None):
305 """
306 sets the trade-off factors for the level-set variation part
307
308 :param mu: new values for the trade-off factors. Values must be positive.
309 :type mu: `float``, ``list`` of ``float`` or ```numpy.array```
310 """
311 numLS=self.getNumLevelSets()
312 if mu is None:
313 if numLS == 1:
314 mu = 1.
315 else:
316 mu = np.ones((numLS,))
317
318 mu=np.asarray(mu)
319 if numLS == 1:
320 if mu.shape == (1,): mu=mu[0]
321 if mu.shape == ():
322 if mu > 0:
323 self.__mu= mu
324 self._new_mu=True
325 else:
326 raise ValueError("Value for trade-off factor must be positive.")
327 else:
328 raise ValueError("Unexpected shape %s for mu."%mu.shape)
329 else:
330 if mu.shape == (numLS,):
331 if min(mu) > 0:
332 self.__mu= mu
333 self._new_mu=True
334 else:
335 raise ValueError("All value for mu must be positive.")
336 else:
337 raise ValueError("Unexpected shape %s for trade-off factor."%mu.shape)
338
339 def setTradeOffFactorsForCrossGradient(self, mu_c=None):
340 """
341 sets the trade-off factors for the cross--gradient terms
342
343 :param mu_c: new values for the trade-off factors for the cross--gradient terms. Values must be positive.
344 if now value is given ones are used. Onky value mu_c[l,k] for l<k are used.
345 :type mu: `float``, ``list`` of ``float`` or ```numpy.array```
346
347 """
348 numLS=self.getNumLevelSets()
349 if mu_c is None or numLS < 2:
350 self.__mu_c = np.ones((numLS,numLS))
351 else:
352 mu_c=np.asarray(mu_c)
353 if mu_c.shape == (numLS,numLS):
354 if not all( [ [ mu_c[l,k] > 0. for l in xrange(k) ] for k in xrange(1,numLS) ]):
355 raise ValueError("All trade-off factors in the lower triangle of mu_c must be positive.")
356 else:
357 self.__mu_c = mu_c
358 self._new_mu=True
359 else:
360 raise ValueError("Unexpected shape %s for mu."%mu_c.shape)
361
362 def getArguments(self, m):
363 """
364 """
365 return ( grad(m),)
366
367
368 def getValue(self, m, grad_m):
369 """
370 returns the value of the cost function J with respect to m.
371
372 :rtype: ``float``
373 """
374 mu=self.__mu
375 mu_c=self.__mu_c
376 DIM=self.getDomain().getDim()
377 numLS=self.getNumLevelSets()
378
379 A=0
380 if self.__w0 is not None:
381 A+=inner(integrate(m**2 * self.__w0), mu)
382
383 if self.__w1 is not None:
384 if numLS == 1:
385 A+=integrate(inner(grad_m**2, self.__w1))*mu
386 else:
387 for k in xrange(numLS):
388 A+=mu[k]*integrate(inner(grad_m[k,:]**2,self.__w1[k,:]))
389
390 if numLS > 1:
391 for k in xrange(numLS):
392 gk=grad_m[k,:]
393 len_gk=length(gk)
394 for l in xrange(k):
395 gl=grad_m[l,:]
396 A+= mu_c[l,k] * integrate( self.__wc[l,k] * ( len_gk * length(gl) )**2 - inner(gk, gl)**2 )
397 return A/2
398
399 def getGradient(self, m, grad_m):
400 """
401 returns the gradient of the cost function J with respect to m.
402 The function returns Y_k=dPsi/dm_k and X_kj=dPsi/dm_kj
403 """
404
405 mu=self.__mu
406 mu_c=self.__mu_c
407 DIM=self.getDomain().getDim()
408 numLS=self.getNumLevelSets()
409
410 if self.__w0 is not None:
411 Y = m * self.__w0 * mu
412 else:
413 if numLS == 1:
414 Y = Scalar(0, grad_m.getFunctionSpace())
415 else:
416 Y = Data(0, (numLS,) , grad_m.getFunctionSpace())
417
418 if self.__w1 is not None:
419 X=grad_m*self.__w1
420 if numLS == 1:
421 X=grad_m* self.__w1*mu
422 else:
423 for k in xrange(numLS):
424 X[k,:]*=mu[k]
425 else:
426 X = Data(0, grad_m.getShape(), grad_m.getFunctionSpace())
427
428 # cross gradient terms:
429 if numLS > 1:
430 for k in xrange(numLS):
431 grad_m_k=grad_m[k,:]
432 l2_grad_m_k = length(grad_m_k)**2
433 for l in xrange(k):
434 grad_m_l=grad_m[l,:]
435 l2_grad_m_l = length(grad_m_l)**2
436 grad_m_lk = inner(grad_m_l, grad_m_k)
437 f= mu_c[l,k]* self.__wc[l,k]
438 X[l,:] += f * ( l2_grad_m_l * grad_m_l - grad_m_lk * grad_m_k )
439 X[k,:] += f * ( l2_grad_m_k * grad_m_k - grad_m_lk * grad_m_l )
440
441 return ArithmeticTuple(Y, X)
442
443
444
445 def getInverseHessianApproximation(self, m, r, grad_m):
446 """
447 """
448 if self._new_mu or self._update_Hessian:
449
450 self._new_mu=False
451 self._update_Hessian=False
452 mu=self.__mu
453 mu_c=self.__mu_c
454
455 DIM=self.getDomain().getDim()
456 numLS=self.getNumLevelSets()
457 if self.__w0 is not None:
458 if numLS == 1:
459 D=self.__w0 * mu
460 else:
461 D=self.getPDE().createCoefficient("D")
462 for k in xrange(numLS): D[k,k]=self.__w0[k] * mu[k]
463 self.getPDE().setValue(D=D)
464
465 A=self.getPDE().createCoefficient("A")
466 if self.__w1 is not None:
467 if numLS == 1:
468 for i in xrange(DIM): A[i,i]=self.__w1[i] * mu
469 else:
470 for k in xrange(numLS):
471 for i in xrange(DIM): A[k,i,k,i]=self.__w1[k,i] * mu[k]
472
473 if numLS > 1:
474 for k in xrange(numLS):
475 grad_m_k=grad_m[k,:]
476 l2_grad_m_k = length(grad_m_k)**2
477 o_kk=outer(grad_m_k, grad_m_k)
478 for l in xrange(k):
479 grad_m_l=grad_m[l,:]
480 l2_grad_m_l = length(grad_m_l)**2
481 i_lk = inner(grad_m_l, grad_m_k)
482 o_lk = outer(grad_m_l, grad_m_k)
483 o_kl = outer(grad_m_k, grad_m_l)
484 o_ll=outer(grad_m_l, grad_m_l)
485 f= mu_c[l,k]* self.__wc[l,k]
486
487 A[l,:,l:] += f * ( l2_grad_m_k * kronecker(DIM) - o_kk )
488 A[l,:,k:] += f * ( 2 * o_lk - o_kl - i_lk * kronecker(DIM) )
489 A[k,:,l:] += f * ( 2 * o_kl - o_lk - i_lk * kronecker(DIM) )
490 A[k,:,k:] += f * ( l2_grad_m_l * kronecker(DIM) - o_ll )
491 self.getPDE().setValue(A=A)
492 #self.getPDE().resetRightHandSideCoefficients()
493 #self.getPDE().setValue(X=r[1])
494 #print "X only: ",self.getPDE().getSolution()
495 #self.getPDE().resetRightHandSideCoefficients()
496 #self.getPDE().setValue(Y=r[0])
497 #print "Y only: ",self.getPDE().getSolution()
498
499 self.getPDE().resetRightHandSideCoefficients()
500 self.getPDE().setValue(X=r[1], Y=r[0])
501 return self.getPDE().getSolution()
502
503 def updateHessian(self):
504 """
505 notify the class to recalculate the Hessian operator
506 """
507 if not self.__useDiagonalHessianApproximation:
508 self._update_Hessian=True
509
510 def getNorm(self, m):
511 """
512 returns the norm of ``m``
513
514 :param m: level set function
515 :type m: `Data`
516 :rtype: ``float``
517 """
518 return sqrt(integrate(length(m)**2)/self.__vol_d)

  ViewVC Help
Powered by ViewVC 1.1.26