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

revision 4099 by gross, Tue Dec 11 04:04:47 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  #  #
# Line 13  Line 13
13  #  #
14  ##############################################################################  ##############################################################################
15
17  http://www.uq.edu.au  http://www.uq.edu.au
# Line 25  __all__ = ['Regularization'] Line 25  __all__ = ['Regularization']
25  from costfunctions import CostFunction  from costfunctions import CostFunction
26
27  import numpy as np  import numpy as np
28  from esys.escript import Data, grad, inner, integrate, interpolate, kronecker, boundingBoxEdgeLengths, vol  from esys.escript import ReducedFunction, outer, Data, Scalar, grad, inner, integrate, interpolate, kronecker, boundingBoxEdgeLengths, vol, sqrt, length
29  from esys.escript.linearPDEs import LinearPDE, IllegalCoefficientValue  from esys.escript.linearPDEs import LinearPDE, IllegalCoefficientValue
30  from esys.escript.pdetools import ArithmeticTuple  from esys.escript.pdetools import ArithmeticTuple
31
# Line 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: ======================================
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*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]*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
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]*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
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()
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,):
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]
300      elif mu.shape == () and numLS ==1:          elif mu.shape == () and numLS ==1:
302      else:          else:
303         raise ValueError("Unexpected shape %s for mu."%(mu.shape,))             raise ValueError("Unexpected shape %s for mu."%(mu.shape,))
304
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
341          """          """
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          """          """
371
372
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
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:
428              if numLS == 1:              if numLS == 1:
430              else:              else:
431                  for k in xrange(numLS):                  for k in xrange(numLS):
432                      X[k,:]*=mu[k]                      X[k,:]*=mu[k]
433          else:          else:
435          if numLS > 1:
436              raise NotImplementedError          # cross gradient terms:
437            if numLS > 1:
438              for  k in xrange(numLS):
441                 for  l in xrange(k):
445                   f=  mu_c[l,k]* self.__wc[l,k]
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):
486                     for  l in xrange(k):
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 484  class Regularization(CostFunction): Line 514  class Regularization(CostFunction):
514          """          """
515          if not self.__useDiagonalHessianApproximation:          if not self.__useDiagonalHessianApproximation:
516              self._update_Hessian=True              self._update_Hessian=True
517
518        def getNorm(self, m):
519            """
520            returns the norm of ``m``
521
522            :param m: level set function
523            :type m: `Data`
524            :rtype: ``float``
525            """
526            return sqrt(integrate(length(m)**2)/self.__vol_d)

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