/[escript]/branches/subworld2/downunder/py_src/regularizations.py
ViewVC logotype

Contents of /branches/subworld2/downunder/py_src/regularizations.py

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.26