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
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)
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.
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 \}$
Dimensions of the underlying domain.
Expected to be the empty tuple, or a tuple of length vdim
.
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
.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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)
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.
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