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

revision 4121 by gross, Wed Dec 12 06:17:03 2012 UTC revision 4122 by gross, Thu Dec 20 05:42:35 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, grad, inner, integrate, interpolate, kronecker, boundingBoxEdgeLengths, vol, sqrt, length  from esys.escript import 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:
# Line 251  class Regularization(CostFunction): Line 248  class Regularization(CostFunction):
248          :rtype: `LinearPDE`          :rtype: `LinearPDE`
249          """          """
250          return self.__pde          return self.__pde
251
252      def getDualProduct(self, m, r):      def getDualProduct(self, m, r):
253          """          """
254          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 410  class Regularization(CostFunction): Line 407  class Regularization(CostFunction):
407          DIM=self.getDomain().getDim()          DIM=self.getDomain().getDim()
408          numLS=self.getNumLevelSets()          numLS=self.getNumLevelSets()
409
n=0

410          if self.__w0 is not None:          if self.__w0 is not None:
411              Y = m * self.__w0 * mu              Y = m * self.__w0 * mu
412          else:          else:
413              Y = Data()          if numLS == 1:
414          n+=numLS            Y = Scalar(0,  grad_m.getFunctionSpace())
415            else:
416                  Y = Data(0, (numLS,) , grad_m.getFunctionSpace())
417
418          if self.__w1 is not None:          if self.__w1 is not None:
420              if numLS == 1:              if numLS == 1:
422              else:              else:
423                  for k in xrange(numLS):                  for k in xrange(numLS):
424                      X[k,:]*=mu[k]                      X[k,:]*=mu[k]
425          else:          else:
426              X = Data()              X = Data(0, grad_m.getShape(), grad_m.getFunctionSpace())
427          if numLS > 1:
428              raise NotImplementedError          # cross gradient terms:
429            if numLS > 1:
430          for  k in xrange(numLS):
433             for  l in xrange(k):
437               f=  mu_c[l,k]* self.__wc[l,k]
438               X[l,:] += f * ( l2_grad_m_l *  grad_m_l - grad_m_lk * grad_m_k )
439               X[k,:] += f * ( l2_grad_m_k *  grad_m_k - grad_m_lk * grad_m_l )
440
441          return ArithmeticTuple(Y, X)          return ArithmeticTuple(Y, X)
442
443
# Line 462  class Regularization(CostFunction): Line 471  class Regularization(CostFunction):
471                          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]
472
473              if numLS > 1:              if numLS > 1:
474                  raise NotImplementedError             for  k in xrange(numLS):
478                 for  l in xrange(k):
485                    f=  mu_c[l,k]* self.__wc[l,k]
486
487                    A[l,:,l:] += f * ( l2_grad_m_k * kronecker(DIM) - o_kk )
488                    A[l,:,k:] += f * ( 2 * o_lk -   o_kl - i_lk * kronecker(DIM) )
489                    A[k,:,l:] += f * ( 2 * o_kl -   o_lk - i_lk * kronecker(DIM) )
490                    A[k,:,k:] += f * ( l2_grad_m_l * kronecker(DIM) - o_ll )
491              self.getPDE().setValue(A=A)              self.getPDE().setValue(A=A)
492          #self.getPDE().resetRightHandSideCoefficients()          #self.getPDE().resetRightHandSideCoefficients()
493          #self.getPDE().setValue(X=r[1])          #self.getPDE().setValue(X=r[1])

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