/[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 4080 - (hide annotations)
Mon Nov 19 01:45:38 2012 UTC (7 years ago) by jfenwick
File MIME type: text/x-python
File size: 15229 byte(s)
Epydoc fixes

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 gross 4079 from esys.escript import Data, grad, inner, integrate, kronecker, boundingBoxEdgeLengths, interpolate, vol
31 gross 4074 from esys.escript.pdetools import ArithmeticTuple
32 caltinay 3946
33 gross 4074 class Regularization(CostFunction):
34 caltinay 3946 """
35 jfenwick 4080 The regularization term for the level set function ``m`` within the cost function J for an inversion:
36 gross 4074
37 jfenwick 4080 *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 gross 4074
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 jfenwick 4080 :type s0: ``Scalar`` if ``numLevelSets``==1 or `Data` object of shape ('numLevelSets`,) if numLevelSets > 1
61 gross 4074 :param s1: scaling factor for the grad(m_i) terms. If not set zero is assumed.
62 jfenwick 4080 :type s1: ``Vector`` if ``numLevelSets``==1 or `Data` object of shape (`numLevelSets`,DIM) if numLevelSets > 1.
63 gross 4074 :param sc: scaling factor for the cross correlation terms. If not set zero is assumed. Used for the case if numLevelSets > 1 only.
64 jfenwick 4080 values ``sc[l,k]`` in the lower triangle (l<k) are used only.
65 gross 4074 :type sc: `Data` object of shape (`numLevelSets`,`numLevelSets`)
66 jfenwick 4080 :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 gross 4074 :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 gross 4079 VV=vol(domain)
109 gross 4074 if not s0 is None:
110     s0 = interpolate(s0,self.__pde.getFunctionSpaceForCoefficient('D'))
111     s=s0.getShape()
112     if numLevelSets == 1 :
113     if s == () :
114     V=integrate(s0)
115     if V > 0:
116     self.__weight_index.append(n)
117 gross 4079 s0*=VV/V
118 gross 4074 else:
119     s0=None
120     else:
121     raise ValueError("Unexpected shape %s for weight s0."%s)
122     else:
123     if s == (numLevelSets,):
124     for k in xrange(numLevelSets):
125     V=integrate(s0[k])
126 gross 4079 if V > 0:
127 gross 4074 self.__weight_index.append(n+k)
128 gross 4079 s0[k]*=VV/V
129 gross 4074 else:
130     raise ValueError("Unexpected shape %s for weight s0."%s)
131     self.__s0=s0
132     n+=numLevelSets
133    
134     # The S1 weighting factor
135     if not s1 is None:
136     s1 = interpolate(s1,self.__pde.getFunctionSpaceForCoefficient('A'))
137     s=s1.getShape()
138    
139     if numLevelSets == 1 :
140     if s == (DIM,) :
141 gross 4079 V=integrate(inner(s1, 1/self.__L2))
142 gross 4074 if V > 0:
143     self.__weight_index.append(n)
144 gross 4079 s1*=VV/V
145     print "REG SCALE = ",s1
146 gross 4074 else:
147     raise ValueError("Unexpected shape %s for weight s1."%s)
148     else:
149     if s == (numLevelSets,DIM):
150     for k in xrange(numLevelSets):
151     for i in xrange(DIM):
152     ww=s1[k,:]
153 gross 4079 V=integrate(inner(ww,1/self.__L2))
154 gross 4074 if V > 0:
155     self.__weight_index.append(n+i)
156 gross 4079 s1[k,:]=ww*(VV/V)
157 gross 4074 else:
158     raise ValueError("Unexpected shape %s for weight s1."%s)
159    
160     self.__s1=s1
161     n+=numLevelSets
162 caltinay 3946
163 gross 4074 # The SC weighting factor
164     if not sc is None:
165     if numLevelSets == 1 :
166     sc=None
167     else:
168     sc = interpolate(sc,self.__pde.getFunctionSpaceForCoefficient('A'))
169     s=sc.getShape()
170     if s == (numLevelSets,numLevelSets):
171     for k in xrange(numLevelSets):
172     sc[k,k]=0.
173     for l in xrange(k):
174     ww=sc[l,k]
175     V=integrate(ww)
176     if V > 0:
177     self.__weight_index.append(n+k*numLevelSets+l)
178 gross 4079 sc[l,k]=ww*VV/V*self.__L4
179 gross 4074 sc[k,l]=0
180     else:
181     raise ValueError("Unexpected shape %s for weight s0."%s)
182    
183     self.__sc=sc
184     self.setWeights()
185    
186     def getDomain(self):
187 caltinay 3946 """
188 gross 4074 return the domain of the regularization term
189     :rtype: ``Domain``
190 caltinay 3946 """
191 gross 4074 return self.__domain
192    
193     def getNumLevelSets(self):
194     """
195     return the number of level set functions
196     :rtype: ``int``
197     """
198     return self.__numLevelSets
199 caltinay 3946
200 gross 4074 def getPDE(self):
201 caltinay 4007 """
202 gross 4074 returns the linear PDE to be solved for the Hessian Operator inverse
203     :rtype: ``LinearPDE``
204 caltinay 4007 """
205 gross 4074 return self.__pde
206    
207     def getDualProduct(self, m, r):
208     """
209     returns the dual product of a gradient represented by X=r[1] and Y=r[0] with a level set function m:
210    
211     *Y_i*m_i + X_ij*m_{i,j}*
212    
213     :type m: ``esys.escript.Data``
214     :type r: ``ArithmeticTuple``
215     :rtype: float
216     """
217     A=0
218     if not r[0].isEmpty(): A+=integrate(inner(r[0], m))
219     if not r[1].isEmpty(): A+=integrate(inner(r[1], grad(m)))
220     return A
221 caltinay 3946
222 gross 4074 def getValue(self, m, grad_m):
223 caltinay 4007 """
224 gross 4074 return the value of the costfunction J
225    
226    
227    
228     :rtype: ``float``
229 caltinay 4007 """
230 gross 4074 mu=self.getWeights( uncompress=True)
231     DIM=self.getDomain().getDim()
232     numLS=self.getNumLevelSets()
233    
234 caltinay 3946 A=0
235 gross 4074 n=0
236    
237     if self.__s0 is not None:
238     A+=inner(integrate(m**2*self.__s0), mu[:numLS])
239     n+=numLS
240    
241     if self.__s1 is not None:
242     if numLS == 1:
243     A+=integrate(inner(grad_m**2, self.__s1))*mu[n]
244     else:
245     for k in xrange(numLS):
246     A+=mu[n+k]*integrate(inner(grad_m[k,:]**2,self.__s1[k,:]))
247     n+=numLS
248    
249     if self.__sc is not None:
250     for k in xrange(numLS):
251     gk=grad_m[k,:]
252     len_gk=length(gk)
253     for l in xrange(k):
254     gl=grad_m[l,:]
255     A+= mu[n+k*numLS+l] * integrate( self.__sc[l,k] * ( len_gk * length(gl) )**2 - inner(gk, gl)**2 )
256     return A/2
257    
258     def getGradient(self, m, grad_m):
259     """
260     returns the gradient of the costfunction J with respect to m. The function returns Y_k=dPsi/dm_k and
261     X_kj=dPsi/dm_kj
262     """
263    
264     mu=self.getWeights( uncompress=True)
265     DIM=self.getDomain().getDim()
266     numLS=self.getNumLevelSets()
267    
268     n=0
269    
270     if self.__s0 is not None:
271     Y = m * self.__s0 * mu[:numLS]
272     else:
273     Y = Data()
274     n+=numLS
275 caltinay 3946
276 gross 4074 if self.__s1 is not None:
277     if numLS == 1:
278     X=grad_m* self.__s1*mu[n]
279     else:
280     X=self.getPDE().createCoefficient("X")
281     for k in xrange(numLS):
282     X[k,:]=mu[n+k]*grad_m[k,:]*self.__s1[k,:]
283     else:
284     X = Data()
285     n+=numLS
286     if self.__sc is not None:
287     raise NotImplementedError
288     return ArithmeticTuple(Y, X)
289    
290     def getArguments(self, m):
291 caltinay 4007 """
292     """
293 gross 4074 return ( grad(m),)
294    
295     def getInverseHessianApproximation(self, m, r, grad_m):
296     """
297     """
298     if self._new_mu or self._update_Hessian:
299    
300     self._new_mu=False
301     self._update_Hessian=False
302    
303     mu=self.getWeights( uncompress=True)
304     DIM=self.getDomain().getDim()
305     numLS=self.getNumLevelSets()
306     n=0
307     if self.__s0 is not None:
308     if numLS == 1:
309     D=self.__s0 * mu[n]
310     else:
311     D=self.getPDE().createCoefficient("D")
312     for k in xrange(numLS): D[k,k]=self.__s0[k] * mu[n+k]
313     self.getPDE().setValue(D=D)
314     n+=numLS
315 caltinay 3946
316 gross 4074 A=self.getPDE().createCoefficient("A")
317     if self.__s1 is not None:
318     if numLS == 1:
319     for i in xrange(DIM): A[i,i]=self.__s1[i] * mu[n]
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 gross 4079 print A
327 gross 4074 self.getPDE().setValue(A=A)
328     self.getPDE().resetRightHandSideCoefficients()
329 gross 4079 self.getPDE().setValue(X=r[1])
330     print "X only: ",self.getPDE().getSolution()
331     self.getPDE().resetRightHandSideCoefficients()
332     self.getPDE().setValue(Y=r[0])
333     print "Y only: ",self.getPDE().getSolution()
334    
335     self.getPDE().resetRightHandSideCoefficients()
336 gross 4074 self.getPDE().setValue(X=r[1], Y=r[0])
337     return self.getPDE().getSolution()
338    
339     def updateHessian(self):
340     """
341     notify the class to recalculate the Hessian operator
342     """
343     if not self.__useDiagonalHessianApproximation:
344     self._update_Hessian=True
345 caltinay 3946
346 gross 4074 # ================ we should factor these out =============================================================
347     def getNumRelevantTerms(self):
348 caltinay 4007 """
349 gross 4074 returns the number of terms in the costfunction of regularization with non-zero scaling factors
350     :rtype: ``int``
351 caltinay 4007 """
352 gross 4074 return len(self.__weight_index)
353    
354     def getNumTerms(self):
355     """
356     returns the number of terms in the costfunction of regularization with non-zero scaling factors
357     :rtype: ``int``
358     """
359     return len(self.__weight_index)
360    
361     def setWeights(self, mu=None):
362     """
363     sets the values for the weighting terms in the costsfunction. Note that only values
364     correspond to entries with non-negative scaling factors are used.
365    
366     :param mu: weights
367     :type mu: ``list`` of ``float`` or ``np,array``
368     """
369     if mu == None:
370     mu = np.ones(self.getNumRelevantTerms())
371     mu=np.asarray(mu)
372    
373    
374     if len(mu) == self.getNumRelevantTerms():
375     if not mu.shape == (self.getNumRelevantTerms(),):
376     raise ValueError("%s values are expected."%self.getNumRelevantTerms())
377     self.__mu=mu
378     else:
379     if not mu.shape == (self.__total_num_weights,):
380     raise ValueError("%s values are expected."%self.__total_num_weights)
381     self.__mu = np.zeros(self.getNumRelevantTerms())
382     for i in xrange(len(self.__weight_index)): self.__mu[i]=mu[self.__weight_index[i]]
383     self._new_mu=True
384 caltinay 3946
385 gross 4074 def setWeightsForS0(self, mu=None):
386     """
387     set the weights for s0-terms
388     """
389     numLS=self.getNumLevelSets()
390     my_mu=self.getWeights(uncompress=True)
391     if mu is None:
392     my_mu[:numLS]=1
393     else:
394     my_mu[:numLS]=mu
395     self.setWeights(my_mu)
396    
397     def setWeightsForS1(self, mu=None):
398     """
399     set the weights for s1-terms
400     """
401     numLS=self.getNumLevelSets()
402     my_mu=self.getWeights(uncompress=True)
403     if mu is None:
404     my_mu[numLS:2*numLS]=1
405     else:
406     my_mu[numLS:2*numLS]=mu
407     self.setWeights(my_mu)
408    
409     def setWeightsForSc(self, mu):
410     """
411     set the weights for s1-terms
412     """
413     numLS=self.getNumLevelSets()
414     my_mu=self.getWeights(uncompress=True)
415     if mu is None:
416     my_mu[2*numLS:]=1
417     else:
418     my_mu[2*numLS:]=mu
419     self.setWeights(my_mu)
420    
421    
422     def getWeights(self, uncompress=False):
423     """
424     Returns the weights for the terms in the costsfunction.
425     The first ``numLevelSets`` values used for the
426     regularization terms and the remaining values for the cross correlation terms.
427    
428     """
429     if uncompress:
430     mu = np.zeros(self.__total_num_weights)
431     for i in xrange(len(self.__weight_index)): mu[self.__weight_index[i]] = self.__mu[i]
432     return mu
433     else:
434     return self.__mu
435    
436     def getWeightIndex(self):
437     """
438     returns an iondex to the contributions of terms with non-zero scaling factor.
439     """
440     return self.__weight_index
441    

  ViewVC Help
Powered by ViewVC 1.1.26