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

revision 4238 by jfenwick, Thu Feb 21 09:58:49 2013 UTC revision 4263 by caltinay, Thu Feb 28 05:59:04 2013 UTC
# Line 29  from esys.escript import Function, outer Line 29  from esys.escript import Function, outer
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
import sys
if sys.version_info.major>2:
xrange=range
32
33  class Regularization(CostFunction):  class Regularization(CostFunction):
34      """      """
# Line 140  class Regularization(CostFunction): Line 137  class Regularization(CostFunction):
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 range(k) ] for k in range(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)
# Line 191  class Regularization(CostFunction): Line 188  class Regularization(CostFunction):
188                 raise ValueError("Non-positive weighting factor detected.")                 raise ValueError("Non-positive weighting factor detected.")
189          else:          else:
190
191               for k in xrange(numLevelSets):               for k in range(numLevelSets):
192                   A=0                   A=0
193                   if w0 is not None:                   if w0 is not None:
194                       A = integrate(w0[k])                       A = integrate(w0[k])
# Line 208  class Regularization(CostFunction): Line 205  class Regularization(CostFunction):
205
207                   if wc is not None:                   if wc is not None:
208                       for l in xrange(k):                       for l in range(k):
209                          A = integrate(wc[l,k])/L4                          A = integrate(wc[l,k])/L4
210                          if A > 0:                          if A > 0:
211                             f = scale_c[l,k]/A                             f = scale_c[l,k]/A
# Line 298  class Regularization(CostFunction): Line 295  class Regularization(CostFunction):
295          if mu.shape == (numTF,):          if mu.shape == (numTF,):
297              mu_c2=np.zeros((numLS,numLS))              mu_c2=np.zeros((numLS,numLS))
298              for k in xrange(numLS):              for k in range(numLS):
299                 for l in xrange(k):                 for l in range(k):
300                     mu_c2[l,k] = mu[numLS+l+((k-1)*k)/2]                     mu_c2[l,k] = mu[numLS+l+((k-1)*k)/2]
302          elif mu.shape == () and numLS ==1:          elif mu.shape == () and numLS ==1:
# Line 360  class Regularization(CostFunction): Line 357  class Regularization(CostFunction):
357          else:          else:
358              mu_c=np.asarray(mu_c)              mu_c=np.asarray(mu_c)
359              if mu_c.shape == (numLS,numLS):              if mu_c.shape == (numLS,numLS):
360                  if not all( [ [ mu_c[l,k] > 0. for l in xrange(k) ] for k in xrange(1,numLS) ]):                  if not all( [ [ mu_c[l,k] > 0. for l in range(k) ] for k in range(1,numLS) ]):
361                       raise ValueError("All trade-off factors in the lower triangle of mu_c must be positive.")                       raise ValueError("All trade-off factors in the lower triangle of mu_c must be positive.")
362                  else:                  else:
363                       self.__mu_c =  mu_c                       self.__mu_c =  mu_c
# Line 392  class Regularization(CostFunction): Line 389  class Regularization(CostFunction):
389              if numLS == 1:              if numLS == 1:
391              else:              else:
392                  for k in xrange(numLS):                  for k in range(numLS):
394
395          if numLS > 1:          if numLS > 1:
396              for k in xrange(numLS):              for k in range(numLS):
398                  len_gk=length(gk)                  len_gk=length(gk)
399                  for l in xrange(k):                  for l in range(k):
401                      A+= mu_c[l,k] * integrate( self.__wc[l,k] * ( len_gk * length(gl) )**2 - inner(gk, gl)**2 )                      A+= mu_c[l,k] * integrate( self.__wc[l,k] * ( len_gk * length(gl) )**2 - inner(gk, gl)**2 )
402          return A/2          return A/2
# Line 429  class Regularization(CostFunction): Line 426  class Regularization(CostFunction):
426              if numLS == 1:              if numLS == 1:
428              else:              else:
429                  for k in xrange(numLS):                  for k in range(numLS):
430                      X[k,:]*=mu[k]                      X[k,:]*=mu[k]
431          else:          else:
433
435          if numLS > 1:          if numLS > 1:
436            for  k in xrange(numLS):            for  k in range(numLS):
439               for  l in xrange(k):               for  l in range(k):
# Line 466  class Regularization(CostFunction): Line 463  class Regularization(CostFunction):
463                       D=self.__w0 * mu                       D=self.__w0 * mu
464                  else:                  else:
465                       D=self.getPDE().createCoefficient("D")                       D=self.getPDE().createCoefficient("D")
466                       for k in xrange(numLS): D[k,k]=self.__w0[k] * mu[k]                       for k in range(numLS): D[k,k]=self.__w0[k] * mu[k]
467                  self.getPDE().setValue(D=D)                  self.getPDE().setValue(D=D)
468
469              A=self.getPDE().createCoefficient("A")              A=self.getPDE().createCoefficient("A")
470              if self.__w1 is not None:              if self.__w1 is not None:
471                 if numLS == 1:                 if numLS == 1:
472                     for i in xrange(DIM): A[i,i]=self.__w1[i] * mu                     for i in range(DIM): A[i,i]=self.__w1[i] * mu
473                 else:                 else:
474                     for k in xrange(numLS):                     for k in range(numLS):
475                          for i in xrange(DIM): A[k,i,k,i]=self.__w1[k,i] * mu[k]                          for i in range(DIM): A[k,i,k,i]=self.__w1[k,i] * mu[k]
476
477              if numLS > 1:              if numLS > 1:
478                 for  k in xrange(numLS):                 for  k in range(numLS):