/[escript]/branches/subworld2/downunder/py_src/regularizations.py
ViewVC logotype

Diff of /branches/subworld2/downunder/py_src/regularizations.py

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 5504 by jfenwick, Wed Mar 4 22:58:13 2015 UTC revision 5505 by jfenwick, Wed Mar 4 23:05:47 2015 UTC
# Line 50  class Regularization(CostFunction): Line 50  class Regularization(CostFunction):
50      """      """
51      def __init__(self, domain, numLevelSets=1,      def __init__(self, domain, numLevelSets=1,
52                         w0=None, w1=None, wc=None,                         w0=None, w1=None, wc=None,
53                         location_of_set_m=Data(),                         location_of_set_m=None,
54                         useDiagonalHessianApproximation=False, tol=1e-8,                         useDiagonalHessianApproximation=False, tol=1e-8,
55                         coordinates=None,                         coordinates=None,
56                         scale=None, scale_c=None):                         scale=None, scale_c=None):
# Line 105  class Regularization(CostFunction): Line 105  class Regularization(CostFunction):
105          if wc is None and numLevelSets>1:          if wc is None and numLevelSets>1:
106              raise ValueError("Values for wc must be given.")              raise ValueError("Values for wc must be given.")
107    
108                
109          self.logger = logging.getLogger('inv.%s'%self.__class__.__name__)          self.logger = logging.getLogger('inv.%s'%self.__class__.__name__)
         self.__domain=domain  
         DIM=self.__domain.getDim()  
110          self.__numLevelSets=numLevelSets          self.__numLevelSets=numLevelSets
111          self.__trafo=makeTranformation(domain, coordinates)          self.__domain=domain
112          self.__pde=LinearPDE(self.__domain, numEquations=self.__numLevelSets, numSolutions=self.__numLevelSets)          
         self.__pde.getSolverOptions().setTolerance(tol)  
         self.__pde.setSymmetryOn()  
         self.__pde.setValue(A=self.__pde.createCoefficient('A'), D=self.__pde.createCoefficient('D'), )  
         try:  
             self.__pde.setValue(q=location_of_set_m)  
         except IllegalCoefficientValue:  
             raise ValueError("Unable to set location of fixed level set function.")  
   
113          # =========== check the shape of the scales: ========================          # =========== check the shape of the scales: ========================
114          if scale is None:          if scale is None:
115              if numLevelSets == 1 :              if numLevelSets == 1 :
# Line 148  class Regularization(CostFunction): Line 139  class Regularization(CostFunction):
139                  if not all( [ [ scale_c[l,k] > 0. for l in range(k) ] for k in range(1,numLevelSets) ]):                  if not all( [ [ scale_c[l,k] > 0. for l in range(k) ] for k in range(1,numLevelSets) ]):
140                      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.")
141              else:              else:
142                  raise ValueError("Unexpected shape %s for scale."%scale_c.shape)                  raise ValueError("Unexpected shape %s for scale."%scale_c.shape)        
143          # ===== check the shape of the weights: =============================          
144          if w0 is not None:      if domain is None:
145              w0 = interpolate(w0,self.__pde.getFunctionSpaceForCoefficient('D'))          if type(dimension)==int:
146              s0=w0.getShape()  #         DIM=dimension
147              if numLevelSets == 1:  #         self.__trafo=None
148                  if not s0 == () :  #         self.__pde=None            
149                      raise ValueError("Unexpected shape %s for weight w0."%(s0,))                    pass
150              else:          else:
151                  if not s0 == (numLevelSets,):            raise ValueError("If domain is not provided, dimension must be given as an integer")
152                      raise ValueError("Unexpected shape %s for weight w0."%(s0,))      else:
153              if not self.__trafo.isCartesian():          self.__domain=None  # bypass the already set check
154                  w0*=self.__trafo.getVolumeFactor()          self.configure(domain, w0, w1, wc, location_of_set_m, useDiagonalHessianApproximation, tol,
155          if not w1 is None:                         coordinates, scale, scale_c)
156              w1 = interpolate(w1,self.__pde.getFunctionSpaceForCoefficient('A'))  #       DIM=self.__domain.getDim()
             s1=w1.getShape()  
             if numLevelSets == 1 :  
                 if not s1 == (DIM,) :  
                     raise ValueError("Unexpected shape %s for weight w1."%(s1,))  
             else:  
                 if not s1 == (numLevelSets,DIM):  
                     raise ValueError("Unexpected shape %s for weight w1."%(s1,))  
             if not self.__trafo.isCartesian():  
                 f=self.__trafo.getScalingFactors()**2*self.__trafo.getVolumeFactor()  
                 if numLevelSets == 1:  
                     w1*=f  
                 else:  
                     for i in range(numLevelSets): w1[i,:]*=f  
   
         if numLevelSets == 1:  
             wc=None  
         else:  
             wc = interpolate(wc,self.__pde.getFunctionSpaceForCoefficient('A'))  
             sc=wc.getShape()  
             if not sc == (numLevelSets, numLevelSets):  
                 raise ValueError("Unexpected shape %s for weight wc."%(sc,))  
             if not self.__trafo.isCartesian():  
                 raise ValueError("Non-cartesian coordinates for cross-gradient term is not supported yet.")  
         # ============= now we rescale weights: =============================  
         L2s=np.asarray(boundingBoxEdgeLengths(domain))**2  
         L4=1/np.sum(1/L2s)**2  
         if numLevelSets == 1:  
             A=0  
             if w0 is not None:  
                 A = integrate(w0)  
             if w1 is not None:  
                 A += integrate(inner(w1, 1/L2s))  
             if A > 0:  
                 f = scale/A  
                 if w0 is not None:  
                     w0*=f  
                 if w1 is not None:  
                     w1*=f  
             else:  
                 raise ValueError("Non-positive weighting factor detected.")  
         else: # numLevelSets > 1  
             for k in range(numLevelSets):  
                 A=0  
                 if w0 is not None:  
                     A = integrate(w0[k])  
                 if w1 is not None:  
                     A += integrate(inner(w1[k,:], 1/L2s))  
                 if A > 0:  
                     f = scale[k]/A  
                     if w0 is not None:  
                         w0[k]*=f  
                     if w1 is not None:  
                         w1[k,:]*=f  
                 else:  
                     raise ValueError("Non-positive weighting factor for level set component %d detected."%k)  
   
                 # and now the cross-gradient:  
                 if wc is not None:  
                     for l in range(k):  
                         A = integrate(wc[l,k])/L4  
                         if A > 0:  
                             f = scale_c[l,k]/A  
                             wc[l,k]*=f  
 #                       else:  
 #                           raise ValueError("Non-positive weighting factor for cross-gradient level set components %d and %d detected."%(l,k))  
   
         self.__w0=w0  
         self.__w1=w1  
         self.__wc=wc  
