agd.Metrics.Seismic.implicit_base

View Source

  1# Copyright 2020 Jean-Marie Mirebeau, University Paris-Sud, CNRS, University Paris-Saclay
  2# Distributed WITHOUT ANY WARRANTY. Licensed under the Apache License, Version 2.0, see http://www.apache.org/licenses/LICENSE-2.0
  3
  4import numpy as np
  5from ... import AutomaticDifferentiation as ad
  6from ... import LinearParallel as lp
  7from ... import FiniteDifferences as fd
  8from ..base import Base
  9import copy
 10
 11class ImplicitBase(Base):
 12	"""
 13	Base class for a metric defined implicitly, 
 14	in terms of a level set function for the unit ball
 15	of the dual metric, and of a linear transformation.
 16
 17	Inputs:
 18	- niter_sqp (int, optional): number of iterations for Sequential Quadratic Programming
 19	- relax_sqp (tuple,optional): relaxation parameter for the first iterations of SQP
 20	- qconv_sqp (real, optional): such that hessian+quasi_sqp*grad^T grad > 0. Used when 
 21	  the constraint is a quasi-convex function, but exp(qconv_sqp*f) is strongly convex
 22	"""
 23
 24	def __init__(self,inverse_transformation=None,
 25		niter_sqp=6,relax_sqp=tuple(),qconv_sqp=0.):
 26		self.inverse_transformation = inverse_transformation
 27		self.niter_sqp = niter_sqp
 28		self.relax_sqp = relax_sqp
 29		self.qconv_sqp = qconv_sqp
 30
 31	def norm(self,v):
 32		v=ad.asarray(v)
 33		a=self.inverse_transformation
 34		if a is not None:  
 35			v,a = fd.common_field((v,a),(1,2))
 36			v = lp.dot_AV(a,v)
 37		# AD case is handled by the envelope theorem
 38		return lp.dot_VV(v,self._gradient(ad.remove_ad(v))) 
 39
 40	def gradient(self,v):
 41		v=ad.asarray(v)
 42		a=self.inverse_transformation
 43		if a is None:
 44			return self._gradient(v)
 45		else:
 46			return lp.dot_AV(lp.transpose(a),self._gradient(lp.dot_AV(a,v)))
 47
 48	def with_cost(self,cost): 
 49		if self.inverse_transformation is None: 
 50			return super(ImplicitBase,self).with_cost(cost)
 51		inv_trans,cost = fd.common_field((self.inverse_transformation,cost),depths=(2,0))
 52		inv_trans = inv_trans*cost
 53		other = copy.copy(self)
 54		other.inverse_transformation = inv_trans
 55		other._to_common_field()
 56		return other
 57
 58	def inv_transform(self,a):
 59		inv_trans,a = fd.common_field((self.inverse_transformation,a),depths=(2,2))
 60		inv_trans = a if inv_trans is None else lp.dot_AA(inv_trans,a)
 61		other = copy.copy(self)
 62		other.inverse_transformation = inv_trans
 63		other._to_common_field()
 64		return other
 65
 66	def is_topographic(self,a=None):
 67		if a is None: a = self.inverse_transformation
 68		if a in None: return True
 69		d = self.vdim
 70		return np.all([a[i,j]==(i==j) for i in range(d) for j in range(d-1)])
 71
 72	def flatten_transform(self,topographic=None):
 73		a = self.inverse_transformation
 74		if a is None: return None
 75
 76		if topographic is None: topographic = self.is_topographic(a)
 77		d=self.vdim
 78		if topographic:
 79			return ad.asarray([a[i,-1] in range(d-1)] + [a[d-1,d-1]-1])
 80		else:
 81			return a.reshape((d*d,)+a.shape[2:])
 82
 83	def _gradient(self,v):
 84		"""
 85		Gradient, ignoring self.a
 86		Note : modifies v where null
 87		"""
 88		v=ad.asarray(v) 
 89		zeros = np.all(v==0.,axis=0)
 90		v = np.where(zeros,np.nan,v) # Silent errors due to non differentiability at origin.
 91		grad = sequential_quadratic(v,self._dual_level,params=self._dual_params(v.shape[1:]),
 92			niter=self.niter_sqp,relax=self.relax_sqp,qconv=self.qconv_sqp)
 93		grad[:,zeros]=0.
 94		return grad
 95
 96
 97	def _dual_level(self,v,params=None,relax=0):
 98		"""
 99		A level set function for the dual unit ball, ignoring self.inverse_transformation.
100		Parameters
101		- v : co-vector
102		- params : Some parameters of the instance can be passed in argument, for AD purposes.
103		- relax : for a relaxation of the level set. 0->exact, np.inf->easy (quadratic).
104		"""
105		raise ValueError('_dual_level is not implemented for this class')
106		
107	def _dual_params(self,*args,**kwargs):
108		"""
109		The parameters to be passed to _dual_level.
110		"""
111		return None
112
113	def __iter__(self):
114		yield self.inverse_transformation
115		yield self.niter_sqp
116		yield self.relax_sqp
117		yield self.qconv_sqp
118
119	def _to_common_field(self,*args,**kwargs):
120		"""Makes compatible the dimensions of the various fields"""
121		raise ValueError("_to_common_field is not implemented for this class")
122
123
124def sequential_quadratic(v,f,niter,x=None,params=tuple(),relax=tuple(),qconv=0.):
125	"""
126	Maximizes <x,v> subject to the constraint f(x,*params)<=0, 
127	using sequential quadratic programming.
128	x : initial guess.
129	relax : relaxation parameters to be used in a preliminary path following phase.
130	params : to be passed to evaluated function. Special treatment if ad types.
131	"""
132	if x is None: x=np.zeros_like(v)
133	x_ad = ad.Dense2.identity(constant=x,shape_free=(len(x),))
134
135	# Fixed point iterations 
136	def step(val,V,D,v):
137		M = lp.inverse(D+qconv*lp.outer_self(V))
138		k = np.sqrt((lp.dot_VAV(V,M,V)-2.*val)/lp.dot_VAV(v,M,v))
139		return lp.dot_AV(M,k*v-V)
140
141	# Initial iterations ignoring AD information in params
142	params_noad = tuple(ad.remove_ad(val) for val in params) 
143	for r in relax + (0.,)*niter:
144		f_ad = f(x_ad,params_noad,relax=r)
145		x_ad = x_ad+step(f_ad.value,f_ad.gradient(),f_ad.hessian(),v)
146
147	x=x_ad.value
148
149	# Terminal iteration to introduce ad information from params
150	adtype = ad.is_ad(params,iterables=(tuple,))
151	if adtype:
152		shape_bound = x.shape[1:]
153		params_dis = tuple(ad.disassociate(value,shape_bound=shape_bound) 
154			if ad.cupy_generic.isndarray(value) else value for value in params)
155		x_ad = ad.Dense2.identity(constant=ad.disassociate(x,shape_bound=shape_bound))
156
157		f_ad = f(x_ad,params_dis,0.)
158		x = x + step(ad.associate(f_ad.value), ad.associate(f_ad.gradient()), ad.associate(f_ad.hessian()), v)
159
160	return x

