/[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 4100 by gross, Tue Dec 11 06:55:20 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, 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 94  class Regularization(CostFunction): Line 94  class Regularization(CostFunction):
94                    
95    
96          """          """
         if numLevelSets>1:  
               raise ValueError("Currently only numLevelSets<=1 is supported.")  
           
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 113  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            print s1, (DIM,)                if numLevelSets is 1 :
158            print w1                     if not s1 == (DIM,) :
159            if numLevelSets is 1 :                        raise ValueError("Unexpected shape %s for weight w1."%s1)
            if not s1 == (DIM,) :  
               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: ======================================
         self.__vol_d=vol(self.__domain)  
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*self.__vol_d/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]*self.__vol_d/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]*self.__vol_d/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 229  class Regularization(CostFunction): Line 224  class Regularization(CostFunction):
224    
225          self.__num_tradeoff_factors=numLevelSets+((numLevelSets-1)*numLevelSets)/2          self.__num_tradeoff_factors=numLevelSets+((numLevelSets-1)*numLevelSets)/2
226          self.setTradeOffFactors()          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 253  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 290  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.")
        else:  
           raise ValueError("Unexpected shape %s for mu."%mu.shape)  
      else:  
        if mu.shape == (numLS,):  
            if min(mu) > 0:  
            self.__mu= mu  
            self._new_mu=True  
            else:  
            raise ValueError("All value for mu must be positive.")  
328             else:             else:
329             raise ValueError("Unexpected shape %s for trade-off factor."%mu.shape)                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):      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 412  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          n=0          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              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.__w1 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.__w1*mu                  X=grad_m* self.__w1*mu
430              else:              else:
                 X=grad_m[k,:]*self.__w1[k,:]  
431                  for k in xrange(numLS):                  for k in xrange(numLS):
432                      X[k,:]*=mu[k]                      X[k,:]*=mu[k]
433          else:          else:
434              X = Data()              X = Data(0, grad_m.getShape(), grad_m.getFunctionSpace())
435          if numLS > 1:          
436              raise NotImplementedError          # cross gradient terms:
437            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    
# Line 464  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                  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 493  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.4100  
changed lines
  Added in v.4154

  ViewVC Help
Powered by ViewVC 1.1.26