/[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 4122 by gross, Thu Dec 20 05:42: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, Scalar, grad, inner, integrate, interpolate, kronecker, boundingBoxEdgeLengths, vol, sqrt, length  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 96  class Regularization(CostFunction): Line 96  class Regularization(CostFunction):
96          """          """
97          if w0 == None and w1==None:          if w0 == None and w1==None:
98                raise ValueError("Values for w0 or for w1 must be given.")                raise ValueError("Values for w0 or for w1 must be given.")
99      if wc == None and  numLevelSets>1:          if wc == None and  numLevelSets>1:
100            raise ValueError("Values for wc 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()
# Line 110  class Regularization(CostFunction): Line 110  class Regularization(CostFunction):
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          # =========== check the shape of the scales: =================================          # =========== check the shape of the scales: =================================
115          if scale is None:          if scale is None:
116          if numLevelSets == 1 :              if numLevelSets == 1 :
117             scale = 1.                 scale = 1.
         else:  
            scale = np.ones((numLevelSets,))  
     else:  
         scale=np.asarray(scale)  
         if numLevelSets == 1 :  
             if scale.shape == ():  
            if not scale > 0 :  
               raise ValueError("Value for scale must be positive.")  
         else:  
            raise ValueError("Unexpected shape %s for scale."%scale.shape)  
118              else:              else:
119               if scale.shape is (numLevelSets,):                 scale = np.ones((numLevelSets,))
120               if not min(scale) > 0:          else:
121                  raise ValueError("All value for scale must be positive.")              scale=np.asarray(scale)
122           else:              if numLevelSets == 1 :
123             raise ValueError("Unexpected shape %s for scale."%scale.shape)                  if scale.shape == ():
124                           if not scale > 0 :
125                          raise ValueError("Value for scale must be positive.")
126                    else:
127                       raise ValueError("Unexpected shape %s for scale."%scale.shape)
128                else:
129                     if scale.shape is (numLevelSets,):
130                         if not min(scale) > 0:
131                            raise ValueError("All value for scale must be positive.")
132                     else:
133                       raise ValueError("Unexpected shape %s for scale."%scale.shape)
134            
135          if scale_c is None or numLevelSets < 2:          if scale_c is None or numLevelSets < 2:
136          scale_c = np.ones((numLevelSets,numLevelSets))              scale_c = np.ones((numLevelSets,numLevelSets))
137      else:          else:
138          scale_c=np.asarray(scale_c)              scale_c=np.asarray(scale_c)
139          if scale_c.shape == (numLevelSets,numLevelSets):              if scale_c.shape == (numLevelSets,numLevelSets):
140              if not all( [ [ scale_c[l,k] > 0. for l in xrange(k) ] for k in xrange(1,numLevelSets) ]):                  if not all( [ [ scale_c[l,k] > 0. for l in xrange(k) ] for k in xrange(1,numLevelSets) ]):
141                  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.")
142              else:              else:
143           raise ValueError("Unexpected shape %s for scale."%scale_c.shape)                   raise ValueError("Unexpected shape %s for scale."%scale_c.shape)
144      # ===== check the shape of the weights: ============================================          # ===== check the shape of the weights: ============================================
145          if w0 is not None:          if w0 is not None:
146            w0 = interpolate(w0,self.__pde.getFunctionSpaceForCoefficient('D'))                w0 = interpolate(w0,self.__pde.getFunctionSpaceForCoefficient('D'))
147            s0=w0.getShape()                s0=w0.getShape()
148            if numLevelSets == 1 :                if numLevelSets == 1 :
149             if  not s0 == () :                     if  not s0 == () :
150                raise ValueError("Unexpected shape %s for weight w0."%s0)                        raise ValueError("Unexpected shape %s for weight w0."%s0)
151                else:                else:
152             if not s0 == (numLevelSets,):                     if not s0 == (numLevelSets,):
153                raise ValueError("Unexpected shape %s for weight w0."%s0)                        raise ValueError("Unexpected shape %s for weight w0."%s0)
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()
157            if numLevelSets is 1 :                if numLevelSets is 1 :
158             if not s1 == (DIM,) :                     if not s1 == (DIM,) :
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 :          if numLevelSets == 1 :
164               wc=None               wc=None
165          else:          else:
166               wc = interpolate(wc,self.__pde.getFunctionSpaceForCoefficient('A'))               wc = interpolate(wc,self.__pde.getFunctionSpaceForCoefficient('A'))
167               sc=wc.getShape()               sc=wc.getShape()
168               raise ValueError("Unexpected shape %s for weight wc."%sc)               if not sc == (numLevelSets, numLevelSets):
169                    raise ValueError("Unexpected shape %s for weight wc."%(sc,))
170          # ============= now we rescale weights: ======================================          # ============= now we rescale weights: ======================================
171          L2s=np.asarray(boundingBoxEdgeLengths(domain))**2          L2s=np.asarray(boundingBoxEdgeLengths(domain))**2
172          L4=1/np.sum(1/L2s)**2          L4=1/np.sum(1/L2s)**2
173          if numLevelSets is 1 :          if numLevelSets == 1 :
174              A=0              A=0
175              if w0 is not None:              if w0 is not None:
176              A = integrate(w0)                  A = integrate(w0)
177          if w1 is not None:              if w1 is not None:
178              A += integrate(inner(w1, 1/L2s))                  A += integrate(inner(w1, 1/L2s))
179          if A > 0:              if A > 0:
180              f = scale/A                  f = scale/A
181              if w0 is not None:                  if w0 is not None:
182                   w0*=f                       w0*=f
183              if w1 is not None:                  if w1 is not None:
184               w1*=f                       w1*=f
185              else:              else:
186             raise ValueError("Non-positive weighting factor detected.")                 raise ValueError("Non-positive weighting factor detected.")
187          else:          else:
188           for k in xrange(numLevelSets):  
189               A=0               for k in xrange(numLevelSets):
190                     A=0
191                   if w0 is not None:                   if w0 is not None:
192                   A = integrate(w0[k])                       A = integrate(w0[k])
193               if w1 is not None:                   if w1 is not None:
194                    A += integrate(inner(w1[k,:], 1/L2s))                        A += integrate(inner(w1[k,:], 1/L2s))
195               if A > 0:                   if A > 0:
196                    f = scale[k]/A                        f = scale[k]/A
197                    if w0 is not None:                        if w0 is not None:
198                       w0[k]*=f                           w0[k]*=f
199                    if w1 is not None:                        if w1 is not None:
200                   w1[k,:]*=f                           w1[k,:]*=f
201                   else:                   else:
202                 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)
203                            
204               # and now the cross-gradient:                   # and now the cross-gradient:
205               for l in xrange(k):                     if wc is not None:
206                   A = integrate(wc[l,k])/L4                       for l in xrange(k):  
207                   if A > 0:                          A = integrate(wc[l,k])/L4
208                      f = scale_c[l,k]/A                          if A > 0:
209                      wc[l,k]*=f                             f = scale_c[l,k]/A
210                       else:                             wc[l,k]*=f
211                      raise ValueError("Non-positive weighting factor for cross-gradient level set components %d and %d detected."%(l,k))  #                        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          self.__w0=w0
215          self.__w1=w1          self.__w1=w1
216          self.__wc=wc          self.__wc=wc
           
   
217    
218          self.__pde_is_set=False                  self.__pde_is_set=False        
219          if self.__numLevelSets > 1:          if self.__numLevelSets > 1:
# Line 285  class Regularization(CostFunction): Line 286  class Regularization(CostFunction):
286          numLS=self.getNumLevelSets()          numLS=self.getNumLevelSets()
287          numTF=self.getNumTradeOffFactors()          numTF=self.getNumTradeOffFactors()
288          if mu is None:          if mu is None:
289         mu = np.ones((numTF,))             mu = np.ones((numTF,))
290      else:          else:
291         mu = np.asarray(mu)             mu = np.asarray(mu)
292    
293      if mu.shape == (numTF,):          if mu.shape == (numTF,):
294          self.setTradeOffFactorsForVariation(mu[:numLS])              self.setTradeOffFactorsForVariation(mu[:numLS])
295          mu_c2=np.zeros((numLS,numLS))              mu_c2=np.zeros((numLS,numLS))
296          for k in xrange(numLS):              for k in xrange(numLS):
297             for l in xrange(k):                 for l in xrange(k):
298             mu_c2[l,k] = mu[numLS+l+((k-1)*k)/2]                     mu_c2[l,k] = mu[numLS+l+((k-1)*k)/2]
299          self.setTradeOffFactorsForCrossGradient(mu_c2)              self.setTradeOffFactorsForCrossGradient(mu_c2)
300      elif mu.shape == () and numLS ==1:          elif mu.shape == () and numLS ==1:
301          self.setTradeOffFactorsForVariation(mu)              self.setTradeOffFactorsForVariation(mu)
302      else:          else:
303         raise ValueError("Unexpected shape %s for mu."%(mu.shape,))             raise ValueError("Unexpected shape %s for mu."%(mu.shape,))
304                    
305      def setTradeOffFactorsForVariation(self, mu=None):      def setTradeOffFactorsForVariation(self, mu=None):
306           """           """
307           sets the trade-off factors for the level-set variation part           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.           :param mu:  new values for the trade-off factors. Values must be positive.
310           :type mu: `float``, ``list`` of ``float`` or ```numpy.array```           :type mu: ``float``, ``list`` of ``float`` or ```numpy.array```
311           """           """
312           numLS=self.getNumLevelSets()           numLS=self.getNumLevelSets()
313           if mu is None:           if mu is None:
314          if numLS == 1:              if numLS == 1:
315             mu = 1.                 mu = 1.
316          else:              else:
317             mu = np.ones((numLS,))                 mu = np.ones((numLS,))
318    
319       mu=np.asarray(mu)           mu=np.asarray(mu)
320       if numLS == 1:           if numLS == 1:
321         if mu.shape == (1,): mu=mu[0]             if mu.shape == (1,): mu=mu[0]
322         if mu.shape == ():             if mu.shape == ():
323            if mu > 0:                if mu > 0:
324           self.__mu= mu                   self.__mu= mu
325           self._new_mu=True                   self._new_mu=True
326            else:                else:
327           raise ValueError("Value for trade-off factor must be positive.")                   raise ValueError("Value for trade-off factor must be positive.")
328         else:             else:
329            raise ValueError("Unexpected shape %s for mu."%mu.shape)                raise ValueError("Unexpected shape %s for mu."%mu.shape)
330       else:           else:
331         if mu.shape == (numLS,):             if mu.shape == (numLS,):
332             if min(mu) > 0:                 if min(mu) > 0:
333             self.__mu= mu                     self.__mu= mu
334             self._new_mu=True                     self._new_mu=True
335             else:                 else:
336             raise ValueError("All value for mu must be positive.")                     raise ValueError("All value for mu must be positive.")
337             else:             else:
338             raise ValueError("Unexpected shape %s for trade-off factor."%mu.shape)                 raise ValueError("Unexpected shape %s for trade-off factor."%mu.shape)
339    
340      def setTradeOffFactorsForCrossGradient(self, mu_c=None):      def setTradeOffFactorsForCrossGradient(self, mu_c=None):
341          """          """
342          sets the trade-off factors for the cross--gradient terms          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 terms. Values must be positive.          :param mu_c: new values for the trade-off factors for the cross-gradient
345                        if now value is given ones are used. Onky value mu_c[l,k] for l<k are used.                       terms. Values must be positive. If no value is given ones
346          :type mu: `float``, ``list`` of ``float`` or ```numpy.array```                       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()          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      else:          if isinstance(mu_c, float) or isinstance(mu_c, int):
354          mu_c=np.asarray(mu_c)              self.__mu_c = np.zeros((numLS,numLS))
355          if mu_c.shape == (numLS,numLS):              self.__mu_c[:,:]=mu_c
356              if not all( [ [ mu_c[l,k] > 0. for l in xrange(k) ] for k in xrange(1,numLS) ]):          else:
357               raise ValueError("All trade-off factors in the lower triangle of mu_c must be positive.")              mu_c=np.asarray(mu_c)
358          else:              if mu_c.shape == (numLS,numLS):
359               self.__mu_c =  mu_c                  if not all( [ [ mu_c[l,k] > 0. for l in xrange(k) ] for k in xrange(1,numLS) ]):
360               self._new_mu=True                       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:              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    
373      def getValue(self, m, grad_m):      def getValue(self, m, grad_m):
374          """          """
# Line 407  class Regularization(CostFunction): Line 412  class Regularization(CostFunction):
412          DIM=self.getDomain().getDim()          DIM=self.getDomain().getDim()
413          numLS=self.getNumLevelSets()          numLS=self.getNumLevelSets()
414    
415            print "WARNING: WRONG FUNCTION SPACE"
416            
417            grad_m=grad(m, ReducedFunction(m.getDomain()))
418          if self.__w0 is not None:          if self.__w0 is not None:
419              Y = m * self.__w0 * mu              Y = m * self.__w0 * mu
420          else:          else:
421          if numLS == 1:              if numLS == 1:
422            Y = Scalar(0,  grad_m.getFunctionSpace())                Y = Scalar(0,  grad_m.getFunctionSpace())
423          else:              else:
424                Y = Data(0, (numLS,) , grad_m.getFunctionSpace())                Y = Data(0, (numLS,) , grad_m.getFunctionSpace())
425    
426          if self.__w1 is not None:          if self.__w1 is not None:
427          X=grad_m*self.__w1              X=grad_m*self.__w1
428              if numLS == 1:              if numLS == 1:
429                  X=grad_m* self.__w1*mu                  X=grad_m* self.__w1*mu
430              else:              else:
# Line 426  class Regularization(CostFunction): Line 434  class Regularization(CostFunction):
434              X = Data(0, grad_m.getShape(), grad_m.getFunctionSpace())              X = Data(0, grad_m.getShape(), grad_m.getFunctionSpace())
435                    
436          # cross gradient terms:          # cross gradient terms:
437          if numLS > 1:              if numLS > 1:    
438        for  k in xrange(numLS):            for  k in xrange(numLS):
439               grad_m_k=grad_m[k,:]               grad_m_k=grad_m[k,:]
440           l2_grad_m_k = length(grad_m_k)**2               l2_grad_m_k = length(grad_m_k)**2
441           for  l in xrange(k):               for  l in xrange(k):
442             grad_m_l=grad_m[l,:]                 grad_m_l=grad_m[l,:]
443             l2_grad_m_l = length(grad_m_l)**2                 l2_grad_m_l = length(grad_m_l)**2
444             grad_m_lk = inner(grad_m_l, grad_m_k)                 grad_m_lk = inner(grad_m_l, grad_m_k)
445             f=  mu_c[l,k]* self.__wc[l,k]                 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 )                 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 )                 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    
# Line 471  class Regularization(CostFunction): Line 479  class Regularization(CostFunction):
479                          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]
480                                                    
481              if numLS > 1:              if numLS > 1:
482             for  k in xrange(numLS):                 for  k in xrange(numLS):
483                   grad_m_k=grad_m[k,:]                   grad_m_k=grad_m[k,:]
484               l2_grad_m_k = length(grad_m_k)**2                   l2_grad_m_k = length(grad_m_k)**2
485               o_kk=outer(grad_m_k, grad_m_k)                   o_kk=outer(grad_m_k, grad_m_k)
486               for  l in xrange(k):                   for  l in xrange(k):
487                  grad_m_l=grad_m[l,:]                      grad_m_l=grad_m[l,:]
488                  l2_grad_m_l = length(grad_m_l)**2                      l2_grad_m_l = length(grad_m_l)**2
489                  i_lk = inner(grad_m_l, grad_m_k)                      i_lk = inner(grad_m_l, grad_m_k)
490                  o_lk = outer(grad_m_l, grad_m_k)                      o_lk = outer(grad_m_l, grad_m_k)
491                  o_kl = outer(grad_m_k, grad_m_l)                      o_kl = outer(grad_m_k, grad_m_l)
492                  o_ll=outer(grad_m_l, grad_m_l)                      o_ll=outer(grad_m_l, grad_m_l)
493                  f=  mu_c[l,k]* self.__wc[l,k]                      f=  mu_c[l,k]* self.__wc[l,k]
494                                        
495                  A[l,:,l:] += f * ( l2_grad_m_k * kronecker(DIM) - o_kk )                      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) )                      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) )                      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 )                      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])
# Line 515  class Regularization(CostFunction): Line 523  class Regularization(CostFunction):
523          :type m: `Data`          :type m: `Data`
524          :rtype: ``float``          :rtype: ``float``
525          """          """
         return sqrt(integrate(length(m)**2)/self.__vol_d)  
526            return sqrt(integrate(length(m)**2)/self.__vol_d)

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

  ViewVC Help
Powered by ViewVC 1.1.26