From 7e1c4c2affa9a3d3e81372f3aeec05fd6d3aba6c Mon Sep 17 00:00:00 2001 From: Pierre Raybaut <1311787+PierreRaybaut@users.noreply.github.com> Date: Mon, 26 Oct 2026 07:47:38 +0100 Subject: [PATCH] fix(signal): keep the product uncertainty defined for many factors Since #86, the first-order uncertainty of a product of signals is computed from the mantissas and exponents of the factors, but the products of the mantissas before and after each factor were computed with np.cumprod. Each mantissa lying in [0.5, 1), these products decreased geometrically with the number of factors and underflowed, even when the product itself was representable: the uncertainty was wrong, then zero, from about 150 float32 or 1074 float64 factors (e.g. 0 instead of 0.33 for 1100 factors equal to 1 with an uncertainty of 0.01 each). The running products of the mantissas are now renormalized at each step (np.frexp of the current product, exponents accumulated as integers), so that they stay in [0.5, 1) whatever the number of factors. Scaling by powers of two is exact: the results are bit-identical wherever no intermediate product underflowed. The result changes for many factors: the contract of sigima.signal.product, published at version 1 with Sigima 1.3.0, goes to version 2, qualified. Closes #105 Assisted-by: Claude Opus 5.5 --- .../LC_MESSAGES/release_notes/release_1.04.po | 8 +-- doc/release_notes/release_1.04.md | 4 +- sigima/proc/signal/arithmetic.py | 44 ++++++++++++----- sigima/tests/signal/operation_unit_test.py | 49 +++++++++++++++++++ 4 files changed, 88 insertions(+), 17 deletions(-) diff --git a/doc/locale/fr/LC_MESSAGES/release_notes/release_1.04.po b/doc/locale/fr/LC_MESSAGES/release_notes/release_1.04.po index c564dbe..0890a1d 100644 --- a/doc/locale/fr/LC_MESSAGES/release_notes/release_1.04.po +++ b/doc/locale/fr/LC_MESSAGES/release_notes/release_1.04.po @@ -76,8 +76,8 @@ msgstr "**Centre de l'image** : `ImageObj.xc` et `ImageObj.yc` sont désormais l msgid "**Ellipse coordinates**: `EllipseCoordinates` (`sigima.objects.shape`) now holds the 5 values `(x, y, a, b, theta)` of ellipse geometry results, which can thus be transformed." msgstr "**Coordonnées d'ellipse** : `EllipseCoordinates` (`sigima.objects.shape`) contient désormais les 5 valeurs `(x, y, a, b, theta)` des résultats géométriques d'ellipse, qui peuvent ainsi être transformés." -msgid "**Operation contract versions**: The operation contracts keep version 1 in this release, although the result of some operations changed since 1.3.0 (for example `resampling`, `extract_roi` and the product of signals): the contracts had just been published with Sigima 1.3.0 and were not used yet. This one-off exception does not apply to later changes: from now on, the version of a contract is increased whenever a change of the computation changes its result, whatever the releases of the applications using Sigima, and calls recorded with an earlier version are refused instead of being replayed (see [Operation contracts and versions](../user_guide/contracts.rst)). The image profiles, changed afterwards (see below), are thus at version 2." -msgstr "**Versions des contrats d'opération** : les contrats d'opération conservent la version 1 dans cette version, bien que le résultat de certaines opérations ait changé depuis la 1.3.0 (par exemple `resampling`, `extract_roi` et le produit de signaux) : les contrats venaient d'être publiés avec Sigima 1.3.0 et n'étaient pas encore exploités. Cette exception ponctuelle ne s'applique pas aux modifications ultérieures : désormais, la version d'un contrat est augmentée dès qu'une modification du calcul change son résultat, indépendamment des publications des applications utilisant Sigima, et les appels enregistrés avec une version antérieure sont refusés au lieu d'être rejoués (voir [Contrats d'opération et versions](../user_guide/contracts.rst)). Les profils d'image, modifiés ensuite (voir ci-dessous), sont donc en version 2." +msgid "**Operation contract versions**: The operation contracts keep version 1 in this release, although the result of some operations changed since 1.3.0 (for example `resampling`, `extract_roi` and the product of signals): the contracts had just been published with Sigima 1.3.0 and were not used yet. This one-off exception does not apply to later changes: from now on, the version of a contract is increased whenever a change of the computation changes its result, whatever the releases of the applications using Sigima, and calls recorded with an earlier version are refused instead of being replayed (see [Operation contracts and versions](../user_guide/contracts.rst)). The image profiles and the product of signals (`product`), changed afterwards (see below), are thus at version 2." +msgstr "**Versions des contrats d'opération** : les contrats d'opération conservent la version 1 dans cette version, bien que le résultat de certaines opérations ait changé depuis la 1.3.0 (par exemple `resampling`, `extract_roi` et le produit de signaux) : les contrats venaient d'être publiés avec Sigima 1.3.0 et n'étaient pas encore exploités. Cette exception ponctuelle ne s'applique pas aux modifications ultérieures : désormais, la version d'un contrat est augmentée dès qu'une modification du calcul change son résultat, indépendamment des publications des applications utilisant Sigima, et les appels enregistrés avec une version antérieure sont refusés au lieu d'être rejoués (voir [Contrats d'opération et versions](../user_guide/contracts.rst)). Les profils d'image et le produit de signaux (`product`), modifiés ensuite (voir ci-dessous), sont donc en version 2." msgid "**Image profiles in physical coordinates**: The line, segment, average and radial profiles of images are now defined and returned in the physical coordinates of the image, like regions of interest and analysis results, instead of pixel indices. Their parameters were renamed, so that scripts using the former names fail explicitly instead of silently using another meaning: `LineProfileParam.row`/`col` became `y`/`x` (the row or column whose pixel centers are the nearest), `SegmentProfileParam.row1`/`col1`/`row2`/`col2` and `AverageProfileParam.row1`/`col1`/`row2`/`col2` became `x1`/`y1`/`x2`/`y2` (for the average profile, `None` stands for the edge of the image instead of `-1`, and the pixels are those whose centers lie in the rectangle, like rectangular ROIs), and `RadialProfileParam.x0`/`y0` became `xc`/`yc`. The X axis of the profiles holds physical coordinates or distances, with the unit and label of the corresponding axis of the image. Segment profiles are sampled at regular steps from their start point to their end point included, with one sample per pixel along the axis where the segment spans the most pixels (as in Sigima 1.3), each sample taking the value of the nearest pixel. Radial profiles are binned in steps of the smallest distance between pixel centers, so that they are circular for anisotropic pixels and non-uniform coordinates. Segment and radial profiles are refused when the X and Y units of the image differ (and radial profiles when the centers of neighboring pixels coincide). Coordinates outside the image are refused instead of being clamped. On images without calibration (origin 0, pixel size 1), line, average and radial profiles are unchanged, and segment profiles take the same pixels for horizontal, vertical and diagonal segments: other oblique segments now take the nearest pixel of each sample instead of truncating its indices. The contracts of `line_profile`, `segment_profile`, `average_profile` and `radial_profile` are at version 2." msgstr "**Profils d'image en coordonnées physiques** : les profils ligne, segment, moyen et radial des images sont désormais définis et renvoyés dans les coordonnées physiques de l'image, comme les régions d'intérêt et les résultats d'analyse, au lieu d'indices de pixels. Leurs paramètres ont été renommés, afin que les scripts utilisant les anciens noms échouent explicitement au lieu d'utiliser silencieusement un autre sens : `LineProfileParam.row`/`col` deviennent `y`/`x` (la ligne ou la colonne dont les centres de pixels sont les plus proches), `SegmentProfileParam.row1`/`col1`/`row2`/`col2` et `AverageProfileParam.row1`/`col1`/`row2`/`col2` deviennent `x1`/`y1`/`x2`/`y2` (pour le profil moyen, `None` désigne le bord de l'image au lieu de `-1`, et les pixels retenus sont ceux dont le centre est dans le rectangle, comme pour les ROI rectangulaires), et `RadialProfileParam.x0`/`y0` deviennent `xc`/`yc`. L'axe X des profils contient des coordonnées ou des distances physiques, avec l'unité et le libellé de l'axe correspondant de l'image. Les profils de segment sont échantillonnés à pas réguliers de leur point de départ à leur point d'arrivée inclus, avec un échantillon par pixel le long de l'axe sur lequel le segment couvre le plus de pixels (comme dans Sigima 1.3), chaque échantillon prenant la valeur du pixel le plus proche. Les profils radiaux sont regroupés par pas égaux à la plus petite distance entre centres de pixels, de sorte qu'ils sont circulaires pour des pixels anisotropes et des coordonnées non uniformes. Les profils de segment et radiaux sont refusés lorsque les unités X et Y de l'image diffèrent (et les profils radiaux lorsque les centres de pixels voisins sont confondus). Les coordonnées hors de l'image sont refusées au lieu d'être ramenées au bord. Sur les images sans calibration (origine 0, taille de pixel 1), les profils ligne, moyen et radial sont inchangés, et les profils de segment prennent les mêmes pixels pour les segments horizontaux, verticaux et diagonaux : les autres segments obliques prennent désormais le pixel le plus proche de chaque échantillon au lieu de tronquer ses indices. Les contrats de `line_profile`, `segment_profile`, `average_profile` et `radial_profile` sont en version 2." @@ -112,8 +112,8 @@ msgstr "**Extraction de ROI d'images non uniformes** : l'extraction de plusieurs msgid "**Arbitrary rotation of non-uniform images**: The arbitrary rotation (`rotate`) now rejects images with non-uniform coordinates, like resizing and resampling: the result kept coordinates that did not match its pixels (nor its shape, when the rotation enlarged the image)." msgstr "**Rotation arbitraire d'images non uniformes** : la rotation arbitraire (`rotate`) refuse désormais les images à coordonnées non uniformes, comme le redimensionnement et le rééchantillonnage : le résultat conservait des coordonnées qui ne correspondaient pas à ses pixels (ni à sa forme, lorsque la rotation agrandissait l'image)." -msgid "**Uncertainty of a product of signals**: The uncertainty of the product of signals is now propagated where a factor is zero: the result had no uncertainty at all there. Complex factors now give a real uncertainty, computed from their moduli. The uncertainty no longer becomes undefined (NaN) when the factors span a very wide range of magnitudes (e.g. 10³⁰⁸ × 1 × 10⁻³⁰⁸) or when the product itself is too small to be represented: it is defined wherever its value can be represented, whatever the order of the signals." -msgstr "**Incertitude d'un produit de signaux** : l'incertitude du produit de signaux est désormais propagée là où un facteur est nul : le résultat n'y avait aucune incertitude. Les facteurs complexes donnent désormais une incertitude réelle, calculée à partir de leurs modules. L'incertitude ne devient plus indéfinie (NaN) lorsque les facteurs couvrent une très grande plage d'ordres de grandeur (par exemple 10³⁰⁸ × 1 × 10⁻³⁰⁸) ou lorsque le produit lui-même est trop petit pour être représenté : elle est définie partout où sa valeur est représentable, quel que soit l'ordre des signaux." +msgid "**Uncertainty of a product of signals**: The uncertainty of the product of signals is now propagated where a factor is zero: the result had no uncertainty at all there. Complex factors now give a real uncertainty, computed from their moduli. The uncertainty no longer becomes undefined (NaN) when the factors span a very wide range of magnitudes (e.g. 10³⁰⁸ × 1 × 10⁻³⁰⁸) or when the product itself is too small to be represented, and no longer vanishes when many signals are multiplied (it was wrong, then zero, from about 150 `float32` or 1 000 `float64` signals): it is defined wherever its value can be represented, whatever the number and the order of the signals. The contract of `product` is at version 2." +msgstr "**Incertitude d'un produit de signaux** : l'incertitude du produit de signaux est désormais propagée là où un facteur est nul : le résultat n'y avait aucune incertitude. Les facteurs complexes donnent désormais une incertitude réelle, calculée à partir de leurs modules. L'incertitude ne devient plus indéfinie (NaN) lorsque les facteurs couvrent une très grande plage d'ordres de grandeur (par exemple 10³⁰⁸ × 1 × 10⁻³⁰⁸) ou lorsque le produit lui-même est trop petit pour être représenté, et ne s'annule plus lorsque de nombreux signaux sont multipliés (elle était fausse, puis nulle, à partir d'environ 150 signaux `float32` ou 1 000 signaux `float64`) : elle est définie partout où sa valeur est représentable, quels que soient le nombre et l'ordre des signaux. Le contrat de `product` est en version 2." msgid "**ROIs and annotations after flips and rotations**: Horizontal and vertical flips and 90° and 270° rotations now move the ROIs and annotations of an image with its pixels: they were mirrored about a point half a pixel away from the center of the pixel grid, hence shifted by one pixel, and non-uniform coordinates were ignored. The move is exact for uniform coordinates. For non-uniform coordinates, reversing the pixels along a non-uniform axis is not an affine mapping: points, bounds and unrotated rectangles are moved exactly, while circles, ellipses and polygons are moved approximately (a circle does not stay a circle). ROIs defined in pixel indices were moved as if they were physical coordinates, and inverse ROIs (which mask the inside of their shape) became normal ROIs after these operations, the transposition and the translation." msgstr "**ROI et annotations après les retournements et les rotations** : les retournements horizontal et vertical et les rotations de 90° et 270° déplacent désormais les ROI et les annotations d'une image avec ses pixels : leur symétrie était centrée sur un point situé à un demi-pixel du centre de la grille de pixels, d'où un décalage d'un pixel, et les coordonnées non uniformes étaient ignorées. Le déplacement est exact pour des coordonnées uniformes. Pour des coordonnées non uniformes, l'inversion des pixels le long d'un axe non uniforme n'est pas une application affine : les points, les bornes et les rectangles non tournés sont déplacés exactement, tandis que les cercles, les ellipses et les polygones le sont de façon approchée (un cercle ne reste pas un cercle). Les ROI définies en indices de pixels étaient déplacées comme des coordonnées physiques, et les ROI inverses (qui masquent l'intérieur de leur forme) devenaient des ROI normales après ces opérations, la transposition et la translation." diff --git a/doc/release_notes/release_1.04.md b/doc/release_notes/release_1.04.md index d2be646..c3fed16 100644 --- a/doc/release_notes/release_1.04.md +++ b/doc/release_notes/release_1.04.md @@ -26,7 +26,7 @@ The [user guide](../user_guide/overview.rst) shows how readers and connectors at * **Replay of processing in image ROIs**: The selection of image pixels by regions of interest is now recorded with each call on images carrying regions of interest, and with each call of the operations taking image regions of interest as parameters (`extract_roi`, `extract_rois`, `erase` and `offset_correction`), whatever their input (rule `sigima.image.roi.pixel_centers`, `sigima.proc.contracts.IMAGE_ROI_CONVENTION`, in the execution context described by `describe_call` and `prepare_execution`). Applications that replay recorded processing (e.g. DataLab capsules) refuse to replay such calls recorded with earlier versions, which selected other pixels, instead of reporting a different result. Other calls on images without regions of interest are not concerned. * **Image center**: `ImageObj.xc` and `ImageObj.yc` are now the center of the pixel grid (halfway between the centers of the first and last pixels), for uniform and non-uniform coordinates: they were half a pixel off for uniform coordinates. `width` and `height` are unchanged. * **Ellipse coordinates**: `EllipseCoordinates` (`sigima.objects.shape`) now holds the 5 values `(x, y, a, b, theta)` of ellipse geometry results, which can thus be transformed. -* **Operation contract versions**: The operation contracts keep version 1 in this release, although the result of some operations changed since 1.3.0 (for example `resampling`, `extract_roi` and the product of signals): the contracts had just been published with Sigima 1.3.0 and were not used yet. This one-off exception does not apply to later changes: from now on, the version of a contract is increased whenever a change of the computation changes its result, whatever the releases of the applications using Sigima, and calls recorded with an earlier version are refused instead of being replayed (see [Operation contracts and versions](../user_guide/contracts.rst)). The image profiles, changed afterwards (see below), are thus at version 2. +* **Operation contract versions**: The operation contracts keep version 1 in this release, although the result of some operations changed since 1.3.0 (for example `resampling`, `extract_roi` and the product of signals): the contracts had just been published with Sigima 1.3.0 and were not used yet. This one-off exception does not apply to later changes: from now on, the version of a contract is increased whenever a change of the computation changes its result, whatever the releases of the applications using Sigima, and calls recorded with an earlier version are refused instead of being replayed (see [Operation contracts and versions](../user_guide/contracts.rst)). The image profiles and the product of signals (`product`), changed afterwards (see below), are thus at version 2. * **Image profiles in physical coordinates**: The line, segment, average and radial profiles of images are now defined and returned in the physical coordinates of the image, like regions of interest and analysis results, instead of pixel indices. Their parameters were renamed, so that scripts using the former names fail explicitly instead of silently using another meaning: `LineProfileParam.row`/`col` became `y`/`x` (the row or column whose pixel centers are the nearest), `SegmentProfileParam.row1`/`col1`/`row2`/`col2` and `AverageProfileParam.row1`/`col1`/`row2`/`col2` became `x1`/`y1`/`x2`/`y2` (for the average profile, `None` stands for the edge of the image instead of `-1`, and the pixels are those whose centers lie in the rectangle, like rectangular ROIs), and `RadialProfileParam.x0`/`y0` became `xc`/`yc`. The X axis of the profiles holds physical coordinates or distances, with the unit and label of the corresponding axis of the image. Segment profiles are sampled at regular steps from their start point to their end point included, with one sample per pixel along the axis where the segment spans the most pixels (as in Sigima 1.3), each sample taking the value of the nearest pixel. Radial profiles are binned in steps of the smallest distance between pixel centers, so that they are circular for anisotropic pixels and non-uniform coordinates. Segment and radial profiles are refused when the X and Y units of the image differ (and radial profiles when the centers of neighboring pixels coincide). Coordinates outside the image are refused instead of being clamped. On images without calibration (origin 0, pixel size 1), line, average and radial profiles are unchanged, and segment profiles take the same pixels for horizontal, vertical and diagonal segments: other oblique segments now take the nearest pixel of each sample instead of truncating its indices. The contracts of `line_profile`, `segment_profile`, `average_profile` and `radial_profile` are at version 2. ### 🛠️ Bug Fixes since version 1.3.0 ### @@ -40,7 +40,7 @@ The [user guide](../user_guide/overview.rst) shows how readers and connectors at * **Independent results**: The results of the transpose, 90° and 270° rotations and horizontal and vertical flips of an image no longer share their data (and their coordinates, for non-uniform images) with the source image: modifying a result in place modified its source too. * **ROI extraction of non-uniform images**: Extracting several regions of interest of an image with non-uniform coordinates into a single image now gives it the coordinates of the extracted pixels: the X coordinates were taken from the extracted rows and the Y coordinates from the extracted columns. * **Arbitrary rotation of non-uniform images**: The arbitrary rotation (`rotate`) now rejects images with non-uniform coordinates, like resizing and resampling: the result kept coordinates that did not match its pixels (nor its shape, when the rotation enlarged the image). -* **Uncertainty of a product of signals**: The uncertainty of the product of signals is now propagated where a factor is zero: the result had no uncertainty at all there. Complex factors now give a real uncertainty, computed from their moduli. The uncertainty no longer becomes undefined (NaN) when the factors span a very wide range of magnitudes (e.g. 10³⁰⁸ × 1 × 10⁻³⁰⁸) or when the product itself is too small to be represented: it is defined wherever its value can be represented, whatever the order of the signals. +* **Uncertainty of a product of signals**: The uncertainty of the product of signals is now propagated where a factor is zero: the result had no uncertainty at all there. Complex factors now give a real uncertainty, computed from their moduli. The uncertainty no longer becomes undefined (NaN) when the factors span a very wide range of magnitudes (e.g. 10³⁰⁸ × 1 × 10⁻³⁰⁸) or when the product itself is too small to be represented, and no longer vanishes when many signals are multiplied (it was wrong, then zero, from about 150 `float32` or 1 000 `float64` signals): it is defined wherever its value can be represented, whatever the number and the order of the signals. The contract of `product` is at version 2. * **ROIs and annotations after flips and rotations**: Horizontal and vertical flips and 90° and 270° rotations now move the ROIs and annotations of an image with its pixels: they were mirrored about a point half a pixel away from the center of the pixel grid, hence shifted by one pixel, and non-uniform coordinates were ignored. The move is exact for uniform coordinates. For non-uniform coordinates, reversing the pixels along a non-uniform axis is not an affine mapping: points, bounds and unrotated rectangles are moved exactly, while circles, ellipses and polygons are moved approximately (a circle does not stay a circle). ROIs defined in pixel indices were moved as if they were physical coordinates, and inverse ROIs (which mask the inside of their shape) became normal ROIs after these operations, the transposition and the translation. * **Coordinates of an extracted ROI**: `extract_roi` now gives the extracted image the coordinates of the extracted pixels. The origin of the source image was counted twice, and the coordinates of non-uniform images were shifted instead of being cut to the extracted pixels. * **Analysis results in ROIs**: The coordinates of the analysis results computed in a ROI (centroid, peaks, contours, etc.) are now offset by the position of the first extracted pixel: they could be off by up to half a pixel when the ROI bounds did not lie on pixel centers. diff --git a/sigima/proc/signal/arithmetic.py b/sigima/proc/signal/arithmetic.py index 420556e..ab452b7 100644 --- a/sigima/proc/signal/arithmetic.py +++ b/sigima/proc/signal/arithmetic.py @@ -66,13 +66,36 @@ def __signals_dy_to_array(signals: list[SignalObj]) -> np.ndarray: return np.array(arrays, dtype=np.result_type(*arrays)) +def __mantissa_products(mant: np.ndarray) -> tuple[np.ndarray, np.ndarray]: + """Return the products of the mantissas preceding each factor. + + Args: + mant: Mantissas of the factors (one row per factor), in [0.5, 1) or zero. + + Returns: + Mantissas and exponents of the products: row i is the product of the rows + before row i (1 for the first row). The running product is renormalized + at each step, so that it never underflows whatever the number of factors. + """ + prod_mant = np.empty_like(mant) + prod_exp = np.empty(mant.shape, dtype=np.int64) + current = np.ones_like(mant[0]) + exponent = np.zeros(mant.shape[1:], dtype=np.int64) + for i, row in enumerate(mant): + prod_mant[i], prod_exp[i] = current, exponent + current, shift = np.frexp(current * row) + exponent = exponent + shift + return prod_mant, prod_exp + + def __product_uncertainty(y_array: np.ndarray, dy_array: np.ndarray) -> np.ndarray: """Return the first-order uncertainty of the product of independent factors. ``σ² = Σ_i σ_i² Π_{j≠i} |y_j|²`` is computed without dividing by the factors, so that it is defined where a factor is zero. Each term is computed from the - mantissas and exponents of the factors (:func:`numpy.frexp`), so that no - intermediate product overflows or underflows whatever the order and the + mantissas and exponents of the factors (:func:`numpy.frexp`), the products of + the mantissas being renormalized at each step, so that no intermediate + product overflows or underflows whatever the number, the order and the magnitudes of the factors, and the terms are scaled by the largest one before being squared. @@ -88,14 +111,13 @@ def __product_uncertainty(y_array: np.ndarray, dy_array: np.ndarray) -> np.ndarr ftype = np.result_type(abs_y.dtype, abs_dy.dtype, np.float32) y_mant, y_exp = np.frexp(abs_y.astype(ftype, copy=False)) dy_mant, dy_exp = np.frexp(abs_dy.astype(ftype, copy=False)) - ones = np.ones_like(y_mant[:1]) - # Products of the mantissas before (prefix) and after (suffix) each factor: - # mantissas lie in [0.5, 1) (or are zero), so these products cannot overflow - prefix = np.cumprod(np.concatenate((ones, y_mant[:-1])), axis=0) - suffix = np.cumprod(np.concatenate((ones, y_mant[:0:-1])), axis=0)[::-1] + # Products of the mantissas before (prefix) and after (suffix) each factor + prefix_mant, prefix_exp = __mantissa_products(y_mant) + suffix_mant, suffix_exp = __mantissa_products(y_mant[::-1]) + suffix_mant, suffix_exp = suffix_mant[::-1], suffix_exp[::-1] # Term i is mant[i] * 2**exp[i] - mant = dy_mant * prefix * suffix - exp = dy_exp + (np.sum(y_exp, axis=0) - y_exp) + mant = dy_mant * prefix_mant * suffix_mant + exp = dy_exp + (np.sum(y_exp, axis=0) - y_exp) + prefix_exp + suffix_exp positive = mant > 0.0 exp_max = np.max(np.where(positive, exp, np.iinfo(exp.dtype).min), axis=0) exp_max = np.where(np.any(positive, axis=0), exp_max, 0) @@ -225,9 +247,9 @@ def standard_deviation(src_list: list[SignalObj]) -> SignalObj: @computation_function( operation_id="sigima.signal.product", - contract_version=1, + contract_version=2, x_alignment=SMALLEST_GRID_LINEAR_HOLD, - qualification=ReplayQualification(contract_version=1, comparison=EXACT), + qualification=ReplayQualification(contract_version=2, comparison=EXACT), ) def product(src_list: list[SignalObj]) -> SignalObj: """Compute the element-wise product of multiple signals. diff --git a/sigima/tests/signal/operation_unit_test.py b/sigima/tests/signal/operation_unit_test.py index 576da69..8845f4d 100644 --- a/sigima/tests/signal/operation_unit_test.py +++ b/sigima/tests/signal/operation_unit_test.py @@ -183,6 +183,55 @@ def test_signal_product_uncertainty(factors, expected_dy) -> None: assert np.allclose(np.real(result.dy), expected_dy, rtol=1e-12, atol=0.0) +def __factor_signals(values, rel: float, dtype=np.float64) -> list: + """Return signals of 3 points, each equal to a factor with a relative + uncertainty *rel*.""" + x = np.arange(3, dtype=dtype) + return [ + sigima.objects.create_signal( + f"s{i}", + x, + np.full(3, value, dtype=dtype), + dy=np.full(3, abs(value) * rel, dtype=dtype), + ) + for i, value in enumerate(values) + ] + + +@pytest.mark.parametrize( + "dtype, count, rtol", + [ + (np.float64, 300, 1e-12), + (np.float64, 1100, 1e-12), + (np.float64, 5000, 1e-12), + (np.float32, 200, 1e-6), + (np.float32, 1000, 1e-6), + ], +) +def test_signal_product_uncertainty_many_factors(dtype, count, rtol) -> None: + """The uncertainty of a product of many factors does not vanish (regression: + the products of the mantissas underflowed, so that the uncertainty was wrong, + then zero, from about 150 float32 or 1074 float64 factors).""" + result = sigima.proc.signal.product(__factor_signals([1.0] * count, 0.01, dtype)) + assert np.array_equal(result.y, np.ones(3)) + expected = 0.01 * np.sqrt(count) + assert np.allclose(result.dy, expected, rtol=rtol, atol=0.0) + + +def test_signal_product_uncertainty_many_unequal_factors() -> None: + """Many factors of various magnitudes: σ = |Π y| sqrt(Σ (σ_i / y_i)²).""" + # A representable product much smaller than each factor: 0.6 ** 1000 + result = sigima.proc.signal.product(__factor_signals([0.6] * 1000, 0.01)) + expected = 0.6**1000 * 0.01 * np.sqrt(1000) + assert np.allclose(result.dy, expected, rtol=1e-11, atol=0.0) + # Alternating large and small factors: the product is close to 1, while the + # product of their mantissas (about 0.5 ** 2000) underflowed + values = [1e3, -1e-3] * 2000 + result = sigima.proc.signal.product(__factor_signals(values, 0.01)) + expected = np.abs(np.prod(values)) * 0.01 * np.sqrt(len(values)) + assert np.allclose(result.dy, expected, rtol=1e-11, atol=0.0) + + def test_signal_product_uncertainty_overflow() -> None: """The uncertainty of a product is NaN where its first-order value overflows.""" x = np.arange(3.0)