/[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 4154 by jfenwick, Tue Jan 22 09:30:23 2013 UTC revision 4213 by caltinay, Tue Feb 19 01:16:29 2013 UTC
# 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 ReducedFunction, outer, Data, Scalar, grad, inner, integrate, interpolate, kronecker, boundingBoxEdgeLengths, vol, sqrt, length  from esys.escript import Function, 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 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.          :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          :type scale: ``Scalar`` if ``numLevelSets`` == 1 or `Data` object of shape
88                    (``numLevelSets`` ,) if ``numLevelSets`` > 1                    (``numLevelSets`` ,) if ``numLevelSets`` > 1
89          :param scale_c: scale for the cross gradient terms. If not set          :param scale_c: scale for the cross gradient terms. If not set
# Line 91  class Regularization(CostFunction): Line 91  class Regularization(CostFunction):
91                     only. Only values ``scale_c[l,k]`` in the lower triangle (l<k)                     only. Only values ``scale_c[l,k]`` in the lower triangle (l<k)
92                     are used.                     are used.
93          :type scale_c: `Data` object of shape (``numLevelSets`` , ``numLevelSets``)          :type scale_c: `Data` object of shape (``numLevelSets`` , ``numLevelSets``)
           
94    
95          """          """
96          if w0 == None and w1==None:          if w0 == None and w1==None:
# Line 102  class Regularization(CostFunction): Line 101  class Regularization(CostFunction):
101          self.__domain=domain          self.__domain=domain
102          DIM=self.__domain.getDim()          DIM=self.__domain.getDim()
103          self.__numLevelSets=numLevelSets          self.__numLevelSets=numLevelSets
104            
105          self.__pde=LinearPDE(self.__domain, numEquations=self.__numLevelSets)          self.__pde=LinearPDE(self.__domain, numEquations=self.__numLevelSets)
106          self.__pde.getSolverOptions().setTolerance(tol)          self.__pde.getSolverOptions().setTolerance(tol)
107          self.__pde.setSymmetryOn()          self.__pde.setSymmetryOn()
# Line 110  class Regularization(CostFunction): Line 109  class Regularization(CostFunction):
109              self.__pde.setValue(q=location_of_set_m)              self.__pde.setValue(q=location_of_set_m)
110          except IllegalCoefficientValue:          except IllegalCoefficientValue:
111              raise ValueError("Unable to set location of fixed level set function.")              raise ValueError("Unable to set location of fixed level set function.")
112              
113          # =========== check the shape of the scales: =================================          # =========== check the shape of the scales: ========================
114          if scale is None:          if scale is None:
115              if numLevelSets == 1 :              if numLevelSets == 1 :
116                 scale = 1.                 scale = 1.
# Line 131  class Regularization(CostFunction): Line 130  class Regularization(CostFunction):
130                          raise ValueError("All value for scale must be positive.")                          raise ValueError("All value for scale must be positive.")
131                   else:                   else:
132                     raise ValueError("Unexpected shape %s for scale."%scale.shape)                     raise ValueError("Unexpected shape %s for scale."%scale.shape)
133            
134          if scale_c is None or numLevelSets < 2:          if scale_c is None or numLevelSets < 2:
135              scale_c = np.ones((numLevelSets,numLevelSets))              scale_c = np.ones((numLevelSets,numLevelSets))
136          else:          else:
# Line 141  class Regularization(CostFunction): Line 140  class Regularization(CostFunction):
140                          raise ValueError("All values in the lower triangle of scale_c must be positive.")                          raise ValueError("All values in the lower triangle of scale_c must be positive.")
141              else:              else:
142                   raise ValueError("Unexpected shape %s for scale."%scale_c.shape)                   raise ValueError("Unexpected shape %s for scale."%scale_c.shape)
143          # ===== check the shape of the weights: ============================================          # ===== check the shape of the weights: =============================
144          if w0 is not None:          if w0 is not None:
145                w0 = interpolate(w0,self.__pde.getFunctionSpaceForCoefficient('D'))                w0 = interpolate(w0,self.__pde.getFunctionSpaceForCoefficient('D'))
146                s0=w0.getShape()                s0=w0.getShape()
# Line 150  class Regularization(CostFunction): Line 149  class Regularization(CostFunction):
149                        raise ValueError("Unexpected shape %s for weight w0."%s0)                        raise ValueError("Unexpected shape %s for weight w0."%s0)
150                else:                else:
151                     if not s0 == (numLevelSets,):                     if not s0 == (numLevelSets,):
152                        raise ValueError("Unexpected shape %s for weight w0."%s0)                        raise ValueError("Unexpected shape %s for weight w0."%s0)
153    
154          if not w1 is None:          if not w1 is None:
155                w1 = interpolate(w1,self.__pde.getFunctionSpaceForCoefficient('A'))                w1 = interpolate(w1,self.__pde.getFunctionSpaceForCoefficient('A'))
156                s1=w1.getShape()                s1=w1.getShape()
# Line 159  class Regularization(CostFunction): Line 159  class Regularization(CostFunction):
159                        raise ValueError("Unexpected shape %s for weight w1."%s1)                        raise ValueError("Unexpected shape %s for weight w1."%s1)
160                else:                else:
161                     if not s1 == (numLevelSets,DIM):                     if not s1 == (numLevelSets,DIM):
162                        raise ValueError("Unexpected shape %s for weight w1."%s1)                        raise ValueError("Unexpected shape %s for weight w1."%s1)
163          if numLevelSets == 1 :  
164            if numLevelSets == 1:
165               wc=None               wc=None
166          else:          else:
167               wc = interpolate(wc,self.__pde.getFunctionSpaceForCoefficient('A'))               wc = interpolate(wc,self.__pde.getFunctionSpaceForCoefficient('A'))
168               sc=wc.getShape()               sc=wc.getShape()
169               if not sc == (numLevelSets, numLevelSets):               if not sc == (numLevelSets, numLevelSets):
170                  raise ValueError("Unexpected shape %s for weight wc."%(sc,))                  raise ValueError("Unexpected shape %s for weight wc."%(sc,))
171          # ============= now we rescale weights: ======================================          # ============= now we rescale weights: =============================
172          L2s=np.asarray(boundingBoxEdgeLengths(domain))**2          L2s=np.asarray(boundingBoxEdgeLengths(domain))**2
173          L4=1/np.sum(1/L2s)**2          L4=1/np.sum(1/L2s)**2
174          if numLevelSets == 1 :          if numLevelSets == 1:
175              A=0              A=0
176              if w0 is not None:              if w0 is not None:
177                  A = integrate(w0)                  A = integrate(w0)
# Line 183  class Regularization(CostFunction): Line 184  class Regularization(CostFunction):
184                  if w1 is not None:                  if w1 is not None:
185                       w1*=f                       w1*=f
186              else:              else:
187                 raise ValueError("Non-positive weighting factor detected.")                 raise ValueError("Non-positive weighting factor detected.")
188          else:          else:
189    
190               for k in xrange(numLevelSets):               for k in xrange(numLevelSets):
# Line 199  class Regularization(CostFunction): Line 200  class Regularization(CostFunction):
200                        if w1 is not None:                        if w1 is not None:
201                           w1[k,:]*=f                           w1[k,:]*=f
202                   else:                   else:
203                     raise ValueError("Non-positive weighting factor for level set component %d detected."%k)                     raise ValueError("Non-positive weighting factor for level set component %d detected."%k)
204                    
205                   # and now the cross-gradient:                   # and now the cross-gradient:
206                   if wc is not None:                   if wc is not None:
207                       for l in xrange(k):                         for l in xrange(k):
208                          A = integrate(wc[l,k])/L4                          A = integrate(wc[l,k])/L4
209                          if A > 0:                          if A > 0:
210                             f = scale_c[l,k]/A                             f = scale_c[l,k]/A
211                             wc[l,k]*=f                             wc[l,k]*=f
212  #                        else:  #                        else:
213  #                          raise ValueError("Non-positive weighting factor for cross-gradient level set components %d and %d detected."%(l,k))  #                          raise ValueError("Non-positive weighting factor for cross-gradient level set components %d and %d detected."%(l,k))
214                        
215          self.__w0=w0          self.__w0=w0
216          self.__w1=w1          self.__w1=w1
217          self.__wc=wc          self.__wc=wc
218    
219          self.__pde_is_set=False                  self.__pde_is_set=False
220          if self.__numLevelSets > 1:          if self.__numLevelSets > 1:
221              self.__useDiagonalHessianApproximation=useDiagonalHessianApproximation              self.__useDiagonalHessianApproximation=useDiagonalHessianApproximation
222          else:          else:
# Line 225  class Regularization(CostFunction): Line 226  class Regularization(CostFunction):
226          self.__num_tradeoff_factors=numLevelSets+((numLevelSets-1)*numLevelSets)/2          self.__num_tradeoff_factors=numLevelSets+((numLevelSets-1)*numLevelSets)/2
227          self.setTradeOffFactors()          self.setTradeOffFactors()
228          self.__vol_d=vol(self.__domain)          self.__vol_d=vol(self.__domain)
229            
230      def getDomain(self):      def getDomain(self):
231          """          """
232          returns the domain of the regularization term          returns the domain of the regularization term
# Line 249  class Regularization(CostFunction): Line 250  class Regularization(CostFunction):
250          :rtype: `LinearPDE`          :rtype: `LinearPDE`
251          """          """
252          return self.__pde          return self.__pde
253        
254      def getDualProduct(self, m, r):      def getDualProduct(self, m, r):
255          """          """
256          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 276  class Regularization(CostFunction): Line 277  class Regularization(CostFunction):
277      def setTradeOffFactors(self, mu=None):      def setTradeOffFactors(self, mu=None):
278          """          """
279          sets the trade-off factors for the level-set variation and the cross-gradient          sets the trade-off factors for the level-set variation and the cross-gradient
280            
281          :param mu: new values for the trade-off factors where values mu[:numLevelSets] are the          :param mu: new values for the trade-off factors where values mu[:numLevelSets] are the
282                     trade-off factors for the level-set variation and the remaining values for                     trade-off factors for the level-set variation and the remaining values for
283                     the cross-gradient part with mu_c[l,k]=mu[numLevelSets+l+((k-1)*k)/2] (l<k).                     the cross-gradient part with mu_c[l,k]=mu[numLevelSets+l+((k-1)*k)/2] (l<k).
284                     If no values for mu is given ones are used. Values must be positive.                     If no values for mu is given ones are used. Values must be positive.
285          :type mu: ``list`` of ``float`` or ```numpy.array```          :type mu: ``list`` of ``float`` or ```numpy.array```
# Line 295  class Regularization(CostFunction): Line 296  class Regularization(CostFunction):
296              mu_c2=np.zeros((numLS,numLS))              mu_c2=np.zeros((numLS,numLS))
297              for k in xrange(numLS):              for k in xrange(numLS):
298                 for l in xrange(k):                 for l in xrange(k):
299                     mu_c2[l,k] = mu[numLS+l+((k-1)*k)/2]                     mu_c2[l,k] = mu[numLS+l+((k-1)*k)/2]
300              self.setTradeOffFactorsForCrossGradient(mu_c2)              self.setTradeOffFactorsForCrossGradient(mu_c2)
301          elif mu.shape == () and numLS ==1:          elif mu.shape == () and numLS ==1:
302              self.setTradeOffFactorsForVariation(mu)              self.setTradeOffFactorsForVariation(mu)
303          else:          else:
304             raise ValueError("Unexpected shape %s for mu."%(mu.shape,))             raise ValueError("Unexpected shape %s for mu."%(mu.shape,))
305              
306      def setTradeOffFactorsForVariation(self, mu=None):      def setTradeOffFactorsForVariation(self, mu=None):
307           """           """
308           sets the trade-off factors for the level-set variation part           sets the trade-off factors for the level-set variation part
309            
310           :param mu:  new values for the trade-off factors. Values must be positive.           :param mu:  new values for the trade-off factors. Values must be positive.
311           :type mu: ``float``, ``list`` of ``float`` or ```numpy.array```           :type mu: ``float``, ``list`` of ``float`` or ```numpy.array```
312           """           """
# Line 318  class Regularization(CostFunction): Line 319  class Regularization(CostFunction):
319    
320           mu=np.asarray(mu)           mu=np.asarray(mu)
321           if numLS == 1:           if numLS == 1:
322             if mu.shape == (1,): mu=mu[0]             if mu.shape == (1,): mu=mu[0]
323             if mu.shape == ():             if mu.shape == ():
324                if mu > 0:                if mu > 0:
325                   self.__mu= mu                   self.__mu= mu
326                   self._new_mu=True                   self._new_mu=True
327                else:                else:
328                   raise ValueError("Value for trade-off factor must be positive.")                   raise ValueError("Value for trade-off factor must be positive.")
329             else:             else:
330                raise ValueError("Unexpected shape %s for mu."%mu.shape)                raise ValueError("Unexpected shape %s for mu."%mu.shape)
331           else:           else:
# Line 335  class Regularization(CostFunction): Line 336  class Regularization(CostFunction):
336                 else:                 else:
337                     raise ValueError("All value for mu must be positive.")                     raise ValueError("All value for mu must be positive.")
338             else:             else:
339                 raise ValueError("Unexpected shape %s for trade-off factor."%mu.shape)                 raise ValueError("Unexpected shape %s for trade-off factor."%mu.shape)
340    
341      def setTradeOffFactorsForCrossGradient(self, mu_c=None):      def setTradeOffFactorsForCrossGradient(self, mu_c=None):
342          """          """
343          sets the trade-off factors for the cross-gradient terms          sets the trade-off factors for the cross-gradient terms
344            
345          :param mu_c: new values for the trade-off factors for the cross-gradient          :param mu_c: new values for the trade-off factors for the cross-gradient
346                       terms. Values must be positive. If no value is given ones                       terms. Values must be positive. If no value is given ones
347                       are used. Onky value mu_c[l,k] for l<k are used.                       are used. Onky value mu_c[l,k] for l<k are used.
348          :type mu_c: ``float``, ``list`` of ``float`` or ``numpy.array``          :type mu_c: ``float``, ``list`` of ``float`` or ``numpy.array``
           
349          """          """
350          numLS=self.getNumLevelSets()          numLS=self.getNumLevelSets()
351          if mu_c is None or numLS < 2:          if mu_c is None or numLS < 2:
352              self.__mu_c = np.ones((numLS,numLS))              self.__mu_c = np.ones((numLS,numLS))
353          if isinstance(mu_c, float) or isinstance(mu_c, int):          if isinstance(mu_c, float) or isinstance(mu_c, int):
354              self.__mu_c = np.zeros((numLS,numLS))              self.__mu_c = np.zeros((numLS,numLS))
355              self.__mu_c[:,:]=mu_c              self.__mu_c[:,:]=mu_c
356          else:          else:
357              mu_c=np.asarray(mu_c)              mu_c=np.asarray(mu_c)
358              if mu_c.shape == (numLS,numLS):              if mu_c.shape == (numLS,numLS):
# Line 363  class Regularization(CostFunction): Line 363  class Regularization(CostFunction):
363                       self._new_mu=True                       self._new_mu=True
364              else:              else:
365                   raise ValueError("Unexpected shape %s for mu."%mu_c.shape)                   raise ValueError("Unexpected shape %s for mu."%mu_c.shape)
366        
367      def getArguments(self, m):      def getArguments(self, m):
368          """          """
369          """          """
370          return ( grad(m),)          return ( grad(m),)
                   
371    
372      def getValue(self, m, grad_m):      def getValue(self, m, grad_m):
373          """          """
# Line 391  class Regularization(CostFunction): Line 390  class Regularization(CostFunction):
390              else:              else:
391                  for k in xrange(numLS):                  for k in xrange(numLS):
392                      A+=mu[k]*integrate(inner(grad_m[k,:]**2,self.__w1[k,:]))                      A+=mu[k]*integrate(inner(grad_m[k,:]**2,self.__w1[k,:]))
393                        
394          if numLS > 1:          if numLS > 1:
395              for k in xrange(numLS):              for k in xrange(numLS):
396                  gk=grad_m[k,:]                  gk=grad_m[k,:]
# Line 412  class Regularization(CostFunction): Line 411  class Regularization(CostFunction):
411          DIM=self.getDomain().getDim()          DIM=self.getDomain().getDim()
412          numLS=self.getNumLevelSets()          numLS=self.getNumLevelSets()
413    
414          print "WARNING: WRONG FUNCTION SPACE"          grad_m=grad(m, Function(m.getDomain()))
           
         grad_m=grad(m, ReducedFunction(m.getDomain()))  
415          if self.__w0 is not None:          if self.__w0 is not None:
416              Y = m * self.__w0 * mu              Y = m * self.__w0 * mu
417          else:          else:
# Line 432  class Regularization(CostFunction): Line 429  class Regularization(CostFunction):
429                      X[k,:]*=mu[k]                      X[k,:]*=mu[k]
430          else:          else:
431              X = Data(0, grad_m.getShape(), grad_m.getFunctionSpace())              X = Data(0, grad_m.getShape(), grad_m.getFunctionSpace())
432            
433          # cross gradient terms:          # cross gradient terms:
434          if numLS > 1:              if numLS > 1:
435            for  k in xrange(numLS):            for  k in xrange(numLS):
436               grad_m_k=grad_m[k,:]               grad_m_k=grad_m[k,:]
437               l2_grad_m_k = length(grad_m_k)**2               l2_grad_m_k = length(grad_m_k)**2
# Line 445  class Regularization(CostFunction): Line 442  class Regularization(CostFunction):
442                 f=  mu_c[l,k]* self.__wc[l,k]                 f=  mu_c[l,k]* self.__wc[l,k]
443                 X[l,:] += f * ( l2_grad_m_l *  grad_m_l - grad_m_lk * grad_m_k )                 X[l,:] += f * ( l2_grad_m_l *  grad_m_l - grad_m_lk * grad_m_k )
444                 X[k,:] += f * ( l2_grad_m_k *  grad_m_k - grad_m_lk * grad_m_l )                 X[k,:] += f * ( l2_grad_m_k *  grad_m_k - grad_m_lk * grad_m_l )
           
         return ArithmeticTuple(Y, X)  
   
445    
446            return ArithmeticTuple(Y, X)
447    
448      def getInverseHessianApproximation(self, m, r, grad_m):      def getInverseHessianApproximation(self, m, r, grad_m):
449          """          """
# Line 477  class Regularization(CostFunction): Line 472  class Regularization(CostFunction):
472                 else:                 else:
473                     for k in xrange(numLS):                     for k in xrange(numLS):
474                          for i in xrange(DIM): A[k,i,k,i]=self.__w1[k,i] * mu[k]                          for i in xrange(DIM): A[k,i,k,i]=self.__w1[k,i] * mu[k]
475                            
476              if numLS > 1:              if numLS > 1:
477                 for  k in xrange(numLS):                 for  k in xrange(numLS):
478                   grad_m_k=grad_m[k,:]                   grad_m_k=grad_m[k,:]
# Line 491  class Regularization(CostFunction): Line 486  class Regularization(CostFunction):
486                      o_kl = outer(grad_m_k, grad_m_l)                      o_kl = outer(grad_m_k, grad_m_l)
487                      o_ll=outer(grad_m_l, grad_m_l)                      o_ll=outer(grad_m_l, grad_m_l)
488                      f=  mu_c[l,k]* self.__wc[l,k]                      f=  mu_c[l,k]* self.__wc[l,k]
489                        
490                      A[l,:,l,:] += f * ( l2_grad_m_k * kronecker(DIM) - o_kk )                      A[l,:,l,:] += f * ( l2_grad_m_k * kronecker(DIM) - o_kk )
491                      A[l,:,k,:] += f * ( 2 * o_lk -   o_kl - i_lk * kronecker(DIM) )                      A[l,:,k,:] += f * ( 2 * o_lk -   o_kl - i_lk * kronecker(DIM) )
492                      A[k,:,l,:] += f * ( 2 * o_kl -   o_lk - i_lk * kronecker(DIM) )                      A[k,:,l,:] += f * ( 2 * o_kl -   o_lk - i_lk * kronecker(DIM) )
# Line 514  class Regularization(CostFunction): Line 509  class Regularization(CostFunction):
509          """          """
510          if not self.__useDiagonalHessianApproximation:          if not self.__useDiagonalHessianApproximation:
511              self._update_Hessian=True              self._update_Hessian=True
512                
513      def getNorm(self, m):      def getNorm(self, m):
514          """          """
515          returns the norm of ``m``          returns the norm of ``m``
# Line 524  class Regularization(CostFunction): Line 519  class Regularization(CostFunction):
519          :rtype: ``float``          :rtype: ``float``
520          """          """
521          return sqrt(integrate(length(m)**2)/self.__vol_d)          return sqrt(integrate(length(m)**2)/self.__vol_d)
522    

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

  ViewVC Help
Powered by ViewVC 1.1.26