class ImplicitBase(agd.Metrics.base.Base): View Source

 12class ImplicitBase(Base):
 13	"""
 14	Base class for a metric defined implicitly, 
 15	in terms of a level set function for the unit ball
 16	of the dual metric, and of a linear transformation.
 17
 18	Inputs:
 19	- niter_sqp (int, optional): number of iterations for Sequential Quadratic Programming
 20	- relax_sqp (tuple,optional): relaxation parameter for the first iterations of SQP
 21	- qconv_sqp (real, optional): such that hessian+quasi_sqp*grad^T grad > 0. Used when 
 22	  the constraint is a quasi-convex function, but exp(qconv_sqp*f) is strongly convex
 23	"""
 24
 25	def __init__(self,inverse_transformation=None,
 26		niter_sqp=6,relax_sqp=tuple(),qconv_sqp=0.):
 27		self.inverse_transformation = inverse_transformation
 28		self.niter_sqp = niter_sqp
 29		self.relax_sqp = relax_sqp
 30		self.qconv_sqp = qconv_sqp
 31
 32	def norm(self,v):
 33		v=ad.asarray(v)
 34		a=self.inverse_transformation
 35		if a is not None:  
 36			v,a = fd.common_field((v,a),(1,2))
 37			v = lp.dot_AV(a,v)
 38		# AD case is handled by the envelope theorem
 39		return lp.dot_VV(v,self._gradient(ad.remove_ad(v))) 
 40
 41	def gradient(self,v):
 42		v=ad.asarray(v)
 43		a=self.inverse_transformation
 44		if a is None:
 45			return self._gradient(v)
 46		else:
 47			return lp.dot_AV(lp.transpose(a),self._gradient(lp.dot_AV(a,v)))
 48
 49	def with_cost(self,cost): 
 50		if self.inverse_transformation is None: 
 51			return super(ImplicitBase,self).with_cost(cost)
 52		inv_trans,cost = fd.common_field((self.inverse_transformation,cost),depths=(2,0))
 53		inv_trans = inv_trans*cost
 54		other = copy.copy(self)
 55		other.inverse_transformation = inv_trans
 56		other._to_common_field()
 57		return other
 58
 59	def inv_transform(self,a):
 60		inv_trans,a = fd.common_field((self.inverse_transformation,a),depths=(2,2))
 61		inv_trans = a if inv_trans is None else lp.dot_AA(inv_trans,a)
 62		other = copy.copy(self)
 63		other.inverse_transformation = inv_trans
 64		other._to_common_field()
 65		return other
 66
 67	def is_topographic(self,a=None):
 68		if a is None: a = self.inverse_transformation
 69		if a in None: return True
 70		d = self.vdim
 71		return np.all([a[i,j]==(i==j) for i in range(d) for j in range(d-1)])
 72
 73	def flatten_transform(self,topographic=None):
 74		a = self.inverse_transformation
 75		if a is None: return None
 76
 77		if topographic is None: topographic = self.is_topographic(a)
 78		d=self.vdim
 79		if topographic:
 80			return ad.asarray([a[i,-1] in range(d-1)] + [a[d-1,d-1]-1])
 81		else:
 82			return a.reshape((d*d,)+a.shape[2:])
 83
 84	def _gradient(self,v):
 85		"""
 86		Gradient, ignoring self.a
 87		Note : modifies v where null
 88		"""
 89		v=ad.asarray(v) 
 90		zeros = np.all(v==0.,axis=0)
 91		v = np.where(zeros,np.nan,v) # Silent errors due to non differentiability at origin.
 92		grad = sequential_quadratic(v,self._dual_level,params=self._dual_params(v.shape[1:]),
 93			niter=self.niter_sqp,relax=self.relax_sqp,qconv=self.qconv_sqp)
 94		grad[:,zeros]=0.
 95		return grad
 96
 97
 98	def _dual_level(self,v,params=None,relax=0):
 99		"""
100		A level set function for the dual unit ball, ignoring self.inverse_transformation.
101		Parameters
102		- v : co-vector
103		- params : Some parameters of the instance can be passed in argument, for AD purposes.
104		- relax : for a relaxation of the level set. 0->exact, np.inf->easy (quadratic).
105		"""
106		raise ValueError('_dual_level is not implemented for this class')
107		
108	def _dual_params(self,*args,**kwargs):
109		"""
110		The parameters to be passed to _dual_level.
111		"""
112		return None
113
114	def __iter__(self):
115		yield self.inverse_transformation
116		yield self.niter_sqp
117		yield self.relax_sqp
118		yield self.qconv_sqp
119
120	def _to_common_field(self,*args,**kwargs):
121		"""Makes compatible the dimensions of the various fields"""
122		raise ValueError("_to_common_field is not implemented for this class")

