/[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 4211 by caltinay, Mon Feb 18 23:54:46 2013 UTC revision 4213 by caltinay, Tue Feb 19 01:16:29 2013 UTC
# 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 430  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 443  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 475  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 489  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 512  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 522  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.4211  
changed lines
  Added in v.4213

  ViewVC Help
Powered by ViewVC 1.1.26