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

Contents of /trunk/downunder/py_src/splitregularizations.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: 25362 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__ = ['SplitRegularization']
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 SplitRegularization(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.__pre_input = None
111 self.__pre_args = None
112 self.logger = logging.getLogger('inv.%s'%self.__class__.__name__)
113 self.__domain=domain
114 DIM=self.__domain.getDim()
115 self.__numLevelSets=numLevelSets
116 self.__trafo=makeTranformation(domain, coordinates)
117 self.__pde=LinearPDE(self.__domain, numEquations=self.__numLevelSets, numSolutions=self.__numLevelSets)
118 self.__pde.getSolverOptions().setTolerance(tol)
119 self.__pde.setSymmetryOn()
120 self.__pde.setValue(A=self.__pde.createCoefficient('A'), D=self.__pde.createCoefficient('D'), )
121 try:
122 self.__pde.setValue(q=location_of_set_m)
123 except IllegalCoefficientValue:
124 raise ValueError("Unable to set location of fixed level set function.")
125
126 # =========== check the shape of the scales: ========================
127 if scale is None:
128 if numLevelSets == 1 :
129 scale = 1.
130 else:
131 scale = np.ones((numLevelSets,))
132 else:
133 scale=np.asarray(scale)
134 if numLevelSets == 1:
135 if scale.shape == ():
136 if not scale > 0 :
137 raise ValueError("Value for scale must be positive.")
138 else:
139 raise ValueError("Unexpected shape %s for scale."%scale.shape)
140 else:
141 if scale.shape is (numLevelSets,):
142 if not min(scale) > 0:
143 raise ValueError("All values for scale must be positive.")
144 else:
145 raise ValueError("Unexpected shape %s for scale."%scale.shape)
146
147 if scale_c is None or numLevelSets < 2:
148 scale_c = np.ones((numLevelSets,numLevelSets))
149 else:
150 scale_c=np.asarray(scale_c)
151 if scale_c.shape == (numLevelSets,numLevelSets):
152 if not all( [ [ scale_c[l,k] > 0. for l in range(k) ] for k in range(1,numLevelSets) ]):
153 raise ValueError("All values in the lower triangle of scale_c must be positive.")
154 else:
155 raise ValueError("Unexpected shape %s for scale."%scale_c.shape)
156 # ===== check the shape of the weights: =============================
157 if w0 is not None:
158 w0 = interpolate(w0,self.__pde.getFunctionSpaceForCoefficient('D'))
159 s0=w0.getShape()
160 if numLevelSets == 1:
161 if not s0 == () :
162 raise ValueError("Unexpected shape %s for weight w0."%(s0,))
163 else:
164 if not s0 == (numLevelSets,):
165 raise ValueError("Unexpected shape %s for weight w0."%(s0,))
166 if not self.__trafo.isCartesian():
167 w0*=self.__trafo.getVolumeFactor()
168 if not w1 is None:
169 w1 = interpolate(w1,self.__pde.getFunctionSpaceForCoefficient('A'))
170 s1=w1.getShape()
171 if numLevelSets == 1 :
172 if not s1 == (DIM,) :
173 raise ValueError("Unexpected shape %s for weight w1."%(s1,))
174 else:
175 if not s1 == (numLevelSets,DIM):
176 raise ValueError("Unexpected shape %s for weight w1."%(s1,))
177 if not self.__trafo.isCartesian():
178 f=self.__trafo.getScalingFactors()**2*self.__trafo.getVolumeFactor()
179 if numLevelSets == 1:
180 w1*=f
181 else:
182 for i in range(numLevelSets): w1[i,:]*=f
183
184 if numLevelSets == 1:
185 wc=None
186 else:
187 wc = interpolate(wc,self.__pde.getFunctionSpaceForCoefficient('A'))
188 sc=wc.getShape()
189 if not sc == (numLevelSets, numLevelSets):
190 raise ValueError("Unexpected shape %s for weight wc."%(sc,))
191 if not self.__trafo.isCartesian():
192 raise ValueError("Non-cartesian coordinates for cross-gradient term is not supported yet.")
193 # ============= now we rescale weights: =============================
194 L2s=np.asarray(boundingBoxEdgeLengths(domain))**2
195 L4=1/np.sum(1/L2s)**2
196 if numLevelSets == 1:
197 A=0
198 if w0 is not None:
199 A = integrate(w0)
200 if w1 is not None:
201 A += integrate(inner(w1, 1/L2s))
202 if A > 0:
203 f = scale/A
204 if w0 is not None:
205 w0*=f
206 if w1 is not None:
207 w1*=f
208 else:
209 raise ValueError("Non-positive weighting factor detected.")
210 else: # numLevelSets > 1
211 for k in range(numLevelSets):
212 A=0
213 if w0 is not None:
214 A = integrate(w0[k])
215 if w1 is not None:
216 A += integrate(inner(w1[k,:], 1/L2s))
217 if A > 0:
218 f = scale[k]/A
219 if w0 is not None:
220 w0[k]*=f
221 if w1 is not None:
222 w1[k,:]*=f
223 else:
224 raise ValueError("Non-positive weighting factor for level set component %d detected."%k)
225
226 # and now the cross-gradient:
227 if wc is not None:
228 for l in range(k):
229 A = integrate(wc[l,k])/L4
230 if A > 0:
231 f = scale_c[l,k]/A
232 wc[l,k]*=f
233 # else:
234 # raise ValueError("Non-positive weighting factor for cross-gradient level set components %d and %d detected."%(l,k))
235
236 self.__w0=w0
237 self.__w1=w1
238 self.__wc=wc
239
240 self.__pde_is_set=False
241 if self.__numLevelSets > 1:
242 self.__useDiagonalHessianApproximation=useDiagonalHessianApproximation
243 else:
244 self.__useDiagonalHessianApproximation=True
245 self._update_Hessian=True
246
247 self.__num_tradeoff_factors=numLevelSets+((numLevelSets-1)*numLevelSets)//2
248 self.setTradeOffFactors()
249 self.__vol_d=vol(self.__domain)
250
251 def getDomain(self):
252 """
253 returns the domain of the regularization term
254
255 :rtype: ``Domain``
256 """
257 return self.__domain
258
259 def getCoordinateTransformation(self):
260 """
261 returns the coordinate transformation being used
262
263 :rtype: `CoordinateTransformation`
264 """
265 return self.__trafo
266
267 def getNumLevelSets(self):
268 """
269 returns the number of level set functions
270
271 :rtype: ``int``
272 """
273 return self.__numLevelSets
274
275 def getPDE(self):
276 """
277 returns the linear PDE to be solved for the Hessian Operator inverse
278
279 :rtype: `LinearPDE`
280 """
281 return self.__pde
282
283 def getDualProduct(self, m, r):
284 """
285 returns the dual product of a gradient represented by X=r[1] and Y=r[0]
286 with a level set function m:
287
288 *Y_i*m_i + X_ij*m_{i,j}*
289
290 :type m: `Data`
291 :type r: `ArithmeticTuple`
292 :rtype: ``float``
293 """
294 A=0
295 if not r[0].isEmpty(): A+=integrate(inner(r[0], m))
296 if not r[1].isEmpty(): A+=integrate(inner(r[1], grad(m)))
297 return A
298
299 def getNumTradeOffFactors(self):
300 """
301 returns the number of trade-off factors being used.
302
303 :rtype: ``int``
304 """
305 return self.__num_tradeoff_factors
306
307 def setTradeOffFactors(self, mu=None):
308 """
309 sets the trade-off factors for the level-set variation and the
310 cross-gradient.
311
312 :param mu: new values for the trade-off factors where values
313 mu[:numLevelSets] are the trade-off factors for the
314 level-set variation and the remaining values for
315 the cross-gradient part with
316 mu_c[l,k]=mu[numLevelSets+l+((k-1)*k)/2] (l<k).
317 If no values for mu are given ones are used.
318 Values must be positive.
319 :type mu: ``list`` of ``float`` or ```numpy.array```
320 """
321 numLS=self.getNumLevelSets()
322 numTF=self.getNumTradeOffFactors()
323 if mu is None:
324 mu = np.ones((numTF,))
325 else:
326 mu = np.asarray(mu)
327
328 if mu.shape == (numTF,):
329 self.setTradeOffFactorsForVariation(mu[:numLS])
330 mu_c2=np.zeros((numLS,numLS))
331 for k in range(numLS):
332 for l in range(k):
333 mu_c2[l,k] = mu[numLS+l+((k-1)*k)//2]
334 self.setTradeOffFactorsForCrossGradient(mu_c2)
335 elif mu.shape == () and numLS ==1:
336 self.setTradeOffFactorsForVariation(mu)
337 else:
338 raise ValueError("Unexpected shape %s for mu."%(mu.shape,))
339
340 def setTradeOffFactorsForVariation(self, mu=None):
341 """
342 sets the trade-off factors for the level-set variation part.
343
344 :param mu: new values for the trade-off factors. Values must be positive.
345 :type mu: ``float``, ``list`` of ``float`` or ```numpy.array```
346 """
347 numLS=self.getNumLevelSets()
348 if mu is None:
349 if numLS == 1:
350 mu = 1.
351 else:
352 mu = np.ones((numLS,))
353 if type(mu) == list:
354 #this is a fix for older versions of numpy where passing in an a list of ints causes
355 #this code to break.
356 mu=np.asarray([float(i) for i in mu])
357 else:
358 mu=np.asarray(mu)
359 if numLS == 1:
360 if mu.shape == (1,): mu=mu[0]
361 if mu.shape == ():
362 if mu > 0:
363 self.__mu= mu
364 self._new_mu=True
365 else:
366 raise ValueError("Value for trade-off factor must be positive.")
367 else:
368 raise ValueError("Unexpected shape %s for mu."%str(mu.shape))
369 else:
370 if mu.shape == (numLS,):
371 if min(mu) > 0:
372 self.__mu= mu
373 self._new_mu=True
374 else:
375 raise ValueError("All values for mu must be positive.")
376 else:
377 raise ValueError("Unexpected shape %s for trade-off factor."%str(mu.shape))
378
379 def setTradeOffFactorsForCrossGradient(self, mu_c=None):
380 """
381 sets the trade-off factors for the cross-gradient terms.
382
383 :param mu_c: new values for the trade-off factors for the cross-gradient
384 terms. Values must be positive. If no value is given ones
385 are used. Only value mu_c[l,k] for l<k are used.
386 :type mu_c: ``float``, ``list`` of ``float`` or ``numpy.array``
387 """
388 numLS=self.getNumLevelSets()
389 if mu_c is None or numLS < 2:
390 self.__mu_c = np.ones((numLS,numLS))
391 elif isinstance(mu_c, float) or isinstance(mu_c, int):
392 self.__mu_c = np.zeros((numLS,numLS))
393 self.__mu_c[:,:]=mu_c
394 else:
395 mu_c=np.asarray(mu_c)
396 if mu_c.shape == (numLS,numLS):
397 if not all( [ [ mu_c[l,k] > 0. for l in range(k) ] for k in range(1,numLS) ]):
398 raise ValueError("All trade-off factors in the lower triangle of mu_c must be positive.")
399 else:
400 self.__mu_c = mu_c
401 self._new_mu=True
402 else:
403 raise ValueError("Unexpected shape %s for mu."%(mu_c.shape,))
404
405 def getArguments(self, m):
406 """
407 """
408 raise RuntimeError("Please use the setPoint interface")
409 self.__pre_args = grad(m)
410 self.__pre_input = m
411 return self.__pre_args,
412
413
414 def getValueAtPoint(self):
415 """
416 returns the value of the cost function J with respect to m.
417 This equation is specified in the inversion cookbook.
418
419 :rtype: ``float``
420 """
421 m=self.__pre_input
422 grad_m=self.__pre_args
423
424 mu=self.__mu
425 mu_c=self.__mu_c
426 DIM=self.getDomain().getDim()
427 numLS=self.getNumLevelSets()
428
429 A=0
430 if self.__w0 is not None:
431 r = inner(integrate(m**2 * self.__w0), mu)
432 self.logger.debug("J_R[m^2] = %e"%r)
433 A += r
434
435 if self.__w1 is not None:
436 if numLS == 1:
437 r = integrate(inner(grad_m**2, self.__w1))*mu
438 self.logger.debug("J_R[grad(m)] = %e"%r)
439 A += r
440 else:
441 for k in range(numLS):
442 r = mu[k]*integrate(inner(grad_m[k,:]**2,self.__w1[k,:]))
443 self.logger.debug("J_R[grad(m)][%d] = %e"%(k,r))
444 A += r
445
446 if numLS > 1:
447 for k in range(numLS):
448 gk=grad_m[k,:]
449 len_gk=length(gk)
450 for l in range(k):
451 gl=grad_m[l,:]
452 r = mu_c[l,k] * integrate( self.__wc[l,k] * ( ( len_gk * length(gl) )**2 - inner(gk, gl)**2 ) )
453 self.logger.debug("J_R[cross][%d,%d] = %e"%(l,k,r))
454 A += r
455 return A/2
456
457
458 def getValue(self, m, grad_m):
459 """
460 returns the value of the cost function J with respect to m.
461 This equation is specified in the inversion cookbook.
462
463 :rtype: ``float``
464 """
465
466 if m!=self.__pre_input:
467 raise RuntimeError("Attempt to change point using getValue")
468 # substituting cached values
469 m=self.__pre_input
470 grad_m=self.__pre_args
471
472
473 mu=self.__mu
474 mu_c=self.__mu_c
475 DIM=self.getDomain().getDim()
476 numLS=self.getNumLevelSets()
477
478 A=0
479 if self.__w0 is not None:
480 r = inner(integrate(m**2 * self.__w0), mu)
481 self.logger.debug("J_R[m^2] = %e"%r)
482 A += r
483
484 if self.__w1 is not None:
485 if numLS == 1:
486 r = integrate(inner(grad_m**2, self.__w1))*mu
487 self.logger.debug("J_R[grad(m)] = %e"%r)
488 A += r
489 else:
490 for k in range(numLS):
491 r = mu[k]*integrate(inner(grad_m[k,:]**2,self.__w1[k,:]))
492 self.logger.debug("J_R[grad(m)][%d] = %e"%(k,r))
493 A += r
494
495 if numLS > 1:
496 for k in range(numLS):
497 gk=grad_m[k,:]
498 len_gk=length(gk)
499 for l in range(k):
500 gl=grad_m[l,:]
501 r = mu_c[l,k] * integrate( self.__wc[l,k] * ( ( len_gk * length(gl) )**2 - inner(gk, gl)**2 ) )
502 self.logger.debug("J_R[cross][%d,%d] = %e"%(l,k,r))
503 A += r
504 return A/2
505
506 def getGradient(self):
507 raise RuntimeError("Split versions do not support getGradient. Use getGradientAtPoint instead.")
508
509 def getGradientAtPoint(self):
510 """
511 returns the gradient of the cost function J with respect to m.
512
513 :note: This implementation returns Y_k=dPsi/dm_k and X_kj=dPsi/dm_kj
514 """
515
516 # Using cached values
517 m=self.__pre_input
518 grad_m=self.__pre_args
519
520 mu=self.__mu
521 mu_c=self.__mu_c
522 DIM=self.getDomain().getDim()
523 numLS=self.getNumLevelSets()
524
525 grad_m=grad(m, Function(m.getDomain()))
526 if self.__w0 is not None:
527 Y = m * self.__w0 * mu
528 else:
529 if numLS == 1:
530 Y = Scalar(0, grad_m.getFunctionSpace())
531 else:
532 Y = Data(0, (numLS,) , grad_m.getFunctionSpace())
533
534 if self.__w1 is not None:
535
536 if numLS == 1:
537 X=grad_m* self.__w1*mu
538 else:
539 X=grad_m*self.__w1
540 for k in range(numLS):
541 X[k,:]*=mu[k]
542 else:
543 X = Data(0, grad_m.getShape(), grad_m.getFunctionSpace())
544
545 # cross gradient terms:
546 if numLS > 1:
547 for k in range(numLS):
548 grad_m_k=grad_m[k,:]
549 l2_grad_m_k = length(grad_m_k)**2
550 for l in range(k):
551 grad_m_l=grad_m[l,:]
552 l2_grad_m_l = length(grad_m_l)**2
553 grad_m_lk = inner(grad_m_l, grad_m_k)
554 f = mu_c[l,k]* self.__wc[l,k]
555 X[l,:] += f * (l2_grad_m_k*grad_m_l - grad_m_lk*grad_m_k)
556 X[k,:] += f * (l2_grad_m_l*grad_m_k - grad_m_lk*grad_m_l)
557
558 return ArithmeticTuple(Y, X)
559
560 def getInverseHessianApproximationAtPoint(self, r, solve=True):
561 """
562 """
563
564 # substituting cached values
565 m=self.__pre_input
566 grad_m=self.__pre_args
567
568 if self._new_mu or self._update_Hessian:
569 self._new_mu=False
570 self._update_Hessian=False
571 mu=self.__mu
572 mu_c=self.__mu_c
573
574 DIM=self.getDomain().getDim()
575 numLS=self.getNumLevelSets()
576 if self.__w0 is not None:
577 if numLS == 1:
578 D=self.__w0 * mu
579 else:
580 D=self.getPDE().getCoefficient("D")
581 D.setToZero()
582 for k in range(numLS): D[k,k]=self.__w0[k] * mu[k]
583 self.getPDE().setValue(D=D)
584
585 A=self.getPDE().getCoefficient("A")
586 A.setToZero()
587 if self.__w1 is not None:
588 if numLS == 1:
589 for i in range(DIM): A[i,i]=self.__w1[i] * mu
590 else:
591 for k in range(numLS):
592 for i in range(DIM): A[k,i,k,i]=self.__w1[k,i] * mu[k]
593
594 if numLS > 1:
595 # this could be make faster by creating caches for grad_m_k, l2_grad_m_k and o_kk
596 for k in range(numLS):
597 grad_m_k=grad_m[k,:]
598 l2_grad_m_k = length(grad_m_k)**2
599 o_kk=outer(grad_m_k, grad_m_k)
600 for l in range(k):
601 grad_m_l=grad_m[l,:]
602 l2_grad_m_l = length(grad_m_l)**2
603 i_lk = inner(grad_m_l, grad_m_k)
604 o_lk = outer(grad_m_l, grad_m_k)
605 o_kl = outer(grad_m_k, grad_m_l)
606 o_ll=outer(grad_m_l, grad_m_l)
607 f= mu_c[l,k]* self.__wc[l,k]
608 Z=f * (2*o_lk - o_kl - i_lk*kronecker(DIM))
609 A[l,:,l,:] += f * (l2_grad_m_k*kronecker(DIM) - o_kk)
610 A[l,:,k,:] += Z
611 A[k,:,l,:] += transpose(Z)
612 A[k,:,k,:] += f * (l2_grad_m_l*kronecker(DIM) - o_ll)
613 self.getPDE().setValue(A=A)
614 #self.getPDE().resetRightHandSideCoefficients()
615 #self.getPDE().setValue(X=r[1])
616 #print "X only: ",self.getPDE().getSolution()
617 #self.getPDE().resetRightHandSideCoefficients()
618 #self.getPDE().setValue(Y=r[0])
619 #print "Y only: ",self.getPDE().getSolution()
620
621 self.getPDE().resetRightHandSideCoefficients()
622 self.getPDE().setValue(X=r[1], Y=r[0])
623 if not solve:
624 return self.getPDE()
625 return self.getPDE().getSolution()
626
627 def updateHessian(self):
628 """
629 notifies the class to recalculate the Hessian operator.
630 """
631 if not self.__useDiagonalHessianApproximation:
632 self._update_Hessian=True
633
634 def getNorm(self, m):
635 """
636 returns the norm of ``m``.
637
638 :param m: level set function
639 :type m: `Data`
640 :rtype: ``float``
641 """
642 return sqrt(integrate(length(m)**2)/self.__vol_d)
643
644 def setPoint(self, m):
645 """
646 sets the point which this function will work with
647
648 :param m: level set function
649 :type m: `Data`
650 """
651 self.__pre_input = m
652 self.__pre_args = grad(m)
653

  ViewVC Help
Powered by ViewVC 1.1.26