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

  ViewVC Help
Powered by ViewVC 1.1.26