Fix phasor_calibrate to handle multi harmonic calibration

[bruno-pannunzio]

Description

This PR fixes the issue of phasor_calibrate being unable to handle correctly multiple harmonic calibration. This is related to issue #68.

Not related, but also fixed a minor issue related to #67, to allow phasor_center to pass other keyword arguments other than axis to numpy functions used in methods.
...

Release note

Summarize the changes in the code block below to be included in the
release notes:

Fix `phasor_calibrate` to handle multiple-harmonic calibration.

Checklist

• The pull request title, summary, and description are concise.
• Related issues are linked in the description.
• New dependencies are explained.
• The source code and documentation can be distributed under the MIT license.
• The source code adheres to code standards.
• New classes, functions, and features are thoroughly tested.
• New, user-facing classes, functions, and features are documented.
• New features are covered in tutorials.
• No files other than source code, documentation, and project settings are added to the repository.

[bruno-pannunzio]

@cgohlke, I implemented some changes to phasor_transform and phasor_calibrate in order to properly handle calibration of multiple harmonics, as stated in issue #69. If you consider this is not necessary, we can close this PR, but I was unable to calibrate multiple harmonics of real data with previous implementations, at least with the pipeline I posted. I also added to the phasorpy_introduction tutorial an example for calibration of multiple harmonics. If you think it is not necessary I can remove it.

Finally, I modified the phasor_center function to allow passing other keyword arguments to numpy functions other than axis. I don't think this is necessary for the moment, but I think it may be useful later, although I can't think of an exact use right now. Also, if you believe this is not necessary, I can remove it.

[cgohlke]

I implemented some changes to phasor_transform and phasor_calibrate in order to properly handle calibration of multiple harmonics, as stated in issue #69. If you consider this is not necessary, we can close this PR

Yes, I think this case should be supported by the phasor_calibrate function if possible. My objection is to "fix" the phasor_transform function. Being a numpy ufunc, it is already capable to broadcast in the most general way. I suggest to move your changes out of phasor_transform into phasor_calibrate. That should work AFAICT. No?

I also added to the phasorpy_introduction tutorial an example for calibration of multiple harmonics. If you think it is not necessary I can remove it.

Right now it looks out of place. However, if you think the tutorial will make use of the second harmonic later, for example, to calculate multiple components, this section should replace the single harmonic calibration.

I modified the phasor_center function to allow passing other keyword arguments to numpy functions other than axis

Good idea. The one issue I can see is passing through the out parameter. To support out in the phasor_center function, it should be a tuple (out_real, out_imag).

[cgohlke]

I suggest to move this into the phasor_calibrate function.

axis=tuple(...) is calculated twice. Anyway, it looks to me that this would not work correctly if skip_axes are not the first axes?

There is some code in phasor_center that looks related. Maybe the checking/normalization of skip_axes can be refactored into a helper function?

phasorpy/src/phasorpy/phasor.py

Lines 1988 to 1996 in 8e9b123

	if skip_axes is None:
	axis = None
	else:
	if not isinstance(skip_axes, Sequence):
	skip_axes = (skip_axes,)
	if any(i >= real.ndim for i in skip_axes):
	raise IndexError(f'{skip_axes=} out of range {real.ndim=}')
	skip_axes = tuple(i % real.ndim for i in skip_axes)
	axis = tuple(i for i in range(real.ndim) if i not in skip_axes)

[bruno-pannunzio]

Thanks @cgohlke ! Your suggestions are excelent, I haven't realized I calculated the axis twice and your way is much more cleaner an efficent.

[bruno-pannunzio]

Also, I tried and it works for different axis, not just the first one.

[bruno-pannunzio]

I totally agree with you @cgohlke! I don't know why i didn't think first to change the phasor_calibrate instead of phasor_transform, it makes more sense. I moved it now with your suggestions also and it works great, also with different axes.

I removed the example from phasorpy_introduction tutorial as suggested. I think it will need to be somewhere, but as you pointed out, maybe in the component analysis, since multiple harmonic calibration will be needed there. Let's see when it is done if it makes sense to add it there.

Regarding the out argument passed to numpy functions, I don't know how we can approach that. Maybe a check or a warning or something in the description of the function? Maybe it is too much, specially if we add other methods and we need to consider more exceptions.

[bruno-pannunzio]

I will refactor the axis determination code so that it can be used for both phasor_calibrate and phasor_center. Should this helper function be in the _utils module or in the phasor module as a hidden function?

[cgohlke]

I would prefer the helper function to also validate and "normalize" the skip_axes parameter.

Should we rename skip_axes to skip_axis to match the numpy axis argument?

The helper function should be tested.

How about something like:

def _parse_skip_axis(
    skip_axis: int | Sequence[int] | None,
    /,
    ndim: int,
) -> tuple[tuple[int, ...], tuple[int, ...]]:
    """Return axes to skip and not to skip.

    This helper function is used to validate and parse `skip_axis`
    parameters.

    Parameters
    ----------
    skip_axis : Sequence of int, or None
        Axes to skip. If None, no axes are skipped.
    ndim : int
        Dimensionality of array in which to skip axes.

    Returns
    -------
    skip_axis
        Ordered, positive values of `skip_axis`.
    other_axis
        Axes indices not included in `skip_axis`.

    Raises
    ------
    IndexError
        If any `skip_axis` value is out of bounds of `ndim`.

    Examples
    --------
    >>> _parse_skip_axis((1, -2), 5)
    (1, 3), (0, 2, 4)

    """
    if ndim < 0:
        raise ValueError(f'invalid {ndim=}')
    if skip_axis is None:
        return (), tuple(range(ndim))
    if not isinstance(skip_axis, Sequence):
        skip_axis = (skip_axis,)
    if any(i >= ndim or i < -ndim for i in skip_axis):
        raise IndexError(f"skip_axis={skip_axis} out of range for {ndim=}")
    skip_axis = tuple(sorted(int(i % ndim) for i in skip_axis))
    other_axis = tuple(i for i in range(ndim) if i not in skip_axis)
    return skip_axis, other_axis

and

def test_parse_skip_axis():
    """Test _parse_skip_axis function."""
    assert _parse_skip_axis(None, 0) == ((), ())
    assert _parse_skip_axis(None, 1) == ((), (0,))
    assert _parse_skip_axis((), 1) == ((), (0,))
    assert _parse_skip_axis(0, 1) == ((0,), ())
    assert _parse_skip_axis(0, 2) == ((0,), (1,))
    assert _parse_skip_axis(-1, 2) == ((1,), (0,))
    assert _parse_skip_axis((1, -2), 5) == ((1, 3), (0, 2, 4))
    with pytest.raises(ValueError):
        _parse_skip_axis(0, -1)
    with pytest.raises(IndexError):
        _parse_skip_axis(0, 0)
    with pytest.raises(IndexError):
        _parse_skip_axis(1, 1)
    with pytest.raises(IndexError):
        _parse_skip_axis(-2, 1)

[bruno-pannunzio]

Thanks @cgohlke, that is a nice improvement and you are right to add the test! I also think it makes sense to change to skip_axis to keep with nomenclature matching numpy, so I modified it along all functions.

I think now it is ready, but please let me know if there is anything worth improving.