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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4122 - (hide annotations)
Thu Dec 20 05:42:35 2012 UTC (6 years, 10 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 caltinay 3946
2 jfenwick 3981 ##############################################################################
3 caltinay 3946 #
4     # Copyright (c) 2003-2012 by University of Queensland
5 jfenwick 3981 # http://www.uq.edu.au
6 caltinay 3946 #
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 jfenwick 3981 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12     # Development since 2012 by School of Earth Sciences
13     #
14     ##############################################################################
15 caltinay 3946
16     __copyright__="""Copyright (c) 2003-2012 by University of Queensland
17 jfenwick 3981 http://www.uq.edu.au
18 caltinay 3946 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 caltinay 3947 __all__ = ['Regularization']
24    
25 gross 4074 from costfunctions import CostFunction
26    
27 caltinay 3946 import numpy as np
28 gross 4122 from esys.escript import Data, Scalar, grad, inner, integrate, interpolate, kronecker, boundingBoxEdgeLengths, vol, sqrt, length
29 gross 4074 from esys.escript.linearPDEs import LinearPDE, IllegalCoefficientValue
30     from esys.escript.pdetools import ArithmeticTuple
31 caltinay 3946
32 gross 4074 class Regularization(CostFunction):
33 caltinay 3946 """
34 caltinay 4097 The regularization term for the level set function ``m`` within the cost
35     function J for an inversion:
36    
37 gross 4099 *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 caltinay 4097
39 gross 4099 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 caltinay 4097 integrals over the domain are constant:
43    
44 gross 4099 *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 caltinay 4097
47 caltinay 3946 """
48 caltinay 4097 def __init__(self, domain, numLevelSets=1,
49 gross 4099 w0=None, w1=None, wc=None,
50 caltinay 4097 location_of_set_m=Data(),
51 gross 4099 useDiagonalHessianApproximation=False, tol=1e-8,
52     scale=None, scale_c=None):
53 caltinay 4007 """
54 caltinay 4095 initialization.
55 caltinay 4097
56     :param domain: domain
57 gross 4074 :type domain: `Domain`
58     :param numLevelSets: number of level sets
59     :type numLevelSets: ``int``
60 gross 4099 :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 caltinay 4097 (``numLevelSets`` ,) if ``numLevelSets`` > 1
63 gross 4099 :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 caltinay 4097 (``numLevelSets`` , DIM) if ``numLevelSets`` > 1
66 gross 4099 :param wc: weighting factor for the cross gradient terms. If not set
67 caltinay 4097 zero is assumed. Used for the case if ``numLevelSets`` > 1
68 gross 4099 only. Only values ``wc[l,k]`` in the lower triangle (l<k)
69 caltinay 4097 are used.
70 gross 4099 :type wc: `Data` object of shape (``numLevelSets`` , ``numLevelSets``)
71 caltinay 4097 :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 gross 4099 :param useDiagonalHessianApproximation: if True cross gradient terms
76 caltinay 4097 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 gross 4074 :type useDiagonalHessianApproximation: ``bool``
82 caltinay 4097 :param tol: tolerance when solving the PDE for the inverse of the
83     Hessian Operator
84 gross 4074 :type tol: positive ``float``
85 gross 4099
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 caltinay 4095
96 caltinay 4007 """
97 gross 4099 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 caltinay 4097
102 caltinay 3946 self.__domain=domain
103 gross 4074 DIM=self.__domain.getDim()
104     self.__numLevelSets=numLevelSets
105 gross 4099
106 caltinay 4095 self.__pde=LinearPDE(self.__domain, numEquations=self.__numLevelSets)
107 gross 4074 self.__pde.getSolverOptions().setTolerance(tol)
108     self.__pde.setSymmetryOn()
109     try:
110 caltinay 4097 self.__pde.setValue(q=location_of_set_m)
111     except IllegalCoefficientValue:
112     raise ValueError("Unable to set location of fixed level set function.")
113 gross 4099
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 gross 4074 else:
129 gross 4099 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 gross 4102 f = scale/A
180 gross 4099 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 gross 4102 f = scale[k]/A
195 gross 4099 if w0 is not None:
196     w0[k]*=f
197     if w1 is not None:
198     w1[k,:]*=f
199 caltinay 4095 else:
200 gross 4099 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 gross 4102 f = scale_c[l,k]/A
207 gross 4099 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 caltinay 4097
216    
217 gross 4099 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 caltinay 4097
224 gross 4099 self.__num_tradeoff_factors=numLevelSets+((numLevelSets-1)*numLevelSets)/2
225     self.setTradeOffFactors()
226 gross 4102 self.__vol_d=vol(self.__domain)
227    
228 gross 4074 def getDomain(self):
229 caltinay 3946 """
230 caltinay 4097 returns the domain of the regularization term
231    
232 gross 4074 :rtype: ``Domain``
233 caltinay 3946 """
234 gross 4074 return self.__domain
235 caltinay 4097
236 gross 4074 def getNumLevelSets(self):
237     """
238 caltinay 4097 returns the number of level set functions
239    
240 gross 4074 :rtype: ``int``
241     """
242     return self.__numLevelSets
243 caltinay 3946
244 gross 4074 def getPDE(self):
245 caltinay 4007 """
246 gross 4074 returns the linear PDE to be solved for the Hessian Operator inverse
247 caltinay 4097
248     :rtype: `LinearPDE`
249 caltinay 4007 """
250 gross 4074 return self.__pde
251 gross 4122
252 gross 4074 def getDualProduct(self, m, r):
253     """
254 caltinay 4097 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 gross 4074 *Y_i*m_i + X_ij*m_{i,j}*
258 caltinay 4097
259     :type m: `Data`
260     :type r: `ArithmeticTuple`
261     :rtype: ``float``
262 gross 4074 """
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 caltinay 4095 return A
267 gross 4099 def getNumTradeOffFactors(self):
268     """
269     returns the number of trade-off factors being used.
270 caltinay 3946
271 gross 4099 :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 gross 4074 def getValue(self, m, grad_m):
369 caltinay 4007 """
370 caltinay 4097 returns the value of the cost function J with respect to m.
371    
372 gross 4074 :rtype: ``float``
373 caltinay 4007 """
374 gross 4099 mu=self.__mu
375     mu_c=self.__mu_c
376 gross 4074 DIM=self.getDomain().getDim()
377     numLS=self.getNumLevelSets()
378 caltinay 4097
379 caltinay 3946 A=0
380 gross 4099 if self.__w0 is not None:
381     A+=inner(integrate(m**2 * self.__w0), mu)
382 caltinay 4097
383 gross 4099 if self.__w1 is not None:
384 caltinay 4095 if numLS == 1:
385 gross 4099 A+=integrate(inner(grad_m**2, self.__w1))*mu
386 caltinay 4095 else:
387     for k in xrange(numLS):
388 gross 4099 A+=mu[k]*integrate(inner(grad_m[k,:]**2,self.__w1[k,:]))
389    
390     if numLS > 1:
391 caltinay 4097 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 gross 4099 A+= mu_c[l,k] * integrate( self.__wc[l,k] * ( len_gk * length(gl) )**2 - inner(gk, gl)**2 )
397 gross 4074 return A/2
398 caltinay 4097
399 gross 4074 def getGradient(self, m, grad_m):
400     """
401 caltinay 4097 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 gross 4074 """
404 caltinay 4097
405 gross 4099 mu=self.__mu
406     mu_c=self.__mu_c
407 gross 4074 DIM=self.getDomain().getDim()
408     numLS=self.getNumLevelSets()
409 caltinay 4097
410 gross 4099 if self.__w0 is not None:
411     Y = m * self.__w0 * mu
412 caltinay 4095 else:
413 gross 4122 if numLS == 1:
414     Y = Scalar(0, grad_m.getFunctionSpace())
415     else:
416     Y = Data(0, (numLS,) , grad_m.getFunctionSpace())
417 caltinay 3946
418 gross 4099 if self.__w1 is not None:
419 gross 4122 X=grad_m*self.__w1
420 caltinay 4095 if numLS == 1:
421 gross 4099 X=grad_m* self.__w1*mu
422 caltinay 4095 else:
423     for k in xrange(numLS):
424 gross 4099 X[k,:]*=mu[k]
425 caltinay 4095 else:
426 gross 4122 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 gross 4074 return ArithmeticTuple(Y, X)
442 caltinay 4097
443    
444 gross 4099
445 gross 4074 def getInverseHessianApproximation(self, m, r, grad_m):
446     """
447     """
448     if self._new_mu or self._update_Hessian:
449 caltinay 4097
450 caltinay 4095 self._new_mu=False
451     self._update_Hessian=False
452 gross 4099 mu=self.__mu
453     mu_c=self.__mu_c
454 caltinay 4097
455 gross 4074 DIM=self.getDomain().getDim()
456     numLS=self.getNumLevelSets()
457 gross 4099 if self.__w0 is not None:
458 caltinay 4095 if numLS == 1:
459 gross 4099 D=self.__w0 * mu
460 caltinay 4097 else:
461 gross 4074 D=self.getPDE().createCoefficient("D")
462 gross 4099 for k in xrange(numLS): D[k,k]=self.__w0[k] * mu[k]
463 caltinay 4095 self.getPDE().setValue(D=D)
464 caltinay 3946
465 caltinay 4095 A=self.getPDE().createCoefficient("A")
466 gross 4099 if self.__w1 is not None:
467 caltinay 4095 if numLS == 1:
468 gross 4099 for i in xrange(DIM): A[i,i]=self.__w1[i] * mu
469 caltinay 4095 else:
470 caltinay 4097 for k in xrange(numLS):
471 gross 4099 for i in xrange(DIM): A[k,i,k,i]=self.__w1[k,i] * mu[k]
472    
473     if numLS > 1:
474 gross 4122 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 gross 4074 self.getPDE().setValue(A=A)
492 gross 4102 #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 caltinay 4097
499 gross 4079 self.getPDE().resetRightHandSideCoefficients()
500 gross 4074 self.getPDE().setValue(X=r[1], Y=r[0])
501     return self.getPDE().getSolution()
502 caltinay 4097
503 gross 4074 def updateHessian(self):
504     """
505     notify the class to recalculate the Hessian operator
506     """
507 caltinay 4095 if not self.__useDiagonalHessianApproximation:
508     self._update_Hessian=True
509 gross 4100
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