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

  ViewVC Help
Powered by ViewVC 1.1.26