/[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 4285 - (show annotations)
Thu Mar 7 01:08:43 2013 UTC (6 years, 8 months ago) by caltinay
File MIME type: text/x-python
File size: 20459 byte(s)
Corrected some spelling.

1
2 ##############################################################################
3 #
4 # Copyright (c) 2003-2013 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-2013 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 import numpy as np
28 from esys.escript import Function, outer, Data, Scalar, grad, inner, integrate, interpolate, kronecker, boundingBoxEdgeLengths, vol, sqrt, length
29 from esys.escript.linearPDEs import LinearPDE, IllegalCoefficientValue
30 from esys.escript.pdetools import ArithmeticTuple
31
32
33 class Regularization(CostFunction):
34 """
35 The regularization term for the level set function ``m`` within the cost
36 function J for an inversion:
37
38 *J(m)=1/2 * sum_k integrate( mu[k] * ( w0[k] * m_k**2 * w1[k,i] * m_{k,i}**2) + sum_l<k mu_c[l,k] wc[l,k] * | curl(m_k) x curl(m_l) |^2*
39
40 where w0[k], w1[k,i] and wc[k,l] are non-negative weighting factors and
41 mu[k] and mu_c[l,k] are trade-off factors which may be altered
42 during the inversion. The weighting factors are normalized such that their
43 integrals over the domain are constant:
44
45 *integrate(w0[k] + inner(w1[k,:],1/L[:]**2))=scale[k]* volume(domain)*
46 *integrate(wc[l,k]*1/L**4)=scale_c[k]* volume(domain) *
47
48 """
49 def __init__(self, domain, numLevelSets=1,
50 w0=None, w1=None, wc=None,
51 location_of_set_m=Data(),
52 useDiagonalHessianApproximation=False, tol=1e-8,
53 scale=None, scale_c=None):
54 """
55 initialization.
56
57 :param domain: domain
58 :type domain: `Domain`
59 :param numLevelSets: number of level sets
60 :type numLevelSets: ``int``
61 :param w0: weighting factor for the m**2 term. If not set zero is assumed.
62 :type w0: ``Scalar`` if ``numLevelSets`` == 1 or `Data` object of shape
63 (``numLevelSets`` ,) if ``numLevelSets`` > 1
64 :param w1: weighting factor for the grad(m_i) terms. If not set zero is assumed
65 :type w1: ``Vector`` if ``numLevelSets`` == 1 or `Data` object of shape
66 (``numLevelSets`` , DIM) if ``numLevelSets`` > 1
67 :param wc: weighting factor for the cross gradient terms. If not set
68 zero is assumed. Used for the case if ``numLevelSets`` > 1
69 only. Only values ``wc[l,k]`` in the lower triangle (l<k)
70 are used.
71 :type wc: `Data` object of shape (``numLevelSets`` , ``numLevelSets``)
72 :param location_of_set_m: marks location of zero values of the level
73 set function ``m`` by a positive entry.
74 :type location_of_set_m: ``Scalar`` if ``numLevelSets`` == 1 or `Data`
75 object of shape (``numLevelSets`` ,) if ``numLevelSets`` > 1
76 :param useDiagonalHessianApproximation: if True cross gradient terms
77 between level set components are ignored when calculating
78 approximations of the inverse of the Hessian Operator.
79 This can speed-up the calculation of the inverse but may
80 lead to an increase of the number of iteration steps in the
81 inversion.
82 :type useDiagonalHessianApproximation: ``bool``
83 :param tol: tolerance when solving the PDE for the inverse of the
84 Hessian Operator
85 :type tol: positive ``float``
86
87 :param scale: weighting factor for level set function variation terms. If not set one is used.
88 :type scale: ``Scalar`` if ``numLevelSets`` == 1 or `Data` object of shape
89 (``numLevelSets`` ,) if ``numLevelSets`` > 1
90 :param scale_c: scale for the cross gradient terms. If not set
91 one is assumed. Used for the case if ``numLevelSets`` > 1
92 only. Only values ``scale_c[l,k]`` in the lower triangle (l<k)
93 are used.
94 :type scale_c: `Data` object of shape (``numLevelSets`` , ``numLevelSets``)
95
96 """
97 if w0 == None and w1==None:
98 raise ValueError("Values for w0 or for w1 must be given.")
99 if wc == None and numLevelSets>1:
100 raise ValueError("Values for wc must be given.")
101
102 self.__domain=domain
103 DIM=self.__domain.getDim()
104 self.__numLevelSets=numLevelSets
105
106 self.__pde=LinearPDE(self.__domain, numEquations=self.__numLevelSets)
107 self.__pde.getSolverOptions().setTolerance(tol)
108 self.__pde.setSymmetryOn()
109 try:
110 self.__pde.setValue(q=location_of_set_m)
111 except IllegalCoefficientValue:
112 raise ValueError("Unable to set location of fixed level set function.")
113
114 # =========== check the shape of the scales: ========================
115 if scale is None:
116 if numLevelSets == 1 :
117 scale = 1.
118 else:
119 scale = np.ones((numLevelSets,))
120 else:
121 scale=np.asarray(scale)
122 if numLevelSets == 1 :
123 if scale.shape == ():
124 if not scale > 0 :
125 raise ValueError("Value for scale must be positive.")
126 else:
127 raise ValueError("Unexpected shape %s for scale."%scale.shape)
128 else:
129 if scale.shape is (numLevelSets,):
130 if not min(scale) > 0:
131 raise ValueError("All value for scale must be positive.")
132 else:
133 raise ValueError("Unexpected shape %s for scale."%scale.shape)
134
135 if scale_c is None or numLevelSets < 2:
136 scale_c = np.ones((numLevelSets,numLevelSets))
137 else:
138 scale_c=np.asarray(scale_c)
139 if scale_c.shape == (numLevelSets,numLevelSets):
140 if not all( [ [ scale_c[l,k] > 0. for l in range(k) ] for k in range(1,numLevelSets) ]):
141 raise ValueError("All values in the lower triangle of scale_c must be positive.")
142 else:
143 raise ValueError("Unexpected shape %s for scale."%scale_c.shape)
144 # ===== check the shape of the weights: =============================
145 if w0 is not None:
146 w0 = interpolate(w0,self.__pde.getFunctionSpaceForCoefficient('D'))
147 s0=w0.getShape()
148 if numLevelSets == 1 :
149 if not s0 == () :
150 raise ValueError("Unexpected shape %s for weight w0."%s0)
151 else:
152 if not s0 == (numLevelSets,):
153 raise ValueError("Unexpected shape %s for weight w0."%s0)
154
155 if not w1 is None:
156 w1 = interpolate(w1,self.__pde.getFunctionSpaceForCoefficient('A'))
157 s1=w1.getShape()
158 if numLevelSets is 1 :
159 if not s1 == (DIM,) :
160 raise ValueError("Unexpected shape %s for weight w1."%s1)
161 else:
162 if not s1 == (numLevelSets,DIM):
163 raise ValueError("Unexpected shape %s for weight w1."%s1)
164
165 if numLevelSets == 1:
166 wc=None
167 else:
168 wc = interpolate(wc,self.__pde.getFunctionSpaceForCoefficient('A'))
169 sc=wc.getShape()
170 if not sc == (numLevelSets, numLevelSets):
171 raise ValueError("Unexpected shape %s for weight wc."%(sc,))
172 # ============= now we rescale weights: =============================
173 L2s=np.asarray(boundingBoxEdgeLengths(domain))**2
174 L4=1/np.sum(1/L2s)**2
175 if numLevelSets == 1:
176 A=0
177 if w0 is not None:
178 A = integrate(w0)
179 if w1 is not None:
180 A += integrate(inner(w1, 1/L2s))
181 if A > 0:
182 f = scale/A
183 if w0 is not None:
184 w0*=f
185 if w1 is not None:
186 w1*=f
187 else:
188 raise ValueError("Non-positive weighting factor detected.")
189 else:
190
191 for k in range(numLevelSets):
192 A=0
193 if w0 is not None:
194 A = integrate(w0[k])
195 if w1 is not None:
196 A += integrate(inner(w1[k,:], 1/L2s))
197 if A > 0:
198 f = scale[k]/A
199 if w0 is not None:
200 w0[k]*=f
201 if w1 is not None:
202 w1[k,:]*=f
203 else:
204 raise ValueError("Non-positive weighting factor for level set component %d detected."%k)
205
206 # and now the cross-gradient:
207 if wc is not None:
208 for l in range(k):
209 A = integrate(wc[l,k])/L4
210 if A > 0:
211 f = scale_c[l,k]/A
212 wc[l,k]*=f
213 # else:
214 # raise ValueError("Non-positive weighting factor for cross-gradient level set components %d and %d detected."%(l,k))
215
216 self.__w0=w0
217 self.__w1=w1
218 self.__wc=wc
219
220 self.__pde_is_set=False
221 if self.__numLevelSets > 1:
222 self.__useDiagonalHessianApproximation=useDiagonalHessianApproximation
223 else:
224 self.__useDiagonalHessianApproximation=True
225 self._update_Hessian=True
226
227 self.__num_tradeoff_factors=numLevelSets+((numLevelSets-1)*numLevelSets)/2
228 self.setTradeOffFactors()
229 self.__vol_d=vol(self.__domain)
230
231 def getDomain(self):
232 """
233 returns the domain of the regularization term
234
235 :rtype: ``Domain``
236 """
237 return self.__domain
238
239 def getNumLevelSets(self):
240 """
241 returns the number of level set functions
242
243 :rtype: ``int``
244 """
245 return self.__numLevelSets
246
247 def getPDE(self):
248 """
249 returns the linear PDE to be solved for the Hessian Operator inverse
250
251 :rtype: `LinearPDE`
252 """
253 return self.__pde
254
255 def getDualProduct(self, m, r):
256 """
257 returns the dual product of a gradient represented by X=r[1] and Y=r[0]
258 with a level set function m:
259
260 *Y_i*m_i + X_ij*m_{i,j}*
261
262 :type m: `Data`
263 :type r: `ArithmeticTuple`
264 :rtype: ``float``
265 """
266 A=0
267 if not r[0].isEmpty(): A+=integrate(inner(r[0], m))
268 if not r[1].isEmpty(): A+=integrate(inner(r[1], grad(m)))
269 return A
270 def getNumTradeOffFactors(self):
271 """
272 returns the number of trade-off factors being used.
273
274 :rtype: ``int``
275 """
276 return self.__num_tradeoff_factors
277
278 def setTradeOffFactors(self, mu=None):
279 """
280 sets the trade-off factors for the level-set variation and the cross-gradient
281
282 :param mu: new values for the trade-off factors where values mu[:numLevelSets] are the
283 trade-off factors for the level-set variation and the remaining values for
284 the cross-gradient part with mu_c[l,k]=mu[numLevelSets+l+((k-1)*k)/2] (l<k).
285 If no values for mu is given ones are used. Values must be positive.
286 :type mu: ``list`` of ``float`` or ```numpy.array```
287 """
288 numLS=self.getNumLevelSets()
289 numTF=self.getNumTradeOffFactors()
290 if mu is None:
291 mu = np.ones((numTF,))
292 else:
293 mu = np.asarray(mu)
294
295 if mu.shape == (numTF,):
296 self.setTradeOffFactorsForVariation(mu[:numLS])
297 mu_c2=np.zeros((numLS,numLS))
298 for k in range(numLS):
299 for l in range(k):
300 mu_c2[l,k] = mu[numLS+l+((k-1)*k)/2]
301 self.setTradeOffFactorsForCrossGradient(mu_c2)
302 elif mu.shape == () and numLS ==1:
303 self.setTradeOffFactorsForVariation(mu)
304 else:
305 raise ValueError("Unexpected shape %s for mu."%(mu.shape,))
306
307 def setTradeOffFactorsForVariation(self, mu=None):
308 """
309 sets the trade-off factors for the level-set variation part
310
311 :param mu: new values for the trade-off factors. Values must be positive.
312 :type mu: ``float``, ``list`` of ``float`` or ```numpy.array```
313 """
314 numLS=self.getNumLevelSets()
315 if mu is None:
316 if numLS == 1:
317 mu = 1.
318 else:
319 mu = np.ones((numLS,))
320
321 mu=np.asarray(mu)
322 if numLS == 1:
323 if mu.shape == (1,): mu=mu[0]
324 if mu.shape == ():
325 if mu > 0:
326 self.__mu= mu
327 self._new_mu=True
328 else:
329 raise ValueError("Value for trade-off factor must be positive.")
330 else:
331 raise ValueError("Unexpected shape %s for mu."%mu.shape)
332 else:
333 if mu.shape == (numLS,):
334 if min(mu) > 0:
335 self.__mu= mu
336 self._new_mu=True
337 else:
338 raise ValueError("All value for mu must be positive.")
339 else:
340 raise ValueError("Unexpected shape %s for trade-off factor."%mu.shape)
341
342 def setTradeOffFactorsForCrossGradient(self, mu_c=None):
343 """
344 sets the trade-off factors for the cross-gradient terms
345
346 :param mu_c: new values for the trade-off factors for the cross-gradient
347 terms. Values must be positive. If no value is given ones
348 are used. Onky value mu_c[l,k] for l<k are used.
349 :type mu_c: ``float``, ``list`` of ``float`` or ``numpy.array``
350 """
351 numLS=self.getNumLevelSets()
352 if mu_c is None or numLS < 2:
353 self.__mu_c = np.ones((numLS,numLS))
354 if isinstance(mu_c, float) or isinstance(mu_c, int):
355 self.__mu_c = np.zeros((numLS,numLS))
356 self.__mu_c[:,:]=mu_c
357 else:
358 mu_c=np.asarray(mu_c)
359 if mu_c.shape == (numLS,numLS):
360 if not all( [ [ mu_c[l,k] > 0. for l in range(k) ] for k in range(1,numLS) ]):
361 raise ValueError("All trade-off factors in the lower triangle of mu_c must be positive.")
362 else:
363 self.__mu_c = mu_c
364 self._new_mu=True
365 else:
366 raise ValueError("Unexpected shape %s for mu."%mu_c.shape)
367
368 def getArguments(self, m):
369 """
370 """
371 return ( grad(m),)
372
373 def getValue(self, m, grad_m):
374 """
375 returns the value of the cost function J with respect to m.
376
377 :rtype: ``float``
378 """
379 mu=self.__mu
380 mu_c=self.__mu_c
381 DIM=self.getDomain().getDim()
382 numLS=self.getNumLevelSets()
383
384 A=0
385 if self.__w0 is not None:
386 A+=inner(integrate(m**2 * self.__w0), mu)
387
388 if self.__w1 is not None:
389 if numLS == 1:
390 A+=integrate(inner(grad_m**2, self.__w1))*mu
391 else:
392 for k in range(numLS):
393 A+=mu[k]*integrate(inner(grad_m[k,:]**2,self.__w1[k,:]))
394
395 if numLS > 1:
396 for k in range(numLS):
397 gk=grad_m[k,:]
398 len_gk=length(gk)
399 for l in range(k):
400 gl=grad_m[l,:]
401 A+= mu_c[l,k] * integrate( self.__wc[l,k] * ( len_gk * length(gl) )**2 - inner(gk, gl)**2 )
402 return A/2
403
404 def getGradient(self, m, grad_m):
405 """
406 returns the gradient of the cost function J with respect to m.
407 The function returns Y_k=dPsi/dm_k and X_kj=dPsi/dm_kj
408 """
409
410 mu=self.__mu
411 mu_c=self.__mu_c
412 DIM=self.getDomain().getDim()
413 numLS=self.getNumLevelSets()
414
415 grad_m=grad(m, Function(m.getDomain()))
416 if self.__w0 is not None:
417 Y = m * self.__w0 * mu
418 else:
419 if numLS == 1:
420 Y = Scalar(0, grad_m.getFunctionSpace())
421 else:
422 Y = Data(0, (numLS,) , grad_m.getFunctionSpace())
423
424 if self.__w1 is not None:
425 X=grad_m*self.__w1
426 if numLS == 1:
427 X=grad_m* self.__w1*mu
428 else:
429 for k in range(numLS):
430 X[k,:]*=mu[k]
431 else:
432 X = Data(0, grad_m.getShape(), grad_m.getFunctionSpace())
433
434 # cross gradient terms:
435 if numLS > 1:
436 for k in range(numLS):
437 grad_m_k=grad_m[k,:]
438 l2_grad_m_k = length(grad_m_k)**2
439 for l in range(k):
440 grad_m_l=grad_m[l,:]
441 l2_grad_m_l = length(grad_m_l)**2
442 grad_m_lk = inner(grad_m_l, grad_m_k)
443 f= mu_c[l,k]* self.__wc[l,k]
444 X[l,:] += f * ( l2_grad_m_l * grad_m_l - grad_m_lk * grad_m_k )
445 X[k,:] += f * ( l2_grad_m_k * grad_m_k - grad_m_lk * grad_m_l )
446
447 return ArithmeticTuple(Y, X)
448
449 def getInverseHessianApproximation(self, m, r, grad_m):
450 """
451 """
452 if self._new_mu or self._update_Hessian:
453
454 self._new_mu=False
455 self._update_Hessian=False
456 mu=self.__mu
457 mu_c=self.__mu_c
458
459 DIM=self.getDomain().getDim()
460 numLS=self.getNumLevelSets()
461 if self.__w0 is not None:
462 if numLS == 1:
463 D=self.__w0 * mu
464 else:
465 D=self.getPDE().createCoefficient("D")
466 for k in range(numLS): D[k,k]=self.__w0[k] * mu[k]
467 self.getPDE().setValue(D=D)
468
469 A=self.getPDE().createCoefficient("A")
470 if self.__w1 is not None:
471 if numLS == 1:
472 for i in range(DIM): A[i,i]=self.__w1[i] * mu
473 else:
474 for k in range(numLS):
475 for i in range(DIM): A[k,i,k,i]=self.__w1[k,i] * mu[k]
476
477 if numLS > 1:
478 for k in range(numLS):
479 grad_m_k=grad_m[k,:]
480 l2_grad_m_k = length(grad_m_k)**2
481 o_kk=outer(grad_m_k, grad_m_k)
482 for l in range(k):
483 grad_m_l=grad_m[l,:]
484 l2_grad_m_l = length(grad_m_l)**2
485 i_lk = inner(grad_m_l, grad_m_k)
486 o_lk = outer(grad_m_l, grad_m_k)
487 o_kl = outer(grad_m_k, grad_m_l)
488 o_ll=outer(grad_m_l, grad_m_l)
489 f= mu_c[l,k]* self.__wc[l,k]
490
491 A[l,:,l,:] += f * ( l2_grad_m_k * kronecker(DIM) - o_kk )
492 A[l,:,k,:] += f * ( 2 * o_lk - o_kl - i_lk * kronecker(DIM) )
493 A[k,:,l,:] += f * ( 2 * o_kl - o_lk - i_lk * kronecker(DIM) )
494 A[k,:,k,:] += f * ( l2_grad_m_l * kronecker(DIM) - o_ll )
495 self.getPDE().setValue(A=A)
496 #self.getPDE().resetRightHandSideCoefficients()
497 #self.getPDE().setValue(X=r[1])
498 #print "X only: ",self.getPDE().getSolution()
499 #self.getPDE().resetRightHandSideCoefficients()
500 #self.getPDE().setValue(Y=r[0])
501 #print "Y only: ",self.getPDE().getSolution()
502
503 self.getPDE().resetRightHandSideCoefficients()
504 self.getPDE().setValue(X=r[1], Y=r[0])
505 return self.getPDE().getSolution()
506
507 def updateHessian(self):
508 """
509 notify the class to recalculate the Hessian operator
510 """
511 if not self.__useDiagonalHessianApproximation:
512 self._update_Hessian=True
513
514 def getNorm(self, m):
515 """
516 returns the norm of ``m``
517
518 :param m: level set function
519 :type m: `Data`
520 :rtype: ``float``
521 """
522 return sqrt(integrate(length(m)**2)/self.__vol_d)
523

  ViewVC Help
Powered by ViewVC 1.1.26