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

A Rander norm takes the form $$ F(x) = \sqrt{< x,m.x>} + < w,x>, $$ where m is a given symmetric positive definite tensor, and w is a given vector subject to the consition $< w,m^{-1} w> < 1$.

Member fields and __init__ arguments :

  • m : an array of shape (vdim,vdim,n1,..,nk) where vdim is the ambient space dimension. The array must be symmetric, a.k.a m[i,j] = m[j,i] for all $0\leq i< j< vdim$.
  • w : an array of shape (vdim,n1,...,nk)
Rander(m, w)
27	def __init__(self,m,w):
28		m,w = (ad.asarray(e) for e in (m,w))
29		self.m,self.w = fd.common_field((m,w),(2,1))
def norm(self, v):
31	def norm(self,v):
32		v,m,w = fd.common_field((ad.asarray(v),self.m,self.w),(1,2,1))
33		return np.sqrt(lp.dot_VAV(v,m,v))+lp.dot_VV(w,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 dual(self):
36	def dual(self):
37		# The dual of a Randers norm is also a Randers norm
38		s = lp.inverse(self.m-lp.outer_self(self.w))
39		omega = lp.dot_AV(s,self.w)
40		return Rander((1+lp.dot_VV(self.w,omega))*s, -omega)

Dual norm, mathematically defined by $N^*(x) = max\{ < x, y> ; N(y)\leq 1 \}$

vdim
42	@property
43	def vdim(self): return len(self.m)

The ambient vector space dimension, often denoted $d$ in mathematical formulas.

shape
45	@property
46	def shape(self): return self.m.shape[2:]	

Dimensions of the underlying domain. Expected to be the empty tuple, or a tuple of length vdim.

def is_definite(self):
48	def is_definite(self):
49		return np.logical_and(Riemann(self.m).is_definite(),
50			lp.dot_VV(lp.solve_AV(self.m,self.w),self.w) <1.)

Attempts to check wether the data defines a mathematically valid norm.

def anisotropy_bound(self):
52	def anisotropy_bound(self):
53		r = np.sqrt(lp.dot_VV(lp.solve_AV(self.m,self.w),self.w))
54		return ((1.+r)/(1.-r))*Riemann(self.m).anisotropy_bound()

An upper bound on the anisotropy of the norm.

def cost_bound(self):
55	def cost_bound(self):
56		return Riemann(self.m).cost_bound() + ad.Optimization.norm(self.w,ord=2,axis=0)

Upper bound on $N(u)$, for any unit vector $u$, where $N$ is the norm defined by the class.

def inv_transform(self, a):
58	def inv_transform(self,a):
59		w,a = fd.common_field((self.w,a),depths=(1,2))
60		return Rander(Riemann(self.m).inv_transform(a).m,lp.dot_VA(w,a))

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

def with_costs(self, costs):
62	def with_costs(self,costs):
63		costs,m,w = fd.common_field((costs,self.m,self.w),depths=(1,2,1))
64		return Rander(m*lp.outer_self(costs),w*costs)

Produces a norm $N'$ defined by $$ N'(x) = N(costs * x) $$ where the multiplication is elementwise.

def flatten(self, inverse_m=False):
66	def flatten(self,inverse_m=False):
67		m,w = self.m,self.w
68		if inverse_m: m = lp.inverse(m)
69		return np.concatenate((misc.flatten_symmetric_matrix(m),w),axis=0)

Flattens and concatenate the member fields into a single array.

@classmethod
def expand(cls, arr):
71	@classmethod
72	def expand(cls,arr):
73		m = misc.expand_symmetric_matrix(arr,extra_length=True)
74		d=len(m)
75		d_sym = (d*(d+1))//2
76		assert(len(arr)==d_sym+d)
77		w=arr[d_sym:]
78		return cls(m,w)

Inverse of the flatten member function. Turns a suitable array into a metric.

def model_HFM(self):
80	def model_HFM(self):
81		return "Rander"+str(self.vdim)

The name of the 'model' for parameter, as input to the HFM library.

@classmethod
def from_cast(cls, metric):
83	@classmethod
84	def from_cast(cls,metric): 
85		if isinstance(metric,cls):	return metric
86		riemann = Riemann.from_cast(metric)
87		return cls(riemann.m,(0.,)*riemann.vdim)

Produces a metric by casting another metric of a compatible type.

@classmethod
def from_Zermelo(cls, metric, drift):
 93	@classmethod
 94	def from_Zermelo(cls,metric,drift):
 95		r"""
 96		Zermelo's navigation problem consists in computing a minimal path for 
 97		whose velocity is unit w.r.t. a Riemannian metric, and which is subject 
 98		to a drift. The accessible velocities take the form 
 99			$x+drift$ where $< x,m.x> \leq 1$
100		This function reformulates it as a shortest path problem 
101		in a Rander manifold.
102
103		Inputs : 
104		- metric : Symmetric positive definite matrix (Riemannian metric)
105		- drift : Vector field, obeying $< drift,metric.drift> < 1$ (Drift)
106
107		Outputs : 
108		- the Rander metric.
109		"""
110		return cls(lp.inverse(metric),drift).dual()

Zermelo's navigation problem consists in computing a minimal path for whose velocity is unit w.r.t. a Riemannian metric, and which is subject to a drift. The accessible velocities take the form $x+drift$ where $< x,m.x> \leq 1$ This function reformulates it as a shortest path problem in a Rander manifold.

Inputs :

  • metric : Symmetric positive definite matrix (Riemannian metric)
  • drift : Vector field, obeying $< drift,metric.drift> < 1$ (Drift)

Outputs :

  • the Rander metric.
def to_Zermelo(self):
112	def to_Zermelo(self):
113		"""
114		Input : Parameters of a Rander metric.
115
116		Output : Parameters of the corresponding Zermelo problem, of motion on a 
117		Riemannian manifold with a drift.
118		"""
119		self_dual = self.dual()
120		return lp.inverse(self_dual.m), self_dual.w

Input : Parameters of a Rander metric.

Output : Parameters of the corresponding Zermelo problem, of motion on a Riemannian manifold with a drift.

def to_Varadhan(eps=1):
122	def to_Varadhan(eps=1):
123		r"""
124		The Rander eikonal equation can be reformulated in an (approximate)
125		linear form, using a logarithmic transformation
126			$u + 2 eps < omega,grad u> - eps**2 Tr(D hess u)$.
127		Then -eps log(u) solves the Rander eikonal equation, 
128		up to a small additional diffusive term.
129
130		Inputs : 
131		- m and w, parameters of the Rander metric
132		- eps (optionnal), relaxation parameter
133
134		Outputs : 
135		- $D$ and $2 * omega$, parameters of the linear PDE. 
136		 ($D * eps^2$ and $2 * omega * eps$ if eps is specified)
137		"""
138		return self.Varadhan_from_Zermelo(self.to_Zermelo(),eps)

The Rander eikonal equation can be reformulated in an (approximate) linear form, using a logarithmic transformation $u + 2 eps < omega,grad u> - eps**2 Tr(D hess u)$. Then -eps log(u) solves the Rander eikonal equation, up to a small additional diffusive term.

Inputs :

  • m and w, parameters of the Rander metric
  • eps (optionnal), relaxation parameter

Outputs :

  • $D$ and $2 * omega$, parameters of the linear PDE. ($D * eps^2$ and $2 * omega * eps$ if eps is specified)
@staticmethod
def Varadhan_from_Zermelo(metric, drift, eps=1):
140	@staticmethod
141	def Varadhan_from_Zermelo(metric,drift,eps=1):
142		"""
143		Zermelo's navigation problem can be turned into a Rander shortest path problem,
144		which itself can be (approximately) expressed in linear form using the logarithmic
145		transformation. This function composes the above two steps.
146		"""
147		return eps**2*(lp.inverse(metric)-lp.outer_self(drift)), 2.*eps*drift

Zermelo's navigation problem can be turned into a Rander shortest path problem, which itself can be (approximately) expressed in linear form using the logarithmic transformation. This function composes the above two steps.

def make_proj_dual(self, **kwargs):
149	def make_proj_dual(self,**kwargs):
150		"""kwargs : passed to Riemann.make_proj_dual"""
151		# Based on the Zermelo transformation, by which the Randers unit ball can be expressed 
152		# as a shifted ellipse
153		proj_riemann = Riemann(self.m).make_proj_dual(**kwargs)
154		def proj(x):
155			x,w = fd.common_field((x,self.w),depths=(1,1))
156			return proj_riemann(x-w)+w
157		return proj

kwargs : passed to Riemann.make_proj_dual