/[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 5148 - (show annotations)
Mon Sep 15 01:25:23 2014 UTC (4 years, 11 months ago) by caltinay
File MIME type: text/x-python
File size: 22832 byte(s)
Merging ripley diagonal storage + CUDA support into trunk.
Options file version has been incremented due to new options
'cuda' and 'nvccflags'.

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

  ViewVC Help
Powered by ViewVC 1.1.26