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

revision 4123 by gross, Thu Dec 20 05:42:35 2012 UTC revision 4124 by gross, Fri Dec 21 01:35:43 2012 UTC
# Line 25  __all__ = ['Regularization'] Line 25  __all__ = ['Regularization']
25  from costfunctions import CostFunction  from costfunctions import CostFunction
26
27  import numpy as np  import numpy as np
28  from esys.escript import 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 165  class Regularization(CostFunction): Line 165  class Regularization(CostFunction):
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)
# Line 184  class Regularization(CostFunction): Line 185  class Regularization(CostFunction):
185              else:              else:
186             raise ValueError("Non-positive weighting factor detected.")             raise ValueError("Non-positive weighting factor detected.")
187          else:          else:
188
189           for k in xrange(numLevelSets):           for k in xrange(numLevelSets):
190               A=0               A=0
191                   if w0 is not None:                   if w0 is not None:
# Line 200  class Regularization(CostFunction): Line 202  class Regularization(CostFunction):
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]/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 348  class Regularization(CostFunction): Line 349  class Regularization(CostFunction):
349          numLS=self.getNumLevelSets()          numLS=self.getNumLevelSets()
350          if mu_c is None or numLS < 2:          if mu_c is None or numLS < 2:
351          self.__mu_c = np.ones((numLS,numLS))          self.__mu_c = np.ones((numLS,numLS))
352        if isinstance(mu_c, float) or isinstance(mu_c, int):
353            self.__mu_c = np.zeros((numLS,numLS))
354            self.__mu_c[:,:]=mu_c
355      else:      else:
356          mu_c=np.asarray(mu_c)          mu_c=np.asarray(mu_c)
357          if mu_c.shape == (numLS,numLS):          if mu_c.shape == (numLS,numLS):
# Line 407  class Regularization(CostFunction): Line 411  class Regularization(CostFunction):
411          DIM=self.getDomain().getDim()          DIM=self.getDomain().getDim()
412          numLS=self.getNumLevelSets()          numLS=self.getNumLevelSets()
413
414            print "WARNING: WRONG FUNCTION SPACE"
415
417          if self.__w0 is not None:          if self.__w0 is not None:
418              Y = m * self.__w0 * mu              Y = m * self.__w0 * mu
419          else:          else:
# Line 484  class Regularization(CostFunction): Line 491  class Regularization(CostFunction):
492                  f=  mu_c[l,k]* self.__wc[l,k]                  f=  mu_c[l,k]* self.__wc[l,k]
493
494                  A[l,:,l:] += f * ( l2_grad_m_k * kronecker(DIM) - o_kk )                  A[l,:,l,:] += f * ( l2_grad_m_k * kronecker(DIM) - o_kk )
495                  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) )
496                  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) )
497                  A[k,:,k:] += f * ( l2_grad_m_l * kronecker(DIM) - o_ll )                  A[k,:,k,:] += f * ( l2_grad_m_l * kronecker(DIM) - o_ll )
498              self.getPDE().setValue(A=A)              self.getPDE().setValue(A=A)
499          #self.getPDE().resetRightHandSideCoefficients()          #self.getPDE().resetRightHandSideCoefficients()
500          #self.getPDE().setValue(X=r[1])          #self.getPDE().setValue(X=r[1])

Legend:
 Removed from v.4123 changed lines Added in v.4124