agd.Metrics.Seismic.tti

  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
  5import copy
  6
  7from ... import LinearParallel as lp
  8from ... import AutomaticDifferentiation as ad
  9from ... import FiniteDifferences as fd
 10from .. import misc
 11from ..riemann import Riemann
 12from .implicit_base import ImplicitBase
 13from .thomsen_data import ThomsenElasticMaterial, ThomsenGeometricMaterial
 14
 15
 16class TTI(ImplicitBase):
 17	"""
 18	A family of reduced models, known as Tilted Transversally Anisotropic,
 19	and arising in seismic tomography.
 20
 21	The *dual* unit ball is defined by an equation of the form
 22	$$
 23	l(X^2+Y^2,Z^2) + (1/2)*q(X^2+Y^2,Z^2) = 1,
 24	$$
 25	where $l$ is linear and $q$ is quadratic, where $X,Y,Z$ are the coefficients 
 26	of the input vector, usually altered by a linear transformation.
 27	In two dimensions, ignore the $Y^2$ term.
 28
 29	The primal norm is obtained implicitly, by solving an optimization problem.
 30
 31	Member fields and __init__ arguments : 
 32	- linear : an array of shape (2,n1,...,nk) encoding the linear part l
 33	- quadratic : an array of shape (2,2,n1,...,nk) encoding the quadratic part q
 34	- vdim (optional) : the ambient space dimension
 35	- *args,**kwargs (optional) : passed to implicit_base
 36	"""
 37
 38	def __init__(self,linear,quadratic,vdim=None,*args,**kwargs): #rotation_angles=None,
 39		super(TTI,self).__init__(*args,**kwargs)
 40		self.linear = ad.asarray(linear)
 41		self.quadratic = ad.asarray(quadratic)
 42		assert len(self.linear) == 2
 43		assert self.quadratic.shape[:2] == (2,2)
 44		self._to_common_field()
 45		
 46		if vdim is None:
 47			if self.inverse_transformation is not None: vdim=len(self.inverse_transformation)
 48			elif self.linear.ndim>1: vdim = self.linear.ndim-1
 49			else: raise ValueError("Unspecified dimension")
 50		self._vdim=vdim
 51		#self.rotation_angles=rotation_angles
 52
 53	@property
 54	def vdim(self): return self._vdim
 55	
 56	@property
 57	def shape(self): return self.linear.shape[1:]
 58
 59	def _dual_level_root(self,v):
 60		"""Level set function defining the dual unit ball in the square root domain"""
 61		l,q = self._dual_params(v.shape[1:])
 62		return lp.dot_VV(l,v) + 0.5*lp.dot_VAV(v,q,v) - 1.
 63
 64	def _dual_level(self,v,params=None,relax=0.):
 65		"""Level set function defining the dual unit ball."""
 66		l,q = self._dual_params(v.shape[1:]) if params is None else params
 67		v2 = v**2
 68		if self.vdim==3: v2 = ad.array([v2[:2].sum(axis=0),v2[2]])
 69		return lp.dot_VV(l,v2) + 0.5*np.exp(-relax)*lp.dot_VAV(v2,q,v2) - 1.
 70	
 71	def cost_bound(self):
 72		return self.Isotropic_approx()[1]
 73		# Ignoring the quadratic term for now.
 74		return self.Riemann_approx().cost_bound()
 75
 76	def _dual_params(self,shape=None):
 77		return fd.common_field((self.linear,self.quadratic),depths=(1,2),shape=shape)
 78
 79	def __iter__(self):
 80		yield self.linear
 81		yield self.quadratic
 82		yield self._vdim
 83#		yield self.rotation_angles
 84		for x in super(TTI,self).__iter__(): yield x
 85
 86	def _to_common_field(self,shape=None):
 87		self.linear,self.quadratic,self.inverse_transformation = fd.common_field(
 88			(self.linear,self.quadratic,self.inverse_transformation),
 89			depths=(1,2,2),shape=shape)
 90
 91	@classmethod
 92	def from_cast(cls,metric):
 93		if isinstance(metric,cls): return metric
 94		else: raise ValueError("No casting operations supported towards the TTI model")
 95		# Even cast from Riemann is non-canonical
 96
 97	def model_HFM(self):
 98		return f"TTI{self.vdim}"
 99
