downunder/py_src/regularizations.py


##############################################################################
#
# Copyright (c) 2003-2012 by University of Queensland
# http://www.uq.edu.au
#
# Primary Business: Queensland, Australia
# Licensed under the Open Software License version 3.0
# http://www.opensource.org/licenses/osl-3.0.php
#
# Development until 2012 by Earth Systems Science Computational Center (ESSCC)
# Development since 2012 by School of Earth Sciences
#
##############################################################################

__copyright__="""Copyright (c) 2003-2012 by University of Queensland
http://www.uq.edu.au
Primary Business: Queensland, Australia"""
__license__="""Licensed under the Open Software License version 3.0
http://www.opensource.org/licenses/osl-3.0.php"""
__url__="https://launchpad.net/escript-finley"

__all__ = ['Regularization']

from costfunctions import CostFunction


import numpy as np
from esys.escript.linearPDEs import LinearPDE, IllegalCoefficientValue
from esys.escript import Data, grad, inner, integrate, kronecker, boundingBoxEdgeLengths, interpolate, vol
from esys.escript.pdetools import ArithmeticTuple

class Regularization(CostFunction):
    """
    The regularization term for the level set function `m` within the cost function J for an inversion:
    
    *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*
    
    where s0[k], s1[k,i] and  sc[k,l] are non-negative scaling factors and mu_0[k], mu_1[k], mu_c[l,k] are weighting factors
    which may be alter during the inversion. The scaling factors are normalized such that their integrals over the
    domain are constant:
   
    *integrate(s0[k])=1*
    *integrate(inner(s1[k,:],L[:])=1*
    *integrate(inner(sc[l,k]*L**4)=1*
    
    """
    def __init__(self, domain, numLevelSets=1,  
                       s0=None, s1=None, sc=None,
                       location_of_set_m=Data(), 
                       useDiagonalHessianApproximation=True, tol=1e-8):
        """
        initialization
        
        :param domain: domain 
        :type domain: `Domain`
        :param numLevelSets: number of level sets
        :type numLevelSets: ``int``
        :param s0: scaling factor for the m**2 term. If not set zero is assumed.
        :type s0: ``Scalar`` if `numLevelSets`==1 or `Data` object of shape ('numLevelSets`,) if numLevelSets > 1
        :param s1: scaling factor for the grad(m_i) terms. If not set zero is assumed.
        :type s1: ``Vector`` if `numLevelSets`==1 or `Data` object of shape (`numLevelSets`,DIM) if numLevelSets > 1.
        :param sc: scaling factor for the cross correlation terms. If not set zero is assumed. Used for the case if numLevelSets > 1 only.
                   values `sc[l,k]``  in the lower triangle (l<k) are used only.
        :type sc: `Data` object of shape (`numLevelSets`,`numLevelSets`)     
        :param location_of_set_m: marks location of zero values of the level set function `m` by a positive entry.
        :type location_of_set_m: ``Scalar`` if `numLevelSets`==1 or `Data` object of shape ('numLevelSets`,) if numLevelSets > 1. 
        :param useDiagonalHessianApproximation: if True cross correllation terms between level set components are ignored when calculating 
                                                approximations of the inverse of the Hessian Operator. This can speep-up the calculation of 
                                                the inverse but may lead to an increase of the number of iteration steps in the inversion.
        :type useDiagonalHessianApproximation: ``bool``
        :param tol: toleramce when solving the PDE for the inverse of the Hessian Operator
        :type tol: positive ``float``
        
        
        """
        if numLevelSets>1:
          raise ValueError("Currently only numLevelSets<=1 is supportered.")
    if s0 == None and s1==None:
          raise ValueError("Values for s0 or s1 must be given.")
    
        self.__domain=domain
        DIM=self.__domain.getDim()
        self.__L2=np.asarray(boundingBoxEdgeLengths(domain))**2
        self.__L4=np.sum(self.__L2)**2
        self.__numLevelSets=numLevelSets
        
        if self.__numLevelSets > 1:
            self.__useDiagonalHessianApproximation=useDiagonalHessianApproximation
        else:
        self.__useDiagonalHessianApproximation=True 
    self._update_Hessian=True
    
    self.__pde=LinearPDE(self.__domain, numEquations=self.__numLevelSets)
        self.__pde.getSolverOptions().setTolerance(tol)
        self.__pde.setSymmetryOn()
        self.__pde_is_set=False
        try:
      self.__pde.setValue(q=location_of_set_m)
    except IllegalCoefficientValue: 
      raise ValueError("Unable to set location of fixed level set function.")
    
    self.__total_num_weights=2*numLevelSets+((numLevelSets-1)*numLevelSets)/2
    self.__weight_index=[]  # this is a mapping from the relevant mu-coefficients to the set of all mu-coefficients
                       # we count s0, then s1, then sc (k<l).
    # THE S0 weighting factor
    n=0
    VV=vol(domain)
        if not s0 is None:
        s0 = interpolate(s0,self.__pde.getFunctionSpaceForCoefficient('D'))     
        s=s0.getShape()
        if numLevelSets == 1 :
             if s == () :
             V=integrate(s0)
             if V > 0:
               self.__weight_index.append(n)
               s0*=VV/V
                 else:
               s0=None
         else:
             raise ValueError("Unexpected shape %s for weight s0."%s)
            else:
             if s == (numLevelSets,):
             for k in xrange(numLevelSets):
                V=integrate(s0[k])
                if V > 0: 
                    self.__weight_index.append(n+k)
                    s0[k]*=VV/V
         else:
             raise ValueError("Unexpected shape %s for weight s0."%s)   
        self.__s0=s0
        n+=numLevelSets
        
        # The S1 weighting factor
        if not s1 is None:
        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)
            else:
             if s == (numLevelSets,DIM):
             for k in xrange(numLevelSets):
                for i in xrange(DIM):
                ww=s1[k,:]
                V=integrate(inner(ww,1/self.__L2))
                    if V > 0: 
                       self.__weight_index.append(n+i)
                       s1[k,:]=ww*(VV/V)
         else:
             raise ValueError("Unexpected shape %s for weight s1."%s)      
           
        self.__s1=s1   
        n+=numLevelSets

        # The SC weighting factor
        if not sc is None:
        if numLevelSets == 1 :
              sc=None
            else:
              sc = interpolate(sc,self.__pde.getFunctionSpaceForCoefficient('A'))
          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)      
           
        self.__sc=sc
    self.setWeights()
        
    def getDomain(self):
        """
        return the domain of the regularization term
        :rtype: ``Domain``
        """
        return self.__domain
        
    def getNumLevelSets(self):
        """
        return the number of level set functions
        :rtype: ``int``
        """
        return self.__numLevelSets

    def getPDE(self):
        """
        returns the linear PDE to be solved for the Hessian Operator inverse
        :rtype: ``LinearPDE``
        """
        return self.__pde
     
    def getDualProduct(self, m, r):
        """
        returns the dual product of a gradient represented by X=r[1] and Y=r[0] with a level set function m:
             
             *Y_i*m_i + X_ij*m_{i,j}*
        
        :type m: ``esys.escript.Data``
        :type r: ``ArithmeticTuple``
        :rtype: float
        """
        A=0
        if not r[0].isEmpty(): A+=integrate(inner(r[0], m))
        if not r[1].isEmpty(): A+=integrate(inner(r[1], grad(m)))
    return A

    def getValue(self, m, grad_m):
        """
        return the value of the costfunction J
        
        
        :rtype: ``float``
        """
        mu=self.getWeights( uncompress=True)
        DIM=self.getDomain().getDim()
        numLS=self.getNumLevelSets()
                
        A=0
        n=0
        
        if self.__s0 is not None:
            A+=inner(integrate(m**2*self.__s0), mu[:numLS])
        n+=numLS
        
        if self.__s1 is not None:
        if numLS == 1:
            A+=integrate(inner(grad_m**2, self.__s1))*mu[n]
        else:
            for k in xrange(numLS):
            A+=mu[n+k]*integrate(inner(grad_m[k,:]**2,self.__s1[k,:]))
        n+=numLS    
        
        if self.__sc is not None:
          for k in xrange(numLS):
           gk=grad_m[k,:]
           len_gk=length(gk)
           for l in xrange(k):
               gl=grad_m[l,:]
               A+= mu[n+k*numLS+l] * integrate( self.__sc[l,k] * ( len_gk * length(gl) )**2 - inner(gk, gl)**2 ) 
        return A/2
        
    def getGradient(self, m,  grad_m):
        """
        returns the gradient of the costfunction J with respect to m. The function returns Y_k=dPsi/dm_k and
        X_kj=dPsi/dm_kj
        """
        
        mu=self.getWeights( uncompress=True)
        DIM=self.getDomain().getDim()
        numLS=self.getNumLevelSets()
        
        n=0
        
        if self.__s0 is not None:
        Y = m * self.__s0 * mu[:numLS]
    else:
        Y = Data()
        n+=numLS

        if self.__s1 is not None:
        if numLS == 1:
            X=grad_m* self.__s1*mu[n]
        else:
            X=self.getPDE().createCoefficient("X")
            for k in xrange(numLS):
            X[k,:]=mu[n+k]*grad_m[k,:]*self.__s1[k,:]
    else:
        X = Data()
        n+=numLS 
        if self.__sc is not None:
           raise NotImplementedError
        return ArithmeticTuple(Y, X)
        
    def getArguments(self, m):
        """
        """
        return ( grad(m),)
         
    def getInverseHessianApproximation(self, m, r, grad_m):
        """
        """
        if self._new_mu or self._update_Hessian:
        
        self._new_mu=False
        self._update_Hessian=False
        
        mu=self.getWeights( uncompress=True)
            DIM=self.getDomain().getDim()
            numLS=self.getNumLevelSets()
            n=0
            if self.__s0 is not None:
            if numLS == 1:
                 D=self.__s0 * mu[n]
                else:   
                     D=self.getPDE().createCoefficient("D")
                 for k in xrange(numLS): D[k,k]=self.__s0[k] * mu[n+k]
            self.getPDE().setValue(D=D)
        n+=numLS

        A=self.getPDE().createCoefficient("A")
            if self.__s1 is not None:
           if numLS == 1:
           for i in xrange(DIM): A[i,i]=self.__s1[i] * mu[n]
           else:
               for k in xrange(numLS): 
                    for i in xrange(DIM): A[k,i,k,i]=self.__s1[k,i] * mu[n+k]
            n+=numLS 
            if self.__sc is not None:
                raise NotImplementedError
            print A
            self.getPDE().setValue(A=A)
        self.getPDE().resetRightHandSideCoefficients()
        self.getPDE().setValue(X=r[1])
        print "X only: ",self.getPDE().getSolution()
        self.getPDE().resetRightHandSideCoefficients()
        self.getPDE().setValue(Y=r[0])
        print "Y only: ",self.getPDE().getSolution()
        
        self.getPDE().resetRightHandSideCoefficients()
        self.getPDE().setValue(X=r[1], Y=r[0])
        return self.getPDE().getSolution()
        
    def updateHessian(self):
        """
        notify the class to recalculate the Hessian operator
        """
    if not self.__useDiagonalHessianApproximation:
        self._update_Hessian=True

   # ================ we should factor these out =============================================================
    def getNumRelevantTerms(self):
        """
        returns the number of terms in the costfunction of regularization with non-zero scaling factors
        :rtype: ``int``
        """
        return len(self.__weight_index)
    
    def getNumTerms(self):
        """
        returns the number of terms in the costfunction 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 costsfunction. Note that only values
        correspond 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):
         """
         set 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):
         """
         set 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): 
         """
         set the weights for s1-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 costsfunction. 
        The first ``numLevelSets`` values used for the 
        regularization terms and the remaining values for the cross correlation terms.
 
       :type mu: ``list`` of ``float``
        """
        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 iondex to the contributions of terms with non-zero scaling factor.
        """
        return self.__weight_index
     
