/[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 4097 by caltinay, Fri Dec 7 01:18:35 2012 UTC revision 4099 by gross, Tue Dec 11 04:04:47 2012 UTC
# Line 34  class Regularization(CostFunction): Line 34  class Regularization(CostFunction):
34      The regularization term for the level set function ``m`` within the cost      The regularization term for the level set function ``m`` within the cost
35      function J for an inversion:      function J for an inversion:
36    
37      *J(m)=1/2 * sum_k imtegrate( mu_0[k] * s0[k] * m_k**2 + mu_1[k] * s1[k,i] * m_{k,i}**2) + sum_l<k mu_c[l,k] sc[l,l] * | curl(m_k) x curl(m_l) |^2*      *J(m)=1/2 * sum_k imtegrate( mu[k] * ( w0[k] * m_k**2 * w1[k,i] * m_{k,i}**2) + sum_l<k mu_c[l,k] wc[l,k] * | curl(m_k) x curl(m_l) |^2*
38    
39      where s0[k], s1[k,i] and  sc[k,l] are non-negative scaling factors and      where w0[k], w1[k,i] and  wc[k,l] are non-negative weighting factors and
40      mu_0[k], mu_1[k], mu_c[l,k] are weighting factors which may be altered      mu[k] and mu_c[l,k] are trade-off factors which may be altered
41      during the inversion. The scaling factors are normalized such that their      during the inversion. The weighting factors are normalized such that their
42      integrals over the domain are constant:      integrals over the domain are constant:
43    
44      *integrate(s0[k])=1*      *integrate(w0[k] + inner(w1[k,:],1/L[:]**2))=scale[k]* volume(domain)*
45      *integrate(inner(s1[k,:],L[:])=1*      *integrate(wc[l,k]*1/L**4)=scale_c[k]* volume(domain) *
     *integrate(inner(sc[l,k]*L**4)=1*  
46    
47      """      """
48      def __init__(self, domain, numLevelSets=1,      def __init__(self, domain, numLevelSets=1,
49                         s0=None, s1=None, sc=None,                         w0=None, w1=None, wc=None,
50                         location_of_set_m=Data(),                         location_of_set_m=Data(),
51                         useDiagonalHessianApproximation=True, tol=1e-8):                         useDiagonalHessianApproximation=False, tol=1e-8,
52                           scale=None, scale_c=None):
53          """          """
54          initialization.          initialization.
55    
# Line 57  class Regularization(CostFunction): Line 57  class Regularization(CostFunction):
57          :type domain: `Domain`          :type domain: `Domain`
58          :param numLevelSets: number of level sets          :param numLevelSets: number of level sets
59          :type numLevelSets: ``int``          :type numLevelSets: ``int``
60          :param s0: scaling factor for the m**2 term. If not set zero is assumed.          :param w0: weighting factor for the m**2 term. If not set zero is assumed.
61          :type s0: ``Scalar`` if ``numLevelSets`` == 1 or `Data` object of shape          :type w0: ``Scalar`` if ``numLevelSets`` == 1 or `Data` object of shape
62                    (``numLevelSets`` ,) if ``numLevelSets`` > 1                    (``numLevelSets`` ,) if ``numLevelSets`` > 1
63          :param s1: scaling factor for the grad(m_i) terms. If not set zero is assumed          :param w1: weighting factor for the grad(m_i) terms. If not set zero is assumed
64          :type s1: ``Vector`` if ``numLevelSets`` == 1 or `Data` object of shape          :type w1: ``Vector`` if ``numLevelSets`` == 1 or `Data` object of shape
65                    (``numLevelSets`` , DIM) if ``numLevelSets`` > 1                    (``numLevelSets`` , DIM) if ``numLevelSets`` > 1
66          :param sc: scaling factor for the cross correlation terms. If not set          :param wc: weighting factor for the cross gradient terms. If not set
67                     zero is assumed. Used for the case if ``numLevelSets`` > 1                     zero is assumed. Used for the case if ``numLevelSets`` > 1
68                     only. Only values ``sc[l,k]`` in the lower triangle (l<k)                     only. Only values ``wc[l,k]`` in the lower triangle (l<k)
69                     are used.                     are used.
70          :type sc: `Data` object of shape (``numLevelSets`` , ``numLevelSets``)          :type wc: `Data` object of shape (``numLevelSets`` , ``numLevelSets``)
71          :param location_of_set_m: marks location of zero values of the level          :param location_of_set_m: marks location of zero values of the level
72                                    set function ``m`` by a positive entry.                                    set function ``m`` by a positive entry.
73          :type location_of_set_m: ``Scalar`` if ``numLevelSets`` == 1 or `Data`          :type location_of_set_m: ``Scalar`` if ``numLevelSets`` == 1 or `Data`
74                  object of shape (``numLevelSets`` ,) if ``numLevelSets`` > 1                  object of shape (``numLevelSets`` ,) if ``numLevelSets`` > 1
75          :param useDiagonalHessianApproximation: if True cross correlation terms          :param useDiagonalHessianApproximation: if True cross gradient terms
76                      between level set components are ignored when calculating                      between level set components are ignored when calculating
77                      approximations of the inverse of the Hessian Operator.                      approximations of the inverse of the Hessian Operator.
78                      This can speed-up the calculation of the inverse but may                      This can speed-up the calculation of the inverse but may
# Line 82  class Regularization(CostFunction): Line 82  class Regularization(CostFunction):
82          :param tol: tolerance when solving the PDE for the inverse of the          :param tol: tolerance when solving the PDE for the inverse of the
83                      Hessian Operator                      Hessian Operator
84          :type tol: positive ``float``          :type tol: positive ``float``
85            
86            :param scale: weighting factor for level set function variation terms. If not set one is used.
87            :type scale: ``Scalar`` if ``numLevelSets`` == 1 or `Data` object of shape
88                      (``numLevelSets`` ,) if ``numLevelSets`` > 1
89            :param scale_c: scale for the cross gradient terms. If not set
90                       one is assumed. Used for the case if ``numLevelSets`` > 1
91                       only. Only values ``scale_c[l,k]`` in the lower triangle (l<k)
92                       are used.
93            :type scale_c: `Data` object of shape (``numLevelSets`` , ``numLevelSets``)
94            
95    
96          """          """
97          if numLevelSets>1:          if numLevelSets>1:
98                raise ValueError("Currently only numLevelSets<=1 is supported.")                raise ValueError("Currently only numLevelSets<=1 is supported.")
99          if s0 == None and s1==None:          
100                raise ValueError("Values for s0 or s1 must be given.")          if w0 == None and w1==None:
101                  raise ValueError("Values for w0 or for w1 must be given.")
102        if wc == None and  numLevelSets>1:
103              raise ValueError("Values for wc must be given.")
104    
105          self.__domain=domain          self.__domain=domain
106          DIM=self.__domain.getDim()          DIM=self.__domain.getDim()
         self.__L2=np.asarray(boundingBoxEdgeLengths(domain))**2  
         self.__L4=np.sum(self.__L2)**2  
107          self.__numLevelSets=numLevelSets          self.__numLevelSets=numLevelSets
108            
         if self.__numLevelSets > 1:  
             self.__useDiagonalHessianApproximation=useDiagonalHessianApproximation  
         else:  
             self.__useDiagonalHessianApproximation=True  
         self._update_Hessian=True  
   
109          self.__pde=LinearPDE(self.__domain, numEquations=self.__numLevelSets)          self.__pde=LinearPDE(self.__domain, numEquations=self.__numLevelSets)
110          self.__pde.getSolverOptions().setTolerance(tol)          self.__pde.getSolverOptions().setTolerance(tol)
111          self.__pde.setSymmetryOn()          self.__pde.setSymmetryOn()
         self.__pde_is_set=False  
112          try:          try:
113              self.__pde.setValue(q=location_of_set_m)              self.__pde.setValue(q=location_of_set_m)
114          except IllegalCoefficientValue:          except IllegalCoefficientValue:
115              raise ValueError("Unable to set location of fixed level set function.")              raise ValueError("Unable to set location of fixed level set function.")
116          
117          self.__total_num_weights=2*numLevelSets+((numLevelSets-1)*numLevelSets)/2          # =========== check the shape of the scales: =================================
118          self.__weight_index=[]  # this is a mapping from the relevant mu-coefficients to the set of all mu-coefficients          if scale is None:
119                             # we count s0, then s1, then sc (k<l).          if numLevelSets == 1 :
120          # THE S0 weighting factor             scale = 1.
121          n=0          else:
122          VV=vol(domain)             scale = np.ones((numLevelSets,))
123          if s0 is not None:      else:
124              s0 = interpolate(s0,self.__pde.getFunctionSpaceForCoefficient('D'))          scale=np.asarray(scale)
125              s=s0.getShape()          if numLevelSets == 1 :
126              if numLevelSets == 1 :              if scale.shape == ():
127                   if s == () :             if not scale > 0 :
128                       V=integrate(s0)                raise ValueError("Value for scale must be positive.")
129                       if V > 0:          else:
130                         self.__weight_index.append(n)             raise ValueError("Unexpected shape %s for scale."%scale.shape)
                        s0*=VV/V  
                      else:  
                        s0=None  
                  else:  
                      raise ValueError("Unexpected shape %s for weight s0."%s)  
131              else:              else:
132                   if s == (numLevelSets,):               if scale.shape is (numLevelSets,):
133                       for k in xrange(numLevelSets):               if not min(scale) > 0:
134                          V=integrate(s0[k])                  raise ValueError("All value for scale must be positive.")
135                          if V > 0:           else:
136                              self.__weight_index.append(n+k)             raise ValueError("Unexpected shape %s for scale."%scale.shape)
137                              s0[k]*=VV/V      
138                   else:          if scale_c is None or numLevelSets < 2:
139                       raise ValueError("Unexpected shape %s for weight s0."%s)          scale_c = np.ones((numLevelSets,numLevelSets))
140          self.__s0=s0      else:
141          n+=numLevelSets          scale_c=np.asarray(scale_c)
142            if scale_c.shape == (numLevelSets,numLevelSets):
143          # The S1 weighting factor              if not all( [ [ scale_c[l,k] > 0. for l in xrange(k) ] for k in xrange(1,numLevelSets) ]):
144          if not s1 is None:                  raise ValueError("All values in the lower triangle of scale_c must be positive.")
             s1 = interpolate(s1,self.__pde.getFunctionSpaceForCoefficient('A'))  
             s=s1.getShape()  
   
             if numLevelSets == 1 :  
                 if s == (DIM,) :  
                     V=integrate(inner(s1, 1/self.__L2))  
                     if V > 0:  
                         self.__weight_index.append(n)  
                         s1*=VV/V  
                     print "REG SCALE = ",s1  
                 else:  
                     raise ValueError("Unexpected shape %s for weight s1."%s)  
145              else:              else:
146                  if s == (numLevelSets,DIM):           raise ValueError("Unexpected shape %s for scale."%scale_c.shape)
147                      for k in xrange(numLevelSets):      # ===== check the shape of the weights: ============================================
148                          for i in xrange(DIM):          if w0 is not None:
149                              ww=s1[k,:]            w0 = interpolate(w0,self.__pde.getFunctionSpaceForCoefficient('D'))
150                              V=integrate(inner(ww,1/self.__L2))            s0=w0.getShape()
151                              if V > 0:            if numLevelSets == 1 :
152                                  self.__weight_index.append(n+i)             if  not s0 == () :
153                                  s1[k,:]=ww*(VV/V)                raise ValueError("Unexpected shape %s for weight w0."%s0)
154                  else:                else:
155                      raise ValueError("Unexpected shape %s for weight s1."%s)             if not s0 == (numLevelSets,):
156                  raise ValueError("Unexpected shape %s for weight w0."%s0)
157        if not w1 is None:
158              w1 = interpolate(w1,self.__pde.getFunctionSpaceForCoefficient('A'))
159              s1=w1.getShape()
160              print s1, (DIM,)
161              print w1
162              if numLevelSets is 1 :
163               if not s1 == (DIM,) :
164                  raise ValueError("Unexpected shape %s for weight w1."%s1)
165                  else:
166               if not s1 == (numLevelSets,DIM):
167                  raise ValueError("Unexpected shape %s for weight w1."%s1)
168            if numLevelSets == 1 :
169                 wc=None
170            else:
171                 wc = interpolate(wc,self.__pde.getFunctionSpaceForCoefficient('A'))
172                 sc=wc.getShape()
173                 raise ValueError("Unexpected shape %s for weight wc."%sc)
174            # ============= now we rescale weights: ======================================
175            vol_d=vol(self.__domain)
176            L2s=np.asarray(boundingBoxEdgeLengths(domain))**2
177            L4=1/np.sum(1/L2s)**2
178            if numLevelSets is 1 :
179                A=0
180                if w0 is not None:
181                A = integrate(w0)
182            if w1 is not None:
183                A += integrate(inner(w1, 1/L2s))
184            if A > 0:
185                f = scale*vol_d/A
186                if w0 is not None:
187                     w0*=f
188                if w1 is not None:
189                 w1*=f
190                else:
191               raise ValueError("Non-positive weighting factor detected.")
192            else:
193             for k in xrange(numLevelSets):
194                 A=0
195                     if w0 is not None:
196                     A = integrate(w0[k])
197                 if w1 is not None:
198                      A += integrate(inner(w1[k,:], 1/L2s))
199                 if A > 0:
200                      f = scale[k]*vol_d/A
201                      if w0 is not None:
202                         w0[k]*=f
203                      if w1 is not None:
204                     w1[k,:]*=f
205                     else:
206                   raise ValueError("Non-positive weighting factor for level set component %d detected."%k)
207            
208                 # and now the cross-gradient:
209                 for l in xrange(k):  
210                     A = integrate(wc[l,k])/L4
211                     if A > 0:
212                        f = scale_c[l,k]*vol_d/A
213                        wc[l,k]*=f
214                         else:
215                        raise ValueError("Non-positive weighting factor for cross-gradient level set components %d and %d detected."%(l,k))
216                    
217            self.__w0=w0
218            self.__w1=w1
219            self.__wc=wc
220            
221    
         self.__s1=s1  
         n+=numLevelSets  
222    
223          # The SC weighting factor          self.__pde_is_set=False        
224          if not sc is None:          if self.__numLevelSets > 1:
225              if numLevelSets == 1 :              self.__useDiagonalHessianApproximation=useDiagonalHessianApproximation
226                  sc=None          else:
227              else:              self.__useDiagonalHessianApproximation=True
228                  sc = interpolate(sc,self.__pde.getFunctionSpaceForCoefficient('A'))          self._update_Hessian=True
                 s=sc.getShape()  
                 if s == (numLevelSets,numLevelSets):  
                     for k in xrange(numLevelSets):  
                         sc[k,k]=0.  
                         for l in xrange(k):  
                             ww=sc[l,k]  
                             V=integrate(ww)  
                             if V > 0:  
                                 self.__weight_index.append(n+k*numLevelSets+l)  
                                 sc[l,k]=ww*VV/V*self.__L4  
                                 sc[k,l]=0  
                 else:  
                     raise ValueError("Unexpected shape %s for weight s0."%s)  
229    
230          self.__sc=sc          self.__num_tradeoff_factors=numLevelSets+((numLevelSets-1)*numLevelSets)/2
231          self.setWeights()          self.setTradeOffFactors()
232    
233      def getDomain(self):      def getDomain(self):
234          """          """
# Line 232  class Regularization(CostFunction): Line 269  class Regularization(CostFunction):
269          if not r[0].isEmpty(): A+=integrate(inner(r[0], m))          if not r[0].isEmpty(): A+=integrate(inner(r[0], m))
270          if not r[1].isEmpty(): A+=integrate(inner(r[1], grad(m)))          if not r[1].isEmpty(): A+=integrate(inner(r[1], grad(m)))
271          return A          return A
272        def getNumTradeOffFactors(self):
273            """
274            returns the number of trade-off factors being used.
275    
276            :rtype: ``int``
277            """
278            return self.__num_tradeoff_factors
279    
280        def setTradeOffFactors(self, mu=None):
281            """
282            sets the trade-off factors for the level-set variation and the cross-gradient
283            
284            :param mu: new values for the trade-off factors where values mu[:numLevelSets] are the
285                       trade-off factors for the level-set variation and the remaining values for
286                       the cross-gradient part with mu_c[l,k]=mu[numLevelSets+l+((k-1)*k)/2] (l<k).
287                       If no values for mu is given ones are used. Values must be positive.
288            :type mu: ``list`` of ``float`` or ```numpy.array```
289            """
290            numLS=self.getNumLevelSets()
291            numTF=self.getNumTradeOffFactors()
292            if mu is None:
293           mu = np.ones((numTF,))
294        else:
295           mu = np.asarray(mu)
296    
297        if mu.shape == (numTF,):
298            self.setTradeOffFactorsForVariation(mu[:numLS])
299            mu_c2=np.zeros((numLS,numLS))
300            for k in xrange(numLS):
301               for l in xrange(k):
302               mu_c2[l,k] = mu[numLS+l+((k-1)*k)/2]
303            self.setTradeOffFactorsForCrossGradient(mu_c2)
304        elif mu.shape == () and numLS ==1:
305            self.setTradeOffFactorsForVariation(mu)
306        else:
307           raise ValueError("Unexpected shape %s for mu."%(mu.shape,))
308          
309        def setTradeOffFactorsForVariation(self, mu=None):
310             """
311             sets the trade-off factors for the level-set variation part
312            
313             :param mu:  new values for the trade-off factors. Values must be positive.
314             :type mu: `float``, ``list`` of ``float`` or ```numpy.array```
315             """
316             numLS=self.getNumLevelSets()
317             if mu is None:
318            if numLS == 1:
319               mu = 1.
320            else:
321               mu = np.ones((numLS,))
322    
323         mu=np.asarray(mu)
324         if numLS == 1:
325           if mu.shape == (1,): mu=mu[0]
326           if mu.shape == ():
327              if mu > 0:
328             self.__mu= mu
329             self._new_mu=True
330              else:
331             raise ValueError("Value for trade-off factor must be positive.")
332           else:
333              raise ValueError("Unexpected shape %s for mu."%mu.shape)
334         else:
335           if mu.shape == (numLS,):
336               if min(mu) > 0:
337               self.__mu= mu
338               self._new_mu=True
339               else:
340               raise ValueError("All value for mu must be positive.")
341               else:
342               raise ValueError("Unexpected shape %s for trade-off factor."%mu.shape)
343    
344        def setTradeOffFactorsForCrossGradient(self, mu_c=None):
345            """
346            sets the trade-off factors for the cross--gradient terms
347            
348            :param mu_c:  new values for the trade-off factors for the cross--gradient terms. Values must be positive.
349                          if now value is given ones are used. Onky value mu_c[l,k] for l<k are used.
350            :type mu: `float``, ``list`` of ``float`` or ```numpy.array```
351            
352            """
353            numLS=self.getNumLevelSets()
354            if mu_c is None or numLS < 2:
355            self.__mu_c = np.ones((numLS,numLS))
356        else:
357            mu_c=np.asarray(mu_c)
358            if mu_c.shape == (numLS,numLS):
359                if not all( [ [ mu_c[l,k] > 0. for l in xrange(k) ] for k in xrange(1,numLS) ]):
360                 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:
365             raise ValueError("Unexpected shape %s for mu."%mu_c.shape)
366        
367        def getArguments(self, m):
368            """
369            """
370            return ( grad(m),)
371            
372    
373      def getValue(self, m, grad_m):      def getValue(self, m, grad_m):
374          """          """
# Line 239  class Regularization(CostFunction): Line 376  class Regularization(CostFunction):
376    
377          :rtype: ``float``          :rtype: ``float``
378          """          """
379          mu=self.getWeights( uncompress=True)          mu=self.__mu
380            mu_c=self.__mu_c
381          DIM=self.getDomain().getDim()          DIM=self.getDomain().getDim()
382          numLS=self.getNumLevelSets()          numLS=self.getNumLevelSets()
383    
384          A=0          A=0
385          n=0          if self.__w0 is not None:
386               A+=inner(integrate(m**2 * self.__w0), mu)
         if self.__s0 is not None:  
             A+=inner(integrate(m**2*self.__s0), mu[:numLS])  
         n+=numLS  
387    
388          if self.__s1 is not None:          if self.__w1 is not None:
389              if numLS == 1:              if numLS == 1:
390                  A+=integrate(inner(grad_m**2, self.__s1))*mu[n]                  A+=integrate(inner(grad_m**2, self.__w1))*mu
391              else:              else:
392                  for k in xrange(numLS):                  for k in xrange(numLS):
393                      A+=mu[n+k]*integrate(inner(grad_m[k,:]**2,self.__s1[k,:]))                      A+=mu[k]*integrate(inner(grad_m[k,:]**2,self.__w1[k,:]))
394          n+=numLS                      
395            if numLS > 1:
         if self.__sc is not None:  
396              for k in xrange(numLS):              for k in xrange(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 xrange(k):
400                      gl=grad_m[l,:]                      gl=grad_m[l,:]
401                      A+= mu[n+k*numLS+l] * integrate( self.__sc[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
403    
404      def getGradient(self, m,  grad_m):      def getGradient(self, m,  grad_m):
# Line 273  class Regularization(CostFunction): Line 407  class Regularization(CostFunction):
407          The function returns Y_k=dPsi/dm_k and X_kj=dPsi/dm_kj          The function returns Y_k=dPsi/dm_k and X_kj=dPsi/dm_kj
408          """          """
409    
410          mu=self.getWeights( uncompress=True)          mu=self.__mu
411            mu_c=self.__mu_c
412          DIM=self.getDomain().getDim()          DIM=self.getDomain().getDim()
413          numLS=self.getNumLevelSets()          numLS=self.getNumLevelSets()
414    
415          n=0          n=0
416    
417          if self.__s0 is not None:          if self.__w0 is not None:
418              Y = m * self.__s0 * mu[:numLS]              Y = m * self.__w0 * mu
419          else:          else:
420              Y = Data()              Y = Data()
421          n+=numLS          n+=numLS
422    
423          if self.__s1 is not None:          if self.__w1 is not None:
424              if numLS == 1:              if numLS == 1:
425                  X=grad_m* self.__s1*mu[n]                  X=grad_m* self.__w1*mu
426              else:              else:
427                  X=self.getPDE().createCoefficient("X")                  X=grad_m[k,:]*self.__w1[k,:]
428                  for k in xrange(numLS):                  for k in xrange(numLS):
429                      X[k,:]=mu[n+k]*grad_m[k,:]*self.__s1[k,:]                      X[k,:]*=mu[k]
430          else:          else:
431              X = Data()              X = Data()
432          n+=numLS          if numLS > 1:
         if self.__sc is not None:  
433              raise NotImplementedError              raise NotImplementedError
434          return ArithmeticTuple(Y, X)          return ArithmeticTuple(Y, X)
435    
436      def getArguments(self, m):  
         """  
         """  
         return ( grad(m),)  
437    
438      def getInverseHessianApproximation(self, m, r, grad_m):      def getInverseHessianApproximation(self, m, r, grad_m):
439          """          """
# Line 311  class Regularization(CostFunction): Line 442  class Regularization(CostFunction):
442    
443              self._new_mu=False              self._new_mu=False
444              self._update_Hessian=False              self._update_Hessian=False
445                mu=self.__mu
446                mu_c=self.__mu_c
447    
             mu=self.getWeights( uncompress=True)  
448              DIM=self.getDomain().getDim()              DIM=self.getDomain().getDim()
449              numLS=self.getNumLevelSets()              numLS=self.getNumLevelSets()
450              n=0              if self.__w0 is not None:
             if self.__s0 is not None:  
451                  if numLS == 1:                  if numLS == 1:
452                       D=self.__s0 * mu[n]                       D=self.__w0 * mu
453                  else:                  else:
454                       D=self.getPDE().createCoefficient("D")                       D=self.getPDE().createCoefficient("D")
455                       for k in xrange(numLS): D[k,k]=self.__s0[k] * mu[n+k]                       for k in xrange(numLS): D[k,k]=self.__w0[k] * mu[k]
456                  self.getPDE().setValue(D=D)                  self.getPDE().setValue(D=D)
             n+=numLS  
457    
458              A=self.getPDE().createCoefficient("A")              A=self.getPDE().createCoefficient("A")
459              if self.__s1 is not None:              if self.__w1 is not None:
460                 if numLS == 1:                 if numLS == 1:
461                     for i in xrange(DIM): A[i,i]=self.__s1[i] * mu[n]                     for i in xrange(DIM): A[i,i]=self.__w1[i] * mu
462                 else:                 else:
463                     for k in xrange(numLS):                     for k in xrange(numLS):
464                          for i in xrange(DIM): A[k,i,k,i]=self.__s1[k,i] * mu[n+k]                          for i in xrange(DIM): A[k,i,k,i]=self.__w1[k,i] * mu[k]
465              n+=numLS                          
466              if self.__sc is not None:              if numLS > 1:
467                  raise NotImplementedError                  raise NotImplementedError
468              print A              print A
469              self.getPDE().setValue(A=A)              self.getPDE().setValue(A=A)
# Line 354  class Regularization(CostFunction): Line 484  class Regularization(CostFunction):
484          """          """
485          if not self.__useDiagonalHessianApproximation:          if not self.__useDiagonalHessianApproximation:
486              self._update_Hessian=True              self._update_Hessian=True
   
    # ================ we should factor these out ==============================  
   
     def getNumRelevantTerms(self):  
         """  
         returns the number of terms in the cost function of regularization  
         with non-zero scaling factors.  
   
         :rtype: ``int``  
         """  
         return len(self.__weight_index)  
   
     def getNumTerms(self):  
         """  
         returns the number of terms in the cost function of regularization  
         with non-zero scaling factors.  
   
         :rtype: ``int``  
         """  
         return len(self.__weight_index)  
   
     def setWeights(self, mu=None):  
         """  
         sets the values for the weighting terms in the cost function.  
         Note that only values corresponding to entries with non-negative  
         scaling factors are used.  
   
         :param mu: weights  
         :type mu: ``list`` of ``float`` or ``np,array``  
         """  
         if mu == None:  
             mu = np.ones(self.getNumRelevantTerms())  
         mu=np.asarray(mu)  
   
         if len(mu) == self.getNumRelevantTerms():  
             if not mu.shape == (self.getNumRelevantTerms(),):  
                 raise ValueError("%s values are expected."%self.getNumRelevantTerms())  
             self.__mu=mu  
         else:  
             if not mu.shape == (self.__total_num_weights,):  
                 raise ValueError("%s values are expected."%self.__total_num_weights)  
             self.__mu = np.zeros(self.getNumRelevantTerms())  
             for i in xrange(len(self.__weight_index)): self.__mu[i]=mu[self.__weight_index[i]]  
         self._new_mu=True  
   
     def setWeightsForS0(self, mu=None):  
          """  
          sets the weights for s0-terms  
          """  
          numLS=self.getNumLevelSets()  
          my_mu=self.getWeights(uncompress=True)  
          if mu is None:  
              my_mu[:numLS]=1  
          else:  
              my_mu[:numLS]=mu  
          self.setWeights(my_mu)  
   
     def setWeightsForS1(self, mu=None):  
          """  
          sets the weights for s1-terms  
          """  
          numLS=self.getNumLevelSets()  
          my_mu=self.getWeights(uncompress=True)  
          if mu is None:  
               my_mu[numLS:2*numLS]=1  
          else:  
               my_mu[numLS:2*numLS]=mu  
          self.setWeights(my_mu)  
   
     def setWeightsForSc(self, mu):  
          """  
          sets the weights for sc-terms  
          """  
          numLS=self.getNumLevelSets()  
          my_mu=self.getWeights(uncompress=True)  
          if mu is None:  
               my_mu[2*numLS:]=1  
          else:  
               my_mu[2*numLS:]=mu  
          self.setWeights(my_mu)  
   
     def getWeights(self, uncompress=False):  
         """  
         Returns the weights for the terms in the cost function.  
         The first ``numLevelSets`` values are used for the regularization terms  
         and the remaining values for the cross correlation terms.  
         """  
         if uncompress:  
             mu = np.zeros(self.__total_num_weights)  
             for i in xrange(len(self.__weight_index)):  
                 mu[self.__weight_index[i]] = self.__mu[i]  
             return mu  
         else:  
             return self.__mu  
   
     def getWeightIndex(self):  
         """  
         returns an index to the contributions of terms with non-zero scaling  
         factor.  
         """  
         return self.__weight_index  
   

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

  ViewVC Help
Powered by ViewVC 1.1.26