Base class for a metric defined implicitly, in terms of a level set function for the unit ball of the dual metric, and of a linear transformation.

Inputs:

niter_sqp (int, optional): number of iterations for Sequential Quadratic Programming
relax_sqp (tuple,optional): relaxation parameter for the first iterations of SQP
qconv_sqp (real, optional): such that hessian+quasi_sqpgrad^T grad > 0. Used when the constraint is a quasi-convex function, but exp(qconv_sqpf) is strongly convex

ImplicitBase( inverse_transformation=None, niter_sqp=6, relax_sqp=(), qconv_sqp=0.0) View Source

25	def __init__(self,inverse_transformation=None,
26		niter_sqp=6,relax_sqp=tuple(),qconv_sqp=0.):
27		self.inverse_transformation = inverse_transformation
28		self.niter_sqp = niter_sqp
29		self.relax_sqp = relax_sqp
30		self.qconv_sqp = qconv_sqp

inverse_transformation

niter_sqp

relax_sqp

qconv_sqp

def norm(self, v): View Source

32	def norm(self,v):
33		v=ad.asarray(v)
34		a=self.inverse_transformation
35		if a is not None:  
36			v,a = fd.common_field((v,a),(1,2))
37			v = lp.dot_AV(a,v)
38		# AD case is handled by the envelope theorem
39		return lp.dot_VV(v,self._gradient(ad.remove_ad(v)))

Norm or quasi-norm defined by the class, often denoted $N$ in mathematical formulas. Unless incorrect data is provided, this member function obeys, for all vectors $u,v\in R^d$ and all $\alpha \geq 0$

$N(u+v) \leq N(u)+N(v)$
$N(\alpha u) = \alpha N(u)$
$N(u)\geq 0$ with equality iff $u=0$.

Broadcasting will occur depending on the shape of $v$ and of the class data.

def gradient(self, v): View Source

