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

  ViewVC Help
Powered by ViewVC 1.1.26