100	def extract_xz(self):
101		"""
102		Extract a two dimensional Hooke tensor from a three dimensional one, 
103		corresponding to a slice through the X and Z axes.
104		Axes transformation information (rotation) is discarded.
105 		"""
106		if len(self.shape)==3: raise ValueError("Three dimensional field")
107		if self.inverse_transformation is not None:
108			raise ValueError("Cannot extract XZ slice from tilted norms")
109		other = copy.copy(self)
110		other._vdim = 2
111		return other
112
113	def flatten(self,transposed_transformation=False,cp_get=False):
114		if self.inverse_transformation is None:
115			xp = ad.cupy_generic.get_array_module(self.linear)
116			trans = fd.as_field(xp.eye(self.vdim,dtype=self.linear.dtype),self.shape,depth=2) 
117		else: trans = self.inverse_transformation
118		if transposed_transformation: trans = lp.transpose(lp.inverse(trans))
119
120		if cp_get: # return a numpy array (memory optimization for GPU intensive)
121			return np.concatenate((self.linear.get(),
122				misc.flatten_symmetric_matrix(self.quadratic.get()),
123				trans.get().reshape((self.vdim**2,)+self.shape)),axis=0)
124		else: 
125			return np.concatenate((self.linear,misc.flatten_symmetric_matrix(self.quadratic),
126				trans.reshape((self.vdim**2,)+self.shape)),axis=0)
127
128	@classmethod
129	def expand(cls,arr):
130		vdim = np.sqrt(len(arr)-(2+3))
131		assert(vdim==int(vdim))
132		vdim = int(vdim)
133		shape = arr.shape[1:]
134
135		linear = arr[0:2]
136		quadratic = misc.expand_symmetric_matrix(arr[2:5])
137		inv_trans = arr[5:].reshape((vdim,vdim)+shape)
138		return cls(linear,quadratic,vdim=vdim,inverse_transformation=inv_trans)
139
140	@classmethod
141	def from_hexagonal(cls,c11,_,c13,c33,c44,vdim=3):
142		linear = [c11+c44,c33+c44]
143		mixed = c13**2-c11*c33+2.*c13*c44
144		quadratic = [[-2.*c11*c44,mixed],[mixed,-2.*c33*c44]]
145		return cls(linear,quadratic,vdim=vdim)
146
147	@classmethod
148	def from_ThomsenElastic(cls,tem,vdim=3):
149		"""Produces a norm from the given Thomsem elasticity parameters."""
150		if not isinstance(tem,ThomsenElasticMaterial): tem = ThomsenElasticMaterial(*tem)
151		hex,ρ = tem.to_hexagonal()
152		return cls.from_hexagonal(*hex,vdim),ρ
153
154	@classmethod
155	def from_ThomsenGeometric(cls,tgm,vdim=3,normalize_Vp=False):
156		"""Produces a norm from the given Thomsen geometric paramters."""
157		if not isinstance(tgm,ThomsenGeometricMaterial): tgm = ThomsenGeometricMaterial(*tgm)
158		if normalize_Vp:
159			Vp,Vs,ϵ,δ = tgm
160			tgm = ThomsenGeometricMaterial(1.,Vs/Vp,ϵ,δ)
161		c11,c13,c33,c44 = tgm.to_c()
162		return cls.from_hexagonal(c11,None,c13,c33,c44,vdim=vdim)
163
164	def α_bounds(self):
165		"""
166		The TTI norm can be written as an enveloppe of ellipses, with eigenvalues 
167		(1-α,α) / μ(α) or (1-α,1-α,α)/μ(α), where α is within the bounds given by this function.
168		Note : another way to obtain the envelope is to use the Riemann_envelope method.
169
170		Returns : αa,αb,mix_is_min
171		"""
172		l,q = self.linear,self._q()
173		a,b = _axes_intersections(l,q)
174		z = np.zeros_like(a)
175		ga = _grad_ratio(l,q,(a,z))
176		gb = _grad_ratio(l,q,(z,b))
177		αa = ga[1]/ga.sum(axis=0)
178		αb = gb[1]/gb.sum(axis=0)
179		return np.minimum(αa,αb),np.maximum(αa,αb),αa<αb
180
181	def μ(self,α):
182		""" See the method α_bounds. """
183		a = ad.array((1-α,α))
184		l = self.linear
185		Q = self.quadratic
186		δ = lp.det(Q)
187		δRef = np.sum(Q**2,axis=(0,1)) # If δ/δRef is small, then degenerate case (two lines)
188		R = ad.array([[Q[1,1],-Q[0,1]],[-Q[1,0],Q[0,0]]]) # Adjugate of Q
189		Rl = lp.dot_AV(R,l)
190		s = 2*δ+lp.dot_VV(l,Rl)
191		ε = np.sign(s)
192		aRa = lp.dot_VAV(a,R,a)
193
194		# For generic parameters, the result is num/den. However, den = num = 0 is possible.
195		# This arises in particular if Q = λ l l^T, in which case the conic consists of two lines.
196		# We must also check wether Q is close to zero, in which case the conic degenerates
197		# to a line. In both these cases, we approximate the conic piece witha line, and
198		# we get μ by using one intersection along an axis : satisfy (1-α,α).(a0,0)/μ(α) = 1.
199		# Another possible criterion, more costly but possibly better, would be to directly 
200		# check the difference between self.α_bounds() - if they are close, then the conic piece
201		# can be well approximated with a line. 
202		tol = 5*np.finfo(Q.dtype).eps;
203		degen = np.sum(Rl**2,axis=0)<=tol*np.sum(l**2,axis=0)**3 # conic piece degenerates to a line
204#		print(np.sum(Rl**2,axis=0)/np.sum(l**2,axis=0)**3,tol)
205#		tol = 100*np.finfo(Q.dtype).eps; degen = np.abs(δ)<=δRef*tol # Degenerate quadratic form => two lines
206#		print(degen,l,R,Rl)
207		sol_degen = (1-α) * _solve2_above(-2,l[0],Q[0,0],0.) # Use intersection along one axis
208		
209		num = lp.det([a,l])**2+2*aRa
210		den = ε*np.sqrt(aRa*s)+lp.dot_VV(a,Rl)
211#		print(δ/δRef,np.finfo(Q.dtype).eps)
212#		print(self.α_bounds(),α)
213#		print(degen,sol_degen,num,den)
214		return np.where(degen,sol_degen, num/den)
215
216	def Isotropic_approx(self):
217		"""
218		Isotropic approximation of the TTI norm.
219		Returns : (cmin,cmax). These costs correspond to the interior and exterior approximations.
220		Assumption : the linear transformation must be a rotation.
221		"""
222		αmin,αmax,mix_is_min = self.α_bounds()
223		l,q = self.linear,self._q()
224		a,b = _axes_intersections(l,q)
225		# Costs corresponding to vertical and horizontal propagation for the Dual norm
226		c0,c1=np.sqrt(1/a),np.sqrt(1/b)
227		pos = (αmin<0.5) & (0.5<αmax) # If pos==True, one extremal speed is not on the axes.
228		αmed = np.where(pos,0.5,(αmin+αmax)/2)
229		# 0.5 corresponds to the isotropy. (αmin+αmax)/2 is a dummy value in [αmin,αmax]
230		ch = np.where(pos,np.sqrt(0.5/self.μ(αmed)),(c0+c1)/2)
231		c = np.sort(ad.asarray([c0,ch,c1]),axis=0)
232		return 1/c[2],1/c[0] # Take inverses, to account for norm duality
233
234	def Riemann_approx(self,avg=False,**kwargs):
235		"""
236		Riemannian approximations of the TTI norm, homothetic to each other, and expectedly good.
237		Returns : 
238		 - Mmin,Mmax. Some interior and exterior approximating Riemannian metrics.
239		 - Mavg, if avg=True. A good approximating Riemannian metric, neither interior nor exterior.
240		 - kwargs : passed to Riemann.dual (typical : avoid_np_linalg_inv=True, for speed)
241		"""
242		import time; top = time.time()
243		l,q = self.linear,self._q()
244		a,b = _axes_intersections(l,q)
245		ai,bi = 1/a,1/b
246		α = bi/(ai+bi)
247		c = 1/(self.μ(α)*(ai+bi))
248		cmin,cmax,cavg = np.minimum(1,c),np.maximum(1,c),((1+np.sqrt(c))/2)**2
249		diag = (ai,bi) if self.vdim==2 else (ai,ai,bi)
250		riem = Riemann.from_diagonal(diag)
251		if self.inverse_transformation is not None: 
252			riem = riem.inv_transform(self.inverse_transformation)
253		m = riem.dual(**kwargs).m # Take inverse, to account for norm duality
254		return m/cavg if avg else (m/cmax,m/cmin)
255
256	def Riemann_envelope(self,nmix,gpu=False):
257		"""
258		Approximation of a TTI norm using an envelope of Riemannian norms.
259		- nmix : number of ellipses used for the approximation.
260		- gpu : same implementation as on the GPU (changes almost nothing)
261		returns
262		- riems : a list of nmix Riemannian norms.
263		- mix_is_min : wether to take the minimum or the maximum of the Riemannian norms.
264		"""
265		# Set the interpolation times 
266		if isinstance(nmix,int):
267			if nmix==1: ts=[0.5]
268			else: ts = np.linspace(0,1,nmix)
269		else: ts = nmix # Hack to get specific interpolation times, must be sorted
270
271		if gpu:
272			l,q = self.linear,self._q()
273			ab = _axes_intersections(l,q)
274			diags = _diags(l,q,ts,ab)
275			mix_is_min = lp.det([diags[0],diags[-1]])>0 if len(diags) else None
276		else:
277			αmin,αmax,mix_is_min = self.α_bounds()
278			diags = [np.array([1-α,α])/self.μ(α) for α in np.linspace(αmin,αmax,nmix)]
279
280		if self.vdim==3: diags = [(a,a,b) for a,b in diags]
281		riems = [Riemann.from_diagonal(1./ad.array(diag)) for diag in diags]
282		if self.inverse_transformation is not None:
283			riems = [riem.inv_transform(self.inverse_transformation) for riem in riems]
284		
285		return mix_is_min, riems
286
287	def _q(self):
288		"""
289		Quadratic part, in format compatible with the C routines adapted below
290		"""
291		return self.quadratic[((0,0,1),(0,1,1))]
292
293# ----- The following code is adapted from agd/Eikonal/HFM_CUDA/cuda/TTI_.h -----
294# It computes the approximation of the TTI norm with an envelope of riemannian norms
295# The chosen collection of ellipses (very slightly) differs from a uniform sampling  
296# within αmin,αmax given by self.α_bounds() 
297
298def _solve2(a,b,c):
299	"""
300	Returns the two roots of a quadratic equation, a + 2 b t + c t^2.
301	The discriminant must be non-negative, but aside from that the equation may be degenerate.
302	"""
303	sdelta = np.sqrt(b*b-a*c);
304	u = -b + sdelta 
305	v = -b - sdelta
306	b0 = np.abs(c)>np.abs(a)
307	xp = ad.cupy_generic.get_array_module(b0)
308	b1 = xp.asarray(a!=0) # Needed with cupy
309#	with np.errstate(divide='ignore'): # Does not silent much
310	return xp.asarray( (np.where(b0,u/c,np.where(b1,a/u,0.)), 
311		np.where(b0,v/c,np.where(b1,a/v,np.inf))) )
312
313def _solve2_above(a, b, c, above):
314	"""
315	Returns the smallest root of the considered quadratic equation above the given threshold.
316	Such a root is assumed to exist.
317	"""
318	r = np.sort(_solve2(a,b,c),axis=0)
319	return np.where(r[0]>=above, r[0], r[1])
320
321def _axes_intersections(l,q):
322	"""
323	Finds the intersections (a,0) and (0,b) of the curve f(x)=0
324	with the axes (the closest intersections to the origin).
325	"""
326	return _solve2_above(-2,l[0],q[0],0.), _solve2_above(-2,l[1],q[2],0.)
327
328def _grad_ratio(l,q,x):
329	"""
330	Returns g(x) := df(x)/<x,df(x)> where f(x):= C + 2 <l,x> + <qx,x> 
331	Note that the curve tangent to the level line of f at x is 
332	<y-x,df(x)> ≥ 0, equivalently <y,g(x)> ≥ 1
333	"""
334	hgrad = ad.array([q[0]*x[0]+q[1]*x[1]+l[0], q[1]*x[0]+q[2]*x[1]+l[1]]) #lp.dot_AV(q,x)+l # df(x)/2
335	return hgrad/lp.dot_VV(x,hgrad)
336
337def _diags(l, q, ts, axes_intersections):
338	"""
339	Samples the curve defined by {x≥0|f(x)=0}, 
340	(the connected component closest to the origin)
341	where f(x):= -2 + 2 <l,x> + <qx,x>,
342	and returns diag(i) := grad f(x)/<x,grad f(x)>.
343	"""
344	a,b=axes_intersections
345	zero = np.zeros_like(a)
346
347	# Change of basis in f, with e0 = {1/2.,1/2.}, e1 = {1/2.,-1/2.}
348	L = ad.array([l[0]+l[1], l[0]-l[1]])/2.
349	Q = ad.array([q[0]+2*q[1]+q[2], q[0]-q[2], q[0]-2*q[1]+q[2]])/4.
350
351	diags=[]
352	for t in ts:
353		if   t==0.: x=(a,zero)
354		elif t==1.: x=(zero,b)
355		else :
356			v = (1.-t)*a - t*b
357			# Solving f(u e0+ v e_1) = 0 w.r.t u
358			u = _solve2_above(-2.+2.*L[1]*v+Q[2]*v*v, L[0]+Q[1]*v, Q[0], np.abs(v))
359			# Inverse change of basis
360			x = ad.array([u+v, u-v])/2.
361		diags.append(_grad_ratio(l,q,x))
362
363	return diags
364
365
366# ---- Some instances of TTI metrics ----
367
368# See Hooke.py file for reference
369TTI.mica = TTI.from_hexagonal(178.,42.4,14.5,54.9,12.2), 2.79
370# Stishovite is tetragonal, not hexagonal, but the P-Wave velocity, in the XZ plane,
371# is equivalent to an hexagonal model.
372TTI.stishovite2 = TTI.from_hexagonal(453,np.nan,203,776,252).extract_xz(), 4.29