/[escript]/trunk/downunder/py_src/regularizations.py
ViewVC logotype

Annotation of /trunk/downunder/py_src/regularizations.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6787 - (hide annotations)
Mon Feb 4 01:41:55 2019 UTC (9 months, 2 weeks ago) by aellery
File MIME type: text/x-python
File size: 23121 byte(s)
eScripts now compiles successfully with the anaconda python interpreter.


1 caltinay 4394
2 jfenwick 3981 ##############################################################################
3 caltinay 3946 #
4 jfenwick 6651 # Copyright (c) 2003-2018 by The University of Queensland
5 jfenwick 3981 # http://www.uq.edu.au
6 caltinay 3946 #
7     # Primary Business: Queensland, Australia
8 jfenwick 6112 # Licensed under the Apache License, version 2.0
9     # http://www.apache.org/licenses/LICENSE-2.0
10 caltinay 3946 #
11 jfenwick 3981 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 jfenwick 4657 # Development 2012-2013 by School of Earth Sciences
13     # Development from 2014 by Centre for Geoscience Computing (GeoComp)
14 jfenwick 3981 #
15     ##############################################################################
16 caltinay 3946
17 sshaw 5706 from __future__ import print_function, division
18    
19 jfenwick 6651 __copyright__="""Copyright (c) 2003-2018 by The University of Queensland
20 jfenwick 3981 http://www.uq.edu.au
21 caltinay 3946 Primary Business: Queensland, Australia"""
22 jfenwick 6112 __license__="""Licensed under the Apache License, version 2.0
23     http://www.apache.org/licenses/LICENSE-2.0"""
24 caltinay 3946 __url__="https://launchpad.net/escript-finley"
25    
26 caltinay 3947 __all__ = ['Regularization']
27    
28 jfenwick 4387
29 caltinay 5112 import logging
30 caltinay 3946 import numpy as np
31 gross 5005 from esys.escript import Function, outer, Data, Scalar, grad, inner, integrate, interpolate, kronecker, boundingBoxEdgeLengths, vol, sqrt, length,Lsup, transpose
32 caltinay 5148 from esys.escript.linearPDEs import LinearPDE, IllegalCoefficientValue,SolverOptions
33 gross 4074 from esys.escript.pdetools import ArithmeticTuple
34 sshaw 5780 from .coordinates import makeTransformation
35 caltinay 4417 from .costfunctions import CostFunction
36 caltinay 3946
37 gross 4074 class Regularization(CostFunction):
38 caltinay 3946 """
39 caltinay 4097 The regularization term for the level set function ``m`` within the cost
40     function J for an inversion:
41    
42 caltinay 4285 *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 caltinay 4097
44 gross 4099 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 caltinay 4097 integrals over the domain are constant:
48    
49 gross 4099 *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 caltinay 4097
52 caltinay 3946 """
53 caltinay 4097 def __init__(self, domain, numLevelSets=1,
54 gross 4099 w0=None, w1=None, wc=None,
55 caltinay 4097 location_of_set_m=Data(),
56 gross 4099 useDiagonalHessianApproximation=False, tol=1e-8,
57 gross 4373 coordinates=None,
58 gross 4099 scale=None, scale_c=None):
59 caltinay 4007 """
60 caltinay 4095 initialization.
61 caltinay 4097
62     :param domain: domain
63 gross 4074 :type domain: `Domain`
64     :param numLevelSets: number of level sets
65     :type numLevelSets: ``int``
66 gross 4099 :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 caltinay 4097 (``numLevelSets`` ,) if ``numLevelSets`` > 1
69 gross 4099 :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 caltinay 4097 (``numLevelSets`` , DIM) if ``numLevelSets`` > 1
72 gross 4099 :param wc: weighting factor for the cross gradient terms. If not set
73 caltinay 4097 zero is assumed. Used for the case if ``numLevelSets`` > 1
74 gross 4099 only. Only values ``wc[l,k]`` in the lower triangle (l<k)
75 caltinay 4097 are used.
76 gross 4099 :type wc: `Data` object of shape (``numLevelSets`` , ``numLevelSets``)
77 caltinay 4097 :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 gross 4099 :param useDiagonalHessianApproximation: if True cross gradient terms
82 caltinay 4097 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 gross 4074 :type useDiagonalHessianApproximation: ``bool``
88 caltinay 4097 :param tol: tolerance when solving the PDE for the inverse of the
89     Hessian Operator
90 gross 4074 :type tol: positive ``float``
91 caltinay 4213
92 gross 4373 :param coordinates: defines coordinate system to be used
93     :type coordinates: ReferenceSystem` or `SpatialCoordinateTransformation`
94 caltinay 4417 :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 gross 4099 :param scale_c: scale for the cross gradient terms. If not set
99     one is assumed. Used for the case if ``numLevelSets`` > 1
100 caltinay 4417 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 caltinay 4095
104 caltinay 4007 """
105 jfenwick 5024 if w0 is None and w1 is None:
106 caltinay 4417 raise ValueError("Values for w0 or for w1 must be given.")
107 jfenwick 5024 if wc is None and numLevelSets>1:
108 caltinay 4417 raise ValueError("Values for wc must be given.")
109 caltinay 4097
110 caltinay 5112 self.logger = logging.getLogger('inv.%s'%self.__class__.__name__)
111 caltinay 3946 self.__domain=domain
112 gross 4074 DIM=self.__domain.getDim()
113     self.__numLevelSets=numLevelSets
114 sshaw 5780 self.__trafo=makeTransformation(domain, coordinates)
115 gross 4729 self.__pde=LinearPDE(self.__domain, numEquations=self.__numLevelSets, numSolutions=self.__numLevelSets)
116 gross 4074 self.__pde.getSolverOptions().setTolerance(tol)
117     self.__pde.setSymmetryOn()
118 gross 5005 self.__pde.setValue(A=self.__pde.createCoefficient('A'), D=self.__pde.createCoefficient('D'), )
119 gross 4074 try:
120 caltinay 4097 self.__pde.setValue(q=location_of_set_m)
121     except IllegalCoefficientValue:
122     raise ValueError("Unable to set location of fixed level set function.")
123 caltinay 4213
124     # =========== check the shape of the scales: ========================
125 gross 4099 if scale is None:
126 caltinay 4132 if numLevelSets == 1 :
127 caltinay 4417 scale = 1.
128 gross 4074 else:
129 caltinay 4417 scale = np.ones((numLevelSets,))
130 caltinay 4132 else:
131     scale=np.asarray(scale)
132 caltinay 4417 if numLevelSets == 1:
133 caltinay 4132 if scale.shape == ():
134 caltinay 4417 if not scale > 0 :
135     raise ValueError("Value for scale must be positive.")
136 caltinay 4132 else:
137 caltinay 4417 raise ValueError("Unexpected shape %s for scale."%scale.shape)
138 caltinay 4132 else:
139 caltinay 4417 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 caltinay 4213
145 gross 4099 if scale_c is None or numLevelSets < 2:
146 caltinay 4132 scale_c = np.ones((numLevelSets,numLevelSets))
147     else:
148     scale_c=np.asarray(scale_c)
149     if scale_c.shape == (numLevelSets,numLevelSets):
150 caltinay 4263 if not all( [ [ scale_c[l,k] > 0. for l in range(k) ] for k in range(1,numLevelSets) ]):
151 caltinay 4417 raise ValueError("All values in the lower triangle of scale_c must be positive.")
152 gross 4099 else:
153 caltinay 4417 raise ValueError("Unexpected shape %s for scale."%scale_c.shape)
154 caltinay 4213 # ===== check the shape of the weights: =============================
155 gross 4099 if w0 is not None:
156 caltinay 4417 w0 = interpolate(w0,self.__pde.getFunctionSpaceForCoefficient('D'))
157     s0=w0.getShape()
158     if numLevelSets == 1:
159     if not s0 == () :
160 gross 4729 raise ValueError("Unexpected shape %s for weight w0."%(s0,))
161 caltinay 4417 else:
162     if not s0 == (numLevelSets,):
163 gross 4729 raise ValueError("Unexpected shape %s for weight w0."%(s0,))
164 caltinay 4417 if not self.__trafo.isCartesian():
165     w0*=self.__trafo.getVolumeFactor()
166 caltinay 4132 if not w1 is None:
167 caltinay 4417 w1 = interpolate(w1,self.__pde.getFunctionSpaceForCoefficient('A'))
168     s1=w1.getShape()
169     if numLevelSets == 1 :
170     if not s1 == (DIM,) :
171 gross 4729 raise ValueError("Unexpected shape %s for weight w1."%(s1,))
172 caltinay 4417 else:
173     if not s1 == (numLevelSets,DIM):
174 gross 4729 raise ValueError("Unexpected shape %s for weight w1."%(s1,))
175 caltinay 4417 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 caltinay 4213 if numLevelSets == 1:
183 caltinay 4417 wc=None
184 gross 4099 else:
185 caltinay 4417 wc = interpolate(wc,self.__pde.getFunctionSpaceForCoefficient('A'))
186     sc=wc.getShape()
187     if not sc == (numLevelSets, numLevelSets):
188 gross 4124 raise ValueError("Unexpected shape %s for weight wc."%(sc,))
189 caltinay 4417 if not self.__trafo.isCartesian():
190     raise ValueError("Non-cartesian coordinates for cross-gradient term is not supported yet.")
191 caltinay 4213 # ============= now we rescale weights: =============================
192 gross 4099 L2s=np.asarray(boundingBoxEdgeLengths(domain))**2
193     L4=1/np.sum(1/L2s)**2
194 caltinay 4213 if numLevelSets == 1:
195 gross 4099 A=0
196     if w0 is not None:
197 caltinay 4132 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 caltinay 4417 w0*=f
204 caltinay 4132 if w1 is not None:
205 caltinay 4417 w1*=f
206 gross 4099 else:
207 caltinay 4417 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 gross 4124
224 caltinay 4417 # and now the cross-gradient:
225     if wc is not None:
226     for l in range(k):
227 caltinay 4132 A = integrate(wc[l,k])/L4
228     if A > 0:
229 caltinay 4417 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 caltinay 4213
234 gross 4099 self.__w0=w0
235     self.__w1=w1
236     self.__wc=wc
237 caltinay 4097
238 caltinay 4213 self.__pde_is_set=False
239 gross 4099 if self.__numLevelSets > 1:
240     self.__useDiagonalHessianApproximation=useDiagonalHessianApproximation
241     else:
242     self.__useDiagonalHessianApproximation=True
243     self._update_Hessian=True
244 caltinay 4097
245 jfenwick 5411 self.__num_tradeoff_factors=numLevelSets+((numLevelSets-1)*numLevelSets)//2
246 gross 4099 self.setTradeOffFactors()
247 gross 4102 self.__vol_d=vol(self.__domain)
248 caltinay 4213
249 gross 4074 def getDomain(self):
250 caltinay 3946 """
251 caltinay 4097 returns the domain of the regularization term
252    
253 gross 4074 :rtype: ``Domain``
254 caltinay 3946 """
255 gross 4074 return self.__domain
256 caltinay 4417
257 gross 4373 def getCoordinateTransformation(self):
258     """
259     returns the coordinate transformation being used
260 caltinay 4097
261 caltinay 4417 :rtype: `CoordinateTransformation`
262 gross 4373 """
263     return self.__trafo
264 caltinay 4417
265 gross 4074 def getNumLevelSets(self):
266     """
267 caltinay 4097 returns the number of level set functions
268    
269 gross 4074 :rtype: ``int``
270     """
271     return self.__numLevelSets
272 caltinay 3946
273 gross 4074 def getPDE(self):
274 caltinay 4007 """
275 gross 4074 returns the linear PDE to be solved for the Hessian Operator inverse
276 caltinay 4097
277     :rtype: `LinearPDE`
278 caltinay 4007 """
279 gross 4074 return self.__pde
280 caltinay 4213
281 gross 4074 def getDualProduct(self, m, r):
282     """
283 caltinay 4097 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 gross 4074 *Y_i*m_i + X_ij*m_{i,j}*
287 caltinay 4097
288     :type m: `Data`
289     :type r: `ArithmeticTuple`
290     :rtype: ``float``
291 gross 4074 """
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 caltinay 4095 return A
296 caltinay 4417
297 gross 4099 def getNumTradeOffFactors(self):
298     """
299     returns the number of trade-off factors being used.
300 caltinay 3946
301 gross 4099 :rtype: ``int``
302     """
303     return self.__num_tradeoff_factors
304    
305     def setTradeOffFactors(self, mu=None):
306     """
307 caltinay 4417 sets the trade-off factors for the level-set variation and the
308     cross-gradient.
309 caltinay 4213
310 caltinay 4417 :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 gross 4099 :type mu: ``list`` of ``float`` or ```numpy.array```
318     """
319     numLS=self.getNumLevelSets()
320     numTF=self.getNumTradeOffFactors()
321     if mu is None:
322 caltinay 4417 mu = np.ones((numTF,))
323 caltinay 4132 else:
324 caltinay 4417 mu = np.asarray(mu)
325 gross 4099
326 caltinay 4132 if mu.shape == (numTF,):
327     self.setTradeOffFactorsForVariation(mu[:numLS])
328     mu_c2=np.zeros((numLS,numLS))
329 caltinay 4263 for k in range(numLS):
330 caltinay 4417 for l in range(k):
331 caltinay 5073 mu_c2[l,k] = mu[numLS+l+((k-1)*k)//2]
332 caltinay 4132 self.setTradeOffFactorsForCrossGradient(mu_c2)
333     elif mu.shape == () and numLS ==1:
334     self.setTradeOffFactorsForVariation(mu)
335     else:
336 caltinay 4417 raise ValueError("Unexpected shape %s for mu."%(mu.shape,))
337 caltinay 4213
338 gross 4099 def setTradeOffFactorsForVariation(self, mu=None):
339 caltinay 4417 """
340     sets the trade-off factors for the level-set variation part.
341 caltinay 4213
342 caltinay 4417 :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 caltinay 4132 if numLS == 1:
348 caltinay 4417 mu = 1.
349 caltinay 4132 else:
350 caltinay 4417 mu = np.ones((numLS,))
351 jduplessis 4931 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 caltinay 4417 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 sshaw 4643 raise ValueError("Unexpected shape %s for mu."%str(mu.shape))
367 caltinay 4417 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 sshaw 4643 raise ValueError("Unexpected shape %s for trade-off factor."%str(mu.shape))
376 gross 4099
377     def setTradeOffFactorsForCrossGradient(self, mu_c=None):
378     """
379 caltinay 4417 sets the trade-off factors for the cross-gradient terms.
380 caltinay 4213
381 caltinay 4132 :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 caltinay 4417 are used. Only value mu_c[l,k] for l<k are used.
384 caltinay 4132 :type mu_c: ``float``, ``list`` of ``float`` or ``numpy.array``
385 gross 4099 """
386     numLS=self.getNumLevelSets()
387     if mu_c is None or numLS < 2:
388 caltinay 4132 self.__mu_c = np.ones((numLS,numLS))
389 gross 4532 elif isinstance(mu_c, float) or isinstance(mu_c, int):
390 caltinay 4132 self.__mu_c = np.zeros((numLS,numLS))
391 caltinay 4213 self.__mu_c[:,:]=mu_c
392 caltinay 4132 else:
393     mu_c=np.asarray(mu_c)
394     if mu_c.shape == (numLS,numLS):
395 caltinay 4263 if not all( [ [ mu_c[l,k] > 0. for l in range(k) ] for k in range(1,numLS) ]):
396 caltinay 4417 raise ValueError("All trade-off factors in the lower triangle of mu_c must be positive.")
397 caltinay 4132 else:
398 caltinay 4417 self.__mu_c = mu_c
399     self._new_mu=True
400 gross 4099 else:
401 gross 4532 raise ValueError("Unexpected shape %s for mu."%(mu_c.shape,))
402 caltinay 4213
403 gross 4099 def getArguments(self, m):
404     """
405     """
406 caltinay 4417 return grad(m),
407 gross 4099
408 gross 4074 def getValue(self, m, grad_m):
409 caltinay 4007 """
410 caltinay 4097 returns the value of the cost function J with respect to m.
411 jfenwick 4416 This equation is specified in the inversion cookbook.
412 caltinay 4097
413 gross 4074 :rtype: ``float``
414 caltinay 4007 """
415 gross 4099 mu=self.__mu
416     mu_c=self.__mu_c
417 gross 4074 DIM=self.getDomain().getDim()
418     numLS=self.getNumLevelSets()
419 caltinay 4097
420 caltinay 3946 A=0
421 gross 4099 if self.__w0 is not None:
422 caltinay 5112 r = inner(integrate(m**2 * self.__w0), mu)
423     self.logger.debug("J_R[m^2] = %e"%r)
424     A += r
425 caltinay 4097
426 gross 4099 if self.__w1 is not None:
427 caltinay 4095 if numLS == 1:
428 caltinay 5112 r = integrate(inner(grad_m**2, self.__w1))*mu
429     self.logger.debug("J_R[grad(m)] = %e"%r)
430     A += r
431 caltinay 4095 else:
432 caltinay 4263 for k in range(numLS):
433 caltinay 5112 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 caltinay 4213
437 gross 4099 if numLS > 1:
438 caltinay 4263 for k in range(numLS):
439 caltinay 4097 gk=grad_m[k,:]
440     len_gk=length(gk)
441 caltinay 4263 for l in range(k):
442 caltinay 4097 gl=grad_m[l,:]
443 caltinay 5112 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 gross 4074 return A/2
447 caltinay 4097
448 gross 4074 def getGradient(self, m, grad_m):
449     """
450 caltinay 4097 returns the gradient of the cost function J with respect to m.
451 caltinay 4417
452     :note: This implementation returns Y_k=dPsi/dm_k and X_kj=dPsi/dm_kj
453 gross 4074 """
454 caltinay 4097
455 gross 4099 mu=self.__mu
456     mu_c=self.__mu_c
457 gross 4074 DIM=self.getDomain().getDim()
458     numLS=self.getNumLevelSets()
459 caltinay 4097
460 caltinay 4211 grad_m=grad(m, Function(m.getDomain()))
461 gross 4099 if self.__w0 is not None:
462     Y = m * self.__w0 * mu
463 caltinay 4095 else:
464 caltinay 4132 if numLS == 1:
465 caltinay 4417 Y = Scalar(0, grad_m.getFunctionSpace())
466 caltinay 4132 else:
467 caltinay 4417 Y = Data(0, (numLS,) , grad_m.getFunctionSpace())
468 caltinay 3946
469 gross 4099 if self.__w1 is not None:
470 gross 5004
471 caltinay 4095 if numLS == 1:
472 gross 4099 X=grad_m* self.__w1*mu
473 caltinay 4095 else:
474 gross 5004 X=grad_m*self.__w1
475 caltinay 4263 for k in range(numLS):
476 gross 4099 X[k,:]*=mu[k]
477 caltinay 4095 else:
478 gross 4122 X = Data(0, grad_m.getShape(), grad_m.getFunctionSpace())
479 caltinay 4213
480 gross 4122 # cross gradient terms:
481 caltinay 4213 if numLS > 1:
482 caltinay 4417 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 gross 5004 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 caltinay 4213
493 gross 4074 return ArithmeticTuple(Y, X)
494 caltinay 4097
495 sshaw 5017 def getInverseHessianApproximation(self, m, r, grad_m, solve=True):
496 gross 4074 """
497     """
498     if self._new_mu or self._update_Hessian:
499 caltinay 4095 self._new_mu=False
500     self._update_Hessian=False
501 gross 4099 mu=self.__mu
502     mu_c=self.__mu_c
503 caltinay 4097
504 gross 4074 DIM=self.getDomain().getDim()
505     numLS=self.getNumLevelSets()
506 gross 4099 if self.__w0 is not None:
507 caltinay 4095 if numLS == 1:
508 caltinay 4417 D=self.__w0 * mu
509 caltinay 4097 else:
510 gross 5005 D=self.getPDE().getCoefficient("D")
511     D.setToZero()
512 caltinay 4417 for k in range(numLS): D[k,k]=self.__w0[k] * mu[k]
513 caltinay 4095 self.getPDE().setValue(D=D)
514 caltinay 3946
515 gross 5005 A=self.getPDE().getCoefficient("A")
516     A.setToZero()
517 gross 4099 if self.__w1 is not None:
518 caltinay 4417 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 caltinay 4263 for i in range(DIM): A[k,i,k,i]=self.__w1[k,i] * mu[k]
523 caltinay 4213
524 gross 4099 if numLS > 1:
525 gross 5005 # this could be make faster by creating caches for grad_m_k, l2_grad_m_k and o_kk
526 caltinay 4417 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 gross 5005 Z=f * (2*o_lk - o_kl - i_lk*kronecker(DIM))
539 caltinay 4417 A[l,:,l,:] += f * (l2_grad_m_k*kronecker(DIM) - o_kk)
540 gross 5005 A[l,:,k,:] += Z
541     A[k,:,l,:] += transpose(Z)
542 caltinay 4417 A[k,:,k,:] += f * (l2_grad_m_l*kronecker(DIM) - o_ll)
543 gross 4074 self.getPDE().setValue(A=A)
544 gross 4102 #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 caltinay 4097
551 gross 4079 self.getPDE().resetRightHandSideCoefficients()
552 gross 4074 self.getPDE().setValue(X=r[1], Y=r[0])
553 sshaw 5017 if not solve:
554     return self.getPDE()
555 aellery 6784 oldsettings = self.getPDE().getSolverOptions().getSolverMethod()
556 aellery 6787 self.getPDE().getSolverOptions().setSolverMethod(SolverOptions.GMRES)
557 aellery 6784 solution = self.getPDE().getSolution()
558     self.getPDE().getSolverOptions().setSolverMethod(oldsettings)
559     return solution
560 caltinay 4097
561 gross 4074 def updateHessian(self):
562     """
563 caltinay 4417 notifies the class to recalculate the Hessian operator.
564 gross 4074 """
565 caltinay 4095 if not self.__useDiagonalHessianApproximation:
566     self._update_Hessian=True
567 caltinay 4213
568 gross 4100 def getNorm(self, m):
569     """
570 caltinay 4417 returns the norm of ``m``.
571 gross 4100
572     :param m: level set function
573     :type m: `Data`
574     :rtype: ``float``
575     """
576 caltinay 4132 return sqrt(integrate(length(m)**2)/self.__vol_d)
577 caltinay 4213

  ViewVC Help
Powered by ViewVC 1.1.26