agd.Selling

  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
  4r"""
  5This file implements Selling's algorithm in dimension $d \in \\{2,3\\}$, which decomposes a 
  6symmetric positive definite matrix $D$, of dimension $d\leq 3$, in the form 
  7$$
  8	D = \sum_{0\leq i < I} a_i e_i e_i^\top,
  9$$
 10where $a_i \geq 0$ and $e_i\in Z^d$ is a vector with integer coordinates, 
 11and where $I = d(d+2)/2$.
 12
 13Selling's decomposition is a central tool in the design of adaptive discretization schemes
 14for anisotropic partial differential equations, on Cartesian grids.
 15"""
 16
 17import numpy as np
 18from itertools import cycle
 19from .LinearParallel import dot_VAV, perp, cross
 20from . import AutomaticDifferentiation as ad
 21samesize_int_t = ad.cupy_generic.samesize_int_t
 22
 23iterMax2 = 100 
 24iterMax3 = 100 
 25
 26# -------- Dimension based dispatch -----
 27
 28def ObtuseSuperbase(m,sb=None):
 29	"""
 30	Input : 
 31	- m : symmetric positive definite matrix, defined as an
 32	 array of shape $(d,d, n_1,...,n_k)$.
 33	- sb (optional) : initial guess for the obtuse superbase.
 34
 35	Ouput : an m-obtuse superbase 
 36	"""
 37	m = ad.remove_ad(m) # The superbase is a discrete object, AD is not applicable.
 38	dim = len(m)
 39	if sb is None: sb = CanonicalSuperbase(m)
 40	if   dim==1: return _ObtuseSuperbase1(m,sb)	
 41	elif dim==2: return _ObtuseSuperbase2(m,sb)
 42	elif dim==3: return _ObtuseSuperbase3(m,sb)
 43	else: raise ValueError("Selling's decomposition only applies in dimension <=3.") 	
 44
 45def Decomposition(m,sb=None):
 46	r"""
 47	Use Selling's algorithm to decompose a symmetric positive definite matrix.
 48	Note : results unspecified if matrices are not symmetric (includes optional AD part) 
 49	or not positive definite.
 50
 51	Input : 
 52	- m : symmetric positive definite matrix, defined as an
 53	 array of shape $(d,d, n_1,...,n_k)$ where $d\leq 3$.
 54	- sb (optional) : superbase to use for the decomposition,
 55	 array of shape $(d,d+1, n_1,...,n_k)$.
 56	Output : the coefficients and offsets of the decomposition.
 57	"""
 58	dim = len(m)
 59
 60	# Strange bugs were observed when running the code below on the GPU with early cupy versions,
 61	# related with the array indexing pattern. Using a custom cuda kernel is faster and safer.
 62	if ad.cupy_generic.from_cupy(m) and dim>1:
 63		from . import Eikonal
 64		return Eikonal.VoronoiDecomposition(m,mode='gpu')
 65
 66	if sb is None: sb = ObtuseSuperbase(m,sb)
 67	if   dim==1: return _Decomposition1(m,sb)	
 68	elif dim==2: return _Decomposition2(m,sb)
 69	elif dim==3: return _Decomposition3(m,sb)
 70	else: raise ValueError("Selling's decomposition only applies in dimension <=3.") 
 71
 72
 73def GatherByOffset(T,Coefs,Offsets):
 74	"""
 75	Get the coefficient of each offset.
 76	This function is essentially used to make nice plots of how the superbase coefficients
 77	and offsets vary as the decomposed tensor varies. Opposite offsets are regarded as identical.
 78	"""
 79	T,Coefs,Offsets = map(ad.cupy_generic.cupy_get,(T,Coefs,Offsets))
 80	TimeCoef = {};
 81	for (i,t) in enumerate(T):
 82		coefs = Coefs[:,i]
 83		offsets = Offsets[:,:,i]
 84		for (j,c) in enumerate(coefs):
 85			offset = tuple(offsets[:,j].astype(int))
 86			offset_m = tuple(-offsets[:,j].astype(int))
 87			if offset<offset_m:
 88				offset=offset_m
 89			if offset in TimeCoef:
 90				TimeCoef[offset][0].append(t)
 91				TimeCoef[offset][1].append(c)
 92			else:
 93				TimeCoef[offset] = ([t],[c])
 94	return TimeCoef
 95
 96
 97def CanonicalSuperbase(m):
 98	"""
 99	Returns a superbase with the same dimensions and array type as m.
100
101	Output : 
102	 - m : array of shape $(d,d, n_1,...,n_k)$
103	"""
104	d=len(m); assert m.shape[1]==d
105	shape=m.shape[2:]
106
107	sb=np.zeros_like(m,shape=(d,d+1,*shape))
108	sb[:,0]=-1
109	for i in range(d):
110		sb[i,i+1]=1
111	return sb
112
113def SuperbasesForConditioning(cond,dim=2):
114	"""
115	Returns a family of superbases. 
116	For any positive definite matrix $M$ with condition number below the given bound,
117	one of these superbases will be $M$-obtuse.
118	(Condition number is the ratio of the largest to the smallest eigenvalue.)
119
120	Input : 
121	 - cond (scalar) : the bound on the condition number.
122	"""
123	if   dim==1: return _SuperbasesForConditioning1(cond)
124	elif dim==2: return _SuperbasesForConditioning2(cond)
125	elif dim==3: return _SuperbasesForConditioning3(cond)
126	else: raise ValueError("Selling's decomposition only applies in dimension <=3.") 
127
128
129# ------- One dimensional variant (trivial) ------
130
131def _ObtuseSuperbase1(m,sb=None):
132	return CanonicalSuperbase(m)
133
134def _Decomposition1(m,sb):
135	shape = m.shape[2:]
136	offsets = np.ones_like(ad.remove_ad(m),shape=(1,1,*shape),dtype=samesize_int_t(m.dtype))
137	coefs = m.reshape((1,*shape))
138	return coefs,offsets
139
140def _SuperbasesForConditioning1(cond):
141	sb = CanonicalSuperbase(np.eye(1))
142	return sb.reshape(sb.shape+(1,))
143
144# ------- Two dimensional variant ------
145
146# We do everyone in parallel, without selection or early abort
147def _ObtuseSuperbase2(m,sb):
148	"""
149		Use Selling's algorithm to compute an obtuse superbase.
150
151		input : symmetric positive definite matrix m, dim=2
152		input/output : superbase b (must be valid at startup)
153		
154		module variable : iterMax2, max number of iterations
155
156		output : wether the algorithm succeeded
157	"""
158	iterReduced = 0
159	for iter,(i,j,k) in zip(range(iterMax2), cycle([(0,1,2),(1,2,0),(2,0,1)]) ):
160		# Test for a positive angle, and modify superbase if necessary
161		acute = dot_VAV(sb[:,i],m,sb[:,j]) > 0
162		if np.any(acute):
163			sb[:,k,acute] = sb[:,i,acute]-sb[:,j,acute]
164			sb[:,i,acute] = -sb[:,i,acute]
165			iterReduced=0
166		elif iterReduced<3: iterReduced+=1
167		else: return sb
168
169	raise ValueError(f"Selling's algorithm did not terminate in iterMax2={iterMax2} iterations")
170	
171# Produce the matrix decomposition
172def _Decomposition2(m,sb):
173	"""
174		Use Selling's algorithm to decompose a tensor
175
176		input : symmetric positive definite tensor 
177		output : coefficients, offsets
178	"""
179	shape = m.shape[2:]
180	coef  = np.zeros_like(m,shape=(3,*shape))
181	for (i,j,k) in [(0,1,2),(1,2,0),(2,0,1)]:
182		coef[i] = -dot_VAV(sb[:,j],m, sb[:,k]) 
183	
184	return coef,perp(sb).astype(samesize_int_t(coef.dtype))
185
186def _SuperbasesForConditioning2(cond):
187	"""
188	Implementation is based on exploring the Stern-Brocot tree, 
189	with a stopping criterion based on the angle between consecutive vectors.
190	"""
191
192	mu = np.sqrt(cond)
193	theta = np.pi/2. - np.arccos( 2/(mu+1./mu))
194
195	u=np.array( (1,0) )
196	l = [np.array( (-1,0) ),np.array( (0,1) )]
197	m = []
198	superbases =[]
199
200	def angle(u,v): return np.arctan2(u[0]*v[1]-u[1]*v[0], u[0]*v[0]+u[1]*v[1])
201
202	while l:
203		v=l[-1]
204		if angle(u,v)<=theta:
205			m.append(u)
206			u=v
207			l.pop()
208		else:
209			l.append(u+v)
210			superbases.append((u,v,-u-v))
211
212
213	return np.array(superbases).transpose((2,1,0))
214
215# ------- Three dimensional variant -------
216
217# We do everyone in parallel, without selection or early abort
218def _ObtuseSuperbase3(m,sb):
219	"""
220		Use Selling's algorithm to compute an obtuse superbase.
221
222		input : symmetric positive definite matrix m, dim=3
223		input/output : superbase b (must be valid at startup)
224		
225		module variable : iterMax3, max number of iterations
226
227		output : wether the algorithm succeeded
228	"""
229#	This energy is reduced at each time step
230#	def energy(): return sum([dot_VAV(sb[:,i],m,sb[:,i]) for i in range(4)])
231#	Introducing a tolerance, for the positivity test, was not needed in our experiments
232#	tol = 100*np.finfo(m.dtype).eps*trace(m) 
233#	if(np.any(tol<=0)): raise ValueError("Matrices must be positive definite (negative trace found)") 
234	iterReduced = 0
235	sigma = cycle([(0,1,2,3),(0,2,1,3),(0,3,1,2),(1,2,0,3),(1,3,0,2),(2,3,0,1)])
236	for iter,(i,j,k,l) in zip(range(iterMax3), sigma):		
237		# Test for a positive angle, and modify superbase if necessary
238		acute = dot_VAV(sb[:,i],m,sb[:,j]) > 0
239		if np.any(acute):
240			sb[:,k,acute] += sb[:,i,acute]
241			sb[:,l,acute] += sb[:,i,acute]
242			sb[:,i,acute] = -sb[:,i,acute]
243			iterReduced=0
244		elif iterReduced<6: iterReduced+=1
245		else: return sb
246	raise ValueError(f"Selling's algorithm did not terminate in iterMax3={iterMax3} iterations")
247	
248def _Decomposition3(m,sb):
249	"""
250		Use Selling's algorithm to decompose a tensor
251
252		input : symmetric positive definite tensor, d=3
253		output : coefficients, offsets
254	"""
255	shape = m.shape[2:]
256	
257	coef   = np.zeros_like(m,shape=(6,*shape))
258	offset = np.zeros_like(ad.remove_ad(m),shape=(3,6,*shape))
259	for iter,(i,j,k,l) in enumerate(
260		[(0,1,2,3),(0,2,1,3),(0,3,1,2),(1,2,0,3),(1,3,0,2),(2,3,0,1)]):
261		coef[iter] = -dot_VAV(sb[:,i], m, sb[:,j]) 
262		offset[:,iter] = cross(sb[:,k], sb[:,l])
263		
264	return coef,offset.astype(samesize_int_t(coef.dtype))
265
266def _SuperbasesForConditioning3(cond):
267	raise ValueError("Sorry, _SuperbasesForConditioning3 is not implemented yet")