1
2	##############################################################################
3	#
4	# Copyright (c) 2003-2012 by University of Queensland
5	# http://www.uq.edu.au
6	#
7	# Primary Business: Queensland, Australia
8	# Licensed under the Open Software License version 3.0
9	# http://www.opensource.org/licenses/osl-3.0.php
10	#
11	# Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12	# Development since 2012 by School of Earth Sciences
13	#
14	##############################################################################
15
16	__copyright__="""Copyright (c) 2003-2012 by University of Queensland
17	http://www.uq.edu.au
18	Primary Business: Queensland, Australia"""
19	__license__="""Licensed under the Open Software License version 3.0
20	http://www.opensource.org/licenses/osl-3.0.php"""
21	__url__="https://launchpad.net/escript-finley"
22
23	__all__ = ['Regularization']
24
25	from costfunctions import CostFunction
26
27
28	import numpy as np
29	from esys.escript.linearPDEs import LinearPDE, IllegalCoefficientValue
30	from esys.escript import Data, grad, inner, integrate, kronecker, boundingBoxEdgeLengths, interpolate, vol
31	from esys.escript.pdetools import ArithmeticTuple
32
33	class Regularization(CostFunction):
34	"""
35	The regularization term for the level set function `m` within the cost 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*
38
39	where s0[k], s1[k,i] and sc[k,l] are non-negative scaling factors and mu_0[k], mu_1[k], mu_c[l,k] are weighting factors
40	which may be alter during the inversion. The scaling factors are normalized such that their integrals over the
41	domain are constant:
42
43	integrate(s0[k])=1
44	integrate(inner(s1[k,:],L[:])=1
45	integrate(inner(sc[l,k]L*4)=1
46
47	"""
48	def __init__(self, domain, numLevelSets=1,
49	s0=None, s1=None, sc=None,
50	location_of_set_m=Data(),
51	useDiagonalHessianApproximation=True, tol=1e-8):
52	"""
53	initialization
54
55	:param domain: domain
56	:type domain: `Domain`
57	:param numLevelSets: number of level sets
58	:type numLevelSets: ``int``
59	: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
61	: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.
63	: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.
65	: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.
67	: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
69	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.
71	:type useDiagonalHessianApproximation: ``bool``
72	:param tol: toleramce when solving the PDE for the inverse of the Hessian Operator
73	:type tol: positive ``float``
74
75
76	"""
77	if numLevelSets>1:
78	raise ValueError("Currently only numLevelSets<=1 is supportered.")
79	if s0 == None and s1==None:
80	raise ValueError("Values for s0 or s1 must be given.")
81
82	self.__domain=domain
83	DIM=self.__domain.getDim()
84	self.__L2=np.asarray(boundingBoxEdgeLengths(domain))**2
85	self.__L4=np.sum(self.__L2)**2
86	self.__numLevelSets=numLevelSets
87
88	if self.__numLevelSets > 1:
89	self.__useDiagonalHessianApproximation=useDiagonalHessianApproximation
90	else:
91	self.__useDiagonalHessianApproximation=True
92	self._update_Hessian=True
93
94	self.__pde=LinearPDE(self.__domain, numEquations=self.__numLevelSets)
95	self.__pde.getSolverOptions().setTolerance(tol)
96	self.__pde.setSymmetryOn()
97	self.__pde_is_set=False
98	try:
99	self.__pde.setValue(q=location_of_set_m)
100	except IllegalCoefficientValue:
101	raise ValueError("Unable to set location of fixed level set function.")
102
103	self.__total_num_weights=2numLevelSets+((numLevelSets-1)numLevelSets)/2
104	self.__weight_index=[] # this is a mapping from the relevant mu-coefficients to the set of all mu-coefficients
105	# we count s0, then s1, then sc (k<l).
106	# THE S0 weighting factor
107	n=0
108	VV=vol(domain)
109	if not s0 is None:
110	s0 = interpolate(s0,self.__pde.getFunctionSpaceForCoefficient('D'))
111	s=s0.getShape()
112	if numLevelSets == 1 :
113	if s == () :
114	V=integrate(s0)
115	if V > 0:
116	self.__weight_index.append(n)
117	s0*=VV/V
118	else:
119	s0=None
120	else:
121	raise ValueError("Unexpected shape %s for weight s0."%s)
122	else:
123	if s == (numLevelSets,):
124	for k in xrange(numLevelSets):
125	V=integrate(s0[k])
126	if V > 0:
127	self.__weight_index.append(n+k)
128	s0[k]*=VV/V
129	else:
130	raise ValueError("Unexpected shape %s for weight s0."%s)
131	self.__s0=s0
132	n+=numLevelSets
133
134	# The S1 weighting factor
135	if not s1 is None:
136	s1 = interpolate(s1,self.__pde.getFunctionSpaceForCoefficient('A'))
137	s=s1.getShape()
138
139	if numLevelSets == 1 :
140	if s == (DIM,) :
141	V=integrate(inner(s1, 1/self.__L2))
142	if V > 0:
143	self.__weight_index.append(n)
144	s1*=VV/V
145	print "REG SCALE = ",s1
146	else:
147	raise ValueError("Unexpected shape %s for weight s1."%s)
148	else:
149	if s == (numLevelSets,DIM):
150	for k in xrange(numLevelSets):
151	for i in xrange(DIM):
152	ww=s1[k,:]
153	V=integrate(inner(ww,1/self.__L2))
154	if V > 0:
155	self.__weight_index.append(n+i)
156	s1[k,:]=ww*(VV/V)
157	else:
158	raise ValueError("Unexpected shape %s for weight s1."%s)
159
160	self.__s1=s1
161	n+=numLevelSets
162
163	# The SC weighting factor
164	if not sc is None:
165	if numLevelSets == 1 :
166	sc=None
167	else:
168	sc = interpolate(sc,self.__pde.getFunctionSpaceForCoefficient('A'))
169	s=sc.getShape()
170	if s == (numLevelSets,numLevelSets):
171	for k in xrange(numLevelSets):
172	sc[k,k]=0.
173	for l in xrange(k):
174	ww=sc[l,k]
175	V=integrate(ww)
176	if V > 0:
177	self.__weight_index.append(n+k*numLevelSets+l)
178	sc[l,k]=wwVV/Vself.__L4
179	sc[k,l]=0
180	else:
181	raise ValueError("Unexpected shape %s for weight s0."%s)
182
183	self.__sc=sc
184	self.setWeights()
185
186	def getDomain(self):
187	"""
188	return the domain of the regularization term
189	:rtype: ``Domain``
190	"""
191	return self.__domain
192
193	def getNumLevelSets(self):
194	"""
195	return the number of level set functions
196	:rtype: ``int``
197	"""
198	return self.__numLevelSets
199
200	def getPDE(self):
201	"""
202	returns the linear PDE to be solved for the Hessian Operator inverse
203	:rtype: ``LinearPDE``
204	"""
205	return self.__pde
206
207	def getDualProduct(self, m, r):
208	"""
209	returns the dual product of a gradient represented by X=r[1] and Y=r[0] with a level set function m:
210
211	Y_im_i + X_ijm_{i,j}
212
213	:type m: ``esys.escript.Data``
214	:type r: ``ArithmeticTuple``
215	:rtype: float
216	"""
217	A=0
218	if not r[0].isEmpty(): A+=integrate(inner(r[0], m))
219	if not r[1].isEmpty(): A+=integrate(inner(r[1], grad(m)))
220	return A
221
222	def getValue(self, m, grad_m):
223	"""
224	return the value of the costfunction J
225
226
227
228	:rtype: ``float``
229	"""
230	mu=self.getWeights( uncompress=True)
231	DIM=self.getDomain().getDim()
232	numLS=self.getNumLevelSets()
233
234	A=0
235	n=0
236
237	if self.__s0 is not None:
238	A+=inner(integrate(m*2self.__s0), mu[:numLS])
239	n+=numLS
240
241	if self.__s1 is not None:
242	if numLS == 1:
243	A+=integrate(inner(grad_m*2, self.__s1))mu[n]
244	else:
245	for k in xrange(numLS):
246	A+=mu[n+k]integrate(inner(grad_m[k,:]*2,self.__s1[k,:]))
247	n+=numLS
248
249	if self.__sc is not None:
250	for k in xrange(numLS):
251	gk=grad_m[k,:]
252	len_gk=length(gk)
253	for l in xrange(k):
254	gl=grad_m[l,:]
255	A+= mu[n+knumLS+l] integrate( self.__sc[l,k] * ( len_gk * length(gl) )2 - inner(gk, gl)2 )
256	return A/2
257
258	def getGradient(self, m, grad_m):
259	"""
260	returns the gradient of the costfunction J with respect to m. The function returns Y_k=dPsi/dm_k and
261	X_kj=dPsi/dm_kj
262	"""
263
264	mu=self.getWeights( uncompress=True)
265	DIM=self.getDomain().getDim()
266	numLS=self.getNumLevelSets()
267
268	n=0
269
270	if self.__s0 is not None:
271	Y = m * self.__s0 * mu[:numLS]
272	else:
273	Y = Data()
274	n+=numLS
275
276	if self.__s1 is not None:
277	if numLS == 1:
278	X=grad_m* self.__s1*mu[n]
279	else:
280	X=self.getPDE().createCoefficient("X")
281	for k in xrange(numLS):
282	X[k,:]=mu[n+k]grad_m[k,:]self.__s1[k,:]
283	else:
284	X = Data()
285	n+=numLS
286	if self.__sc is not None:
287	raise NotImplementedError
288	return ArithmeticTuple(Y, X)
289
290	def getArguments(self, m):
291	"""
292	"""
293	return ( grad(m),)
294
295	def getInverseHessianApproximation(self, m, r, grad_m):
296	"""
297	"""
298	if self._new_mu or self._update_Hessian:
299
300	self._new_mu=False
301	self._update_Hessian=False
302
303	mu=self.getWeights( uncompress=True)
304	DIM=self.getDomain().getDim()
305	numLS=self.getNumLevelSets()
306	n=0
307	if self.__s0 is not None:
308	if numLS == 1:
309	D=self.__s0 * mu[n]
310	else:
311	D=self.getPDE().createCoefficient("D")
312	for k in xrange(numLS): D[k,k]=self.__s0[k] * mu[n+k]
313	self.getPDE().setValue(D=D)
314	n+=numLS
315
316	A=self.getPDE().createCoefficient("A")
317	if self.__s1 is not None:
318	if numLS == 1:
319	for i in xrange(DIM): A[i,i]=self.__s1[i] * mu[n]
320	else:
321	for k in xrange(numLS):
322	for i in xrange(DIM): A[k,i,k,i]=self.__s1[k,i] * mu[n+k]
323	n+=numLS
324	if self.__sc is not None:
325	raise NotImplementedError
326	print A
327	self.getPDE().setValue(A=A)
328	self.getPDE().resetRightHandSideCoefficients()
329	self.getPDE().setValue(X=r[1])
330	print "X only: ",self.getPDE().getSolution()
331	self.getPDE().resetRightHandSideCoefficients()
332	self.getPDE().setValue(Y=r[0])
333	print "Y only: ",self.getPDE().getSolution()
334
335	self.getPDE().resetRightHandSideCoefficients()
336	self.getPDE().setValue(X=r[1], Y=r[0])
337	return self.getPDE().getSolution()
338
339	def updateHessian(self):
340	"""
341	notify the class to recalculate the Hessian operator
342	"""
343	if not self.__useDiagonalHessianApproximation:
344	self._update_Hessian=True
345
346	# ================ we should factor these out =============================================================
347	def getNumRelevantTerms(self):
348	"""
349	returns the number of terms in the costfunction of regularization with non-zero scaling factors
350	:rtype: ``int``
351	"""
352	return len(self.__weight_index)
353
354	def getNumTerms(self):
355	"""
356	returns the number of terms in the costfunction of regularization with non-zero scaling factors
357	:rtype: ``int``
358	"""
359	return len(self.__weight_index)
360
361	def setWeights(self, mu=None):
362	"""
363	sets the values for the weighting terms in the costsfunction. Note that only values
364	correspond to entries with non-negative scaling factors are used.
365
366	:param mu: weights
367	:type mu: ``list`` of ``float`` or ``np,array``
368	"""
369	if mu == None:
370	mu = np.ones(self.getNumRelevantTerms())
371	mu=np.asarray(mu)
372
373
374	if len(mu) == self.getNumRelevantTerms():
375	if not mu.shape == (self.getNumRelevantTerms(),):
376	raise ValueError("%s values are expected."%self.getNumRelevantTerms())
377	self.__mu=mu
378	else:
379	if not mu.shape == (self.__total_num_weights,):
380	raise ValueError("%s values are expected."%self.__total_num_weights)
381	self.__mu = np.zeros(self.getNumRelevantTerms())
382	for i in xrange(len(self.__weight_index)): self.__mu[i]=mu[self.__weight_index[i]]
383	self._new_mu=True
384
385	def setWeightsForS0(self, mu=None):
386	"""
387	set the weights for s0-terms
388	"""
389	numLS=self.getNumLevelSets()
390	my_mu=self.getWeights(uncompress=True)
391	if mu is None:
392	my_mu[:numLS]=1
393	else:
394	my_mu[:numLS]=mu
395	self.setWeights(my_mu)
396
397	def setWeightsForS1(self, mu=None):
398	"""
399	set the weights for s1-terms
400	"""
401	numLS=self.getNumLevelSets()
402	my_mu=self.getWeights(uncompress=True)
403	if mu is None:
404	my_mu[numLS:2*numLS]=1
405	else:
406	my_mu[numLS:2*numLS]=mu
407	self.setWeights(my_mu)
408
409	def setWeightsForSc(self, mu):
410	"""
411	set the weights for s1-terms
412	"""
413	numLS=self.getNumLevelSets()
414	my_mu=self.getWeights(uncompress=True)
415	if mu is None:
416	my_mu[2*numLS:]=1
417	else:
418	my_mu[2*numLS:]=mu
419	self.setWeights(my_mu)
420
421
422	def getWeights(self, uncompress=False):
423	"""
424	Returns the weights for the terms in the costsfunction.
425	The first ``numLevelSets`` values used for the
426	regularization terms and the remaining values for the cross correlation terms.
427
428	:type mu: ``list`` of ``float``
429	"""
430	if uncompress:
431	mu = np.zeros(self.__total_num_weights)
432	for i in xrange(len(self.__weight_index)): mu[self.__weight_index[i]] = self.__mu[i]
433	return mu
434	else:
435	return self.__mu
436
437	def getWeightIndex(self):
438	"""
439	returns an iondex to the contributions of terms with non-zero scaling factor.
440	"""
441	return self.__weight_index
442