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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 4097 by caltinay, Fri Dec 7 01:18:35 2012 UTC revision 4154 by jfenwick, Tue Jan 22 09:30:23 2013 UTC
# Line 1  Line 1 
1    
2  ##############################################################################  ##############################################################################
3  #  #
4  # Copyright (c) 2003-2012 by University of Queensland  # Copyright (c) 2003-2013 by University of Queensland
5  # http://www.uq.edu.au  # http://www.uq.edu.au
6  #  #
7  # Primary Business: Queensland, Australia  # Primary Business: Queensland, Australia
# Line 13  Line 13 
13  #  #
14  ##############################################################################  ##############################################################################
15    
16  __copyright__="""Copyright (c) 2003-2012 by University of Queensland  __copyright__="""Copyright (c) 2003-2013 by University of Queensland
17  http://www.uq.edu.au  http://www.uq.edu.au
18  Primary Business: Queensland, Australia"""  Primary Business: Queensland, Australia"""
19  __license__="""Licensed under the Open Software License version 3.0  __license__="""Licensed under the Open Software License version 3.0
# Line 25  __all__ = ['Regularization'] Line 25  __all__ = ['Regularization']
25  from costfunctions import CostFunction  from costfunctions import CostFunction
26    
27  import numpy as np  import numpy as np
28  from esys.escript import Data, grad, inner, integrate, interpolate, kronecker, boundingBoxEdgeLengths, vol  from esys.escript import ReducedFunction, outer, Data, Scalar, grad, inner, integrate, interpolate, kronecker, boundingBoxEdgeLengths, vol, sqrt, length
29  from esys.escript.linearPDEs import LinearPDE, IllegalCoefficientValue  from esys.escript.linearPDEs import LinearPDE, IllegalCoefficientValue
30  from esys.escript.pdetools import ArithmeticTuple  from esys.escript.pdetools import ArithmeticTuple
31    
# Line 34  class Regularization(CostFunction): Line 34  class Regularization(CostFunction):
34      The regularization term for the level set function ``m`` within the cost      The regularization term for the level set function ``m`` within the cost
35      function J for an inversion:      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*      *J(m)=1/2 * sum_k imtegrate( 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*
38    
39      where s0[k], s1[k,i] and  sc[k,l] are non-negative scaling factors and      where w0[k], w1[k,i] and  wc[k,l] are non-negative weighting factors and
40      mu_0[k], mu_1[k], mu_c[l,k] are weighting factors which may be altered      mu[k] and mu_c[l,k] are trade-off factors which may be altered
41      during the inversion. The scaling factors are normalized such that their      during the inversion. The weighting factors are normalized such that their
42      integrals over the domain are constant:      integrals over the domain are constant:
43    
44      *integrate(s0[k])=1*      *integrate(w0[k] + inner(w1[k,:],1/L[:]**2))=scale[k]* volume(domain)*
45      *integrate(inner(s1[k,:],L[:])=1*      *integrate(wc[l,k]*1/L**4)=scale_c[k]* volume(domain) *
     *integrate(inner(sc[l,k]*L**4)=1*  
46    
47      """      """
48      def __init__(self, domain, numLevelSets=1,      def __init__(self, domain, numLevelSets=1,
49                         s0=None, s1=None, sc=None,                         w0=None, w1=None, wc=None,
50                         location_of_set_m=Data(),                         location_of_set_m=Data(),
51                         useDiagonalHessianApproximation=True, tol=1e-8):                         useDiagonalHessianApproximation=False, tol=1e-8,
52                           scale=None, scale_c=None):
53          """          """
54          initialization.          initialization.
55    
# Line 57  class Regularization(CostFunction): Line 57  class Regularization(CostFunction):
57          :type domain: `Domain`          :type domain: `Domain`
58          :param numLevelSets: number of level sets          :param numLevelSets: number of level sets
59          :type numLevelSets: ``int``          :type numLevelSets: ``int``
60          :param s0: scaling factor for the m**2 term. If not set zero is assumed.          :param w0: weighting factor for the m**2 term. If not set zero is assumed.
61          :type s0: ``Scalar`` if ``numLevelSets`` == 1 or `Data` object of shape          :type w0: ``Scalar`` if ``numLevelSets`` == 1 or `Data` object of shape
62                    (``numLevelSets`` ,) if ``numLevelSets`` > 1                    (``numLevelSets`` ,) if ``numLevelSets`` > 1
63          :param s1: scaling factor for the grad(m_i) terms. If not set zero is assumed          :param w1: weighting factor for the grad(m_i) terms. If not set zero is assumed
64          :type s1: ``Vector`` if ``numLevelSets`` == 1 or `Data` object of shape          :type w1: ``Vector`` if ``numLevelSets`` == 1 or `Data` object of shape
65                    (``numLevelSets`` , DIM) if ``numLevelSets`` > 1                    (``numLevelSets`` , DIM) if ``numLevelSets`` > 1
66          :param sc: scaling factor for the cross correlation terms. If not set          :param wc: weighting factor for the cross gradient terms. If not set
67                     zero is assumed. Used for the case if ``numLevelSets`` > 1                     zero is assumed. Used for the case if ``numLevelSets`` > 1
68                     only. Only values ``sc[l,k]`` in the lower triangle (l<k)                     only. Only values ``wc[l,k]`` in the lower triangle (l<k)
69                     are used.                     are used.
70          :type sc: `Data` object of shape (``numLevelSets`` , ``numLevelSets``)          :type wc: `Data` object of shape (``numLevelSets`` , ``numLevelSets``)
71          :param location_of_set_m: marks location of zero values of the level          :param location_of_set_m: marks location of zero values of the level
72                                    set function ``m`` by a positive entry.                                    set function ``m`` by a positive entry.
73          :type location_of_set_m: ``Scalar`` if ``numLevelSets`` == 1 or `Data`          :type location_of_set_m: ``Scalar`` if ``numLevelSets`` == 1 or `Data`
74                  object of shape (``numLevelSets`` ,) if ``numLevelSets`` > 1                  object of shape (``numLevelSets`` ,) if ``numLevelSets`` > 1
75          :param useDiagonalHessianApproximation: if True cross correlation terms          :param useDiagonalHessianApproximation: if True cross gradient terms
76                      between level set components are ignored when calculating                      between level set components are ignored when calculating
77                      approximations of the inverse of the Hessian Operator.                      approximations of the inverse of the Hessian Operator.
78                      This can speed-up the calculation of the inverse but may                      This can speed-up the calculation of the inverse but may
# Line 82  class Regularization(CostFunction): Line 82  class Regularization(CostFunction):
82          :param tol: tolerance when solving the PDE for the inverse of the          :param tol: tolerance when solving the PDE for the inverse of the
83                      Hessian Operator                      Hessian Operator
84          :type tol: positive ``float``          :type tol: positive ``float``
85            
86            :param scale: weighting factor for level set function variation terms. If not set one is used.
87            :type scale: ``Scalar`` if ``numLevelSets`` == 1 or `Data` object of shape
88                      (``numLevelSets`` ,) if ``numLevelSets`` > 1
89            :param scale_c: scale for the cross gradient terms. If not set
90                       one is assumed. Used for the case if ``numLevelSets`` > 1
91                       only. Only values ``scale_c[l,k]`` in the lower triangle (l<k)
92                       are used.
93            :type scale_c: `Data` object of shape (``numLevelSets`` , ``numLevelSets``)
94            
95    
96          """          """
97          if numLevelSets>1:          if w0 == None and w1==None:
98                raise ValueError("Currently only numLevelSets<=1 is supported.")                raise ValueError("Values for w0 or for w1 must be given.")
99          if s0 == None and s1==None:          if wc == None and  numLevelSets>1:
100                raise ValueError("Values for s0 or s1 must be given.")                raise ValueError("Values for wc must be given.")
101    
102          self.__domain=domain          self.__domain=domain
103          DIM=self.__domain.getDim()          DIM=self.__domain.getDim()
         self.__L2=np.asarray(boundingBoxEdgeLengths(domain))**2  
         self.__L4=np.sum(self.__L2)**2  
104          self.__numLevelSets=numLevelSets          self.__numLevelSets=numLevelSets
105            
         if self.__numLevelSets > 1:  
             self.__useDiagonalHessianApproximation=useDiagonalHessianApproximation  
         else:  
             self.__useDiagonalHessianApproximation=True  
         self._update_Hessian=True  
   
106          self.__pde=LinearPDE(self.__domain, numEquations=self.__numLevelSets)          self.__pde=LinearPDE(self.__domain, numEquations=self.__numLevelSets)
107          self.__pde.getSolverOptions().setTolerance(tol)          self.__pde.getSolverOptions().setTolerance(tol)
108          self.__pde.setSymmetryOn()          self.__pde.setSymmetryOn()
         self.__pde_is_set=False  
109          try:          try:
110              self.__pde.setValue(q=location_of_set_m)              self.__pde.setValue(q=location_of_set_m)
111          except IllegalCoefficientValue:          except IllegalCoefficientValue:
112              raise ValueError("Unable to set location of fixed level set function.")              raise ValueError("Unable to set location of fixed level set function.")
113              
114          self.__total_num_weights=2*numLevelSets+((numLevelSets-1)*numLevelSets)/2          # =========== check the shape of the scales: =================================
115          self.__weight_index=[]  # this is a mapping from the relevant mu-coefficients to the set of all mu-coefficients          if scale is None:
                            # we count s0, then s1, then sc (k<l).  
         # THE S0 weighting factor  
         n=0  
         VV=vol(domain)  
         if s0 is not None:  
             s0 = interpolate(s0,self.__pde.getFunctionSpaceForCoefficient('D'))  
             s=s0.getShape()  
116              if numLevelSets == 1 :              if numLevelSets == 1 :
117                   if s == () :                 scale = 1.
                      V=integrate(s0)  
                      if V > 0:  
                        self.__weight_index.append(n)  
                        s0*=VV/V  
                      else:  
                        s0=None  
                  else:  
                      raise ValueError("Unexpected shape %s for weight s0."%s)  
118              else:              else:
119                   if s == (numLevelSets,):                 scale = np.ones((numLevelSets,))
120                       for k in xrange(numLevelSets):          else:
121                          V=integrate(s0[k])              scale=np.asarray(scale)
                         if V > 0:  
                             self.__weight_index.append(n+k)  
                             s0[k]*=VV/V  
                  else:  
                      raise ValueError("Unexpected shape %s for weight s0."%s)  
         self.__s0=s0  
         n+=numLevelSets  
   
         # The S1 weighting factor  
         if not s1 is None:  
             s1 = interpolate(s1,self.__pde.getFunctionSpaceForCoefficient('A'))  
             s=s1.getShape()  
   
122              if numLevelSets == 1 :              if numLevelSets == 1 :
123                  if s == (DIM,) :                  if scale.shape == ():
124                      V=integrate(inner(s1, 1/self.__L2))                     if not scale > 0 :
125                      if V > 0:                        raise ValueError("Value for scale must be positive.")
                         self.__weight_index.append(n)  
                         s1*=VV/V  
                     print "REG SCALE = ",s1  
126                  else:                  else:
127                      raise ValueError("Unexpected shape %s for weight s1."%s)                     raise ValueError("Unexpected shape %s for scale."%scale.shape)
128              else:              else:
129                  if s == (numLevelSets,DIM):                   if scale.shape is (numLevelSets,):
130                      for k in xrange(numLevelSets):                       if not min(scale) > 0:
131                          for i in xrange(DIM):                          raise ValueError("All value for scale must be positive.")
132                              ww=s1[k,:]                   else:
133                              V=integrate(inner(ww,1/self.__L2))                     raise ValueError("Unexpected shape %s for scale."%scale.shape)
134                              if V > 0:          
135                                  self.__weight_index.append(n+i)          if scale_c is None or numLevelSets < 2:
136                                  s1[k,:]=ww*(VV/V)              scale_c = np.ones((numLevelSets,numLevelSets))
137                  else:          else:
138                      raise ValueError("Unexpected shape %s for weight s1."%s)              scale_c=np.asarray(scale_c)
139                if scale_c.shape == (numLevelSets,numLevelSets):
140          self.__s1=s1                  if not all( [ [ scale_c[l,k] > 0. for l in xrange(k) ] for k in xrange(1,numLevelSets) ]):
141          n+=numLevelSets                          raise ValueError("All values in the lower triangle of scale_c must be positive.")
   
         # The SC weighting factor  
         if not sc is None:  
             if numLevelSets == 1 :  
                 sc=None  
142              else:              else:
143                  sc = interpolate(sc,self.__pde.getFunctionSpaceForCoefficient('A'))                   raise ValueError("Unexpected shape %s for scale."%scale_c.shape)
144                  s=sc.getShape()          # ===== check the shape of the weights: ============================================
145                  if s == (numLevelSets,numLevelSets):          if w0 is not None:
146                      for k in xrange(numLevelSets):                w0 = interpolate(w0,self.__pde.getFunctionSpaceForCoefficient('D'))
147                          sc[k,k]=0.                s0=w0.getShape()
148                          for l in xrange(k):                if numLevelSets == 1 :
149                              ww=sc[l,k]                     if  not s0 == () :
150                              V=integrate(ww)                        raise ValueError("Unexpected shape %s for weight w0."%s0)
151                              if V > 0:                else:
152                                  self.__weight_index.append(n+k*numLevelSets+l)                     if not s0 == (numLevelSets,):
153                                  sc[l,k]=ww*VV/V*self.__L4                        raise ValueError("Unexpected shape %s for weight w0."%s0)
154                                  sc[k,l]=0          if not w1 is None:
155                  else:                w1 = interpolate(w1,self.__pde.getFunctionSpaceForCoefficient('A'))
156                      raise ValueError("Unexpected shape %s for weight s0."%s)                s1=w1.getShape()
157                  if numLevelSets is 1 :
158                       if not s1 == (DIM,) :
159                          raise ValueError("Unexpected shape %s for weight w1."%s1)
160                  else:
161                       if not s1 == (numLevelSets,DIM):
162                          raise ValueError("Unexpected shape %s for weight w1."%s1)
163            if numLevelSets == 1 :
164                 wc=None
165            else:
166                 wc = interpolate(wc,self.__pde.getFunctionSpaceForCoefficient('A'))
167                 sc=wc.getShape()
168                 if not sc == (numLevelSets, numLevelSets):
169                    raise ValueError("Unexpected shape %s for weight wc."%(sc,))
170            # ============= now we rescale weights: ======================================
171            L2s=np.asarray(boundingBoxEdgeLengths(domain))**2
172            L4=1/np.sum(1/L2s)**2
173            if numLevelSets == 1 :
174                A=0
175                if w0 is not None:
176                    A = integrate(w0)
177                if w1 is not None:
178                    A += integrate(inner(w1, 1/L2s))
179                if A > 0:
180                    f = scale/A
181                    if w0 is not None:
182                         w0*=f
183                    if w1 is not None:
184                         w1*=f
185                else:
186                   raise ValueError("Non-positive weighting factor detected.")
187            else:
188    
189                 for k in xrange(numLevelSets):
190                     A=0
191                     if w0 is not None:
192                         A = integrate(w0[k])
193                     if w1 is not None:
194                          A += integrate(inner(w1[k,:], 1/L2s))
195                     if A > 0:
196                          f = scale[k]/A
197                          if w0 is not None:
198                             w0[k]*=f
199                          if w1 is not None:
200                             w1[k,:]*=f
201                     else:
202                       raise ValueError("Non-positive weighting factor for level set component %d detected."%k)
203                    
204                     # and now the cross-gradient:
205                     if wc is not None:
206                         for l in xrange(k):  
207                            A = integrate(wc[l,k])/L4
208                            if A > 0:
209                               f = scale_c[l,k]/A
210                               wc[l,k]*=f
211    #                        else:
212    #                          raise ValueError("Non-positive weighting factor for cross-gradient level set components %d and %d detected."%(l,k))
213                        
214            self.__w0=w0
215            self.__w1=w1
216            self.__wc=wc
217    
218          self.__sc=sc          self.__pde_is_set=False        
219          self.setWeights()          if self.__numLevelSets > 1:
220                self.__useDiagonalHessianApproximation=useDiagonalHessianApproximation
221            else:
222                self.__useDiagonalHessianApproximation=True
223            self._update_Hessian=True
224    
225            self.__num_tradeoff_factors=numLevelSets+((numLevelSets-1)*numLevelSets)/2
226            self.setTradeOffFactors()
227            self.__vol_d=vol(self.__domain)
228            
229      def getDomain(self):      def getDomain(self):
230          """          """
231          returns the domain of the regularization term          returns the domain of the regularization term
# Line 216  class Regularization(CostFunction): Line 249  class Regularization(CostFunction):
249          :rtype: `LinearPDE`          :rtype: `LinearPDE`
250          """          """
251          return self.__pde          return self.__pde
252        
253      def getDualProduct(self, m, r):      def getDualProduct(self, m, r):
254          """          """
255          returns the dual product of a gradient represented by X=r[1] and Y=r[0]          returns the dual product of a gradient represented by X=r[1] and Y=r[0]
# Line 232  class Regularization(CostFunction): Line 265  class Regularization(CostFunction):
265          if not r[0].isEmpty(): A+=integrate(inner(r[0], m))          if not r[0].isEmpty(): A+=integrate(inner(r[0], m))
266          if not r[1].isEmpty(): A+=integrate(inner(r[1], grad(m)))          if not r[1].isEmpty(): A+=integrate(inner(r[1], grad(m)))
267          return A          return A
268        def getNumTradeOffFactors(self):
269            """
270            returns the number of trade-off factors being used.
271    
272            :rtype: ``int``
273            """
274            return self.__num_tradeoff_factors
275    
276        def setTradeOffFactors(self, mu=None):
277            """
278            sets the trade-off factors for the level-set variation and the cross-gradient
279            
280            :param mu: new values for the trade-off factors where values mu[:numLevelSets] are the
281                       trade-off factors for the level-set variation and the remaining values for
282                       the cross-gradient part with mu_c[l,k]=mu[numLevelSets+l+((k-1)*k)/2] (l<k).
283                       If no values for mu is given ones are used. Values must be positive.
284            :type mu: ``list`` of ``float`` or ```numpy.array```
285            """
286            numLS=self.getNumLevelSets()
287            numTF=self.getNumTradeOffFactors()
288            if mu is None:
289               mu = np.ones((numTF,))
290            else:
291               mu = np.asarray(mu)
292    
293            if mu.shape == (numTF,):
294                self.setTradeOffFactorsForVariation(mu[:numLS])
295                mu_c2=np.zeros((numLS,numLS))
296                for k in xrange(numLS):
297                   for l in xrange(k):
298                       mu_c2[l,k] = mu[numLS+l+((k-1)*k)/2]
299                self.setTradeOffFactorsForCrossGradient(mu_c2)
300            elif mu.shape == () and numLS ==1:
301                self.setTradeOffFactorsForVariation(mu)
302            else:
303               raise ValueError("Unexpected shape %s for mu."%(mu.shape,))
304              
305        def setTradeOffFactorsForVariation(self, mu=None):
306             """
307             sets the trade-off factors for the level-set variation part
308            
309             :param mu:  new values for the trade-off factors. Values must be positive.
310             :type mu: ``float``, ``list`` of ``float`` or ```numpy.array```
311             """
312             numLS=self.getNumLevelSets()
313             if mu is None:
314                if numLS == 1:
315                   mu = 1.
316                else:
317                   mu = np.ones((numLS,))
318    
319             mu=np.asarray(mu)
320             if numLS == 1:
321               if mu.shape == (1,): mu=mu[0]
322               if mu.shape == ():
323                  if mu > 0:
324                     self.__mu= mu
325                     self._new_mu=True
326                  else:
327                     raise ValueError("Value for trade-off factor must be positive.")
328               else:
329                  raise ValueError("Unexpected shape %s for mu."%mu.shape)
330             else:
331               if mu.shape == (numLS,):
332                   if min(mu) > 0:
333                       self.__mu= mu
334                       self._new_mu=True
335                   else:
336                       raise ValueError("All value for mu must be positive.")
337               else:
338                   raise ValueError("Unexpected shape %s for trade-off factor."%mu.shape)
339    
340        def setTradeOffFactorsForCrossGradient(self, mu_c=None):
341            """
342            sets the trade-off factors for the cross-gradient terms
343            
344            :param mu_c: new values for the trade-off factors for the cross-gradient
345                         terms. Values must be positive. If no value is given ones
346                         are used. Onky value mu_c[l,k] for l<k are used.
347            :type mu_c: ``float``, ``list`` of ``float`` or ``numpy.array``
348            
349            """
350            numLS=self.getNumLevelSets()
351            if mu_c is None or numLS < 2:
352                self.__mu_c = np.ones((numLS,numLS))
353            if isinstance(mu_c, float) or isinstance(mu_c, int):
354                self.__mu_c = np.zeros((numLS,numLS))
355                self.__mu_c[:,:]=mu_c
356            else:
357                mu_c=np.asarray(mu_c)
358                if mu_c.shape == (numLS,numLS):
359                    if not all( [ [ mu_c[l,k] > 0. for l in xrange(k) ] for k in xrange(1,numLS) ]):
360                         raise ValueError("All trade-off factors in the lower triangle of mu_c must be positive.")
361                    else:
362                         self.__mu_c =  mu_c
363                         self._new_mu=True
364                else:
365                     raise ValueError("Unexpected shape %s for mu."%mu_c.shape)
366        
367        def getArguments(self, m):
368            """
369            """
370            return ( grad(m),)
371                    
372    
373      def getValue(self, m, grad_m):      def getValue(self, m, grad_m):
374          """          """
# Line 239  class Regularization(CostFunction): Line 376  class Regularization(CostFunction):
376    
377          :rtype: ``float``          :rtype: ``float``
378          """          """
379          mu=self.getWeights( uncompress=True)          mu=self.__mu
380            mu_c=self.__mu_c
381          DIM=self.getDomain().getDim()          DIM=self.getDomain().getDim()
382          numLS=self.getNumLevelSets()          numLS=self.getNumLevelSets()
383    
384          A=0          A=0
385          n=0          if self.__w0 is not None:
386               A+=inner(integrate(m**2 * self.__w0), mu)
387    
388          if self.__s0 is not None:          if self.__w1 is not None:
             A+=inner(integrate(m**2*self.__s0), mu[:numLS])  
         n+=numLS  
   
         if self.__s1 is not None:  
389              if numLS == 1:              if numLS == 1:
390                  A+=integrate(inner(grad_m**2, self.__s1))*mu[n]                  A+=integrate(inner(grad_m**2, self.__w1))*mu
391              else:              else:
392                  for k in xrange(numLS):                  for k in xrange(numLS):
393                      A+=mu[n+k]*integrate(inner(grad_m[k,:]**2,self.__s1[k,:]))                      A+=mu[k]*integrate(inner(grad_m[k,:]**2,self.__w1[k,:]))
394          n+=numLS                      
395            if numLS > 1:
         if self.__sc is not None:  
396              for k in xrange(numLS):              for k in xrange(numLS):
397                  gk=grad_m[k,:]                  gk=grad_m[k,:]
398                  len_gk=length(gk)                  len_gk=length(gk)
399                  for l in xrange(k):                  for l in xrange(k):
400                      gl=grad_m[l,:]                      gl=grad_m[l,:]
401                      A+= mu[n+k*numLS+l] * integrate( self.__sc[l,k] * ( len_gk * length(gl) )**2 - inner(gk, gl)**2 )                      A+= mu_c[l,k] * integrate( self.__wc[l,k] * ( len_gk * length(gl) )**2 - inner(gk, gl)**2 )
402          return A/2          return A/2
403    
404      def getGradient(self, m,  grad_m):      def getGradient(self, m,  grad_m):
# Line 273  class Regularization(CostFunction): Line 407  class Regularization(CostFunction):
407          The function returns Y_k=dPsi/dm_k and X_kj=dPsi/dm_kj          The function returns Y_k=dPsi/dm_k and X_kj=dPsi/dm_kj
408          """          """
409    
410          mu=self.getWeights( uncompress=True)          mu=self.__mu
411            mu_c=self.__mu_c
412          DIM=self.getDomain().getDim()          DIM=self.getDomain().getDim()
413          numLS=self.getNumLevelSets()          numLS=self.getNumLevelSets()
414    
415          n=0          print "WARNING: WRONG FUNCTION SPACE"
416            
417          if self.__s0 is not None:          grad_m=grad(m, ReducedFunction(m.getDomain()))
418              Y = m * self.__s0 * mu[:numLS]          if self.__w0 is not None:
419                Y = m * self.__w0 * mu
420          else:          else:
421              Y = Data()              if numLS == 1:
422          n+=numLS                Y = Scalar(0,  grad_m.getFunctionSpace())
423                else:
424                  Y = Data(0, (numLS,) , grad_m.getFunctionSpace())
425    
426          if self.__s1 is not None:          if self.__w1 is not None:
427                X=grad_m*self.__w1
428              if numLS == 1:              if numLS == 1:
429                  X=grad_m* self.__s1*mu[n]                  X=grad_m* self.__w1*mu
430              else:              else:
                 X=self.getPDE().createCoefficient("X")  
431                  for k in xrange(numLS):                  for k in xrange(numLS):
432                      X[k,:]=mu[n+k]*grad_m[k,:]*self.__s1[k,:]                      X[k,:]*=mu[k]
433          else:          else:
434              X = Data()              X = Data(0, grad_m.getShape(), grad_m.getFunctionSpace())
435          n+=numLS          
436          if self.__sc is not None:          # cross gradient terms:
437              raise NotImplementedError          if numLS > 1:    
438              for  k in xrange(numLS):
439                 grad_m_k=grad_m[k,:]
440                 l2_grad_m_k = length(grad_m_k)**2
441                 for  l in xrange(k):
442                   grad_m_l=grad_m[l,:]
443                   l2_grad_m_l = length(grad_m_l)**2
444                   grad_m_lk = inner(grad_m_l, grad_m_k)
445                   f=  mu_c[l,k]* self.__wc[l,k]
446                   X[l,:] += f * ( l2_grad_m_l *  grad_m_l - grad_m_lk * grad_m_k )
447                   X[k,:] += f * ( l2_grad_m_k *  grad_m_k - grad_m_lk * grad_m_l )
448            
449          return ArithmeticTuple(Y, X)          return ArithmeticTuple(Y, X)
450    
451      def getArguments(self, m):  
         """  
         """  
         return ( grad(m),)  
452    
453      def getInverseHessianApproximation(self, m, r, grad_m):      def getInverseHessianApproximation(self, m, r, grad_m):
454          """          """
# Line 311  class Regularization(CostFunction): Line 457  class Regularization(CostFunction):
457    
458              self._new_mu=False              self._new_mu=False
459              self._update_Hessian=False              self._update_Hessian=False
460                mu=self.__mu
461                mu_c=self.__mu_c
462    
             mu=self.getWeights( uncompress=True)  
463              DIM=self.getDomain().getDim()              DIM=self.getDomain().getDim()
464              numLS=self.getNumLevelSets()              numLS=self.getNumLevelSets()
465              n=0              if self.__w0 is not None:
             if self.__s0 is not None:  
466                  if numLS == 1:                  if numLS == 1:
467                       D=self.__s0 * mu[n]                       D=self.__w0 * mu
468                  else:                  else:
469                       D=self.getPDE().createCoefficient("D")                       D=self.getPDE().createCoefficient("D")
470                       for k in xrange(numLS): D[k,k]=self.__s0[k] * mu[n+k]                       for k in xrange(numLS): D[k,k]=self.__w0[k] * mu[k]
471                  self.getPDE().setValue(D=D)                  self.getPDE().setValue(D=D)
             n+=numLS  
472    
473              A=self.getPDE().createCoefficient("A")              A=self.getPDE().createCoefficient("A")
474              if self.__s1 is not None:              if self.__w1 is not None:
475                 if numLS == 1:                 if numLS == 1:
476                     for i in xrange(DIM): A[i,i]=self.__s1[i] * mu[n]                     for i in xrange(DIM): A[i,i]=self.__w1[i] * mu
477                 else:                 else:
478                     for k in xrange(numLS):                     for k in xrange(numLS):
479                          for i in xrange(DIM): A[k,i,k,i]=self.__s1[k,i] * mu[n+k]                          for i in xrange(DIM): A[k,i,k,i]=self.__w1[k,i] * mu[k]
480              n+=numLS                          
481              if self.__sc is not None:              if numLS > 1:
482                  raise NotImplementedError                 for  k in xrange(numLS):
483              print A                   grad_m_k=grad_m[k,:]
484                     l2_grad_m_k = length(grad_m_k)**2
485                     o_kk=outer(grad_m_k, grad_m_k)
486                     for  l in xrange(k):
487                        grad_m_l=grad_m[l,:]
488                        l2_grad_m_l = length(grad_m_l)**2
489                        i_lk = inner(grad_m_l, grad_m_k)
490                        o_lk = outer(grad_m_l, grad_m_k)
491                        o_kl = outer(grad_m_k, grad_m_l)
492                        o_ll=outer(grad_m_l, grad_m_l)
493                        f=  mu_c[l,k]* self.__wc[l,k]
494                        
495                        A[l,:,l,:] += f * ( l2_grad_m_k * kronecker(DIM) - o_kk )
496                        A[l,:,k,:] += f * ( 2 * o_lk -   o_kl - i_lk * kronecker(DIM) )
497                        A[k,:,l,:] += f * ( 2 * o_kl -   o_lk - i_lk * kronecker(DIM) )
498                        A[k,:,k,:] += f * ( l2_grad_m_l * kronecker(DIM) - o_ll )
499              self.getPDE().setValue(A=A)              self.getPDE().setValue(A=A)
500          self.getPDE().resetRightHandSideCoefficients()          #self.getPDE().resetRightHandSideCoefficients()
501          self.getPDE().setValue(X=r[1])          #self.getPDE().setValue(X=r[1])
502          print "X only: ",self.getPDE().getSolution()          #print "X only: ",self.getPDE().getSolution()
503          self.getPDE().resetRightHandSideCoefficients()          #self.getPDE().resetRightHandSideCoefficients()
504          self.getPDE().setValue(Y=r[0])          #self.getPDE().setValue(Y=r[0])
505          print "Y only: ",self.getPDE().getSolution()          #print "Y only: ",self.getPDE().getSolution()
506    
507          self.getPDE().resetRightHandSideCoefficients()          self.getPDE().resetRightHandSideCoefficients()
508          self.getPDE().setValue(X=r[1], Y=r[0])          self.getPDE().setValue(X=r[1], Y=r[0])
# Line 354  class Regularization(CostFunction): Line 514  class Regularization(CostFunction):
514          """          """
515          if not self.__useDiagonalHessianApproximation:          if not self.__useDiagonalHessianApproximation:
516              self._update_Hessian=True              self._update_Hessian=True
517                
518     # ================ we should factor these out ==============================      def getNorm(self, m):
   
     def getNumRelevantTerms(self):  
         """  
         returns the number of terms in the cost function of regularization  
         with non-zero scaling factors.  
   
         :rtype: ``int``  
         """  
         return len(self.__weight_index)  
   
     def getNumTerms(self):  
         """  
         returns the number of terms in the cost function of regularization  
         with non-zero scaling factors.  
   
         :rtype: ``int``  
519          """          """
520          return len(self.__weight_index)          returns the norm of ``m``
521    
522      def setWeights(self, mu=None):          :param m: level set function
523          """          :type m: `Data`
524          sets the values for the weighting terms in the cost function.          :rtype: ``float``
         Note that only values corresponding to entries with non-negative  
         scaling factors are used.  
   
         :param mu: weights  
         :type mu: ``list`` of ``float`` or ``np,array``  
         """  
         if mu == None:  
             mu = np.ones(self.getNumRelevantTerms())  
         mu=np.asarray(mu)  
   
         if len(mu) == self.getNumRelevantTerms():  
             if not mu.shape == (self.getNumRelevantTerms(),):  
                 raise ValueError("%s values are expected."%self.getNumRelevantTerms())  
             self.__mu=mu  
         else:  
             if not mu.shape == (self.__total_num_weights,):  
                 raise ValueError("%s values are expected."%self.__total_num_weights)  
             self.__mu = np.zeros(self.getNumRelevantTerms())  
             for i in xrange(len(self.__weight_index)): self.__mu[i]=mu[self.__weight_index[i]]  
         self._new_mu=True  
   
     def setWeightsForS0(self, mu=None):  
          """  
          sets the weights for s0-terms  
          """  
          numLS=self.getNumLevelSets()  
          my_mu=self.getWeights(uncompress=True)  
          if mu is None:  
              my_mu[:numLS]=1  
          else:  
              my_mu[:numLS]=mu  
          self.setWeights(my_mu)  
   
     def setWeightsForS1(self, mu=None):  
          """  
          sets the weights for s1-terms  
          """  
          numLS=self.getNumLevelSets()  
          my_mu=self.getWeights(uncompress=True)  
          if mu is None:  
               my_mu[numLS:2*numLS]=1  
          else:  
               my_mu[numLS:2*numLS]=mu  
          self.setWeights(my_mu)  
   
     def setWeightsForSc(self, mu):  
          """  
          sets the weights for sc-terms  
          """  
          numLS=self.getNumLevelSets()  
          my_mu=self.getWeights(uncompress=True)  
          if mu is None:  
               my_mu[2*numLS:]=1  
          else:  
               my_mu[2*numLS:]=mu  
          self.setWeights(my_mu)  
   
     def getWeights(self, uncompress=False):  
         """  
         Returns the weights for the terms in the cost function.  
         The first ``numLevelSets`` values are used for the regularization terms  
         and the remaining values for the cross correlation terms.  
         """  
         if uncompress:  
             mu = np.zeros(self.__total_num_weights)  
             for i in xrange(len(self.__weight_index)):  
                 mu[self.__weight_index[i]] = self.__mu[i]  
             return mu  
         else:  
             return self.__mu  
   
     def getWeightIndex(self):  
         """  
         returns an index to the contributions of terms with non-zero scaling  
         factor.  
525          """          """
526          return self.__weight_index          return sqrt(integrate(length(m)**2)/self.__vol_d)
   

Legend:
Removed from v.4097  
changed lines
  Added in v.4154

  ViewVC Help
Powered by ViewVC 1.1.26