/[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 4079 - (show annotations)
Fri Nov 16 07:59:01 2012 UTC (7 years ago) by gross
File MIME type: text/x-python
File size: 15248 byte(s)
some modifications to scaling in downunder. still not perfect.
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 """
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 VV=vol(domain)
109 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 s0*=VV/V
118 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 if V > 0:
127 self.__weight_index.append(n+k)
128 s0[k]*=VV/V
129 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 V=integrate(inner(s1, 1/self.__L2))
142 if V > 0:
143 self.__weight_index.append(n)
144 s1*=VV/V
145 print "REG SCALE = ",s1
146 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 V=integrate(inner(ww,1/self.__L2))
154 if V > 0:
155 self.__weight_index.append(n+i)
156 s1[k,:]=ww*(VV/V)
157 else:
158 raise ValueError("Unexpected shape %s for weight s1."%s)
159
160 self.__s1=s1
161 n+=numLevelSets
162
163 # 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 sc[l,k]=ww*VV/V*self.__L4
179 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 """
188 return the domain of the regularization term
189 :rtype: ``Domain``
190 """
191 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
200 def getPDE(self):
201 """
202 returns the linear PDE to be solved for the Hessian Operator inverse
203 :rtype: ``LinearPDE``
204 """
205 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
222 def getValue(self, m, grad_m):
223 """
224 return the value of the costfunction J
225
226
227
228 :rtype: ``float``
229 """
230 mu=self.getWeights( uncompress=True)
231 DIM=self.getDomain().getDim()
232 numLS=self.getNumLevelSets()
233
234 A=0
235 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
276 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 """
292 """
293 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
316 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 print A
327 self.getPDE().setValue(A=A)
328 self.getPDE().resetRightHandSideCoefficients()
329 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 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
346 # ================ we should factor these out =============================================================
347 def getNumRelevantTerms(self):
348 """
349 returns the number of terms in the costfunction of regularization with non-zero scaling factors
350 :rtype: ``int``
351 """
352 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
385 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 :type mu: ``list`` of ``float``
429 """
430 if uncompress:
431 mu = np.zeros(self.__total_num_weights)
432 for i in xrange(len(self.__weight_index)): mu[self.__weight_index[i]] = self.__mu[i]
433 return mu
434 else:
435 return self.__mu
436
437 def getWeightIndex(self):
438 """
439 returns an iondex to the contributions of terms with non-zero scaling factor.
440 """
441 return self.__weight_index
442

  ViewVC Help
Powered by ViewVC 1.1.26