# Diff of /trunk/downunder/py_src/regularizations.py

revision 4073 by caltinay, Tue Oct 30 03:35:17 2012 UTC revision 4074 by gross, Thu Nov 15 03:30:59 2012 UTC
# Line 22  __url__="https://launchpad.net/escript-f Line 22  __url__="https://launchpad.net/escript-f
22
23  __all__ = ['Regularization']  __all__ = ['Regularization']
24
25    from costfunctions import CostFunction
26
27
28  import numpy as np  import numpy as np
29  from esys.escript.linearPDEs import LinearSinglePDE  from esys.escript.linearPDEs import LinearPDE, IllegalCoefficientValue
30  from esys.escript import Data, grad, inner, integrate, kronecker, vol  from esys.escript import Data, grad, inner, integrate, kronecker, boundingBoxEdgeLengths, interpolate
31    from esys.escript.pdetools import ArithmeticTuple
32
33  class Regularization(object):  class Regularization(CostFunction):
34      """      """
35        The regularization term for the level set function `m` within the cost function J for an inversion:
36
37        *J(m)=1/2 * sum_k imtegrate( mu_0[k]*s0[k] * m_k**2 + mu_1[k]*s1[k,i] * m_{k,i}**2) + sum_l<k mu_c[l,k] sc[l,l]* |curl(m_k) x curl(m_l)|^2*
38
39        where s0[k], s1[k,i] and  sc[k,l] are non-negative scaling factors and mu_0[k], mu_1[k], mu_c[l,k] are weighting factors
40        which may be alter during the inversion. The scaling factors are normalized such that their integrals over the
41        domain are constant:
42
43        *integrate(s0[k])=1*
44        *integrate(inner(s1[k,:],L[:])=1*
45        *integrate(inner(sc[l,k]*L**4)=1*
46
47      """      """
48      def __init__(self, domain, m_ref=0, w0=None, w=None, location_of_set_m=Data(), tol=1e-8):      def __init__(self, domain, numLevelSets=1,
49          """                         s0=None, s1=None, sc=None,
50          """                         location_of_set_m=Data(),
51                           useDiagonalHessianApproximation=True, tol=1e-8):
52            """
53            initialization
54
55            :param domain: domain
56            :type domain: `Domain`
57            :param numLevelSets: number of level sets
58            :type numLevelSets: ``int``
59            :param s0: scaling factor for the m**2 term. If not set zero is assumed.
60            :type s0: ``Scalar`` if `numLevelSets`==1 or `Data` object of shape ('numLevelSets`,) if numLevelSets > 1
61            :param s1: scaling factor for the grad(m_i) terms. If not set zero is assumed.
62            :type s1: ``Vector`` if `numLevelSets`==1 or `Data` object of shape (`numLevelSets`,DIM) if numLevelSets > 1.
63            :param sc: scaling factor for the cross correlation terms. If not set zero is assumed. Used for the case if numLevelSets > 1 only.
64                       values `sc[l,k]``  in the lower triangle (l<k) are used only.
65            :type sc: `Data` object of shape (`numLevelSets`,`numLevelSets`)
66            :param location_of_set_m: marks location of zero values of the level set function `m` by a positive entry.
67            :type location_of_set_m: ``Scalar`` if `numLevelSets`==1 or `Data` object of shape ('numLevelSets`,) if numLevelSets > 1.
68            :param useDiagonalHessianApproximation: if True cross correllation terms between level set components are ignored when calculating
69                                                    approximations of the inverse of the Hessian Operator. This can speep-up the calculation of
70                                                    the inverse but may lead to an increase of the number of iteration steps in the inversion.
71            :type useDiagonalHessianApproximation: ``bool``
72            :param tol: toleramce when solving the PDE for the inverse of the Hessian Operator
73            :type tol: positive ``float``
74
75
76            """
77            if numLevelSets>1:
78              raise ValueError("Currently only numLevelSets<=1 is supportered.")
79        if s0 == None and s1==None:
80              raise ValueError("Values for s0 or s1 must be given.")
81
82          self.__domain=domain          self.__domain=domain
83          self.__m_ref=m_ref          DIM=self.__domain.getDim()
84          self.location_of_set_m=location_of_set_m          self.__L2=np.asarray(boundingBoxEdgeLengths(domain))**2
85          if w0 is None:          self.__L4=np.sum(self.__L2)**2
86              self._w0 = None          self.__numLevelSets=numLevelSets
87          else:
88              self._w0=np.asarray(w0) / vol(self.__domain)          if self.__numLevelSets > 1:
89          if  w is None:              self.__useDiagonalHessianApproximation=useDiagonalHessianApproximation
self._w = None
90          else:          else:
91              self._w=np.asarray(w)          self.__useDiagonalHessianApproximation=True
92        self._update_Hessian=True
93          self.__projector=LinearSinglePDE(domain)
94          self.__projector.getSolverOptions().setTolerance(tol)      self.__pde=LinearPDE(self.__domain, numEquations=self.__numLevelSets)
95          self.__projector.setValue(A=kronecker(domain), q=location_of_set_m, D=0.)          self.__pde.getSolverOptions().setTolerance(tol)
96            self.__pde.setSymmetryOn()
97      def getInner(self, f0, f1):          self.__pde_is_set=False
98          """          try:
99          returns the inner product of two gradients.        self.__pde.setValue(q=location_of_set_m)
100          """      except IllegalCoefficientValue:
101          return integrate(inner(grad(f0),grad(f1)))        raise ValueError("Unable to set location of fixed level set function.")
102
103      def project(self, Y=Data(), X=Data()):      self.__total_num_weights=2*numLevelSets+((numLevelSets-1)*numLevelSets)/2
104          """      self.__weight_index=[]  # this is a mapping from the relevant mu-coefficients to the set of all mu-coefficients
105          """                         # we count s0, then s1, then sc (k<l).
106          self.__projector.setValue(Y=Y, X=X)      # THE S0 weighting factor
107          return  self.__projector.getSolution()      n=0
108            if not s0 is None:
109      def getValue(self, m):          s0 = interpolate(s0,self.__pde.getFunctionSpaceForCoefficient('D'))
110          """          s=s0.getShape()
111            if numLevelSets == 1 :
112                 if s == () :
113                 V=integrate(s0)
114                 if V > 0:
115                   self.__weight_index.append(n)
116                   s0*=1./V
117                     else:
118                   s0=None
119             else:
120                 raise ValueError("Unexpected shape %s for weight s0."%s)
121                else:
122                 if s == (numLevelSets,):
123                 for k in xrange(numLevelSets):
124                    V=integrate(s0[k])
125                    if Lsup(V) > 0:
126                        self.__weight_index.append(n+k)
127                        s0[k]*=1/V
128             else:
129                 raise ValueError("Unexpected shape %s for weight s0."%s)
130            self.__s0=s0
131            n+=numLevelSets
132
133            # The S1 weighting factor
134            if not s1 is None:
135            s1 = interpolate(s1,self.__pde.getFunctionSpaceForCoefficient('A'))
136            s=s1.getShape()
137
138            if numLevelSets == 1 :
139                 if s == (DIM,) :
140                   V=integrate(inner(s1, self.__L2))
141                   if V > 0:
142                  self.__weight_index.append(n)
143                      s1*=1/V
144             else:
145                 raise ValueError("Unexpected shape %s for weight s1."%s)
146                else:
147                 if s == (numLevelSets,DIM):
148                 for k in xrange(numLevelSets):
149                    for i in xrange(DIM):
150                    ww=s1[k,:]
151                    V=integrate(inner(ww,self.__L2))
152                        if V > 0:
153                           self.__weight_index.append(n+i)
154                           s1[k,:]=ww*(1./V)
155             else:
156                 raise ValueError("Unexpected shape %s for weight s1."%s)
157
158            self.__s1=s1
159            n+=numLevelSets
160
161            # The SC weighting factor
162            if not sc is None:
163            if numLevelSets == 1 :
164                  sc=None
165                else:
166                  sc = interpolate(sc,self.__pde.getFunctionSpaceForCoefficient('A'))
167              s=sc.getShape()
168              if s == (numLevelSets,numLevelSets):
169                 for k in xrange(numLevelSets):
170                    sc[k,k]=0.
171                    for l in xrange(k):
172                    ww=sc[l,k]
173                    V=integrate(ww)
174                        if V > 0:
175                           self.__weight_index.append(n+k*numLevelSets+l)
176                           sc[l,k]=ww/V*self.__L4
177                           sc[k,l]=0
178                  else:
179                 raise ValueError("Unexpected shape %s for weight s0."%s)
180
181            self.__sc=sc
182        self.setWeights()
183
184        def getDomain(self):
185            """
186            return the domain of the regularization term
187            :rtype: ``Domain``
188            """
189            return self.__domain
190
191        def getNumLevelSets(self):
192            """
193            return the number of level set functions
194            :rtype: ``int``
195            """
196            return self.__numLevelSets
197
198        def getPDE(self):
199            """
200            returns the linear PDE to be solved for the Hessian Operator inverse
201            :rtype: ``LinearPDE``
202            """
203            return self.__pde
204
205        def getDualProduct(self, m, r):
206            """
207            returns the dual product of a gradient represented by X=r[1] and Y=r[0] with a level set function m:
208
209                 *Y_i*m_i + X_ij*m_{i,j}*
210
211            :type m: ``esys.escript.Data``
212            :type r: ``ArithmeticTuple``
213            :rtype: float
214          """          """
215          A=0          A=0
216          if self._w0 is not None:          if not r[0].isEmpty(): A+=integrate(inner(r[0], m))
217              A=(m-self.__m_ref)**2 * self._w0          if not r[1].isEmpty(): A+=integrate(inner(r[1], grad(m)))
218          if self._w is not None:      return A
219              A=inner(self._w, grad(m-self.__m_ref)**2) + A
220          return integrate(A)      def getValue(self, m, grad_m):
221            """
222      def getGradient(self, m):          return the value of the costfunction J
223
224
225
226            :rtype: ``float``
227            """
228            mu=self.getWeights( uncompress=True)
229            DIM=self.getDomain().getDim()
230            numLS=self.getNumLevelSets()
231
232            A=0
233            n=0
234
235            if self.__s0 is not None:
236                A+=inner(integrate(m**2*self.__s0), mu[:numLS])
237            n+=numLS
238
239            if self.__s1 is not None:
240            if numLS == 1:
241                A+=integrate(inner(grad_m**2, self.__s1))*mu[n]
242            else:
243                for k in xrange(numLS):
244                A+=mu[n+k]*integrate(inner(grad_m[k,:]**2,self.__s1[k,:]))
245            n+=numLS
246
247            if self.__sc is not None:
248              for k in xrange(numLS):
249               gk=grad_m[k,:]
250               len_gk=length(gk)
251               for l in xrange(k):
252                   gl=grad_m[l,:]
253                   A+= mu[n+k*numLS+l] * integrate( self.__sc[l,k] * ( len_gk * length(gl) )**2 - inner(gk, gl)**2 )
254            return A/2
255
256        def getDirectionalDerivative(self, m, d, grad_m):
257            """
258            returns the directional derivative
259            """
260            mu=self.getWeights( uncompress=True)
261            DIM=self.getDomain().getDim()
262            numLS=self.getNumLevelSets()
263            grad_q=grad(q)
264
265            A=0
266            n=0
267
268            if self.__s0 is not None:
269                A+=inner(integrate(m*q*self.__s0), mu[:numLS])
270            n+=numLS
271
272            if self.__s1 is not None:
273            if numLS == 1:
274                A+=integrate(inner(grad_m* self.__s1,grad_q))*mu[n]
275            else:
276                for k in xrange(numLS):
277                A+=mu[n+k]*integrate(inner(grad_m[k,:]*self.__s1[k,:],grad_q[k,:]))
278            n+=numLS*DIM
279
280            # this needs to be checked tested:
281            if self.__sc is not None:
282              for k in xrange(numLS):
283               gm_k=grad_m[k,:]
284               gq_k=grad_q[k,:]
285               len2_gm_k=length(gm_k)**2
286               gmq_k = inner(gm_k,gq_k)
287               for l in xrange(k):
288                   gm_l=grad_m[l,:]
289                   gq_l=grad_q[l,:]
290                   A+= mu[n+k*numLS+l] * integrate( self.__sc[l,k] * (
291                      ( gmq_k * length(gm_l)**2 + len2_gm_k * inner(gm_l,gq_l)
292                               - inner(gm_k, gm_l)*( inner(gm_k, gq_l) + inner(gq_k, gm_l) ) )
293                      ) )
294
295            return A
296
297        def getGradient(self, m,  grad_m):
298            """
299            returns the gradient of the costfunction J with respect to m. The function returns Y_k=dPsi/dm_k and
300            X_kj=dPsi/dm_kj
301            """
302
303            mu=self.getWeights( uncompress=True)
304            DIM=self.getDomain().getDim()
305            numLS=self.getNumLevelSets()
306
307            n=0
308
309            if self.__s0 is not None:
310            Y = m * self.__s0 * mu[:numLS]
311        else:
312            Y = Data()
313            n+=numLS
314
315            if self.__s1 is not None:
316            if numLS == 1:
317                X=grad_m* self.__s1*mu[n]
318            else:
319                X=self.getPDE().createCoefficient("X")
320                for k in xrange(numLS):
321                X[k,:]=mu[n+k]*grad_m[k,:]*self.__s1[k,:]
322        else:
323            X = Data()
324            n+=numLS
325            if self.__sc is not None:
326               raise NotImplementedError
327            return ArithmeticTuple(Y, X)
328
329        def getArguments(self, m):
330          """          """
331          """          """
332          if not self._w0 == None:          return ( grad(m),)
333              Y=2. * (m-self.__m_ref) * self._w0
334          else:      def getInverseHessianApproximation(self, m, r, grad_m):
335              Y=0.          """
336          if not self._w == None:          """
337              X=2. * grad(m-self.__m_ref) * self._w          if self._new_mu or self._update_Hessian:
338          else:
339              X=np.zeros((self.__domain.getDim(),))          self._new_mu=False
340            self._update_Hessian=False
341          return Y, X
342            mu=self.getWeights( uncompress=True)
343      def getArguments(self, m):              DIM=self.getDomain().getDim()
344                numLS=self.getNumLevelSets()
345                n=0
346                if self.__s0 is not None:
347                if numLS == 1:
348                     D=self.__s0 * mu[n]
349                     print "D =",D
350                    else:
351                         D=self.getPDE().createCoefficient("D")
352                     for k in xrange(numLS): D[k,k]=self.__s0[k] * mu[n+k]
353                self.getPDE().setValue(D=D)
354            n+=numLS
355
356            A=self.getPDE().createCoefficient("A")
357                if self.__s1 is not None:
358               if numLS == 1:
359               for i in xrange(DIM): A[i,i]=self.__s1[i] * mu[n]
360               print "A =",A
361               else:
362                   for k in xrange(numLS):
363                        for i in xrange(DIM): A[k,i,k,i]=self.__s1[k,i] * mu[n+k]
364                n+=numLS
365                if self.__sc is not None:
366                    raise NotImplementedError
367                self.getPDE().setValue(A=A)
368            self.getPDE().resetRightHandSideCoefficients()
369            print A
370            print r[1]
371            print r[0]
372            self.getPDE().setValue(X=r[1], Y=r[0])
373            return self.getPDE().getSolution()
374
375        def updateHessian(self):
376            """
377            notify the class to recalculate the Hessian operator
378            """
379        if not self.__useDiagonalHessianApproximation:
380            self._update_Hessian=True
381
382       # ================ we should factor these out =============================================================
383        def getNumRelevantTerms(self):
384            """
385            returns the number of terms in the costfunction of regularization with non-zero scaling factors
386            :rtype: ``int``
387            """
388            return len(self.__weight_index)
389
390        def getNumTerms(self):
391            """
392            returns the number of terms in the costfunction of regularization with non-zero scaling factors
393            :rtype: ``int``
394            """
395            return len(self.__weight_index)
396
397        def setWeights(self, mu=None):
398            """
399            sets the values for the weighting terms in the costsfunction. Note that only values
400            correspond to entries with non-negative scaling factors are used.
401
402            :param mu: weights
403            :type mu: ``list`` of ``float`` or ``np,array``
404            """
405            if mu == None:
406            mu = np.ones(self.getNumRelevantTerms())
407        mu=np.asarray(mu)
408
409
410        if len(mu) == self.getNumRelevantTerms():
411             if not mu.shape == (self.getNumRelevantTerms(),):
412            raise ValueError("%s values are expected."%self.getNumRelevantTerms())
413             self.__mu=mu
414        else:
415             if not mu.shape == (self.__total_num_weights,):
416            raise ValueError("%s values are expected."%self.__total_num_weights)
417             self.__mu = np.zeros(self.getNumRelevantTerms())
418             for i in xrange(len(self.__weight_index)): self.__mu[i]=mu[self.__weight_index[i]]
419        self._new_mu=True
420
421        def setWeightsForS0(self, mu=None):
422             """
423             set the weights for s0-terms
424             """
425             numLS=self.getNumLevelSets()
426             my_mu=self.getWeights(uncompress=True)
427             if mu is None:
428                 my_mu[:numLS]=1
429             else:
430             my_mu[:numLS]=mu
431             self.setWeights(my_mu)
432
433        def setWeightsForS1(self, mu=None):
434             """
435             set the weights for s1-terms
436             """
437             numLS=self.getNumLevelSets()
438             my_mu=self.getWeights(uncompress=True)
439             if mu is None:
440                  my_mu[numLS:2*numLS]=1
441             else:
442                  my_mu[numLS:2*numLS]=mu
443             self.setWeights(my_mu)
444
445        def setWeightsForSc(self, mu):
446             """
447             set the weights for s1-terms
448             """
449             numLS=self.getNumLevelSets()
450             my_mu=self.getWeights(uncompress=True)
451             if mu is None:
452                  my_mu[2*numLS:]=1
453             else:
454                  my_mu[2*numLS:]=mu
455             self.setWeights(my_mu)
456
457
458        def getWeights(self, uncompress=False):
459            """
460            Returns the weights for the terms in the costsfunction.
461            The first ``numLevelSets`` values used for the
462            regularization terms and the remaining values for the cross correlation terms.
463
464           :type mu: ``list`` of ``float``
465            """
466            if uncompress:
467             mu = np.zeros(self.__total_num_weights)
468             for i in xrange(len(self.__weight_index)): mu[self.__weight_index[i]] = self.__mu[i]
469             return mu
470        else:
471               return self.__mu
472
473        def getWeightIndex(self):
474          """          """
475            returns an iondex to the contributions of terms with non-zero scaling factor.
476          """          """
477          return ()          return self.__weight_index
478

Legend:
 Removed from v.4073 changed lines Added in v.4074

 ViewVC Help Powered by ViewVC 1.1.26