/[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 4094 by jfenwick, Mon Nov 19 01:45:38 2012 UTC revision 4095 by caltinay, Wed Dec 5 05:32:22 2012 UTC
# Line 50  class Regularization(CostFunction): Line 50  class Regularization(CostFunction):
50                         location_of_set_m=Data(),                         location_of_set_m=Data(),
51                         useDiagonalHessianApproximation=True, tol=1e-8):                         useDiagonalHessianApproximation=True, tol=1e-8):
52          """          """
53          initialization          initialization.
54                    
55          :param domain: domain          :param domain: domain
56          :type domain: `Domain`          :type domain: `Domain`
57          :param numLevelSets: number of level sets          :param numLevelSets: number of level sets
58          :type numLevelSets: ``int``          :type numLevelSets: ``int``
59          :param s0: scaling factor for the m**2 term. If not set zero is assumed.          :param s0: scaling factor for the m**2 term. If not set zero is assumed.
60          :type s0: ``Scalar`` if ``numLevelSets``==1 or `Data` object of shape ('numLevelSets`,) if numLevelSets > 1          :type s0: ``Scalar`` if ``numLevelSets`` == 1 or `Data` object of shape (``numLevelSets`` ,) if ``numLevelSets`` > 1
61          :param s1: scaling factor for the grad(m_i) terms. If not set zero is assumed.          :param s1: scaling factor for the grad(m_i) terms. If not set zero is assumed.
62          :type s1: ``Vector`` if ``numLevelSets``==1 or `Data` object of shape (`numLevelSets`,DIM) if numLevelSets > 1.          :type s1: ``Vector`` if ``numLevelSets`` == 1 or `Data` object of shape (``numLevelSets`` , DIM) if ``numLevelSets`` > 1.
63          :param sc: scaling factor for the cross correlation terms. If not set zero is assumed. Used for the case if numLevelSets > 1 only.          :param sc: scaling factor for the cross correlation terms. If not set zero is assumed. Used for the case if ``numLevelSets`` > 1 only.
64                     values ``sc[l,k]``  in the lower triangle (l<k) are used only.                     values ``sc[l,k]``  in the lower triangle (l<k) are used only.
65          :type sc: `Data` object of shape (`numLevelSets`,`numLevelSets`)              :type sc: `Data` object of shape (``numLevelSets`` , ``numLevelSets``)
66          :param location_of_set_m: marks location of zero values of the level set function ``m`` by a positive entry.          :param location_of_set_m: marks location of zero values of the level set function ``m`` by a positive entry.
67          :type location_of_set_m: ``Scalar`` if ``numLevelSets``==1 or `Data` object of shape (``numLevelSets``,) if numLevelSets > 1.          :type location_of_set_m: ``Scalar`` if ``numLevelSets`` == 1 or `Data` object of shape (``numLevelSets`` ,) if ``numLevelSets`` > 1.
68          :param useDiagonalHessianApproximation: if True cross correllation terms between level set components are ignored when calculating          :param useDiagonalHessianApproximation: if True cross correllation terms between level set components are ignored when calculating
69                                                  approximations of the inverse of the Hessian Operator. This can speep-up the calculation of                                                  approximations of the inverse of the Hessian Operator. This can speep-up the calculation of
70                                                  the inverse but may lead to an increase of the number of iteration steps in the inversion.                                                  the inverse but may lead to an increase of the number of iteration steps in the inversion.
71          :type useDiagonalHessianApproximation: ``bool``          :type useDiagonalHessianApproximation: ``bool``
72          :param tol: toleramce when solving the PDE for the inverse of the Hessian Operator          :param tol: toleramce when solving the PDE for the inverse of the Hessian Operator
73          :type tol: positive ``float``          :type tol: positive ``float``
74            
           
75          """          """
76          if numLevelSets>1:          if numLevelSets>1:
77            raise ValueError("Currently only numLevelSets<=1 is supportered.")                raise ValueError("Currently only numLevelSets<=1 is supportered.")
78      if s0 == None and s1==None:          if s0 == None and s1==None:
79            raise ValueError("Values for s0 or s1 must be given.")                raise ValueError("Values for s0 or s1 must be given.")
80                
81          self.__domain=domain          self.__domain=domain
82          DIM=self.__domain.getDim()          DIM=self.__domain.getDim()
83          self.__L2=np.asarray(boundingBoxEdgeLengths(domain))**2          self.__L2=np.asarray(boundingBoxEdgeLengths(domain))**2
# Line 88  class Regularization(CostFunction): Line 87  class Regularization(CostFunction):
87          if self.__numLevelSets > 1:          if self.__numLevelSets > 1:
88              self.__useDiagonalHessianApproximation=useDiagonalHessianApproximation              self.__useDiagonalHessianApproximation=useDiagonalHessianApproximation
89          else:          else:
90          self.__useDiagonalHessianApproximation=True              self.__useDiagonalHessianApproximation=True
91      self._update_Hessian=True          self._update_Hessian=True
92                
93      self.__pde=LinearPDE(self.__domain, numEquations=self.__numLevelSets)          self.__pde=LinearPDE(self.__domain, numEquations=self.__numLevelSets)
94          self.__pde.getSolverOptions().setTolerance(tol)          self.__pde.getSolverOptions().setTolerance(tol)
95          self.__pde.setSymmetryOn()          self.__pde.setSymmetryOn()
96          self.__pde_is_set=False          self.__pde_is_set=False
97          try:          try:
98        self.__pde.setValue(q=location_of_set_m)            self.__pde.setValue(q=location_of_set_m)
99      except IllegalCoefficientValue:          except IllegalCoefficientValue:
100        raise ValueError("Unable to set location of fixed level set function.")            raise ValueError("Unable to set location of fixed level set function.")
101                
102      self.__total_num_weights=2*numLevelSets+((numLevelSets-1)*numLevelSets)/2          self.__total_num_weights=2*numLevelSets+((numLevelSets-1)*numLevelSets)/2
103      self.__weight_index=[]  # this is a mapping from the relevant mu-coefficients to the set of all mu-coefficients          self.__weight_index=[]  # this is a mapping from the relevant mu-coefficients to the set of all mu-coefficients
104                         # we count s0, then s1, then sc (k<l).                             # we count s0, then s1, then sc (k<l).
105      # THE S0 weighting factor          # THE S0 weighting factor
106      n=0          n=0
107      VV=vol(domain)          VV=vol(domain)
108          if not s0 is None:          if not s0 is None:
109          s0 = interpolate(s0,self.__pde.getFunctionSpaceForCoefficient('D'))                  s0 = interpolate(s0,self.__pde.getFunctionSpaceForCoefficient('D'))    
110          s=s0.getShape()              s=s0.getShape()
111          if numLevelSets == 1 :              if numLevelSets == 1 :
112               if s == () :                   if s == () :
113               V=integrate(s0)                       V=integrate(s0)
114               if V > 0:                       if V > 0:
115                 self.__weight_index.append(n)                         self.__weight_index.append(n)
116                 s0*=VV/V                         s0*=VV/V
117                   else:                       else:
118                 s0=None                         s0=None
119           else:                   else:
120               raise ValueError("Unexpected shape %s for weight s0."%s)                       raise ValueError("Unexpected shape %s for weight s0."%s)
121              else:              else:
122               if s == (numLevelSets,):                   if s == (numLevelSets,):
123               for k in xrange(numLevelSets):                       for k in xrange(numLevelSets):
124                  V=integrate(s0[k])                          V=integrate(s0[k])
125                  if V > 0:                          if V > 0:
126                      self.__weight_index.append(n+k)                              self.__weight_index.append(n+k)
127                      s0[k]*=VV/V                              s0[k]*=VV/V
128           else:                   else:
129               raise ValueError("Unexpected shape %s for weight s0."%s)                         raise ValueError("Unexpected shape %s for weight s0."%s)  
130          self.__s0=s0          self.__s0=s0
131          n+=numLevelSets          n+=numLevelSets
132                    
133          # The S1 weighting factor          # The S1 weighting factor
134          if not s1 is None:          if not s1 is None:
135          s1 = interpolate(s1,self.__pde.getFunctionSpaceForCoefficient('A'))              s1 = interpolate(s1,self.__pde.getFunctionSpaceForCoefficient('A'))
136          s=s1.getShape()              s=s1.getShape()
137                        
138          if numLevelSets == 1 :              if numLevelSets == 1 :
139               if s == (DIM,) :                   if s == (DIM,) :
140                 V=integrate(inner(s1, 1/self.__L2))                         V=integrate(inner(s1, 1/self.__L2))
141                 if V > 0:                         if V > 0:
142                self.__weight_index.append(n)                            self.__weight_index.append(n)
143                    s1*=VV/V                            s1*=VV/V
144                 print "REG SCALE = ",s1                         print "REG SCALE = ",s1
145           else:                   else:
146               raise ValueError("Unexpected shape %s for weight s1."%s)                       raise ValueError("Unexpected shape %s for weight s1."%s)
147              else:              else:
148               if s == (numLevelSets,DIM):                   if s == (numLevelSets,DIM):
149               for k in xrange(numLevelSets):                       for k in xrange(numLevelSets):
150                  for i in xrange(DIM):                          for i in xrange(DIM):
151                  ww=s1[k,:]                              ww=s1[k,:]
152                  V=integrate(inner(ww,1/self.__L2))                              V=integrate(inner(ww,1/self.__L2))
153                      if V > 0:                              if V > 0:
154                         self.__weight_index.append(n+i)                                 self.__weight_index.append(n+i)
155                         s1[k,:]=ww*(VV/V)                                 s1[k,:]=ww*(VV/V)
156           else:                   else:
157               raise ValueError("Unexpected shape %s for weight s1."%s)                             raise ValueError("Unexpected shape %s for weight s1."%s)      
158                                
159          self.__s1=s1            self.__s1=s1  
160          n+=numLevelSets          n+=numLevelSets
161    
162          # The SC weighting factor          # The SC weighting factor
163          if not sc is None:          if not sc is None:
164          if numLevelSets == 1 :              if numLevelSets == 1 :
165                sc=None                sc=None
166              else:              else:
167                sc = interpolate(sc,self.__pde.getFunctionSpaceForCoefficient('A'))                sc = interpolate(sc,self.__pde.getFunctionSpaceForCoefficient('A'))
168            s=sc.getShape()                s=sc.getShape()
169            if s == (numLevelSets,numLevelSets):                if s == (numLevelSets,numLevelSets):
170               for k in xrange(numLevelSets):                       for k in xrange(numLevelSets):
171                  sc[k,k]=0.                          sc[k,k]=0.
172                  for l in xrange(k):                          for l in xrange(k):
173                  ww=sc[l,k]                              ww=sc[l,k]
174                  V=integrate(ww)                              V=integrate(ww)
175                      if V > 0:                              if V > 0:
176                         self.__weight_index.append(n+k*numLevelSets+l)                                 self.__weight_index.append(n+k*numLevelSets+l)
177                         sc[l,k]=ww*VV/V*self.__L4                                 sc[l,k]=ww*VV/V*self.__L4
178                         sc[k,l]=0                                 sc[k,l]=0
179                else:                else:
180               raise ValueError("Unexpected shape %s for weight s0."%s)                             raise ValueError("Unexpected shape %s for weight s0."%s)      
181                                
182          self.__sc=sc          self.__sc=sc
183      self.setWeights()          self.setWeights()
184                        
185      def getDomain(self):      def getDomain(self):
186          """          """
187          return the domain of the regularization term          return the domain of the regularization term
# Line 217  class Regularization(CostFunction): Line 216  class Regularization(CostFunction):
216          A=0          A=0
217          if not r[0].isEmpty(): A+=integrate(inner(r[0], m))          if not r[0].isEmpty(): A+=integrate(inner(r[0], m))
218          if not r[1].isEmpty(): A+=integrate(inner(r[1], grad(m)))          if not r[1].isEmpty(): A+=integrate(inner(r[1], grad(m)))
219      return A          return A
220    
221      def getValue(self, m, grad_m):      def getValue(self, m, grad_m):
222          """          """
# Line 239  class Regularization(CostFunction): Line 238  class Regularization(CostFunction):
238          n+=numLS          n+=numLS
239                    
240          if self.__s1 is not None:          if self.__s1 is not None:
241          if numLS == 1:              if numLS == 1:
242              A+=integrate(inner(grad_m**2, self.__s1))*mu[n]                  A+=integrate(inner(grad_m**2, self.__s1))*mu[n]
243          else:              else:
244              for k in xrange(numLS):                  for k in xrange(numLS):
245              A+=mu[n+k]*integrate(inner(grad_m[k,:]**2,self.__s1[k,:]))                      A+=mu[n+k]*integrate(inner(grad_m[k,:]**2,self.__s1[k,:]))
246          n+=numLS              n+=numLS    
247                    
248          if self.__sc is not None:          if self.__sc is not None:
249            for k in xrange(numLS):                for k in xrange(numLS):
250             gk=grad_m[k,:]                     gk=grad_m[k,:]
251             len_gk=length(gk)                     len_gk=length(gk)
252             for l in xrange(k):                     for l in xrange(k):
253                 gl=grad_m[l,:]                         gl=grad_m[l,:]
254                 A+= mu[n+k*numLS+l] * integrate( self.__sc[l,k] * ( len_gk * length(gl) )**2 - inner(gk, gl)**2 )                         A+= mu[n+k*numLS+l] * integrate( self.__sc[l,k] * ( len_gk * length(gl) )**2 - inner(gk, gl)**2 )
255          return A/2          return A/2
256                    
257      def getGradient(self, m,  grad_m):      def getGradient(self, m,  grad_m):
# Line 268  class Regularization(CostFunction): Line 267  class Regularization(CostFunction):
267          n=0          n=0
268                    
269          if self.__s0 is not None:          if self.__s0 is not None:
270          Y = m * self.__s0 * mu[:numLS]              Y = m * self.__s0 * mu[:numLS]
271      else:          else:
272          Y = Data()              Y = Data()
273          n+=numLS          n+=numLS
274    
275          if self.__s1 is not None:          if self.__s1 is not None:
276          if numLS == 1:              if numLS == 1:
277              X=grad_m* self.__s1*mu[n]                  X=grad_m* self.__s1*mu[n]
278          else:              else:
279              X=self.getPDE().createCoefficient("X")                  X=self.getPDE().createCoefficient("X")
280              for k in xrange(numLS):                  for k in xrange(numLS):
281              X[k,:]=mu[n+k]*grad_m[k,:]*self.__s1[k,:]                      X[k,:]=mu[n+k]*grad_m[k,:]*self.__s1[k,:]
282      else:          else:
283          X = Data()              X = Data()
284          n+=numLS          n+=numLS
285          if self.__sc is not None:          if self.__sc is not None:
286             raise NotImplementedError             raise NotImplementedError
# Line 296  class Regularization(CostFunction): Line 295  class Regularization(CostFunction):
295          """          """
296          """          """
297          if self._new_mu or self._update_Hessian:          if self._new_mu or self._update_Hessian:
298                        
299          self._new_mu=False              self._new_mu=False
300          self._update_Hessian=False              self._update_Hessian=False
301                        
302          mu=self.getWeights( uncompress=True)              mu=self.getWeights( uncompress=True)
303              DIM=self.getDomain().getDim()              DIM=self.getDomain().getDim()
304              numLS=self.getNumLevelSets()              numLS=self.getNumLevelSets()
305              n=0              n=0
306              if self.__s0 is not None:              if self.__s0 is not None:
307              if numLS == 1:                  if numLS == 1:
308                   D=self.__s0 * mu[n]                       D=self.__s0 * mu[n]
309                  else:                    else:  
310                       D=self.getPDE().createCoefficient("D")                       D=self.getPDE().createCoefficient("D")
311                   for k in xrange(numLS): D[k,k]=self.__s0[k] * mu[n+k]                       for k in xrange(numLS): D[k,k]=self.__s0[k] * mu[n+k]
312              self.getPDE().setValue(D=D)                  self.getPDE().setValue(D=D)
313          n+=numLS              n+=numLS
314    
315          A=self.getPDE().createCoefficient("A")              A=self.getPDE().createCoefficient("A")
316              if self.__s1 is not None:              if self.__s1 is not None:
317             if numLS == 1:                 if numLS == 1:
318             for i in xrange(DIM): A[i,i]=self.__s1[i] * mu[n]                     for i in xrange(DIM): A[i,i]=self.__s1[i] * mu[n]
319             else:                 else:
320                 for k in xrange(numLS):                     for k in xrange(numLS):
321                      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.__s1[k,i] * mu[n+k]
322              n+=numLS              n+=numLS
323              if self.__sc is not None:              if self.__sc is not None:
324                  raise NotImplementedError                  raise NotImplementedError
# Line 340  class Regularization(CostFunction): Line 339  class Regularization(CostFunction):
339          """          """
340          notify the class to recalculate the Hessian operator          notify the class to recalculate the Hessian operator
341          """          """
342      if not self.__useDiagonalHessianApproximation:          if not self.__useDiagonalHessianApproximation:
343          self._update_Hessian=True              self._update_Hessian=True
344    
345     # ================ we should factor these out =============================================================     # ================ we should factor these out =============================================================
346      def getNumRelevantTerms(self):      def getNumRelevantTerms(self):
# Line 367  class Regularization(CostFunction): Line 366  class Regularization(CostFunction):
366          :type mu: ``list`` of ``float`` or ``np,array``          :type mu: ``list`` of ``float`` or ``np,array``
367          """          """
368          if mu == None:          if mu == None:
369          mu = np.ones(self.getNumRelevantTerms())              mu = np.ones(self.getNumRelevantTerms())
370      mu=np.asarray(mu)          mu=np.asarray(mu)
371                
372                
373      if len(mu) == self.getNumRelevantTerms():          if len(mu) == self.getNumRelevantTerms():
374           if not mu.shape == (self.getNumRelevantTerms(),):               if not mu.shape == (self.getNumRelevantTerms(),):
375          raise ValueError("%s values are expected."%self.getNumRelevantTerms())                  raise ValueError("%s values are expected."%self.getNumRelevantTerms())
376           self.__mu=mu               self.__mu=mu
377      else:          else:
378           if not mu.shape == (self.__total_num_weights,):               if not mu.shape == (self.__total_num_weights,):
379          raise ValueError("%s values are expected."%self.__total_num_weights)                  raise ValueError("%s values are expected."%self.__total_num_weights)
380           self.__mu = np.zeros(self.getNumRelevantTerms())               self.__mu = np.zeros(self.getNumRelevantTerms())
381           for i in xrange(len(self.__weight_index)): self.__mu[i]=mu[self.__weight_index[i]]               for i in xrange(len(self.__weight_index)): self.__mu[i]=mu[self.__weight_index[i]]
382      self._new_mu=True          self._new_mu=True
383    
384      def setWeightsForS0(self, mu=None):      def setWeightsForS0(self, mu=None):
385           """           """
# Line 391  class Regularization(CostFunction): Line 390  class Regularization(CostFunction):
390           if mu is None:           if mu is None:
391               my_mu[:numLS]=1               my_mu[:numLS]=1
392           else:           else:
393           my_mu[:numLS]=mu               my_mu[:numLS]=mu
394           self.setWeights(my_mu)           self.setWeights(my_mu)
395                    
396      def setWeightsForS1(self, mu=None):      def setWeightsForS1(self, mu=None):
397           """           """
398           set the weights for s1-terms           set the weights for s1-terms
# Line 427  class Regularization(CostFunction): Line 426  class Regularization(CostFunction):
426    
427          """          """
428          if uncompress:          if uncompress:
429           mu = np.zeros(self.__total_num_weights)               mu = np.zeros(self.__total_num_weights)
430           for i in xrange(len(self.__weight_index)): mu[self.__weight_index[i]] = self.__mu[i]               for i in xrange(len(self.__weight_index)): mu[self.__weight_index[i]] = self.__mu[i]
431           return mu               return mu
432      else:          else:
433             return self.__mu             return self.__mu
434            
435      def getWeightIndex(self):      def getWeightIndex(self):
# Line 438  class Regularization(CostFunction): Line 437  class Regularization(CostFunction):
437          returns an iondex to the contributions of terms with non-zero scaling factor.          returns an iondex to the contributions of terms with non-zero scaling factor.
438          """          """
439          return self.__weight_index          return self.__weight_index
440                
Legend:
Removed from v.4094  
changed lines
  Added in v.4095
  ViewVC Help
Powered by ViewVC 1.1.26