Coverage for /opt/conda/envs/apienv/lib/python3.10/site-packages/daiquiri/core/transform/maps.py: 91%
216 statements
« prev ^ index » next coverage.py v7.6.5, created at 2024-11-15 02:12 +0000
1# -*- coding: utf-8 -*-
2"""
3Coordinate transformations on any domain (unbound or bound)
4can be derived from one of these base classes:
6 Map: f: D -> C
7 D is the domain to the map
8 C is the codomain to the map
9 f(D) ⊆ C is the image of the map
10 Injection: inverse f^(-1) not defined on the entire codomain
11 Surjection: inverse f^(-1) not unique (needs boundary conditions)
12 Bijection: inverse f^(-1) is defined on the entire codomain and unique
13 LinearMap: coordinate transformation can be represented
14 by a matrix M (m x n) multiplication
15 Isomorphism: linear and bijective, matrix M (n x n) is invertible
16 LinearCombination: linear and surjective with matrix M (1 x n)
17 AxisLinearCombination: linear combination with specific inverse
19An example
21 .. code-block:: python
23 class StepperPiezo(AxisLinearCombination):
24 \"""y(nm) = samy(mm) + sampy(um) + off(nm)
25 \"""
27 def __init__(self, offset=0):
28 coefficients = [1e6, 1e3]
29 limits = [[-1, 1], [0, 100]]
30 super().__init__(coefficients=coefficients, offset=offset, limits=limits)
31"""
33from abc import ABC, abstractmethod, abstractproperty
34import numpy
35import warnings
36from daiquiri.core.transform import sets
38FLOAT_ROUND_PRECISION = 5
41class TransformError(ValueError):
42 pass
45class DimensionError(TransformError):
46 pass
49class Map(ABC):
50 """f: D -> C: a==b then f(a)==f(b)"""
52 def __init__(self, domain=None, codomain=None, **kwargs):
53 """
54 :param Set domain:
55 :param Set codomain:
56 """
57 self._domain = domain
58 self._codomain = codomain
60 @property
61 def domain(self):
62 """Instance of Set"""
63 return self._domain
65 @property
66 def codomain(self):
67 """Instance of Set"""
68 return self._codomain
70 @property
71 def ndomain(self):
72 return self.domain.ndim
74 @property
75 def ncodomain(self):
76 return self.codomain.ndim
78 @abstractproperty
79 def image(self):
80 """Subspace of the codomain reached by this map
82 :returns Set:
83 """
84 pass
86 def __repr__(self):
87 return "{}:{}->{}".format(
88 self.__class__.__name__, str(self.domain), str(self.codomain)
89 )
91 def __call__(self, elements):
92 """
93 :param num or sequence elements:
94 :returns sequence:
95 """
96 return self.forward(elements)
98 def forward(self, elements):
99 """
100 :param num or sequence elements:
101 :returns sequence:
102 """
103 elements = self.domain.as_sequence(elements)
104 self.domain.raiseIfNotIsElement(elements)
105 return self._forward(elements)
107 @abstractmethod
108 def _forward(self, elements):
109 """
110 :param num or sequence elements:
111 :returns sequence:
112 """
113 pass
115 def inverse(self, elements):
116 """
117 :param num or sequence elements:
118 :returns sequence:
119 """
120 elements = self.codomain.as_sequence(elements)
121 self.image.raiseIfNotIsElement(elements)
122 return self._inverse(elements)
124 @abstractmethod
125 def _inverse(self, elements):
126 """
127 :param num or sequence elements:
128 :returns sequence:
129 """
130 pass
133class Surjection(Map):
134 """f: D -> C: x->f(x) and forall y in C there is
135 one (or more) x in R^n for which f(x) == y.
137 The inverse is not unique and needs boundary conditions.
138 """
140 @property
141 def image(self):
142 """Subspace of the codomain reached by this map.
143 For a surjection, image and codomain are identical.
145 :returns Set:
146 """
147 return self.codomain
149 def inverse(self, elements, **boundary_conditions):
150 """
151 :param num or sequence elements:
152 :param boundary_conditions: needed to select one solution
153 :returns sequence:
154 """
155 elements = self.codomain.as_sequence(elements)
156 self.codomain.raiseIfNotIsElement(elements)
157 return self._inverse(elements, **boundary_conditions)
159 def _inverse(self, elements, **boundary_conditions):
160 """
161 :param num or sequence elements:
162 param boundary_conditions: needed to select one solution
163 :returns sequence:
164 """
165 return self._right_inverse(elements, **boundary_conditions)
167 @abstractmethod
168 def _right_inverse(self, elements, **boundary_conditions):
169 """
170 :param num or sequence elements:
171 :param boundary_conditions: needed to select one solution
172 :returns sequence:
173 """
174 pass
177class Injection(Map):
178 """f: D -> C : x->f(x) and f(a) != f(b)
180 The inverse is not defined on the entire codomain.
181 """
183 def _inverse(self, elements):
184 """
185 :param num or sequence elements:
186 :returns sequence:
187 """
188 return self._left_inverse(elements)
190 @abstractmethod
191 def _left_inverse(self, elements):
192 """
193 :param num or sequence elements:
194 :returns sequence:
195 """
196 pass
199class Bijection(Surjection, Injection):
200 """f: D ⊆ R^n -> D ⊆ R^n: x->f(x) and the inverse exists and is unique"""
202 @property
203 def image(self):
204 return self.codomain
206 @abstractmethod
207 def _inverse(self, elements):
208 """
209 :param num or sequence elements:
210 :returns sequence:
211 """
212 pass
214 def _right_inverse(self, elements):
215 """
216 :param num or sequence elements:
217 :returns sequence:
218 """
219 return self._inverse(elements)
221 def _left_inverse(self, elements):
222 """
223 :param num or sequence elements:
224 :returns sequence:
225 """
226 return self._inverse(elements)
229class LinearMap(Map):
230 """This is an abstract class because the inverse transformation
231 does not always exist.
232 """
234 def __init__(self, matrix=None, **kwargs):
235 """
236 :param array-like matrix: shape (ncodomain, ndomain)
237 """
238 super().__init__(**kwargs)
239 self._matrix = self._parse_matrix(matrix)
241 @property
242 def matrix(self):
243 """
244 :returns array: shape is (ncodomain, ndomain)
245 """
246 return self._matrix
248 def _parse_matrix(self, value):
249 """
250 :param array-like matrix: shape (ncodomain, ndomain)
251 :returns array: shape (ncodomain, ndomain)
252 """
253 value = numpy.atleast_2d(value)
254 shape = self.ncodomain, self.ndomain
255 if value.shape != shape:
256 raise DimensionError(
257 "Transformation matrix should have shape {}".format(shape)
258 )
259 return value
261 def _forward(self, coord):
262 """
263 :param array coord: shape is (npoints, ndomain)
264 :returns array: shape is (npoints, ncodomain)
265 """
266 return self.matrix.dot(coord.T).T
269class InvertibleLinearMap(LinearMap):
270 """Linear map which has an invertible matrix (not necessarily bijective)"""
272 def _parse_matrix(self, value):
273 """
274 :param array-like matrix: shape (ncodomain, ndomain)
275 :returns array: shape (ncodomain, ndomain)
276 """
277 value = numpy.atleast_2d(value)
278 # Raises exception when not invertible
279 self.__imatrix = numpy.linalg.inv(value)
280 return value
282 @property
283 def inverse_matrix(self):
284 """
285 :returns array: shape is (ncodomain, ndomain)
286 """
287 return self.__imatrix
289 def _inverse(self, icoord):
290 """
291 :param num or array-like icoord: shape is (npoints, ncodomain)
292 :returns array-like: shape is (npoints, ndomain)
293 """
294 return self.inverse_matrix.dot(icoord.T).T
297class Isomorphism(InvertibleLinearMap, Bijection):
298 """f: D ⊆ R^n -> D ⊆ R^n: x->f(x) and f bijective and linear"""
300 pass
303class LinearInvertibleInjection(InvertibleLinearMap, Injection):
304 """f: D ⊂ R^n -> C ⊂ R^n: x->f(x) and f injective and linear"""
306 def _left_inverse(self, icoord):
307 """
308 :param num or array-like icoord: shape is (npoints, ncodomain)
309 :returns array-like: shape is (npoints, ndomain)
310 """
311 # TODO: This can never have worked
312 # return self._inverse(elements)
314 @property
315 def image(self):
316 # TODO: image of linear transformation is not a RealInterval
317 warnings.warn(
318 "Image is not implemented: use codomain (allows inverse transformation for coordinates outside the image)"
319 )
320 return self.codomain
323class LinearCombination(LinearMap, Surjection):
324 """f: D ⊆ R^n -> C ⊆ R: x->c[0]*x[0] + ... + c[-1]*x[-1] + offset"""
326 def __init__(self, coefficients=None, offset=0, **kwargs):
327 """
328 :param array-like coefficients: ndomain
329 :param num offset:
330 """
331 self.offset = offset
332 super().__init__(codomain=None, matrix=coefficients, **kwargs)
334 @property
335 def ncodomain(self):
336 return 1
338 @property
339 def codomain(self):
340 if isinstance(self.domain, sets.RealInterval):
341 coord = self.domain.border_inlim
342 icoord = self.forward(coord)
343 limits = [icoord.min(), icoord.max()]
344 return sets.RealInterval(ndim=1, limits=limits)
345 else:
346 warnings.warn(
347 "Codomain is not implemented: use R (allows inverse transformation for coordinates outside the image)"
348 )
349 return sets.RealVectorSpace(ndim=1)
351 def _forward(self, coord):
352 """
353 :param array coord: shape is (npoints, ndomain)
354 :returns array: shape is (npoints, ncodomain)
355 """
356 return super()._forward(coord) + self.offset
358 @property
359 def coefficients(self):
360 """
361 :returns array: ndomain
362 """
363 return self.matrix[0]
365 def _right_inverse(self, icoord, **boundary_conditions):
366 """
367 :param array icoord: shape is (npoints, 1)
368 :param ctype: boundary condition type
369 :param boundary_conditions: needed to select one solution
370 :returns array-like: shape is (npoints, ndomain)
371 """
372 icoord2 = icoord - self.offset
373 if self.ndomain == 1:
374 # This is actually a bijection
375 return icoord2 / self.coefficients[0]
376 else:
377 return self._lc_inverse(icoord2, **boundary_conditions)
379 @abstractmethod
380 def _lc_inverse(self, icoord, **boundary_conditions):
381 """Offset is already subtracted
383 :param array icoord: shape is (npoints, 1)
384 :param ctype: boundary condition type
385 :param boundary_conditions: needed to select one solution
386 :returns array-like: shape is (npoints, ndomain)
387 """
388 pass
391class AxisLinearCombination(LinearCombination):
392 """A linear combination with sets.RealInterval as domain and
393 specific boundary conditions for the inverse transformation
394 appropriate for parallel axes with different travel ranges.
396 For example:
397 y(nm) = samy(mm) + sampy(um) + offset(nm)
398 """
400 _smallest_fixed = False
401 _largest_fixed = False
403 def __init__(self, coefficients=None, offset=0, limits=None):
404 super().__init__(
405 domain=sets.RealInterval(limits=limits),
406 coefficients=coefficients,
407 offset=offset,
408 )
410 @property
411 def smallest_fixed(self) -> bool:
412 """Fix the smallest axis and center each time"""
413 return self._smallest_fixed
415 @smallest_fixed.setter
416 def smallest_fixed(self, value: bool):
417 self._smallest_fixed = value
419 @property
420 def largest_fixed(self) -> bool:
421 """Fix the largest axis"""
422 return self._largest_fixed
424 @largest_fixed.setter
425 def largest_fixed(self, value: bool):
426 self._largest_fixed = value
428 def codomain_limits(self, axis, reference=None):
429 """Axis limits projected on the codomain
431 :param num axis:
432 :param array reference: shape is (1, ndomain)
433 :returns array-like: shape is (2,)
434 """
435 reference = self._prepare_reference(reference)
436 low = self.domain.limits[axis].low
437 high = self.domain.limits[axis].high
438 pts = numpy.tile(reference, (2, 1))
439 pts[:, axis] = [low, high]
440 icoord = self.forward(pts)[:, 0]
441 icoord.sort()
442 return icoord
444 def _lc_inverse(self, icoord, ctype="single", **boundary_conditions):
445 """Offset is already subtracted
447 :param array icoord: shape is (npoints, 1)
448 :param ctype: boundary condition type
449 :param boundary_conditions: needed to select one solution
450 :returns array-like: shape is (npoints, ndomain)
451 """
452 if ctype == "single":
453 return self._inverse_solution_single(icoord, **boundary_conditions)
454 else:
455 raise ValueError("Unknown boundary condition")
457 @property
458 def default_reference(self):
459 """Default reference position to calculate the inverse"""
460 return self.domain.center
462 def _prepare_reference(self, reference=None):
463 """
464 :param array reference: shape is (1, ndomain)
465 :returns array-like: shape is (1, ndomain)
466 """
467 if reference is None:
468 reference = self.default_reference
469 reference = self.domain.as_sequence(reference).copy()
470 # fixed = self.domain.length == 0
471 # if fixed.any():
472 # reference[0][fixed] = self.domain.center[fixed]
473 return reference
475 def _inverse_solution_single(self, icoord, reference=None):
476 """Try to reach all points with and axis with the smallest
477 possible range. Center axes with a lower range and modify
478 axes with higher range when needed.
480 :param array icoord: shape is (npoints, 1)
481 :param array reference: shape is (1, ndomain)
482 :returns array-like: shape is (npoints, ndomain)
483 """
484 reference = self._prepare_reference(reference=reference)
485 # Range in the domain reached by each axis
486 ranges = self.domain.length_inlim
487 # Range in the domain covered by the points
488 coeff = self.coefficients
489 pcoord = icoord / coeff
490 pranges = pcoord.max(axis=0) - pcoord.min(axis=0)
491 # Potential for float rounding error here
492 # [0, 100.0000000001], [1, 100]
493 pranges = numpy.round(pranges, FLOAT_ROUND_PRECISION)
494 # Can all points be reached by one axis?
495 inrange = pranges <= ranges
496 if not inrange.any():
497 raise ValueError("Points are not within the limits of a single motor")
498 # Loop over all axes (smallest range first)
499 center = self.domain.center
500 idx = numpy.argsort(self.domain.length_inlim * coeff)
501 for i, axis in enumerate(idx):
502 # Center the axis
503 reference[0][axis] = center[axis]
504 if inrange[axis]:
505 if not (i == 0 and self._smallest_fixed and len(idx) > 1):
506 # Try to reach all points with one axis at the
507 # current reference position
508 coord = self._calc_inverse_single_axis(
509 icoord, axis, reference, coeff
510 )
511 if coord in self.domain:
512 return coord
514 if i == 0 and self._largest_fixed and len(idx) > 1:
515 raise ValueError("Cannot reach all points with largest axis fixed")
516 # Try to reach the centroid of all points by moving
517 # the reference position. Only move axes with a longer range
518 # than the current axis.
519 self._calc_inverse_single_point(
520 icoord.mean(), idx[i + 1 :], reference, coeff, center
521 )
522 # Try to reach all points by moving one axis at the
523 # current reference position
524 coord = self._calc_inverse_single_axis(icoord, axis, reference, coeff)
525 if coord in self.domain:
526 return coord
527 raise ValueError("Cannot reach all points with the given boundary conditions")
529 def _calc_inverse_single_axis(self, y, j, x, coeff):
530 """Solve y = sum_i(ci.xi) by fixing all but xj and calculate xj
532 :param array y: shape is (npoints, 1)
533 :param num axis:
534 :param array x: shape is (1, ndomain)
535 :param array coeff: shape is (ndomain, )
536 :returns array: shape is (npoints, ndomain)
537 """
538 x = numpy.tile(x, (y.size, 1))
539 yarr = x[0] * coeff
540 yarr[j] = 0
541 x[:, j] = (y[:, 0] - yarr.sum()) / coeff[j]
542 # To avoid potential float round issues (see above)
543 return numpy.round(x, FLOAT_ROUND_PRECISION)
545 def _calc_inverse_single_point(self, y, idx, x, coeff, center):
546 """Solve y = sum_i(ci.xi) by moving xj with j in `idx`
548 :param num y:
549 :param array idx:
550 :param array x: shape is (1, ndomain)
551 :param array coeff: shape is (ndomain, )
552 :param array center: shape is (ndomain, )
553 """
554 # TODO: this is not enough
555 for j in idx:
556 yarr = x[0] * coeff
557 yarr[j] = 0
558 x[0, j] = (y - yarr.sum()) / coeff[j]
559 if x in self.domain:
560 break
561 x[0, j] = center[j]