geometry.manifolds.differentiable

1 from .. import assert_allclose, formatm, printm, contract, new_contract 2 from abc import ABCMeta, abstractmethod 3 from collections import namedtuple

4 5 6 -class DifferentiableManifold(object):

7 ''' This is the base class for differentiable manifolds. ''' 8 __metaclass__ = ABCMeta 9

10 - def __init__(self, dimension):

11 self._isomorphisms = {} 12 self._embedding = {} 13 self._projection = {} 14 self.dimension = dimension 15 16 self.atol_distance = 1e-8 17 18 # Save reference, but do not create straight away. 19 self._tangent_bundle = None

20 21 @abstractmethod 22 @new_contract

23 - def belongs(self, x):

24 ''' 25 Raises an Exception if the point does not belong to this manifold. 26 27 This function wraps some checks around :py:func:`belongs_`, 28 which is implemented by the subclasses. 29 '''

30 31 @new_contract 32 @contract(bv='tuple(belongs, *)')

33 - def belongs_ts(self, bv):

34 ''' 35 Checks that a vector *vx* belongs to the tangent space 36 at the given point *base*. 37 38 ''' 39 bvp = self.project_ts(bv) 40 assert_allclose(bv[1], bvp[1], atol=self.atol_distance)

41 42 @abstractmethod 43 @contract(bv='tuple(belongs, *)')

44 - def project_ts(self, bv): # TODO: test

45 ''' 46 Projects a vector *v_ambient* in the ambient space 47 to the tangent space at point *base*. 48 '''

49 50 @abstractmethod 51 @contract(a='belongs', b='belongs', returns='>=0')

52 - def distance(self, a, b):

53 ''' 54 Computes the geodesic distance between two points. 55 '''

56 57 # @contract(returns='DifferentiableManifold') # Circular ref

58 - def tangent_bundle(self):

59 ''' Returns the manifold corresponding to the tangent bundle. 60 The default gives a generic implementation. 61 MatrixLieGroup have a different one. 62 ''' 63 if self._tangent_bundle is None: 64 from . import TangentBundle 65 self._tangent_bundle = TangentBundle(self) 66 return self._tangent_bundle

67 68 @abstractmethod 69 @contract(base='belongs', p='belongs', returns='belongs_ts')

70 - def logmap(self, base, p):

71 ''' 72 Computes the logarithmic map from base point *base* to target *b*. 73 # XXX: what should we do in the case there is more than one logmap? 74 '''

75 76 @abstractmethod 77 @contract(bv='belongs_ts', returns='belongs')

78 - def expmap(self, bv):

79 ''' 80 Computes the exponential map from *base* for the velocity 81 vector *v*. 82 83 This function wraps some checks around :py:func:`expmap_`, 84 which is implemented by the subclasses. 85 86 '''

87 88 @contract(returns='list(belongs)')

89 - def interesting_points(self):

90 ''' 91 Returns a list of "interesting points" on this manifold that 92 should be used for testing various properties. 93 ''' 94 return []

95 96 @contract(a='belongs', b='belongs', t='>=0,<=1', returns='belongs')

97 - def geodesic(self, a, b, t):

98 ''' Returns the point interpolated along the geodesic. ''' 99 bv = self.logmap(a, b) 100 return self.expmap((bv[0], bv[1] * t))

101 102 @contract(a='belongs')

103 - def friendly(self, a):

104 ''' 105 Returns a friendly description string for a point on the manifold. 106 ''' 107 return a.__str__()

108 109 @contract(a='belongs', b='belongs')

110 - def assert_close(self, a, b, atol=1e-8, msg=None):

111 ''' 112 Asserts that two points on the manifold are close to the given 113 tolerance. 114 ''' 115 distance = self.distance(a, b) 116 if msg is None: 117 msg = "" 118 if distance > atol: 119 msg += "\nThe two points should be the same:\n" 120 msg += "- a: %s\n" % self.friendly(a) 121 msg += "- b: %s\n" % self.friendly(b) 122 msg += formatm('a', a, 'b', b) 123 assert_allclose(distance, 0, atol=atol, err_msg=msg) 124 return distance

125 126 Isomorphism = namedtuple('Isomorphism', 127 'A B A_to_B B_to_A steps type desc') 128 Embedding = namedtuple('Embedding', 129 'A B A_to_B B_to_A steps type desc') 130 131 @staticmethod

132 - def isomorphism(A, B, A_to_B, B_to_A, itype='user', steps=None, desc=None):

133 if A.dimension != B.dimension: 134 msg = ('You are trying to define an isomorphism' 135 ' between manifolds of different dimension:\n' 136 '- %s has dimension %d;\n' 137 '- %s has dimension %d.\n' % (A, A.dimension, 138 B, B.dimension)) 139 raise ValueError(msg) 140 141 Iso = DifferentiableManifold.Isomorphism 142 if steps is None: 143 steps = [(A, '~', B)] 144 A._isomorphisms[B] = Iso(A, B, A_to_B, B_to_A, steps, itype, desc) 145 B._isomorphisms[A] = Iso(B, A, B_to_A, A_to_B, steps, itype, desc)

