Cubic Spline Interpolator #719
Conversation
Pull Request Checklist
|
|
@MothNik thanks a lot for taking the time to turn the long stand Github issue into a PR ❤️ I am going to try to reply to all the points you have raised here and in your last comment in the Github issue, I will later do a proper review of the code and provide inline comments:
and seems to fail
probably nothing that will be ever perceived, but just curious 😉 1) not sure what you mean, if I run the code I gave you above with |
|
@mrava87 |
|
Perfect 😄 Do not worry about some of the CI failing, these are problems unrelated to your PR (as long as the doc and main Github action CI passes it is all good...) |
- stripped down `matvec` + `matmat` methods to solely rely on `matvec` - stripped down `reshape`-logic
…ssin.interpspline`
- moved utilities for interpolation to `_interp_utils.py` - fixed floating point rounding error of penultimate sample clipping behaviour for linear interpolation - transitioned `CubicSplineInterpolator` to `InterpCubicSpline` and enriched its docs (1 TODO, Notes missing)
- exposed `InterpCubicSpline` - made clipping logic for penultimate sample more natural by just going via the size of `iava`
- exposed `InterpCubicSpline`
- added `InterpCubicSpline`
- fixed missing naming update for `iava` clipping
- implemented matrix solver for when `dims[axis] == 2` - updated docs for minimum number of samples
- thoroughly tests `InterpCubicSpline` against `scipy.interpolate.CubicSpline`
- added `bc_type` argument to be forwards-compatible with the potential addition for further boundary conditions
- fixed typo
- fixed failing checks introduced with `bc_type`
- made casting to complex more explicit
- used `pytest.params` with `id` to make tests self-documenting without relying on comments - transitioned test specifications to a `dataclass` to avoid typos in key-based access with hard-coded strings - added typing
- removed duplicated code for non-unique `iava`-checks
- check `with pytest.raises` for actual error message - this uncovered a bug in the tests for `Bilinear` that only tested for a `ValueError` being raised in the `iava`-duplicate logic, but this was already raised during the attempt of element duplication in `iava` itself
|
@mrava87 The lines with pytest.raises(ValueError):
iava_rep = iava.copy()
iava_rep[-2] = [0, 0]
iava_rep[-1] = [0, 0]
_, _ = Bilinear(
iava_rep, dims=(par.y_num, par.x_num, par.t_num), dtype=par.dtype
)passed the assertion, but never hit the with pytest.raises(ValueError, match="repeated"):
iava_rep = iava.copy()
iava_rep[-2] = [0, 0]
iava_rep[-1] = [0, 0]
_, _ = Bilinear(
iava_rep, dims=(par.y_num, par.x_num, par.t_num), dtype=par.dtype
)revealed the following > with pytest.raises(ValueError, match="repeated"):
E AssertionError: Regex pattern did not match.
E Expected regex: 'repeated'
E Actual message: 'could not broadcast input array from shape (2,) into shape (10,)'The error happened already in Checking against common exceptions like Altneratively, and probably more readable, one can use something more verbose like with pytest.raises(ValueError) as exception_info:
iava_rep = iava.copy()
iava_rep[-2] = [0, 0]
iava_rep[-1] = [0, 0]
_, _ = Bilinear(
iava_rep, dims=(par.y_num, par.x_num, par.t_num), dtype=par.dtype
)
error_message = str(exception_info.value)
assert "iava" in error_message
assert "repeated" in error message |
- added `References` and `Notes` to `InterpCubicSpline`
- fixed typos and phrasing in documentation of `InterpCubicSpline`
Oh yeah, good catch! I agree, sometimes our with... tests have too generic assertments leading to unexpected behaviors like the one you spotted. I like the second one so we can check multiple words or even partial sentences in the error message 😉 |
| from pylops.utils.typing import NumericNDArray | ||
|
|
||
|
|
||
| def ensure_iava_is_unique( |
There was a problem hiding this comment.
Can you please add _ to both methods to re-enforce the fact they are private methods :)
| def ensure_iava_is_unique( | ||
| iava: NumericNDArray, | ||
| axis: int | None = None, | ||
| ) -> None: |
There was a problem hiding this comment.
Add docstring """Ensure that iava has only unique values"""... since it is a private method it is enough to just state what the method does.
| iava: NumericNDArray, | ||
| sample_size: int, | ||
| ) -> None: | ||
| # ensure that samples are not beyond the last sample, in that case set to |
There was a problem hiding this comment.
convert this comment into docstring
|
|
||
| assert np.allclose(y_eval_pylops, y_eval_scipy) | ||
|
|
||
| return |
| pytest.param(5.0, id="upsampling"), | ||
| ], | ||
| ) | ||
| def test_natural_cubic_spline_against_scipy( |
|
|
||
| Parameters | ||
| ---------- | ||
| dims : :obj:`int` |
There was a problem hiding this comment.
I think this should be the same as in Interp
dims : :obj:`list` or :obj:`int`
Number of samples for each dimension.
You can of course add A cubic spline requires ``dims[axis] > 2``. as it makes sense here.
| name : :obj:`str`, optional | ||
| Name of operator (to be used by :func:`pylops.utils.describe.describe`). | ||
|
|
||
| Returns |
There was a problem hiding this comment.
No Returns on a class.... instead you should use Attributes with public facing class attributes (basically we report all but those named self._*
Attributes
----------
dims : :obj:`tuple`
Shape of the array after the adjoint, but before flattening.
For example, ``x_reshaped = (Op.H * y.ravel()).reshape(Op.dims)``.
dimsd : :obj:`tuple`
Shape of the array after the forward, but before flattening.
For example, ``y_reshaped = (Op * x.ravel()).reshape(Op.dimsd)``.
shape : :obj:`tuple`
Operator shape.
lhs_matrix_lu **ADD type and description**
lhs_matrix_transposed_lu **ADD type and description**
|
|
||
| The adjoint operator :math:`\mathbf{S}^{H}` can be derived by rearranging the | ||
| involved operators. All of them are purely real and consequently, a transpose is | ||
| sufficient. This yields |
There was a problem hiding this comment.
This doesnt render well, I think you need an empty line here
| - :math:`p_{j,3}\left(t\right) = \frac{1}{6}\cdot\left(\left(t - \mathbf{k}_{i}\right)^{3} - \left(t - \mathbf{k}_{i}\right)\right)` | ||
|
|
||
| These base polynomials then need to be linearly combined using the coefficients | ||
| :math:`\mathbf{c} = \mathbf{F}\mathbf{y}`. Here, :math:`\mathbf{F}` is a |
There was a problem hiding this comment.
Is this y or x? Above you say the forward is Sx=PFx so I would assume c=Fx?
| .. [1] Wikipedia (German), Spline Interpolation | ||
| https://de.wikipedia.org/wiki/Spline-Interpolation#Kubisch_(Polynome_3._Grades) | ||
|
|
||
| Notes |
2 participants
This pull request aims to close #632 by adding the
kind="cubic_spline"-option topylops.signalprocessing.Interp.It is implemented in a way that can handle 1D-Arrays and 2D-Arrays along a single axis (0 or 1) in a vectorized fashion that heavily reuses pre-computations (optimised to match speed of
kind="linear"while providing similar accuracy askind="sinc"with less oscillations).Currently, only a
numpy-backend is available.