Skip to content
Snippets Groups Projects
Commit 8423b4d7 authored by Wolfgang Kastaun's avatar Wolfgang Kastaun
Browse files

Include newly written shift inversion modules in manual 

Change math markup from mkdocs to sphinx syntax
parent db14caf0
No related branches found
No related tags found
No related merge requests found
Pipeline #386814 passed
Stage: test
Stage: build
Hide whitespace changes
Inline Side-by-side
docs/index.rst
+ 1
− 0
View file @ 8423b4d7
......@@ -35,6 +35,7 @@ artifacts), etc.
hexagon
dsp
containers
shift_inversion
.. toctree::
:maxdepth: 1
......
docs/shift_inversion.rst 0 → 100644
+ 3
− 0
View file @ 8423b4d7
.. automodule:: lisainstrument.shift_inversion_numpy
.. automodule:: lisainstrument.shift_inversion_dask
lisainstrument/shift_inversion_dask.py
+ 27
− 8
View file @ 8423b4d7
"""Functions for applying dynamic real-valued shifts to dask arrays using Lagrange interpolation
"""
Dask-based time coordinate inversion
====================================
Functions for inverting a 1D (time)coordinate transform given as dask array.
This is the same functionality as provided by module shift_inversion_numpy for
numpy arrays.
.. autoclass:: ShiftInverseDask
:members:
:private-members:
:special-members:
There are two choices for the interpolator used internally. One is the
newly written Lagrange interpolator and the other is the Lagrange interpolator
from dsp.timeshift.
.. autofunction:: make_shift_inverse_lagrange_dask
.. autofunction:: make_shift_inverse_dsp_dask
Use make_dynamic_shift_lagrange_dask to create a Lagrange interpolator for dask arrays.
"""
from __future__ import annotations
......@@ -40,8 +59,8 @@ class ShiftInverseDask:
See shift_inversion_numpy.ShiftInverseNumpy for details
Arguments:
fsample: Sample rate $f_s > 0$ [s]
max_abs_shift: Upper limit $S_\mathrm{max} \ge 0 $ [s] for coordinate shift
fsample: Sample rate :math:`f_s > 0` [s]
max_abs_shift: Upper limit :math:`S_\mathrm{max} \ge 0` [s] for coordinate shift
interp: Interpolation operator
max_iter: Maximum iterations before fail
tolerance: Maximum absolute error [s] of result
......@@ -167,8 +186,8 @@ def make_shift_inverse_lagrange_dask(
Arguments:
order: Order of the Lagrange polynomials
fsample: Sample rate $f_s > 0$ [s]
max_abs_shift: Upper limit $S_\mathrm{max} \ge 0 $ [s] for coordinate shift
fsample: Sample rate :math:`f_s > 0` [s]
max_abs_shift: Upper limit :math:`S_\mathrm{max} \ge 0` [s] for coordinate shift
max_iter: Maximum iterations before fail
tolerance: Maximum absolute error of result
Returns:
......@@ -195,8 +214,8 @@ def make_shift_inverse_dsp_dask(
Arguments:
order: Order of the Lagrange polynomials
fsample: Sample rate $f_s > 0$ [s]
max_abs_shift: Upper limit $S_\mathrm{max} \ge 0 $ [s] for coordinate shift
fsample: Sample rate :math:`f_s > 0` [s]
max_abs_shift: Upper limit :math:`S_\mathrm{max} \ge 0` [s] for coordinate shift
max_iter: Maximum iterations before fail
tolerance: Maximum absolute error of result
......
lisainstrument/shift_inversion_numpy.py
+ 64
− 43
View file @ 8423b4d7
"""Functions for inverting a 1D (time)coordinate transform given as numpy array
"""
Time coordinate inversion
=========================
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.
.. autoclass:: ShiftInverseNumpy
:members:
:private-members:
:special-members:
.. autofunction:: make_shift_inverse_lagrange_numpy
.. autofunction:: fixed_point_iter
"""
from __future__ import annotations
......@@ -29,17 +47,20 @@ def fixed_point_iter(
r"""Perform a fixed point iteration for functions operating on a 1D array.
This uses fixed-point iteration to find a solution for
$$ x = f(x) $$
where $x$ is a 1D array and $f$ returns an array of the same size.
.. math::
x = f(x)
where :math:`x` is a 1D array and :math:`f` returns an array of the same size.
The convergence criterion is provided by the user
via a function $r(x)$ that returns a scalar error measure. The iteration
is performed until $r(x) < \epsilon$. If convergence is not achieved
via a function :math:`r(x)` that returns a scalar error measure. The iteration
is performed until :math:`r(x) < \epsilon`. If convergence is not achieved
after a given number of iterations, an exception is raised.
Arguments:
f: The function $f(x)$
ferr: The error measure $r(x)$
f: The function :math:`f(x)`
ferr: The error measure :math:`r(x)`
x: The initial data for the iteration.
tolerance: The tolerance $\epsilon$
max_iter: Maximum number of iterations
......@@ -64,25 +85,25 @@ class ShiftInverseNumpy: # pylint: disable=too-few-public-methods
r"""Invert time coordinate transformation given as numpy array
The purpose of this class is the inversion of a 1D coordinate transform
between some coordinates $u$ and $v$.
between some coordinates :math:`u` and :math:`v`.
The transform is expressed as a
shift
$ S(v) = u(v) - v $.
:math:`S(v) = u(v) - v`.
It will be provided as samples
$ S_k = S(v_k) $
:math:`S_k = S(v_k)`
at regularly spaced sample locations
$ v_k = v_0 + k \Delta v $.
:math:`v_k = v_0 + k \Delta v`.
The inverse will be expressed as
$\hat{S}(u) = u - v(u)$.
:math:`\hat{S}(u) = u - v(u)`.
It will be computed for locations
$ u_l = u_0 + l \Delta u $ regularly spaced with respect to the
$u$-coordinate, i.e. the output will be the sequence
$ \hat{S}_l = u_l - v(u_l) $.
:math:`u_l = u_0 + l \Delta u` regularly spaced with respect to the
:math:`u`-coordinate, i.e. the output will be the sequence
:math:`\hat{S}_l = u_l - v(u_l)`.
Currently, we restrict to the special case where $u_k = v_k$,
i.e. $v_0 = u_0$ and $\Delta u = \Delta v$.
Currently, we restrict to the special case where :math:`u_k = v_k`,
i.e. :math:`v_0 = u_0` and :math:`\Delta u = \Delta v`.
By convention, the coordinates refer to times, and coordinate shift,
tolerance, and sample rate are all given in SI units.
......@@ -98,12 +119,12 @@ class ShiftInverseNumpy: # pylint: disable=too-few-public-methods
):
r"""Set up the coordinate inversion operator.
One needs to provide the sample rate $f_s$ used
One needs to provide the sample rate :math:`f_s` used
both for input and output sequences. All calls to the operator will
assume regularly sampled sequences with $\Delta u = \Delta v = 1 / f_s$.
assume regularly sampled sequences with :math:`\Delta u = \Delta v = 1 / f_s`.
For technical reasons, one also needs to provide a limit for the maximum
shift between the coordinates, i.e. $|S(v)| < S_\mathrm{max}$.
shift between the coordinates, i.e. :math:`|S(v)| < S_\mathrm{max}`.
Another parameter is an interpolation operator to be used in the internal
fixed-point iteration algorithm. One also needs to provide the desired
......@@ -111,8 +132,8 @@ class ShiftInverseNumpy: # pylint: disable=too-few-public-methods
allowed number of iterations before aborting with an exception.
Arguments:
fsample: Sample rate $f_s > 0$ [s]
max_abs_shift: Upper limit $S_\mathrm{max} \ge 0 $ [s] for coordinate shift
fsample: Sample rate :math:`f_s > 0` [s]
max_abs_shift: Upper limit :math:`S_\mathrm{max} \ge 0` [s] for coordinate shift
interp: Interpolation operator
max_iter: Maximum iterations before fail
tolerance: Maximum absolute error [s] of result
......@@ -157,44 +178,44 @@ class ShiftInverseNumpy: # pylint: disable=too-few-public-methods
r"""Compute the inverse coordinate transform
The coordinate transform is given an array with the sequence
$ S_k = S(v_k) $. The sample locations $v_k$ do not have to be provided
:math:`S_k = S(v_k)`. The sample locations :math:`v_k` do not have to be provided
but are assumed to be regularly spaced, i.e.
$ v_k = v_0 + k \Delta v $, and the first array element corresponds to
$k=0$.
:math:`v_k = v_0 + k \Delta v`, and the first array element corresponds to
:math:`k=0`.
The output is given as an array with the sequence
$ \hat{S}_l $ (see constructor docstring) providing the shift
:math:`\hat{S}_l` (see constructor docstring) providing the shift
at sample locations regularly spaced in the transformed coordinate
$u$, that is, at points $ u_l = u_0 + l \Delta u $. The first
element of the output corresponds to $l=0$.
:math:`u`, that is, at points :math:`u_l = u_0 + l \Delta u`. The first
element of the output corresponds to :math:`l=0`.
The output sample locations are not returned, but are implicitly
equal to the input ones, meaning $ v_l = u_l $.
equal to the input ones, meaning :math:`v_l = u_l`.
The algorithm works using fixed point iteration
$ \hat{S}_l^{n+1} = f(l,\hat{S}_l^{n}) $,
:math:`\hat{S}_l^{n+1} = f(l,\hat{S}_l^{n})`,
where
$ f(l,d) = I[v_k, S_k](v_l - d) $, and $I[v_k, S_k]$ is a
function obtained by interpolating the samples $v_k, S_k$,
approximating $I[v_k, S_k](v) \approx S(v) $.
:math:`f(l,d) = I[v_k, S_k](v_l - d)`, and :math:`I[v_k, S_k]` is a
function obtained by interpolating the samples :math:`v_k, S_k`,
approximating :math:`I[v_k, S_k](v) \approx S(v)`.
On the technical level, the interpolation operator is implemented
using shifts directly, with an interface of the form
$I[S_k](d_l) \approx S(v_l + d_l)$.
:math:`I[S_k](d_l) \approx S(v_l + d_l)`.
If the iteration converges, it converges to a solution
$ \bar{S}_l = f(l,\bar{S}_l) \approx S(v_l - \bar{S}_l) = S(u_l - \bar{S}_l)$,
where we used the implicit convention that $u_l = v_l$.
This equation fulfilled for $\bar{S}_l $ is indeed the one that
needs to be fulfilled for the desired quantity $\hat{S}_l$, which
:math:`\bar{S}_l = f(l,\bar{S}_l) \approx S(v_l - \bar{S}_l) = S(u_l - \bar{S}_l)`,
where we used the implicit convention that :math:`u_l = v_l`.
This equation fulfilled for :math:`\bar{S}_l` is indeed the one that
needs to be fulfilled for the desired quantity :math:`\hat{S}_l`, which
can be shown as follows:
$u_l = u(v(u_l)) = S(v(u_l)) + v(u_l)$ and hence
$\hat{S}_l = u_l - v(u_l) = S(v(u_l)) = S(u_l - \hat{S}_l).
:math:`u_l = u(v(u_l)) = S(v(u_l)) + v(u_l)` and hence
:math:`\hat{S}_l = u_l - v(u_l) = S(v(u_l)) = S(u_l - \hat{S}_l)`.
This shows that the iteration converges to the correct solution if
it converges.
The initial value is $ \hat{S}_l^0 = S_l $.
The initial value is :math:`\hat{S}_l^0 = S_l`.
The iteration is repeated until
$ \max_l |\hat{S}_l^{n+1} - \hat{S}_l^{n} | < \epsilon $,
where $\epsilon$ is the tolerance specified when constructing
:math:`\max_l |\hat{S}_l^{n+1} - \hat{S}_l^{n} | < \epsilon`,
where :math:`\epsilon` is the tolerance specified when constructing
the operator.
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment