/[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 4095 - (show annotations)
Wed Dec 5 05:32:22 2012 UTC (6 years, 11 months ago) by caltinay
File MIME type: text/x-python
File size: 16614 byte(s)
A bit of doco cleanup.

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, vol
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 if numLevelSets>1:
77 raise ValueError("Currently only numLevelSets<=1 is supportered.")
78 if s0 == None and s1==None:
79 raise ValueError("Values for s0 or s1 must be given.")
80
81 self.__domain=domain
82 DIM=self.__domain.getDim()
83 self.__L2=np.asarray(boundingBoxEdgeLengths(domain))**2
84 self.__L4=np.sum(self.__L2)**2
85 self.__numLevelSets=numLevelSets
86
87 if self.__numLevelSets > 1:
88 self.__useDiagonalHessianApproximation=useDiagonalHessianApproximation
89 else:
90 self.__useDiagonalHessianApproximation=True
91 self._update_Hessian=True
92
93 self.__pde=LinearPDE(self.__domain, numEquations=self.__numLevelSets)
94 self.__pde.getSolverOptions().setTolerance(tol)
95 self.__pde.setSymmetryOn()
96 self.__pde_is_set=False
97 try:
98 self.__pde.setValue(q=location_of_set_m)
99 except IllegalCoefficientValue:
100 raise ValueError("Unable to set location of fixed level set function.")
101
102 self.__total_num_weights=2*numLevelSets+((numLevelSets-1)*numLevelSets)/2
103 self.__weight_index=[] # this is a mapping from the relevant mu-coefficients to the set of all mu-coefficients
104 # we count s0, then s1, then sc (k<l).
105 # THE S0 weighting factor
106 n=0
107 VV=vol(domain)
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*=VV/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 V > 0:
126 self.__weight_index.append(n+k)
127 s0[k]*=VV/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, 1/self.__L2))
141 if V > 0:
142 self.__weight_index.append(n)
143 s1*=VV/V
144 print "REG SCALE = ",s1
145 else:
146 raise ValueError("Unexpected shape %s for weight s1."%s)
147 else:
148 if s == (numLevelSets,DIM):
149 for k in xrange(numLevelSets):
150 for i in xrange(DIM):
151 ww=s1[k,:]
152 V=integrate(inner(ww,1/self.__L2))
153 if V > 0:
154 self.__weight_index.append(n+i)
155 s1[k,:]=ww*(VV/V)
156 else:
157 raise ValueError("Unexpected shape %s for weight s1."%s)
158
159 self.__s1=s1
160 n+=numLevelSets
161
162 # The SC weighting factor
163 if not sc is None:
164 if numLevelSets == 1 :
165 sc=None
166 else:
167 sc = interpolate(sc,self.__pde.getFunctionSpaceForCoefficient('A'))
168 s=sc.getShape()
169 if s == (numLevelSets,numLevelSets):
170 for k in xrange(numLevelSets):
171 sc[k,k]=0.
172 for l in xrange(k):
173 ww=sc[l,k]
174 V=integrate(ww)
175 if V > 0:
176 self.__weight_index.append(n+k*numLevelSets+l)
177 sc[l,k]=ww*VV/V*self.__L4
178 sc[k,l]=0
179 else:
180 raise ValueError("Unexpected shape %s for weight s0."%s)
181
182 self.__sc=sc
183 self.setWeights()
184
185 def getDomain(self):
186 """
187 return the domain of the regularization term
188 :rtype: ``Domain``
189 """
190 return self.__domain
191
192 def getNumLevelSets(self):
193 """
194 return the number of level set functions
195 :rtype: ``int``
196 """
197 return self.__numLevelSets
198
199 def getPDE(self):
200 """
201 returns the linear PDE to be solved for the Hessian Operator inverse
202 :rtype: ``LinearPDE``
203 """
204 return self.__pde
205
206 def getDualProduct(self, m, r):
207 """
208 returns the dual product of a gradient represented by X=r[1] and Y=r[0] with a level set function m:
209
210 *Y_i*m_i + X_ij*m_{i,j}*
211
212 :type m: ``esys.escript.Data``
213 :type r: ``ArithmeticTuple``
214 :rtype: float
215 """
216 A=0
217 if not r[0].isEmpty(): A+=integrate(inner(r[0], m))
218 if not r[1].isEmpty(): A+=integrate(inner(r[1], grad(m)))
219 return A
220
221 def getValue(self, m, grad_m):
222 """
223 return the value of the costfunction J
224
225
226
227 :rtype: ``float``
228 """
229 mu=self.getWeights( uncompress=True)
230 DIM=self.getDomain().getDim()
231 numLS=self.getNumLevelSets()
232
233 A=0
234 n=0
235
236 if self.__s0 is not None:
237 A+=inner(integrate(m**2*self.__s0), mu[:numLS])
238 n+=numLS
239
240 if self.__s1 is not None:
241 if numLS == 1:
242 A+=integrate(inner(grad_m**2, self.__s1))*mu[n]
243 else:
244 for k in xrange(numLS):
245 A+=mu[n+k]*integrate(inner(grad_m[k,:]**2,self.__s1[k,:]))
246 n+=numLS
247
248 if self.__sc is not None:
249 for k in xrange(numLS):
250 gk=grad_m[k,:]
251 len_gk=length(gk)
252 for l in xrange(k):
253 gl=grad_m[l,:]
254 A+= mu[n+k*numLS+l] * integrate( self.__sc[l,k] * ( len_gk * length(gl) )**2 - inner(gk, gl)**2 )
255 return A/2
256
257 def getGradient(self, m, grad_m):
258 """
259 returns the gradient of the costfunction J with respect to m. The function returns Y_k=dPsi/dm_k and
260 X_kj=dPsi/dm_kj
261 """
262
263 mu=self.getWeights( uncompress=True)
264 DIM=self.getDomain().getDim()
265 numLS=self.getNumLevelSets()
266
267 n=0
268
269 if self.__s0 is not None:
270 Y = m * self.__s0 * mu[:numLS]
271 else:
272 Y = Data()
273 n+=numLS
274
275 if self.__s1 is not None:
276 if numLS == 1:
277 X=grad_m* self.__s1*mu[n]
278 else:
279 X=self.getPDE().createCoefficient("X")
280 for k in xrange(numLS):
281 X[k,:]=mu[n+k]*grad_m[k,:]*self.__s1[k,:]
282 else:
283 X = Data()
284 n+=numLS
285 if self.__sc is not None:
286 raise NotImplementedError
287 return ArithmeticTuple(Y, X)
288
289 def getArguments(self, m):
290 """
291 """
292 return ( grad(m),)
293
294 def getInverseHessianApproximation(self, m, r, grad_m):
295 """
296 """
297 if self._new_mu or self._update_Hessian:
298
299 self._new_mu=False
300 self._update_Hessian=False
301
302 mu=self.getWeights( uncompress=True)
303 DIM=self.getDomain().getDim()
304 numLS=self.getNumLevelSets()
305 n=0
306 if self.__s0 is not None:
307 if numLS == 1:
308 D=self.__s0 * mu[n]
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 else:
320 for k in xrange(numLS):
321 for i in xrange(DIM): A[k,i,k,i]=self.__s1[k,i] * mu[n+k]
322 n+=numLS
323 if self.__sc is not None:
324 raise NotImplementedError
325 print A
326 self.getPDE().setValue(A=A)
327 self.getPDE().resetRightHandSideCoefficients()
328 self.getPDE().setValue(X=r[1])
329 print "X only: ",self.getPDE().getSolution()
330 self.getPDE().resetRightHandSideCoefficients()
331 self.getPDE().setValue(Y=r[0])
332 print "Y only: ",self.getPDE().getSolution()
333
334 self.getPDE().resetRightHandSideCoefficients()
335 self.getPDE().setValue(X=r[1], Y=r[0])
336 return self.getPDE().getSolution()
337
338 def updateHessian(self):
339 """
340 notify the class to recalculate the Hessian operator
341 """
342 if not self.__useDiagonalHessianApproximation:
343 self._update_Hessian=True
344
345 # ================ we should factor these out =============================================================
346 def getNumRelevantTerms(self):
347 """
348 returns the number of terms in the costfunction of regularization with non-zero scaling factors
349 :rtype: ``int``
350 """
351 return len(self.__weight_index)
352
353 def getNumTerms(self):
354 """
355 returns the number of terms in the costfunction of regularization with non-zero scaling factors
356 :rtype: ``int``
357 """
358 return len(self.__weight_index)
359
360 def setWeights(self, mu=None):
361 """
362 sets the values for the weighting terms in the costsfunction. Note that only values
363 correspond to entries with non-negative scaling factors are used.
364
365 :param mu: weights
366 :type mu: ``list`` of ``float`` or ``np,array``
367 """
368 if mu == None:
369 mu = np.ones(self.getNumRelevantTerms())
370 mu=np.asarray(mu)
371
372
373 if len(mu) == self.getNumRelevantTerms():
374 if not mu.shape == (self.getNumRelevantTerms(),):
375 raise ValueError("%s values are expected."%self.getNumRelevantTerms())
376 self.__mu=mu
377 else:
378 if not mu.shape == (self.__total_num_weights,):
379 raise ValueError("%s values are expected."%self.__total_num_weights)
380 self.__mu = np.zeros(self.getNumRelevantTerms())
381 for i in xrange(len(self.__weight_index)): self.__mu[i]=mu[self.__weight_index[i]]
382 self._new_mu=True
383
384 def setWeightsForS0(self, mu=None):
385 """
386 set the weights for s0-terms
387 """
388 numLS=self.getNumLevelSets()
389 my_mu=self.getWeights(uncompress=True)
390 if mu is None:
391 my_mu[:numLS]=1
392 else:
393 my_mu[:numLS]=mu
394 self.setWeights(my_mu)
395
396 def setWeightsForS1(self, mu=None):
397 """
398 set the weights for s1-terms
399 """
400 numLS=self.getNumLevelSets()
401 my_mu=self.getWeights(uncompress=True)
402 if mu is None:
403 my_mu[numLS:2*numLS]=1
404 else:
405 my_mu[numLS:2*numLS]=mu
406 self.setWeights(my_mu)
407
408 def setWeightsForSc(self, mu):
409 """
410 set the weights for s1-terms
411 """
412 numLS=self.getNumLevelSets()
413 my_mu=self.getWeights(uncompress=True)
414 if mu is None:
415 my_mu[2*numLS:]=1
416 else:
417 my_mu[2*numLS:]=mu
418 self.setWeights(my_mu)
419
420
421 def getWeights(self, uncompress=False):
422 """
423 Returns the weights for the terms in the costsfunction.
424 The first ``numLevelSets`` values used for the
425 regularization terms and the remaining values for the cross correlation terms.
426
427 """
428 if uncompress:
429 mu = np.zeros(self.__total_num_weights)
430 for i in xrange(len(self.__weight_index)): mu[self.__weight_index[i]] = self.__mu[i]
431 return mu
432 else:
433 return self.__mu
434
435 def getWeightIndex(self):
436 """
437 returns an iondex to the contributions of terms with non-zero scaling factor.
438 """
439 return self.__weight_index
440

  ViewVC Help
Powered by ViewVC 1.1.26