41	def gradient(self,v):
42		v=ad.asarray(v)
43		a=self.inverse_transformation
44		if a is None:
45			return self._gradient(v)
46		else:
47			return lp.dot_AV(lp.transpose(a),self._gradient(lp.dot_AV(a,v)))

Gradient of the norm defined by the class.

def with_cost(self, cost): View Source

49	def with_cost(self,cost): 
50		if self.inverse_transformation is None: 
51			return super(ImplicitBase,self).with_cost(cost)
52		inv_trans,cost = fd.common_field((self.inverse_transformation,cost),depths=(2,0))
53		inv_trans = inv_trans*cost
54		other = copy.copy(self)
55		other.inverse_transformation = inv_trans
56		other._to_common_field()
57		return other

Produces a norm $N'$ obeying $N'(x) = N(cost*x)$.

def inv_transform(self, a): View Source

59	def inv_transform(self,a):
60		inv_trans,a = fd.common_field((self.inverse_transformation,a),depths=(2,2))
61		inv_trans = a if inv_trans is None else lp.dot_AA(inv_trans,a)
62		other = copy.copy(self)
63		other.inverse_transformation = inv_trans
64		other._to_common_field()
65		return other

Affine transformation of the norm. The new unit ball is the inverse image of the previous one.

def is_topographic(self, a=None): View Source

67	def is_topographic(self,a=None):
68		if a is None: a = self.inverse_transformation
69		if a in None: return True
70		d = self.vdim
71		return np.all([a[i,j]==(i==j) for i in range(d) for j in range(d-1)])

def flatten_transform(self, topographic=None): View Source

73	def flatten_transform(self,topographic=None):
74		a = self.inverse_transformation
75		if a is None: return None
76
77		if topographic is None: topographic = self.is_topographic(a)
78		d=self.vdim
79		if topographic:
80			return ad.asarray([a[i,-1] in range(d-1)] + [a[d-1,d-1]-1])
81		else:
82			return a.reshape((d*d,)+a.shape[2:])

Inherited Members

agd.Metrics.base.Base: dual; vdim; shape; disassociate; is_definite; anisotropy; anisotropy_bound; cost_bound; cos_asym; cos; angle; transform; rotate; rotate_by; with_costs; with_speeds; with_speed; flatten; expand; to_HFM; model_HFM; from_HFM; from_generator; from_cast; array_float_caster; norm2; gradient2; set_interpolation; at; make_proj_dual; make_proj; make_prox

def sequential_quadratic(v, f, niter, x=None, params=(), relax=(), qconv=0.0): View Source

125def sequential_quadratic(v,f,niter,x=None,params=tuple(),relax=tuple(),qconv=0.):
126	"""
127	Maximizes <x,v> subject to the constraint f(x,*params)<=0, 
128	using sequential quadratic programming.
129	x : initial guess.
130	relax : relaxation parameters to be used in a preliminary path following phase.
131	params : to be passed to evaluated function. Special treatment if ad types.
132	"""
133	if x is None: x=np.zeros_like(v)
134	x_ad = ad.Dense2.identity(constant=x,shape_free=(len(x),))
135
136	# Fixed point iterations 
137	def step(val,V,D,v):
138		M = lp.inverse(D+qconv*lp.outer_self(V))
139		k = np.sqrt((lp.dot_VAV(V,M,V)-2.*val)/lp.dot_VAV(v,M,v))
140		return lp.dot_AV(M,k*v-V)
141
142	# Initial iterations ignoring AD information in params
143	params_noad = tuple(ad.remove_ad(val) for val in params) 
144	for r in relax + (0.,)*niter:
145		f_ad = f(x_ad,params_noad,relax=r)
146		x_ad = x_ad+step(f_ad.value,f_ad.gradient(),f_ad.hessian(),v)
147
148	x=x_ad.value
149
150	# Terminal iteration to introduce ad information from params
151	adtype = ad.is_ad(params,iterables=(tuple,))
152	if adtype:
153		shape_bound = x.shape[1:]
154		params_dis = tuple(ad.disassociate(value,shape_bound=shape_bound) 
155			if ad.cupy_generic.isndarray(value) else value for value in params)
156		x_ad = ad.Dense2.identity(constant=ad.disassociate(x,shape_bound=shape_bound))
157
158		f_ad = f(x_ad,params_dis,0.)
159		x = x + step(ad.associate(f_ad.value), ad.associate(f_ad.gradient()), ad.associate(f_ad.hessian()), v)
160
161	return x

Maximizes subject to the constraint f(x,*params)<=0, using sequential quadratic programming. x : initial guess. relax : relaxation parameters to be used in a preliminary path following phase. params : to be passed to evaluated function. Special treatment if ad types.