/[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 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    
206                   # and now the cross-gradient:                   # and now the cross-gradient:
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,):
296              self.setTradeOffFactorsForVariation(mu[:numLS])              self.setTradeOffFactorsForVariation(mu[:numLS])
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]
301              self.setTradeOffFactorsForCrossGradient(mu_c2)              self.setTradeOffFactorsForCrossGradient(mu_c2)
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:
390                  A+=integrate(inner(grad_m**2, self.__w1))*mu                  A+=integrate(inner(grad_m**2, self.__w1))*mu
391              else:              else:
392                  for k in xrange(numLS):                  for k in range(numLS):
393                      A+=mu[k]*integrate(inner(grad_m[k,:]**2,self.__w1[k,:]))                      A+=mu[k]*integrate(inner(grad_m[k,:]**2,self.__w1[k,:]))
394    
395          if numLS > 1:          if numLS > 1:
396              for k in xrange(numLS):              for k in range(numLS):
397                  gk=grad_m[k,:]                  gk=grad_m[k,:]
398                  len_gk=length(gk)                  len_gk=length(gk)
399                  for l in xrange(k):                  for l in range(k):
400                      gl=grad_m[l,:]                      gl=grad_m[l,:]
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:
427                  X=grad_m* self.__w1*mu                  X=grad_m* self.__w1*mu
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:
432              X = Data(0, grad_m.getShape(), grad_m.getFunctionSpace())              X = Data(0, grad_m.getShape(), grad_m.getFunctionSpace())
433    
434          # cross gradient terms:          # cross gradient terms:
435          if numLS > 1:          if numLS > 1:
436            for  k in xrange(numLS):            for  k in range(numLS):
437               grad_m_k=grad_m[k,:]               grad_m_k=grad_m[k,:]
438               l2_grad_m_k = length(grad_m_k)**2               l2_grad_m_k = length(grad_m_k)**2
439               for  l in xrange(k):               for  l in range(k):
440                 grad_m_l=grad_m[l,:]                 grad_m_l=grad_m[l,:]
441                 l2_grad_m_l = length(grad_m_l)**2                 l2_grad_m_l = length(grad_m_l)**2
442                 grad_m_lk = inner(grad_m_l, grad_m_k)                 grad_m_lk = inner(grad_m_l, grad_m_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):
479                   grad_m_k=grad_m[k,:]                   grad_m_k=grad_m[k,:]
480                   l2_grad_m_k = length(grad_m_k)**2                   l2_grad_m_k = length(grad_m_k)**2
481                   o_kk=outer(grad_m_k, grad_m_k)                   o_kk=outer(grad_m_k, grad_m_k)
482                   for  l in xrange(k):                   for  l in range(k):
483                      grad_m_l=grad_m[l,:]                      grad_m_l=grad_m[l,:]
484                      l2_grad_m_l = length(grad_m_l)**2                      l2_grad_m_l = length(grad_m_l)**2
485                      i_lk = inner(grad_m_l, grad_m_k)                      i_lk = inner(grad_m_l, grad_m_k)
Legend:
Removed from v.4238  
changed lines
  Added in v.4263
  ViewVC Help
Powered by ViewVC 1.1.26