157    
158          self.__pde_is_set=False          self.__pde_is_set=False
159          if self.__numLevelSets > 1:          if self.__numLevelSets > 1:
# Line 242  class Regularization(CostFunction): Line 164  class Regularization(CostFunction):
164    
165          self.__num_tradeoff_factors=numLevelSets+((numLevelSets-1)*numLevelSets)//2          self.__num_tradeoff_factors=numLevelSets+((numLevelSets-1)*numLevelSets)//2
166          self.setTradeOffFactors()          self.setTradeOffFactors()
167          self.__vol_d=vol(self.__domain)  
168        def configure(self, domain,
169                           w0=None, w1=None, wc=None,
170                           location_of_set_m=None,
171                           useDiagonalHessianApproximation=False, tol=1e-8,
172                           coordinates=None,
173                           scale=None, scale_c=None):
174        if self.__domain is not None:
175            raise ValueError("Attempt to configure a regularization which was already configured.")
176            self.__domain=domain
177        DIM=self.__domain.getDim()
178        if location_of_set_m is None:
179              location_of_set_m=Data()
180        self.__trafo=makeTranformation(domain, coordinates)
181        self.__pde=LinearPDE(self.__domain, numEquations=self.__numLevelSets, numSolutions=self.__numLevelSets)
182        self.__pde.getSolverOptions().setTolerance(tol)
183        self.__pde.setSymmetryOn()
184        self.__pde.setValue(A=self.__pde.createCoefficient('A'), D=self.__pde.createCoefficient('D'), )
185        try:
186          self.__pde.setValue(q=location_of_set_m)
187        except IllegalCoefficientValue:
188          raise ValueError("Unable to set location of fixed level set function.")
189        # ===== check the shape of the weights: =============================
190        if w0 is not None:
191          w0 = interpolate(w0,self.__pde.getFunctionSpaceForCoefficient('D'))
192          s0=w0.getShape()
193          if numLevelSets == 1:
194            if not s0 == () :
195            raise ValueError("Unexpected shape %s for weight w0."%(s0,))
196            else:
197            if not s0 == (numLevelSets,):
198                raise ValueError("Unexpected shape %s for weight w0."%(s0,))
199          if not self.__trafo.isCartesian():
200                w0*=self.__trafo.getVolumeFactor()
201        if not w1 is None:
202              w1 = interpolate(w1,self.__pde.getFunctionSpaceForCoefficient('A'))
203              s1=w1.getShape()
204              if self.__numLevelSets == 1 :
205                if not s1 == (DIM,) :
206                  raise ValueError("Unexpected shape %s for weight w1."%(s1,))
207              else:
208                if not s1 == (self.__numLevelSets,DIM):
209                  raise ValueError("Unexpected shape %s for weight w1."%(s1,))
210              if not self.__trafo.isCartesian():
211                f=self.__trafo.getScalingFactors()**2*self.__trafo.getVolumeFactor()
212                if self.__numLevelSets == 1:
213                  w1*=f
214                else:
215                  for i in range(self.__numLevelSets): w1[i,:]*=f
216    
217            if self.__numLevelSets == 1:
218              wc=None
219            else:
220              wc = interpolate(wc,self.__pde.getFunctionSpaceForCoefficient('A'))
221              sc=wc.getShape()
222              if not sc == (self.__numLevelSets, self.__numLevelSets):
223                raise ValueError("Unexpected shape %s for weight wc."%(sc,))
224              if not self.__trafo.isCartesian():
225                raise ValueError("Non-cartesian coordinates for cross-gradient term is not supported yet.")
226            # ============= now we rescale weights: =============================
227        L2s=np.asarray(boundingBoxEdgeLengths(domain))**2
228        L4=1/np.sum(1/L2s)**2
229        if self.__numLevelSets == 1:
230          A=0
231          if w0 is not None:
232            A = integrate(w0)
233          if w1 is not None:
234            A += integrate(inner(w1, 1/L2s))
235          if A > 0:
236            f = scale/A
237            if w0 is not None:
238              w0*=f
239            if w1 is not None:
240              w1*=f
241          else:
242            raise ValueError("Non-positive weighting factor detected.")
243        else: # numLevelSets > 1
244          for k in range(self.__numLevelSets):
245            A=0
246            if w0 is not None:
247              A = integrate(w0[k])
248            if w1 is not None:
249              A += integrate(inner(w1[k,:], 1/L2s))
250            if A > 0:
251              f = scale[k]/A
252              if w0 is not None:
253                w0[k]*=f
254              if w1 is not None:
255                w1[k,:]*=f
256            else:
257                    raise ValueError("Non-positive weighting factor for level set component %d detected."%k)
258    
259            # and now the cross-gradient:
260            if wc is not None:
261              for l in range(k):
262                A = integrate(wc[l,k])/L4
263                    if A > 0:
264              f = scale_c[l,k]/A
265              wc[l,k]*=f
266    #                       else:
267    #                           raise ValueError("Non-positive weighting factor for cross-gradient level set components %d and %d detected."%(l,k))
268    
269            self.__w0=w0
270        self.__w1=w1
271        self.__wc=wc
272        self.__vol_d=vol(self.__domain)
273    
274      def getDomain(self):      def getDomain(self):
275          """          """
# Line 250  class Regularization(CostFunction): Line 277  class Regularization(CostFunction):
277    
278          :rtype: ``Domain``          :rtype: ``Domain``
279          """          """
280            if self.__domain is None:
281          raise ValueError("This regularization object has not been fully configured yet.")
282          return self.__domain          return self.__domain
283    
284      def getCoordinateTransformation(self):      def getCoordinateTransformation(self):
# Line 258  class Regularization(CostFunction): Line 287  class Regularization(CostFunction):
287    
288          :rtype: `CoordinateTransformation`          :rtype: `CoordinateTransformation`
289          """          """
290            if self.__domain is None:
291          raise ValueError("This regularization object has not been fully configured yet.")        
292          return self.__trafo          return self.__trafo
293    
294      def getNumLevelSets(self):      def getNumLevelSets(self):
# Line 274  class Regularization(CostFunction): Line 305  class Regularization(CostFunction):
305    
306          :rtype: `LinearPDE`          :rtype: `LinearPDE`
307          """          """
308            if self.__domain is None:
309          raise ValueError("This regularization object has not been fully configured yet.")        
310          return self.__pde          return self.__pde
311    
312      def getDualProduct(self, m, r):      def getDualProduct(self, m, r):
# Line 298  class Regularization(CostFunction): Line 331  class Regularization(CostFunction):
331    
332          :rtype: ``int``          :rtype: ``int``
333          """          """
334            if self.__domain is None:
335          raise ValueError("This regularization object has not been fully configured yet.")        
336          return self.__num_tradeoff_factors          return self.__num_tradeoff_factors
337    
338      def setTradeOffFactors(self, mu=None):      def setTradeOffFactors(self, mu=None):
# Line 314  class Regularization(CostFunction): Line 349  class Regularization(CostFunction):
349                     Values must be positive.                     Values must be positive.
350          :type mu: ``list`` of ``float`` or ```numpy.array```          :type mu: ``list`` of ``float`` or ```numpy.array```
351          """          """
352            if self.__domain is None:
353          raise ValueError("This regularization object has not been fully configured yet.")        
354          numLS=self.getNumLevelSets()          numLS=self.getNumLevelSets()
355          numTF=self.getNumTradeOffFactors()          numTF=self.getNumTradeOffFactors()
356          if mu is None:          if mu is None:
# Line 340  class Regularization(CostFunction): Line 377  class Regularization(CostFunction):
377          :param mu: new values for the trade-off factors. Values must be positive.          :param mu: new values for the trade-off factors. Values must be positive.
378          :type mu: ``float``, ``list`` of ``float`` or ```numpy.array```          :type mu: ``float``, ``list`` of ``float`` or ```numpy.array```
379          """          """
380            if self.__domain is None:
381          raise ValueError("This regularization object has not been fully configured yet.")        
382          numLS=self.getNumLevelSets()          numLS=self.getNumLevelSets()
383          if mu is None:          if mu is None:
384              if numLS == 1:              if numLS == 1:
# Line 381  class Regularization(CostFunction): Line 420  class Regularization(CostFunction):
420                       are used. Only value mu_c[l,k] for l<k are used.                       are used. Only value mu_c[l,k] for l<k are used.
421          :type mu_c: ``float``, ``list`` of ``float`` or ``numpy.array``          :type mu_c: ``float``, ``list`` of ``float`` or ``numpy.array``
422          """          """
423            if self.__domain is None:
424          raise ValueError("This regularization object has not been fully configured yet.")        
425          numLS=self.getNumLevelSets()          numLS=self.getNumLevelSets()
426          if mu_c is None or numLS < 2:          if mu_c is None or numLS < 2:
427              self.__mu_c = np.ones((numLS,numLS))              self.__mu_c = np.ones((numLS,numLS))
# Line 401  class Regularization(CostFunction): Line 442  class Regularization(CostFunction):
442      def getArguments(self, m):      def getArguments(self, m):
443          """          """
444          """          """
445            if self.__domain is None:
446          raise ValueError("This regularization object has not been fully configured yet.")        
447          return grad(m),          return grad(m),
448    
449      def getValue(self, m, grad_m):      def getValue(self, m, grad_m):
# Line 410  class Regularization(CostFunction): Line 453  class Regularization(CostFunction):
453    
454          :rtype: ``float``          :rtype: ``float``
455          """          """
456            if self.__domain is None:
457          raise ValueError("This regularization object has not been fully configured yet.")        
458          mu=self.__mu          mu=self.__mu
459          mu_c=self.__mu_c          mu_c=self.__mu_c
460          DIM=self.getDomain().getDim()          DIM=self.getDomain().getDim()
# Line 449  class Regularization(CostFunction): Line 494  class Regularization(CostFunction):
494    
495          :note: This implementation returns Y_k=dPsi/dm_k and X_kj=dPsi/dm_kj          :note: This implementation returns Y_k=dPsi/dm_k and X_kj=dPsi/dm_kj
496          """          """
497            if self.__domain is None:
498          raise ValueError("This regularization object has not been fully configured yet.")
499          mu=self.__mu          mu=self.__mu
500          mu_c=self.__mu_c          mu_c=self.__mu_c
501          DIM=self.getDomain().getDim()          DIM=self.getDomain().getDim()
# Line 493  class Regularization(CostFunction): Line 539  class Regularization(CostFunction):
539      def getInverseHessianApproximation(self, m, r, grad_m, solve=True):      def getInverseHessianApproximation(self, m, r, grad_m, solve=True):
540          """          """
541          """          """
542            if self.__domain is None:
543          raise ValueError("This regularization object has not been fully configured yet.")        
544          if self._new_mu or self._update_Hessian:          if self._new_mu or self._update_Hessian:
545              self._new_mu=False              self._new_mu=False
546              self._update_Hessian=False              self._update_Hessian=False
# Line 556  class Regularization(CostFunction): Line 604  class Regularization(CostFunction):
604          """          """
605          notifies the class to recalculate the Hessian operator.          notifies the class to recalculate the Hessian operator.
606          """          """
607            if self.__domain is None:
608          raise ValueError("This regularization object has not been fully configured yet.")        
609          if not self.__useDiagonalHessianApproximation:          if not self.__useDiagonalHessianApproximation:
610              self._update_Hessian=True              self._update_Hessian=True
611    
# Line 567  class Regularization(CostFunction): Line 617  class Regularization(CostFunction):
617          :type m: `Data`          :type m: `Data`
618          :rtype: ``float``          :rtype: ``float``
619          """          """
620            if self.__domain is None:
621          raise ValueError("This regularization object has not been fully configured yet.")        
622          return sqrt(integrate(length(m)**2)/self.__vol_d)          return sqrt(integrate(length(m)**2)/self.__vol_d)
623    

Legend:
Removed from v.5504  
changed lines
  Added in v.5505

  ViewVC Help
Powered by ViewVC 1.1.26