/[escript]/trunk/downunder/py_src/regularizations.py
ViewVC logotype

Contents of /trunk/downunder/py_src/regularizations.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4076 - (show annotations)
Thu Nov 15 03:45:24 2012 UTC (7 years ago) by gross
File MIME type: text/x-python
File size: 14969 byte(s)
    def getDirectionalDerivative(self, m, d, grad_m):
        
1
2 ##############################################################################
3 #
4 # Copyright (c) 2003-2012 by University of Queensland
5 # http://www.uq.edu.au
6 #
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 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 # Development since 2012 by School of Earth Sciences
13 #
14 ##############################################################################
15
16 __copyright__="""Copyright (c) 2003-2012 by University of Queensland
17 http://www.uq.edu.au
18 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 __all__ = ['Regularization']
24
25 from costfunctions import CostFunction
26
27
28 import numpy as np
29 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
33 class Regularization(CostFunction):
34 """
35 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 """
48 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 """
53 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 """
77 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 self.__domain=domain
83 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 else:
91 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
161 # 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 """
186 return the domain of the regularization term
187 :rtype: ``Domain``
188 """
189 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
198 def getPDE(self):
199 """
200 returns the linear PDE to be solved for the Hessian Operator inverse
201 :rtype: ``LinearPDE``
202 """
203 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
220 def getValue(self, m, grad_m):
221 """
222 return the value of the costfunction J
223
224
225
226 :rtype: ``float``
227 """
228 mu=self.getWeights( uncompress=True)
229 DIM=self.getDomain().getDim()
230 numLS=self.getNumLevelSets()
231
232 A=0
233 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
274 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 """
290 """
291 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
315 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
341 # ================ we should factor these out =============================================================
342 def getNumRelevantTerms(self):
343 """
344 returns the number of terms in the costfunction of regularization with non-zero scaling factors
345 :rtype: ``int``
346 """
347 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
380 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