agd.Metrics.base
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 .. import AutomaticDifferentiation as ad 6from ..AutomaticDifferentiation import cupy_support as cps 7from .. import LinearParallel as lp 8from .. import FiniteDifferences as fd 9from .. import Interpolation 10 11class Base: 12 """ 13 Base class for a metric, in other words a family of norms. 14 """ 15 16 def __repr__(self): 17 return type(self).__name__ + repr(tuple(self)) 18 19 def __str__(self): 20 return type(self).__name__ + str(tuple(self)) 21 22 def norm(self,v): 23 r""" 24 Norm or quasi-norm defined by the class, often denoted $N$ in mathematical 25 formulas. Unless incorrect data is provided, this member function obeys, 26 for all vectors $u,v\in R^d$ and all $\alpha \geq 0$ 27 - $N(u+v) \leq N(u)+N(v)$ 28 - $N(\alpha u) = \alpha N(u)$ 29 - $N(u)\geq 0$ with equality iff $u=0$. 30 31 Broadcasting will occur depending on the shape of $v$ and of the class data. 32 """ 33 raise NotImplementedError("""Error : norm must be specialized in subclass""") 34 35 def gradient(self,v): 36 """ 37 Gradient of the `norm` defined by the class. 38 """ 39 if ad.is_ad(v) or ad.is_ad(self,iterables=(type(self),)): 40 v_dis = ad.disassociate(v,shape_bound=v.shape[1:]) 41 grad_dis = self.disassociate().gradient(v_dis) 42 return ad.associate(grad_dis) 43 44 v_ad = ad.Dense.identity(constant=v,shape_free=(len(v),)) 45 return np.moveaxis(self.norm(v_ad).coef,-1,0) 46 47 def dual(self): 48 r""" 49 Dual `norm`, mathematically defined by 50 $N^*(x) = max\\{ < x, y> ; N(y)\leq 1 \\}$ 51 """ 52 raise NotImplementedError("dual is not implemented for this norm") 53 54 @property 55 def vdim(self): 56 """ 57 The ambient vector space dimension, often denoted $d$ in mathematical formulas. 58 """ 59 raise NotImplementedError("vdim is not implemented for this norm") 60 61 @property 62 def shape(self): 63 """ 64 Dimensions of the underlying domain. 65 Expected to be the empty tuple, or a tuple of length `vdim`. 66 """ 67 raise NotImplementedError("Shape not implemented for this norm") 68 69 def disassociate(self): 70 """ 71 Hide the automatic differentiation (AD) information of the member fields. 72 See AutomaticDifferentiation.disassociate 73 """ 74 def dis(x): 75 if ad.isndarray(x) and x.shape[-self.vdim:]==self.shape: 76 return ad.disassociate(x,shape_bound=self.shape) 77 return x 78 return self.from_generator(dis(x) for x in self) 79 80# ---- Well posedness related methods ----- 81 def is_definite(self): 82 """ 83 Attempts to check wether the data defines a mathematically valid `norm`. 84 """ 85 raise NotImplementedError("is_definite is not implemented for this norm") 86 87 def anisotropy(self): 88 r""" 89 Anisotropy ratio of the `norm` denoted $N$. 90 Defined as 91 $$ 92 \max_{|u| = |v| = 1} \frac {N(u)}{N(v)}. 93 $$ 94 """ 95 raise NotImplementedError("anisotropy is not implemented for this norm") 96 97 def anisotropy_bound(self): 98 """ 99 An upper bound on the `anisotropy` of the norm. 100 """ 101 return self.anisotropy() 102 def cost_bound(self): 103 """ 104 Upper bound on $N(u)$, for any unit vector $u$, where $N$ is the `norm` 105 defined by the class. 106 """ 107 raise NotImplementedError("cost_bound is not implemented for this norm") 108 109# ---- Causality and acuteness related methods ---- 110 111 def cos_asym(self,u,v): 112 r""" 113 Generalized asymmetric cosine defined by the norm $N$, of the angle between two vectors $u,v$. 114 Defined as 115 $\cos^a_N(u,v) := < \nabla N(u), v> / N(v)$. 116 """ 117 u,v=(ad.asarray(e) for e in (u,v)) 118 return lp.dot_VV(self.gradient(u),v)/self.norm(v) 119 120 def cos(self,u,v): 121 r""" 122 Generalized cosine defined by the norm $N$. Expression 123 $\cos_N(u,v) := \min( \cos^a_N (u,v), \cos^a_N(v,u))$, 124 where $\cos^a_N=$`cos_asym`. 125 """ 126 u,v=(ad.asarray(e) for e in (u,v)) 127 gu,gv=self.gradient(u),self.gradient(v) 128 guu,guv = lp.dot_VV(gu,u),lp.dot_VV(gu,v) 129 gvu,gvv = lp.dot_VV(gv,u),lp.dot_VV(gv,v) 130 return np.minimum(guv/gvv,gvu/guu) 131 132 def angle(self,u,v): 133 r""" 134 Generalized unoriented angle defined by the norm, 135 defined as $\measuredangle_N(u,v) := \arccos(\cos_N(u,v))$. 136 See `cos` member functions. 137 """ 138 c = ad.asarray(self.cos(u,v)) 139 mask=c < -1. 140 c[mask]=0. 141 result = ad.asarray(np.arccos(c)) 142 result[mask]=np.inf 143 return result 144 145# ---- Geometric transformations ---- 146 147 def inv_transform(self,a): 148 """ 149 Affine transformation of the norm. 150 The new unit ball is the inverse image of the previous one. 151 """ 152 raise NotImplementedError("Affine transformation not implemented for this norm") 153 154 def transform(self,a): 155 """ 156 Affine transformation of the norm. 157 The new unit ball is the direct image of the previous one. 158 """ 159 return self.inv_transform(lp.inverse(a)) 160 161 def rotate(self,r): 162 """ 163 Rotation of the norm, by a given rotation matrix. 164 The new unit ball is the direct image of the previous one. 165 """ 166 # For orthogonal matrices, inversion is transposition 167 return self.inv_transform(lp.transpose(r)) 168 169 def rotate_by(self,*args,**kwargs): 170 """ 171 Rotation of the norm, based on rotation parameters : angle (and axis in 3D). 172 """ 173 return self.rotate(lp.rotation(*args,**kwargs)) 174 175 def with_costs(self,costs): 176 """ 177 Produces a norm $N'$ defined by 178 $$ 179 N'(x) = N(costs * x) 180 $$ 181 where the multiplication is elementwise. 182 """ 183 a = np.zeros_like(costs,shape=(len(costs),)+costs.shape) 184 for i,cost in enumerate(costs): a[i,i] = cost 185 return self.inv_transform(a) 186 187 def with_speeds(self,speeds): 188 """ 189 Produces a norm $N'$ obeying 190 $$ 191 N'(x) = N(x/speeds) 192 $$ 193 where the division is elementwise. 194 """ 195 return self.with_costs(1./speeds) 196 197 def with_cost(self,cost): 198 """ 199 Produces a norm $N'$ obeying $N'(x) = N(cost*x)$. 200 """ 201 cost = ad.asarray(cost) 202 costs = np.broadcast_to(cost,(self.vdim,)+cost.shape) 203 return self.with_costs(costs) 204 205 206 def with_speed(self,speed): 207 """ 208 Produces a norm $N'$ obeying $N'(x) = N(x/speed)$. 209 """ 210 return self.with_cost(1/speed) 211 212# ---- Import and export ---- 213 214 def flatten(self): 215 """ 216 Flattens and concatenate the member fields into a single array. 217 """ 218 raise NotImplementedError("Flattening not implemented for this norm") 219 220 @classmethod 221 def expand(cls,arr): 222 """ 223 Inverse of the flatten member function. 224 Turns a suitable array into a metric. 225 """ 226 raise NotImplementedError("Expansion not implemented for this norm") 227 228 def to_HFM(self): 229 """ 230 Formats a metric for the HFM library. 231 This may include flattening some symmetric matrices, 232 concatenating with vector fields, and moving the first axis last. 233 """ 234 return np.moveaxis(self.flatten(),0,-1) 235 236 def model_HFM(self): 237 """ 238 The name of the 'model' for parameter, as input to the HFM library. 239 """ 240 raise NotImplementedError("HFM name is not specified for this norm") 241 242 @classmethod 243 def from_HFM(cls,arr): 244 """ 245 Inverse of the to_HFM member function. 246 Turns a suitable array into a metric. 247 """ 248 return cls.expand(np.moveaxis(arr,-1,0)) 249 250 def __iter__(self): 251 """ 252 Iterate over the member fields. 253 """ 254 raise NotImplementedError("__iter__ not implemented for this norm") 255 256 @classmethod 257 def from_generator(cls,gen): 258 """ 259 Produce a metric from a suitable generator expression. 260 """ 261 return cls(*gen) 262 263 @classmethod 264 def from_cast(cls,metric): 265 """ 266 Produces a metric by casting another metric of a compatible type. 267 """ 268 raise NotImplementedError("from_cast not implemented for this norm") 269 270 @property 271 def array_float_caster(self): 272 """ 273 Returns a caster function, which can be used to turn lists, etc, into 274 arrays with the suitable floating point type, from the suitable library 275 (numpy or cupy), depending on the member fields. 276 """ 277 return ad.cupy_generic.array_float_caster(self,iterables=type(self)) 278 279# ---- Related with Lagrandian and Hamiltonian interpretation ---- 280 281 def norm2(self,v): 282 r""" 283 Half square of the `norm`, defined by $F(v) := \frac 1 2 N(v)^2$. 284 """ 285 n = self.norm(v) 286 return 0.5*n**2 287 288 def gradient2(self,v): 289 """ 290 Gradient of `norm2`the half squared norm. 291 """ 292 g = self.gradient(v) 293 return lp.dot_VV(g,v)*g 294 295 def set_interpolation(self,grid,**kwargs): 296 """ 297 Sets interpolation_data, required to specialize the norm 298 at a given position. 299 300 Inputs: 301 - grid (optional). Coordinate system (required on first call). 302 - kwargs. Passed to UniformGridInterpolation (includes order) 303 """ 304 vdim = len(grid) 305 try: assert self.vdim == vdim 306 except ValueError: pass # Constant isotropic metrics have no dimension 307 308 def make_interp(value): 309 if hasattr(value,'shape') and value.shape[-vdim:]==grid.shape[1:]: 310 return Interpolation.UniformGridInterpolation(grid,value,**kwargs) 311 return value 312 313 self.interpolation_data = tuple(make_interp(value) for value in self) 314 315 def at(self,x): 316 """ 317 Interpolates the metric to a given position, on a grid given beforehand. 318 319 Inputs : 320 - x. Place where interpolation is needed. 321 """ 322 return self.from_generator( 323 field(x) if callable(field) else field 324 for field in self.interpolation_data) 325 # isinstance(field,Interpolation.UniformGridInterpolation) 326 327 328# ---- Proximal methods ---- 329 def make_proj_dual(self,**kwargs): 330 """ 331 Returns the orthogonal projection operator onto the unit ball of the dual metric. 332 """ 333 raise NotImplementedError("Sorry, make_proj_dual is not implemented for this metric") 334 335 def make_proj(self,**kwargs): 336 """ 337 Returns the orthogonal projection operator onto the unit ball. 338 - **kwargs : passed to self.make_proj_dual 339 """ 340 return self.dual().make_proj_dual(**kwargs) 341 342 def make_prox(self,**kwargs): 343 """ 344 Returns the proximal operator associated with the metric. 345 argmin_x (1/2) |x-y|^2 + τ F(y) 346 - **kwargs : passed to self.make_proj_dual 347 """ 348 # Implementation is based on the Moreau formula x = prox_{τ f}(x) + τ prox_{f^*/τ} (x/τ) 349 # noting that f^* is the characteristic function of the dual unit ball, hence 350 # prox_{f^*/τ} is the orthogonal projection onto the dual unit ball. 351 proj = self.make_proj_dual(**kwargs) 352 return lambda x,τ=1.: x - τ*proj(x/τ) 353
class
Base:
12class Base: 13 """ 14 Base class for a metric, in other words a family of norms. 15 """ 16 17 def __repr__(self): 18 return type(self).__name__ + repr(tuple(self)) 19 20 def __str__(self): 21 return type(self).__name__ + str(tuple(self)) 22 23 def norm(self,v): 24 r""" 25 Norm or quasi-norm defined by the class, often denoted $N$ in mathematical 26 formulas. Unless incorrect data is provided, this member function obeys, 27 for all vectors $u,v\in R^d$ and all $\alpha \geq 0$ 28 - $N(u+v) \leq N(u)+N(v)$ 29 - $N(\alpha u) = \alpha N(u)$ 30 - $N(u)\geq 0$ with equality iff $u=0$. 31 32 Broadcasting will occur depending on the shape of $v$ and of the class data. 33 """ 34 raise NotImplementedError("""Error : norm must be specialized in subclass""") 35 36 def gradient(self,v): 37 """ 38 Gradient of the `norm` defined by the class. 39 """ 40 if ad.is_ad(v) or ad.is_ad(self,iterables=(type(self),)): 41 v_dis = ad.disassociate(v,shape_bound=v.shape[1:]) 42 grad_dis = self.disassociate().gradient(v_dis) 43 return ad.associate(grad_dis) 44 45 v_ad = ad.Dense.identity(constant=v,shape_free=(len(v),)) 46 return np.moveaxis(self.norm(v_ad).coef,-1,0) 47 48 def dual(self): 49 r""" 50 Dual `norm`, mathematically defined by 51 $N^*(x) = max\\{ < x, y> ; N(y)\leq 1 \\}$ 52 """ 53 raise NotImplementedError("dual is not implemented for this norm") 54 55 @property 56 def vdim(self): 57 """ 58 The ambient vector space dimension, often denoted $d$ in mathematical formulas. 59 """ 60 raise NotImplementedError("vdim is not implemented for this norm") 61 62 @property 63 def shape(self): 64 """ 65 Dimensions of the underlying domain. 66 Expected to be the empty tuple, or a tuple of length `vdim`. 67 """ 68 raise NotImplementedError("Shape not implemented for this norm") 69 70 def disassociate(self): 71 """ 72 Hide the automatic differentiation (AD) information of the member fields. 73 See AutomaticDifferentiation.disassociate 74 """ 75 def dis(x): 76 if ad.isndarray(x) and x.shape[-self.vdim:]==self.shape: 77 return ad.disassociate(x,shape_bound=self.shape) 78 return x 79 return self.from_generator(dis(x) for x in self) 80 81# ---- Well posedness related methods ----- 82 def is_definite(self): 83 """ 84 Attempts to check wether the data defines a mathematically valid `norm`. 85 """ 86 raise NotImplementedError("is_definite is not implemented for this norm") 87 88 def anisotropy(self): 89 r""" 90 Anisotropy ratio of the `norm` denoted $N$. 91 Defined as 92 $$ 93 \max_{|u| = |v| = 1} \frac {N(u)}{N(v)}. 94 $$ 95 """ 96 raise NotImplementedError("anisotropy is not implemented for this norm") 97 98 def anisotropy_bound(self): 99 """ 100 An upper bound on the `anisotropy` of the norm. 101 """ 102 return self.anisotropy() 103 def cost_bound(self): 104 """ 105 Upper bound on $N(u)$, for any unit vector $u$, where $N$ is the `norm` 106 defined by the class. 107 """ 108 raise NotImplementedError("cost_bound is not implemented for this norm") 109 110# ---- Causality and acuteness related methods ---- 111 112 def cos_asym(self,u,v): 113 r""" 114 Generalized asymmetric cosine defined by the norm $N$, of the angle between two vectors $u,v$. 115 Defined as 116 $\cos^a_N(u,v) := < \nabla N(u), v> / N(v)$. 117 """ 118 u,v=(ad.asarray(e) for e in (u,v)) 119 return lp.dot_VV(self.gradient(u),v)/self.norm(v) 120 121 def cos(self,u,v): 122 r""" 123 Generalized cosine defined by the norm $N$. Expression 124 $\cos_N(u,v) := \min( \cos^a_N (u,v), \cos^a_N(v,u))$, 125 where $\cos^a_N=$`cos_asym`. 126 """ 127 u,v=(ad.asarray(e) for e in (u,v)) 128 gu,gv=self.gradient(u),self.gradient(v) 129 guu,guv = lp.dot_VV(gu,u),lp.dot_VV(gu,v) 130 gvu,gvv = lp.dot_VV(gv,u),lp.dot_VV(gv,v) 131 return np.minimum(guv/gvv,gvu/guu) 132 133 def angle(self,u,v): 134 r""" 135 Generalized unoriented angle defined by the norm, 136 defined as $\measuredangle_N(u,v) := \arccos(\cos_N(u,v))$. 137 See `cos` member functions. 138 """ 139 c = ad.asarray(self.cos(u,v)) 140 mask=c < -1. 141 c[mask]=0. 142 result = ad.asarray(np.arccos(c)) 143 result[mask]=np.inf 144 return result 145 146# ---- Geometric transformations ---- 147 148 def inv_transform(self,a): 149 """ 150 Affine transformation of the norm. 151 The new unit ball is the inverse image of the previous one. 152 """ 153 raise NotImplementedError("Affine transformation not implemented for this norm") 154 155 def transform(self,a): 156 """ 157 Affine transformation of the norm. 158 The new unit ball is the direct image of the previous one. 159 """ 160 return self.inv_transform(lp.inverse(a)) 161 162 def rotate(self,r): 163 """ 164 Rotation of the norm, by a given rotation matrix. 165 The new unit ball is the direct image of the previous one. 166 """ 167 # For orthogonal matrices, inversion is transposition 168 return self.inv_transform(lp.transpose(r)) 169 170 def rotate_by(self,*args,**kwargs): 171 """ 172 Rotation of the norm, based on rotation parameters : angle (and axis in 3D). 173 """ 174 return self.rotate(lp.rotation(*args,**kwargs)) 175 176 def with_costs(self,costs): 177 """ 178 Produces a norm $N'$ defined by 179 $$ 180 N'(x) = N(costs * x) 181 $$ 182 where the multiplication is elementwise. 183 """ 184 a = np.zeros_like(costs,shape=(len(costs),)+costs.shape) 185 for i,cost in enumerate(costs): a[i,i] = cost 186 return self.inv_transform(a) 187 188 def with_speeds(self,speeds): 189 """ 190 Produces a norm $N'$ obeying 191 $$ 192 N'(x) = N(x/speeds) 193 $$ 194 where the division is elementwise. 195 """ 196 return self.with_costs(1./speeds) 197 198 def with_cost(self,cost): 199 """ 200 Produces a norm $N'$ obeying $N'(x) = N(cost*x)$. 201 """ 202 cost = ad.asarray(cost) 203 costs = np.broadcast_to(cost,(self.vdim,)+cost.shape) 204 return self.with_costs(costs) 205 206 207 def with_speed(self,speed): 208 """ 209 Produces a norm $N'$ obeying $N'(x) = N(x/speed)$. 210 """ 211 return self.with_cost(1/speed) 212 213# ---- Import and export ---- 214 215 def flatten(self): 216 """ 217 Flattens and concatenate the member fields into a single array. 218 """ 219 raise NotImplementedError("Flattening not implemented for this norm") 220 221 @classmethod 222 def expand(cls,arr): 223 """ 224 Inverse of the flatten member function. 225 Turns a suitable array into a metric. 226 """ 227 raise NotImplementedError("Expansion not implemented for this norm") 228 229 def to_HFM(self): 230 """ 231 Formats a metric for the HFM library. 232 This may include flattening some symmetric matrices, 233 concatenating with vector fields, and moving the first axis last. 234 """ 235 return np.moveaxis(self.flatten(),0,-1) 236 237 def model_HFM(self): 238 """ 239 The name of the 'model' for parameter, as input to the HFM library. 240 """ 241 raise NotImplementedError("HFM name is not specified for this norm") 242 243 @classmethod 244 def from_HFM(cls,arr): 245 """ 246 Inverse of the to_HFM member function. 247 Turns a suitable array into a metric. 248 """ 249 return cls.expand(np.moveaxis(arr,-1,0)) 250 251 def __iter__(self): 252 """ 253 Iterate over the member fields. 254 """ 255 raise NotImplementedError("__iter__ not implemented for this norm") 256 257 @classmethod 258 def from_generator(cls,gen): 259 """ 260 Produce a metric from a suitable generator expression. 261 """ 262 return cls(*gen) 263 264 @classmethod 265 def from_cast(cls,metric): 266 """ 267 Produces a metric by casting another metric of a compatible type. 268 """ 269 raise NotImplementedError("from_cast not implemented for this norm") 270 271 @property 272 def array_float_caster(self): 273 """ 274 Returns a caster function, which can be used to turn lists, etc, into 275 arrays with the suitable floating point type, from the suitable library 276 (numpy or cupy), depending on the member fields. 277 """ 278 return ad.cupy_generic.array_float_caster(self,iterables=type(self)) 279 280# ---- Related with Lagrandian and Hamiltonian interpretation ---- 281 282 def norm2(self,v): 283 r""" 284 Half square of the `norm`, defined by $F(v) := \frac 1 2 N(v)^2$. 285 """ 286 n = self.norm(v) 287 return 0.5*n**2 288 289 def gradient2(self,v): 290 """ 291 Gradient of `norm2`the half squared norm. 292 """ 293 g = self.gradient(v) 294 return lp.dot_VV(g,v)*g 295 296 def set_interpolation(self,grid,**kwargs): 297 """ 298 Sets interpolation_data, required to specialize the norm 299 at a given position. 300 301 Inputs: 302 - grid (optional). Coordinate system (required on first call). 303 - kwargs. Passed to UniformGridInterpolation (includes order) 304 """ 305 vdim = len(grid) 306 try: assert self.vdim == vdim 307 except ValueError: pass # Constant isotropic metrics have no dimension 308 309 def make_interp(value): 310 if hasattr(value,'shape') and value.shape[-vdim:]==grid.shape[1:]: 311 return Interpolation.UniformGridInterpolation(grid,value,**kwargs) 312 return value 313 314 self.interpolation_data = tuple(make_interp(value) for value in self) 315 316 def at(self,x): 317 """ 318 Interpolates the metric to a given position, on a grid given beforehand. 319 320 Inputs : 321 - x. Place where interpolation is needed. 322 """ 323 return self.from_generator( 324 field(x) if callable(field) else field 325 for field in self.interpolation_data) 326 # isinstance(field,Interpolation.UniformGridInterpolation) 327 328 329# ---- Proximal methods ---- 330 def make_proj_dual(self,**kwargs): 331 """ 332 Returns the orthogonal projection operator onto the unit ball of the dual metric. 333 """ 334 raise NotImplementedError("Sorry, make_proj_dual is not implemented for this metric") 335 336 def make_proj(self,**kwargs): 337 """ 338 Returns the orthogonal projection operator onto the unit ball. 339 - **kwargs : passed to self.make_proj_dual 340 """ 341 return self.dual().make_proj_dual(**kwargs) 342 343 def make_prox(self,**kwargs): 344 """ 345 Returns the proximal operator associated with the metric. 346 argmin_x (1/2) |x-y|^2 + τ F(y) 347 - **kwargs : passed to self.make_proj_dual 348 """ 349 # Implementation is based on the Moreau formula x = prox_{τ f}(x) + τ prox_{f^*/τ} (x/τ) 350 # noting that f^* is the characteristic function of the dual unit ball, hence 351 # prox_{f^*/τ} is the orthogonal projection onto the dual unit ball. 352 proj = self.make_proj_dual(**kwargs) 353 return lambda x,τ=1.: x - τ*proj(x/τ)
Base class for a metric, in other words a family of norms.
def
norm(self, v):
23 def norm(self,v): 24 r""" 25 Norm or quasi-norm defined by the class, often denoted $N$ in mathematical 26 formulas. Unless incorrect data is provided, this member function obeys, 27 for all vectors $u,v\in R^d$ and all $\alpha \geq 0$ 28 - $N(u+v) \leq N(u)+N(v)$ 29 - $N(\alpha u) = \alpha N(u)$ 30 - $N(u)\geq 0$ with equality iff $u=0$. 31 32 Broadcasting will occur depending on the shape of $v$ and of the class data. 33 """ 34 raise NotImplementedError("""Error : norm must be specialized in subclass""")
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
gradient(self, v):
36 def gradient(self,v): 37 """ 38 Gradient of the `norm` defined by the class. 39 """ 40 if ad.is_ad(v) or ad.is_ad(self,iterables=(type(self),)): 41 v_dis = ad.disassociate(v,shape_bound=v.shape[1:]) 42 grad_dis = self.disassociate().gradient(v_dis) 43 return ad.associate(grad_dis) 44 45 v_ad = ad.Dense.identity(constant=v,shape_free=(len(v),)) 46 return np.moveaxis(self.norm(v_ad).coef,-1,0)
Gradient of the norm defined by the class.
def
dual(self):
48 def dual(self): 49 r""" 50 Dual `norm`, mathematically defined by 51 $N^*(x) = max\\{ < x, y> ; N(y)\leq 1 \\}$ 52 """ 53 raise NotImplementedError("dual is not implemented for this norm")
Dual norm, mathematically defined by
$N^*(x) = max\{ < x, y> ; N(y)\leq 1 \}$
vdim
55 @property 56 def vdim(self): 57 """ 58 The ambient vector space dimension, often denoted $d$ in mathematical formulas. 59 """ 60 raise NotImplementedError("vdim is not implemented for this norm")
The ambient vector space dimension, often denoted $d$ in mathematical formulas.
shape
62 @property 63 def shape(self): 64 """ 65 Dimensions of the underlying domain. 66 Expected to be the empty tuple, or a tuple of length `vdim`. 67 """ 68 raise NotImplementedError("Shape not implemented for this norm")
Dimensions of the underlying domain.
Expected to be the empty tuple, or a tuple of length vdim.
def
disassociate(self):
70 def disassociate(self): 71 """ 72 Hide the automatic differentiation (AD) information of the member fields. 73 See AutomaticDifferentiation.disassociate 74 """ 75 def dis(x): 76 if ad.isndarray(x) and x.shape[-self.vdim:]==self.shape: 77 return ad.disassociate(x,shape_bound=self.shape) 78 return x 79 return self.from_generator(dis(x) for x in self)
Hide the automatic differentiation (AD) information of the member fields. See AutomaticDifferentiation.disassociate
def
is_definite(self):
82 def is_definite(self): 83 """ 84 Attempts to check wether the data defines a mathematically valid `norm`. 85 """ 86 raise NotImplementedError("is_definite is not implemented for this norm")
Attempts to check wether the data defines a mathematically valid norm.
def
anisotropy(self):
88 def anisotropy(self): 89 r""" 90 Anisotropy ratio of the `norm` denoted $N$. 91 Defined as 92 $$ 93 \max_{|u| = |v| = 1} \frac {N(u)}{N(v)}. 94 $$ 95 """ 96 raise NotImplementedError("anisotropy is not implemented for this norm")
Anisotropy ratio of the norm denoted $N$.
Defined as
$$
\max_{|u| = |v| = 1} \frac {N(u)}{N(v)}.
$$
def
anisotropy_bound(self):
98 def anisotropy_bound(self): 99 """ 100 An upper bound on the `anisotropy` of the norm. 101 """ 102 return self.anisotropy()
An upper bound on the anisotropy of the norm.
def
cost_bound(self):
103 def cost_bound(self): 104 """ 105 Upper bound on $N(u)$, for any unit vector $u$, where $N$ is the `norm` 106 defined by the class. 107 """ 108 raise NotImplementedError("cost_bound is not implemented for this norm")
Upper bound on $N(u)$, for any unit vector $u$, where $N$ is the norm
defined by the class.
def
cos_asym(self, u, v):
112 def cos_asym(self,u,v): 113 r""" 114 Generalized asymmetric cosine defined by the norm $N$, of the angle between two vectors $u,v$. 115 Defined as 116 $\cos^a_N(u,v) := < \nabla N(u), v> / N(v)$. 117 """ 118 u,v=(ad.asarray(e) for e in (u,v)) 119 return lp.dot_VV(self.gradient(u),v)/self.norm(v)
Generalized asymmetric cosine defined by the norm $N$, of the angle between two vectors $u,v$. Defined as $\cos^a_N(u,v) := < \nabla N(u), v> / N(v)$.
def
cos(self, u, v):
121 def cos(self,u,v): 122 r""" 123 Generalized cosine defined by the norm $N$. Expression 124 $\cos_N(u,v) := \min( \cos^a_N (u,v), \cos^a_N(v,u))$, 125 where $\cos^a_N=$`cos_asym`. 126 """ 127 u,v=(ad.asarray(e) for e in (u,v)) 128 gu,gv=self.gradient(u),self.gradient(v) 129 guu,guv = lp.dot_VV(gu,u),lp.dot_VV(gu,v) 130 gvu,gvv = lp.dot_VV(gv,u),lp.dot_VV(gv,v) 131 return np.minimum(guv/gvv,gvu/guu)
Generalized cosine defined by the norm $N$. Expression
$\cos_N(u,v) := \min( \cos^a_N (u,v), \cos^a_N(v,u))$,
where $\cos^a_N=$cos_asym.
def
angle(self, u, v):
133 def angle(self,u,v): 134 r""" 135 Generalized unoriented angle defined by the norm, 136 defined as $\measuredangle_N(u,v) := \arccos(\cos_N(u,v))$. 137 See `cos` member functions. 138 """ 139 c = ad.asarray(self.cos(u,v)) 140 mask=c < -1. 141 c[mask]=0. 142 result = ad.asarray(np.arccos(c)) 143 result[mask]=np.inf 144 return result
Generalized unoriented angle defined by the norm,
defined as $\measuredangle_N(u,v) := \arccos(\cos_N(u,v))$.
See cos member functions.
def
inv_transform(self, a):
148 def inv_transform(self,a): 149 """ 150 Affine transformation of the norm. 151 The new unit ball is the inverse image of the previous one. 152 """ 153 raise NotImplementedError("Affine transformation not implemented for this norm")
Affine transformation of the norm. The new unit ball is the inverse image of the previous one.
def
transform(self, a):
155 def transform(self,a): 156 """ 157 Affine transformation of the norm. 158 The new unit ball is the direct image of the previous one. 159 """ 160 return self.inv_transform(lp.inverse(a))
Affine transformation of the norm. The new unit ball is the direct image of the previous one.
def
rotate(self, r):
162 def rotate(self,r): 163 """ 164 Rotation of the norm, by a given rotation matrix. 165 The new unit ball is the direct image of the previous one. 166 """ 167 # For orthogonal matrices, inversion is transposition 168 return self.inv_transform(lp.transpose(r))
Rotation of the norm, by a given rotation matrix. The new unit ball is the direct image of the previous one.
def
rotate_by(self, *args, **kwargs):
170 def rotate_by(self,*args,**kwargs): 171 """ 172 Rotation of the norm, based on rotation parameters : angle (and axis in 3D). 173 """ 174 return self.rotate(lp.rotation(*args,**kwargs))
Rotation of the norm, based on rotation parameters : angle (and axis in 3D).
def
with_costs(self, costs):
176 def with_costs(self,costs): 177 """ 178 Produces a norm $N'$ defined by 179 $$ 180 N'(x) = N(costs * x) 181 $$ 182 where the multiplication is elementwise. 183 """ 184 a = np.zeros_like(costs,shape=(len(costs),)+costs.shape) 185 for i,cost in enumerate(costs): a[i,i] = cost 186 return self.inv_transform(a)
Produces a norm $N'$ defined by $$ N'(x) = N(costs * x) $$ where the multiplication is elementwise.
def
with_speeds(self, speeds):
188 def with_speeds(self,speeds): 189 """ 190 Produces a norm $N'$ obeying 191 $$ 192 N'(x) = N(x/speeds) 193 $$ 194 where the division is elementwise. 195 """ 196 return self.with_costs(1./speeds)
Produces a norm $N'$ obeying $$ N'(x) = N(x/speeds) $$ where the division is elementwise.
def
with_cost(self, cost):
198 def with_cost(self,cost): 199 """ 200 Produces a norm $N'$ obeying $N'(x) = N(cost*x)$. 201 """ 202 cost = ad.asarray(cost) 203 costs = np.broadcast_to(cost,(self.vdim,)+cost.shape) 204 return self.with_costs(costs)
Produces a norm $N'$ obeying $N'(x) = N(cost*x)$.
def
with_speed(self, speed):
207 def with_speed(self,speed): 208 """ 209 Produces a norm $N'$ obeying $N'(x) = N(x/speed)$. 210 """ 211 return self.with_cost(1/speed)
Produces a norm $N'$ obeying $N'(x) = N(x/speed)$.
def
flatten(self):
215 def flatten(self): 216 """ 217 Flattens and concatenate the member fields into a single array. 218 """ 219 raise NotImplementedError("Flattening not implemented for this norm")
Flattens and concatenate the member fields into a single array.
@classmethod
def
expand(cls, arr):
221 @classmethod 222 def expand(cls,arr): 223 """ 224 Inverse of the flatten member function. 225 Turns a suitable array into a metric. 226 """ 227 raise NotImplementedError("Expansion not implemented for this norm")
Inverse of the flatten member function. Turns a suitable array into a metric.
def
to_HFM(self):
229 def to_HFM(self): 230 """ 231 Formats a metric for the HFM library. 232 This may include flattening some symmetric matrices, 233 concatenating with vector fields, and moving the first axis last. 234 """ 235 return np.moveaxis(self.flatten(),0,-1)
Formats a metric for the HFM library. This may include flattening some symmetric matrices, concatenating with vector fields, and moving the first axis last.
def
model_HFM(self):
237 def model_HFM(self): 238 """ 239 The name of the 'model' for parameter, as input to the HFM library. 240 """ 241 raise NotImplementedError("HFM name is not specified for this norm")
The name of the 'model' for parameter, as input to the HFM library.
@classmethod
def
from_HFM(cls, arr):
243 @classmethod 244 def from_HFM(cls,arr): 245 """ 246 Inverse of the to_HFM member function. 247 Turns a suitable array into a metric. 248 """ 249 return cls.expand(np.moveaxis(arr,-1,0))
Inverse of the to_HFM member function. Turns a suitable array into a metric.
@classmethod
def
from_generator(cls, gen):
257 @classmethod 258 def from_generator(cls,gen): 259 """ 260 Produce a metric from a suitable generator expression. 261 """ 262 return cls(*gen)
Produce a metric from a suitable generator expression.
@classmethod
def
from_cast(cls, metric):
264 @classmethod 265 def from_cast(cls,metric): 266 """ 267 Produces a metric by casting another metric of a compatible type. 268 """ 269 raise NotImplementedError("from_cast not implemented for this norm")
Produces a metric by casting another metric of a compatible type.
array_float_caster
271 @property 272 def array_float_caster(self): 273 """ 274 Returns a caster function, which can be used to turn lists, etc, into 275 arrays with the suitable floating point type, from the suitable library 276 (numpy or cupy), depending on the member fields. 277 """ 278 return ad.cupy_generic.array_float_caster(self,iterables=type(self))
Returns a caster function, which can be used to turn lists, etc, into arrays with the suitable floating point type, from the suitable library (numpy or cupy), depending on the member fields.
def
norm2(self, v):
282 def norm2(self,v): 283 r""" 284 Half square of the `norm`, defined by $F(v) := \frac 1 2 N(v)^2$. 285 """ 286 n = self.norm(v) 287 return 0.5*n**2
Half square of the norm, defined by $F(v) := \frac 1 2 N(v)^2$.
def
gradient2(self, v):
289 def gradient2(self,v): 290 """ 291 Gradient of `norm2`the half squared norm. 292 """ 293 g = self.gradient(v) 294 return lp.dot_VV(g,v)*g
Gradient of norm2the half squared norm.
def
set_interpolation(self, grid, **kwargs):
296 def set_interpolation(self,grid,**kwargs): 297 """ 298 Sets interpolation_data, required to specialize the norm 299 at a given position. 300 301 Inputs: 302 - grid (optional). Coordinate system (required on first call). 303 - kwargs. Passed to UniformGridInterpolation (includes order) 304 """ 305 vdim = len(grid) 306 try: assert self.vdim == vdim 307 except ValueError: pass # Constant isotropic metrics have no dimension 308 309 def make_interp(value): 310 if hasattr(value,'shape') and value.shape[-vdim:]==grid.shape[1:]: 311 return Interpolation.UniformGridInterpolation(grid,value,**kwargs) 312 return value 313 314 self.interpolation_data = tuple(make_interp(value) for value in self)
Sets interpolation_data, required to specialize the norm at a given position.
Inputs:
- grid (optional). Coordinate system (required on first call).
- kwargs. Passed to UniformGridInterpolation (includes order)
def
at(self, x):
316 def at(self,x): 317 """ 318 Interpolates the metric to a given position, on a grid given beforehand. 319 320 Inputs : 321 - x. Place where interpolation is needed. 322 """ 323 return self.from_generator( 324 field(x) if callable(field) else field 325 for field in self.interpolation_data) 326 # isinstance(field,Interpolation.UniformGridInterpolation)
Interpolates the metric to a given position, on a grid given beforehand.
Inputs :
- x. Place where interpolation is needed.
def
make_proj_dual(self, **kwargs):
330 def make_proj_dual(self,**kwargs): 331 """ 332 Returns the orthogonal projection operator onto the unit ball of the dual metric. 333 """ 334 raise NotImplementedError("Sorry, make_proj_dual is not implemented for this metric")
Returns the orthogonal projection operator onto the unit ball of the dual metric.
def
make_proj(self, **kwargs):
336 def make_proj(self,**kwargs): 337 """ 338 Returns the orthogonal projection operator onto the unit ball. 339 - **kwargs : passed to self.make_proj_dual 340 """ 341 return self.dual().make_proj_dual(**kwargs)
Returns the orthogonal projection operator onto the unit ball.
- **kwargs : passed to self.make_proj_dual
def
make_prox(self, **kwargs):
343 def make_prox(self,**kwargs): 344 """ 345 Returns the proximal operator associated with the metric. 346 argmin_x (1/2) |x-y|^2 + τ F(y) 347 - **kwargs : passed to self.make_proj_dual 348 """ 349 # Implementation is based on the Moreau formula x = prox_{τ f}(x) + τ prox_{f^*/τ} (x/τ) 350 # noting that f^* is the characteristic function of the dual unit ball, hence 351 # prox_{f^*/τ} is the orthogonal projection onto the dual unit ball. 352 proj = self.make_proj_dual(**kwargs) 353 return lambda x,τ=1.: x - τ*proj(x/τ)
Returns the proximal operator associated with the metric. argmin_x (1/2) |x-y|^2 + τ F(y)
- **kwargs : passed to self.make_proj_dual