146 147 @staticmethod

148 - def embedding(A, B, A_to_B, B_to_A, itype='user', steps=None, desc=None):

149 if A.dimension > B.dimension: 150 msg = ('You are trying to define an embedding' 151 ' from a large to a smaller manifold:\n' 152 '- %s has dimension %d;\n' 153 '- %s has dimension %d.\n' % (A, A.dimension, 154 B, B.dimension)) 155 raise ValueError(msg) 156 157 Embed = DifferentiableManifold.Embedding 158 if steps is None: 159 steps = [(A, '=', B)] 160 A._embedding[B] = Embed(A, B, A_to_B, B_to_A, steps, itype, desc) 161 B._projection[A] = Embed(B, A, B_to_A, A_to_B, steps, itype, desc) 162 163 # TODO: move somewhere 164 if False: 165 #if development: 166 try: 167 for a in A.interesting_points(): 168 A.belongs(a) 169 b = A_to_B(a) 170 B.belongs(b) 171 except: 172 print('Invalid embedding:\n %s -> %s using %s' % (A, B, 173 A_to_B)) 174 printm('a', a) 175 raise 176 177 try: 178 for b in B.interesting_points(): 179 B.belongs(b) 180 a = B_to_A(b) 181 A.belongs(a) 182 except: 183 printm('b', b) 184 print('Invalid embedding:\n %s <- %s using %s' % 185 (A, B, B_to_A)) 186 raise

187

188 - def relations_descriptions(self):

189 s = ('[= %s >= %s <= %s]' % 190 (" ".join([str(a) for a in self._isomorphisms]), 191 " ".join([str(a) for a in self._projection]), 192 " ".join([str(a) for a in self._embedding])) 193 ) 194 return s

195 196 @contract(my_point='belongs')

197 - def embed_in(self, M, my_point):

198 ''' Embeds a point on this manifold to the target manifold M. ''' 199 #self.belongs(my_point) 200 if not self.embeddable_in(M): 201 msg = ('%s is not embeddable in %s; %s' % 202 (self, M, self.relations_descriptions())) 203 raise ValueError(msg) 204 x = self._embedding[M].A_to_B(my_point) 205 #M.belongs(x, msg='Error while embedding %s < %s point %s' % 206 # (self, M, my_point)) 207 return x

208 209 @contract(returns='belongs')

210 - def project_from(self, M, his_point):

211 ''' Projects a point on a bigger manifold to this manifold. ''' 212 if not self.embeddable_in(M): 213 msg = ('Cannot project from %s to %s; %s' % 214 (self, M, self.relations_descriptions())) 215 raise ValueError(msg) 216 x = self._embedding[M].B_to_A(his_point) 217 return x

218 219 @contract(my_point='belongs')

220 - def project_to(self, m, my_point):

221 if not self.can_represent(m): 222 msg = ('%s does not contain %s; %s' % 223 (self, m, self.relations_descriptions())) 224 raise ValueError(msg) 225 x = self._projection[m].A_to_B(my_point) 226 return x

227 228 @contract(my_point='belongs')

229 - def convert_to(self, m, my_point):

230 if not self.can_convert_to(m): 231 msg = ('%s cannot be converted to %s; %s' % 232 (self, m, self.relations_descriptions())) 233 raise ValueError(msg) 234 x = self._isomorphisms[m].A_to_B(my_point) 235 return x

236

237 - def can_convert_to(self, manifold):

238 return manifold in self._isomorphisms

239

240 - def can_represent(self, manifold): # XXX: change name

241 return manifold in self._projection 242

243 - def embeddable_in(self, manifold):

244 return manifold in self._embedding

245

246 - def get_dimension(self):

247 ''' Returns the intrinsic dimension of this manifold. ''' 248 return self.dimension

249 250 @contract(returns='belongs', yaml_structure='list|dict')

251 - def from_yaml(self, yaml_structure):

252 ''' Recovers a value from a Yaml structure. ''' 253 # TODO: explicit check 254 # TODO: add testing 255 from geometry.yaml import from_yaml 256 return from_yaml(yaml_structure)

257 258 @contract(x='belongs')

259 - def to_yaml(self, x):

260 # TODO: add testing 261 from geometry.yaml import to_yaml 262 return to_yaml('%s' % self, x)

263

264 265 -class RandomManifold(DifferentiableManifold):

266 ''' This is the base class for manifolds that have the ability 267 to sample random points. ''' 268 269 @abstractmethod

270 - def sample_uniform(self):

271 ''' 272 Samples a random point in this manifold according to the Haar 273 measure. Raises exception if the measure is improper (e.g., R^n). 274 '''

275 276 @abstractmethod 277 @contract(a='belongs')

278 - def sample_velocity(self, a):

279 ''' Samples a random velocity with length 1 at the base point a'''

280

Source Code for Module geometry.manifolds.differentiable_manifold