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

  ViewVC Help
Powered by ViewVC 1.1.26