Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
"""Functions for inverting a 1D (time)coordinate transform given as numpy array
The inversion function internally requires an interpolation operator implementing the
RegularInterpolator interface, and which is provided by the user. Use
make_shift_inverse_lagrange_numpy to create an inversion operator employing Lagrange
interpolation.
"""
from __future__ import annotations
# ~ from dataclasses import dataclass
# ~ from enum import Enum
from typing import Final
import numpy as np
# ~ from typing_extensions import assert_never
from lisainstrument.fir_filters_numpy import NumpyArray1D
from lisainstrument.regular_interpolators import (
RegularInterpolator,
make_regular_interpolator_lagrange,
# ~ make_regular_interpolator_linear,
)
class ShiftInverseNumpy: # pylint: disable=too-few-public-methods
"""Invert coordinate transformation given as shift"""
def __init__(
self,
max_abs_shift: float,
interp: RegularInterpolator,
max_iter: int,
tolerance: float,
):
"""Set up interpolator.
Arguments:
max_abs_shift: Upper limit for absolute difference between coordinate frames w.r.t index space
interp: Interpolation method
max_iter: Maximum iterations before fail
tolerance: Maximum absolute error of result
"""
self._max_abs_shift: Final = float(max_abs_shift)
self._interp_np: Final = interp
self._max_iter = int(max_iter)
self._tolerance = float(tolerance)
@property
def _margin_left(self) -> int:
"""Left margin size.
Specifies how many samples on the left have to be added by boundary conditions.
"""
return self._interp_np.margin_left + self._max_abs_shift
@property
def _margin_right(self) -> int:
"""Right margin size.
Specifies how many samples on the right have to be added by boundary conditions.
"""
return self._interp_np.margin_right + self._max_abs_shift
def _fixed_point_iter(self, f, x):
for _ in range(self._max_iter):
x_next = f(x)
error = np.max(np.abs(x - x_next))
if error < self._tolerance:
return x_next
x = x_next
msg = (
f"ShiftInverseNumpy: iteration did not converge ({error=}, "
f"iterations={self._max_iter}, tolerance={self._tolerance})"
)
raise RuntimeError(msg)
def __call__(self, shift: np.ndarray) -> NumpyArray1D:
"""Find the shift for the inverse transformation given by a shift.
Arguments:
shift: 1D numpy array with shifts of the coordinate transform
Returns:
1D numpy array with shift for inverse coordinate transform
"""
shift_pad = np.pad(
shift,
(self._margin_left, self._margin_right),
mode="constant",
constant_values=(shift[0], shift[-1]),
)
def f_iter(x: np.ndarray) -> np.ndarray:
return self._interp_np.apply_shift(shift_pad, -x, self._margin_left)
return self._fixed_point_iter(f_iter, shift)
def make_shift_inverse_lagrange_numpy(
order: int,
max_abs_shift: float,
max_iter: int,
tolerance: float,
) -> ShiftInverseNumpy:
"""Set up ShiftInverseNumpy instance with Lagrange interpolation method.
Arguments:
order: Order of the Lagrange polynomials
max_abs_shift: Upper limit for absolute difference between coordinate frames w.r.t index space
max_iter: Maximum iterations before fail
tolerance: Maximum absolute error of result
Returns:
Inversion function of type ShiftInverseNumpy
"""
interp = make_regular_interpolator_lagrange(order)
return ShiftInverseNumpy(max_abs_shift, interp, max_iter, tolerance)