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

  ViewVC Help
Powered by ViewVC 1.1.26