/[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 4074 - (hide annotations)
Thu Nov 15 03:30:59 2012 UTC (7 years ago) by gross
File MIME type: text/x-python
File size: 16237 byte(s)
this will break the downundermodule: major revision of the classes.
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 getDirectionalDerivative(self, m, d, grad_m):
257     """
258     returns the directional derivative
259     """
260     mu=self.getWeights( uncompress=True)
261     DIM=self.getDomain().getDim()
262     numLS=self.getNumLevelSets()
263     grad_q=grad(q)
264    
265     A=0
266     n=0
267    
268     if self.__s0 is not None:
269     A+=inner(integrate(m*q*self.__s0), mu[:numLS])
270     n+=numLS
271    
272     if self.__s1 is not None:
273     if numLS == 1:
274     A+=integrate(inner(grad_m* self.__s1,grad_q))*mu[n]
275     else:
276     for k in xrange(numLS):
277     A+=mu[n+k]*integrate(inner(grad_m[k,:]*self.__s1[k,:],grad_q[k,:]))
278     n+=numLS*DIM
279    
280     # this needs to be checked tested:
281     if self.__sc is not None:
282     for k in xrange(numLS):
283     gm_k=grad_m[k,:]
284     gq_k=grad_q[k,:]
285     len2_gm_k=length(gm_k)**2
286     gmq_k = inner(gm_k,gq_k)
287     for l in xrange(k):
288     gm_l=grad_m[l,:]
289     gq_l=grad_q[l,:]
290     A+= mu[n+k*numLS+l] * integrate( self.__sc[l,k] * (
291     ( gmq_k * length(gm_l)**2 + len2_gm_k * inner(gm_l,gq_l)
292     - inner(gm_k, gm_l)*( inner(gm_k, gq_l) + inner(gq_k, gm_l) ) )
293     ) )
294    
295     return A
296    
297     def getGradient(self, m, grad_m):
298     """
299     returns the gradient of the costfunction J with respect to m. The function returns Y_k=dPsi/dm_k and
300     X_kj=dPsi/dm_kj
301     """
302    
303     mu=self.getWeights( uncompress=True)
304     DIM=self.getDomain().getDim()
305     numLS=self.getNumLevelSets()
306    
307     n=0
308    
309     if self.__s0 is not None:
310     Y = m * self.__s0 * mu[:numLS]
311     else:
312     Y = Data()
313     n+=numLS
314 caltinay 3946
315 gross 4074 if self.__s1 is not None:
316     if numLS == 1:
317     X=grad_m* self.__s1*mu[n]
318     else:
319     X=self.getPDE().createCoefficient("X")
320     for k in xrange(numLS):
321     X[k,:]=mu[n+k]*grad_m[k,:]*self.__s1[k,:]
322     else:
323     X = Data()
324     n+=numLS
325     if self.__sc is not None:
326     raise NotImplementedError
327     return ArithmeticTuple(Y, X)
328    
329     def getArguments(self, m):
330 caltinay 4007 """
331     """
332 gross 4074 return ( grad(m),)
333    
334     def getInverseHessianApproximation(self, m, r, grad_m):
335     """
336     """
337     if self._new_mu or self._update_Hessian:
338    
339     self._new_mu=False
340     self._update_Hessian=False
341    
342     mu=self.getWeights( uncompress=True)
343     DIM=self.getDomain().getDim()
344     numLS=self.getNumLevelSets()
345     n=0
346     if self.__s0 is not None:
347     if numLS == 1:
348     D=self.__s0 * mu[n]
349     print "D =",D
350     else:
351     D=self.getPDE().createCoefficient("D")
352     for k in xrange(numLS): D[k,k]=self.__s0[k] * mu[n+k]
353     self.getPDE().setValue(D=D)
354     n+=numLS
355 caltinay 3946
356 gross 4074 A=self.getPDE().createCoefficient("A")
357     if self.__s1 is not None:
358     if numLS == 1:
359     for i in xrange(DIM): A[i,i]=self.__s1[i] * mu[n]
360     print "A =",A
361     else:
362     for k in xrange(numLS):
363     for i in xrange(DIM): A[k,i,k,i]=self.__s1[k,i] * mu[n+k]
364     n+=numLS
365     if self.__sc is not None:
366     raise NotImplementedError
367     self.getPDE().setValue(A=A)
368     self.getPDE().resetRightHandSideCoefficients()
369     print A
370     print r[1]
371     print r[0]
372     self.getPDE().setValue(X=r[1], Y=r[0])
373     return self.getPDE().getSolution()
374    
375     def updateHessian(self):
376     """
377     notify the class to recalculate the Hessian operator
378     """
379     if not self.__useDiagonalHessianApproximation:
380     self._update_Hessian=True
381 caltinay 3946
382 gross 4074 # ================ we should factor these out =============================================================
383     def getNumRelevantTerms(self):
384 caltinay 4007 """
385 gross 4074 returns the number of terms in the costfunction of regularization with non-zero scaling factors
386     :rtype: ``int``
387 caltinay 4007 """
388 gross 4074 return len(self.__weight_index)
389    
390     def getNumTerms(self):
391     """
392     returns the number of terms in the costfunction of regularization with non-zero scaling factors
393     :rtype: ``int``
394     """
395     return len(self.__weight_index)
396    
397     def setWeights(self, mu=None):
398     """
399     sets the values for the weighting terms in the costsfunction. Note that only values
400     correspond to entries with non-negative scaling factors are used.
401    
402     :param mu: weights
403     :type mu: ``list`` of ``float`` or ``np,array``
404     """
405     if mu == None:
406     mu = np.ones(self.getNumRelevantTerms())
407     mu=np.asarray(mu)
408    
409    
410     if len(mu) == self.getNumRelevantTerms():
411     if not mu.shape == (self.getNumRelevantTerms(),):
412     raise ValueError("%s values are expected."%self.getNumRelevantTerms())
413     self.__mu=mu
414     else:
415     if not mu.shape == (self.__total_num_weights,):
416     raise ValueError("%s values are expected."%self.__total_num_weights)
417     self.__mu = np.zeros(self.getNumRelevantTerms())
418     for i in xrange(len(self.__weight_index)): self.__mu[i]=mu[self.__weight_index[i]]
419     self._new_mu=True
420 caltinay 3946
421 gross 4074 def setWeightsForS0(self, mu=None):
422     """
423     set the weights for s0-terms
424     """
425     numLS=self.getNumLevelSets()
426     my_mu=self.getWeights(uncompress=True)
427     if mu is None:
428     my_mu[:numLS]=1
429     else:
430     my_mu[:numLS]=mu
431     self.setWeights(my_mu)
432    
433     def setWeightsForS1(self, mu=None):
434     """
435     set the weights for s1-terms
436     """
437     numLS=self.getNumLevelSets()
438     my_mu=self.getWeights(uncompress=True)
439     if mu is None:
440     my_mu[numLS:2*numLS]=1
441     else:
442     my_mu[numLS:2*numLS]=mu
443     self.setWeights(my_mu)
444    
445     def setWeightsForSc(self, mu):
446     """
447     set the weights for s1-terms
448     """
449     numLS=self.getNumLevelSets()
450     my_mu=self.getWeights(uncompress=True)
451     if mu is None:
452     my_mu[2*numLS:]=1
453     else:
454     my_mu[2*numLS:]=mu
455     self.setWeights(my_mu)
456    
457    
458     def getWeights(self, uncompress=False):
459     """
460     Returns the weights for the terms in the costsfunction.
461     The first ``numLevelSets`` values used for the
462     regularization terms and the remaining values for the cross correlation terms.
463    
464     :type mu: ``list`` of ``float``
465     """
466     if uncompress:
467     mu = np.zeros(self.__total_num_weights)
468     for i in xrange(len(self.__weight_index)): mu[self.__weight_index[i]] = self.__mu[i]
469     return mu
470     else:
471     return self.__mu
472    
473     def getWeightIndex(self):
474     """
475     returns an iondex to the contributions of terms with non-zero scaling factor.
476     """
477     return self.__weight_index
478    

  ViewVC Help
Powered by ViewVC 1.1.26