/[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 4074 - (show annotations)
Thu Nov 15 03:30:59 2012 UTC (6 years, 11 months 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
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 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
315 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 """
331 """
332 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
356 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
382 # ================ we should factor these out =============================================================
383 def getNumRelevantTerms(self):
384 """
385 returns the number of terms in the costfunction of regularization with non-zero scaling factors
386 :rtype: ``int``
387 """
388 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
421 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