agd.Metrics.rander

  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 .base import Base
  6from .riemann import Riemann
  7from . import misc
  8from .. import LinearParallel as lp
  9from .. import AutomaticDifferentiation as ad
 10from .. import FiniteDifferences as fd
 11
 12class Rander(Base):
 13	r"""
 14	A Rander norm takes the form 
 15	$$
 16	F(x) = \sqrt{< x,m.x>} + < w,x>,
 17	$$
 18	where m is a given symmetric positive definite tensor, 
 19	and w is a given vector subject to the consition $< w,m^{-1} w> < 1$.
 20
 21	Member fields and __init__ arguments : 
 22	- m : an array of shape (vdim,vdim,n1,..,nk) where vdim is the ambient space dimension.
 23	The array must be symmetric, a.k.a m[i,j] = m[j,i] for all $0\leq i< j< vdim$.
 24	- w : an array of shape (vdim,n1,...,nk)
 25	"""
 26	def __init__(self,m,w):
 27		m,w = (ad.asarray(e) for e in (m,w))
 28		self.m,self.w = fd.common_field((m,w),(2,1))
 29
 30	def norm(self,v):
 31		v,m,w = fd.common_field((ad.asarray(v),self.m,self.w),(1,2,1))
 32		return np.sqrt(lp.dot_VAV(v,m,v))+lp.dot_VV(w,v)
 33
 34
 35	def dual(self):
 36		# The dual of a Randers norm is also a Randers norm
 37		s = lp.inverse(self.m-lp.outer_self(self.w))
 38		omega = lp.dot_AV(s,self.w)
 39		return Rander((1+lp.dot_VV(self.w,omega))*s, -omega)
 40
 41	@property
 42	def vdim(self): return len(self.m)
 43
 44	@property
 45	def shape(self): return self.m.shape[2:]	
 46	
 47	def is_definite(self):
 48		return np.logical_and(Riemann(self.m).is_definite(),
 49			lp.dot_VV(lp.solve_AV(self.m,self.w),self.w) <1.)
 50
 51	def anisotropy_bound(self):
 52		r = np.sqrt(lp.dot_VV(lp.solve_AV(self.m,self.w),self.w))
 53		return ((1.+r)/(1.-r))*Riemann(self.m).anisotropy_bound()
 54	def cost_bound(self):
 55		return Riemann(self.m).cost_bound() + ad.Optimization.norm(self.w,ord=2,axis=0)
 56
 57	def inv_transform(self,a):
 58		w,a = fd.common_field((self.w,a),depths=(1,2))
 59		return Rander(Riemann(self.m).inv_transform(a).m,lp.dot_VA(w,a))
 60
 61	def with_costs(self,costs):
 62		costs,m,w = fd.common_field((costs,self.m,self.w),depths=(1,2,1))
 63		return Rander(m*lp.outer_self(costs),w*costs)
 64
 65	def flatten(self,inverse_m=False):
 66		m,w = self.m,self.w
 67		if inverse_m: m = lp.inverse(m)
 68		return np.concatenate((misc.flatten_symmetric_matrix(m),w),axis=0)
 69	
 70	@classmethod
 71	def expand(cls,arr):
 72		m = misc.expand_symmetric_matrix(arr,extra_length=True)
 73		d=len(m)
 74		d_sym = (d*(d+1))//2
 75		assert(len(arr)==d_sym+d)
 76		w=arr[d_sym:]
 77		return cls(m,w)
 78
 79	def model_HFM(self):
 80		return "Rander"+str(self.vdim)
 81
 82	@classmethod
 83	def from_cast(cls,metric): 
 84		if isinstance(metric,cls):	return metric
 85		riemann = Riemann.from_cast(metric)
 86		return cls(riemann.m,(0.,)*riemann.vdim)
 87
 88	def __iter__(self):
 89		yield self.m
 90		yield self.w
 91
 92	@classmethod
 93	def from_Zermelo(cls,metric,drift):
 94		r"""
 95		Zermelo's navigation problem consists in computing a minimal path for 
 96		whose velocity is unit w.r.t. a Riemannian metric, and which is subject 
 97		to a drift. The accessible velocities take the form 
 98			$x+drift$ where $< x,m.x> \leq 1$
 99		This function reformulates it as a shortest path problem 
100		in a Rander manifold.
101
102		Inputs : 
103		- metric : Symmetric positive definite matrix (Riemannian metric)
104		- drift : Vector field, obeying $< drift,metric.drift> < 1$ (Drift)
105
106		Outputs : 
107		- the Rander metric.
108		"""
109		return cls(lp.inverse(metric),drift).dual()
110
111	def to_Zermelo(self):
112		"""
113		Input : Parameters of a Rander metric.
114
115		Output : Parameters of the corresponding Zermelo problem, of motion on a 
116		Riemannian manifold with a drift.
117		"""
118		self_dual = self.dual()
119		return lp.inverse(self_dual.m), self_dual.w
120
121	def to_Varadhan(eps=1):
122		r"""
123		The Rander eikonal equation can be reformulated in an (approximate)
124		linear form, using a logarithmic transformation
125			$u + 2 eps < omega,grad u> - eps**2 Tr(D hess u)$.
126		Then -eps log(u) solves the Rander eikonal equation, 
127		up to a small additional diffusive term.
128
129		Inputs : 
130		- m and w, parameters of the Rander metric
131		- eps (optionnal), relaxation parameter
132
133		Outputs : 
134		- $D$ and $2 * omega$, parameters of the linear PDE. 
135		 ($D * eps^2$ and $2 * omega * eps$ if eps is specified)
136		"""
137		return self.Varadhan_from_Zermelo(self.to_Zermelo(),eps)
138
139	@staticmethod
140	def Varadhan_from_Zermelo(metric,drift,eps=1):
141		"""
142		Zermelo's navigation problem can be turned into a Rander shortest path problem,
143		which itself can be (approximately) expressed in linear form using the logarithmic
144		transformation. This function composes the above two steps.
145		"""
146		return eps**2*(lp.inverse(metric)-lp.outer_self(drift)), 2.*eps*drift
147
148	def make_proj_dual(self,**kwargs):
149		"""kwargs : passed to Riemann.make_proj_dual"""
150		# Based on the Zermelo transformation, by which the Randers unit ball can be expressed 
151		# as a shifted ellipse
152		proj_riemann = Riemann(self.m).make_proj_dual(**kwargs)
153		def proj(x):
154			x,w = fd.common_field((x,self.w),depths=(1,1))
155			return proj_riemann(x-w)+w
156		return proj