diff --git a/docs/docs/api-reference/index.md b/docs/docs/api-reference/index.md index 4bd6f561..d530e187 100644 --- a/docs/docs/api-reference/index.md +++ b/docs/docs/api-reference/index.md @@ -18,6 +18,9 @@ available in EasyScience. persisting EasyScience objects (serializers and components). - [job](job.md) – Job and experiment abstractions for running and organizing analyses. +- [math](math.md) – Maths functions on `Parameter` objects (for example + `es.math.sin`, `es.math.log`, `es.math.arctan2`) that propagate units, + uncertainties and fitting bounds. - [models](models.md) – Predefined model implementations (e.g. polynomial models) used in fitting workflows. - [utils](utils.md) – Miscellaneous utility functions and helpers (class diff --git a/docs/docs/api-reference/math.md b/docs/docs/api-reference/math.md new file mode 100644 index 00000000..8fff14c2 --- /dev/null +++ b/docs/docs/api-reference/math.md @@ -0,0 +1 @@ +::: easyscience.math diff --git a/docs/mkdocs.yml b/docs/mkdocs.yml index 91165ebc..47bc6fef 100644 --- a/docs/mkdocs.yml +++ b/docs/mkdocs.yml @@ -206,6 +206,7 @@ nav: - global_object: api-reference/global_object.md - io: api-reference/io.md - job: api-reference/job.md + - math: api-reference/math.md - models: api-reference/models.md - utils: api-reference/utils.md - variable: api-reference/variable.md diff --git a/pixi.lock b/pixi.lock index eb033db0..0c633805 100644 --- a/pixi.lock +++ b/pixi.lock @@ -1,13 +1,13 @@ version: 7 platforms: -- name: p1 +- name: linux-64-glibc-2-35 subdir: linux-64 virtual-packages: - __glibc=2.35 - __unix=0=0 - __linux=4.18 - __archspec=0=x86_64 -- name: p2 +- name: osx-arm64-macos-14-0 subdir: osx-arm64 virtual-packages: - __osx=14.0 @@ -25,7 +25,7 @@ environments: indexes: - https://pypi.org/simple packages: - p1: + linux-64-glibc-2-35: - conda: https://conda.anaconda.org/conda-forge/linux-64/_openmp_mutex-4.5-20_gnu.conda - conda: https://conda.anaconda.org/conda-forge/linux-64/argon2-cffi-bindings-25.1.0-py314h5bd0f2a_2.conda - conda: https://conda.anaconda.org/conda-forge/linux-64/bzip2-1.0.8-hda65f42_9.conda @@ -298,7 +298,7 @@ environments: - pypi: https://files.pythonhosted.org/packages/fd/7b/122376b1fd3c62c1ed9dc80c931ace4844b3c55407b6fb2d199377c9736f/pydantic-2.13.4-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/fe/3b/8ec5074bcfc450fe84273713b4b0a0dd47c0249358f5d82eb8104ffe2520/multidict-6.7.1-cp314-cp314-manylinux2014_x86_64.manylinux_2_17_x86_64.manylinux_2_28_x86_64.whl - pypi: https://files.pythonhosted.org/packages/fe/d7/a71b12b625c44cfcb95d41397752edccc0a731903081609cffd6720f671d/asteval-1.0.9-py3-none-any.whl - p2: + osx-arm64-macos-14-0: - conda: https://conda.anaconda.org/conda-forge/noarch/_python_abi3_support-1.0-hd8ed1ab_3.conda - conda: https://conda.anaconda.org/conda-forge/noarch/anyio-4.14.2-pyhcf101f3_0.conda - conda: https://conda.anaconda.org/conda-forge/noarch/appnope-0.1.4-pyhd8ed1ab_1.conda @@ -832,7 +832,7 @@ environments: indexes: - https://pypi.org/simple packages: - p1: + linux-64-glibc-2-35: - conda: https://conda.anaconda.org/conda-forge/linux-64/_openmp_mutex-4.5-20_gnu.conda - conda: https://conda.anaconda.org/conda-forge/linux-64/argon2-cffi-bindings-25.1.0-py311h49ec1c0_2.conda - conda: https://conda.anaconda.org/conda-forge/linux-64/bzip2-1.0.8-hda65f42_9.conda @@ -1103,7 +1103,7 @@ environments: - pypi: https://files.pythonhosted.org/packages/fb/e2/79c688af8b210d232694e31e59da9f6ec747bae31c3f5946e4e9b98860d5/click-8.4.2-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/fd/7b/122376b1fd3c62c1ed9dc80c931ace4844b3c55407b6fb2d199377c9736f/pydantic-2.13.4-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/fe/d7/a71b12b625c44cfcb95d41397752edccc0a731903081609cffd6720f671d/asteval-1.0.9-py3-none-any.whl - p2: + osx-arm64-macos-14-0: - conda: https://conda.anaconda.org/conda-forge/noarch/anyio-4.14.2-pyhcf101f3_0.conda - conda: https://conda.anaconda.org/conda-forge/noarch/appnope-0.1.4-pyhd8ed1ab_1.conda - conda: https://conda.anaconda.org/conda-forge/noarch/argon2-cffi-25.1.0-pyhd8ed1ab_0.conda @@ -1628,7 +1628,7 @@ environments: indexes: - https://pypi.org/simple packages: - p1: + linux-64-glibc-2-35: - conda: https://conda.anaconda.org/conda-forge/linux-64/_openmp_mutex-4.5-20_gnu.conda - conda: https://conda.anaconda.org/conda-forge/linux-64/argon2-cffi-bindings-25.1.0-py314h5bd0f2a_2.conda - conda: https://conda.anaconda.org/conda-forge/linux-64/bzip2-1.0.8-hda65f42_9.conda @@ -1901,7 +1901,7 @@ environments: - pypi: https://files.pythonhosted.org/packages/fd/7b/122376b1fd3c62c1ed9dc80c931ace4844b3c55407b6fb2d199377c9736f/pydantic-2.13.4-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/fe/3b/8ec5074bcfc450fe84273713b4b0a0dd47c0249358f5d82eb8104ffe2520/multidict-6.7.1-cp314-cp314-manylinux2014_x86_64.manylinux_2_17_x86_64.manylinux_2_28_x86_64.whl - pypi: https://files.pythonhosted.org/packages/fe/d7/a71b12b625c44cfcb95d41397752edccc0a731903081609cffd6720f671d/asteval-1.0.9-py3-none-any.whl - p2: + osx-arm64-macos-14-0: - conda: https://conda.anaconda.org/conda-forge/noarch/_python_abi3_support-1.0-hd8ed1ab_3.conda - conda: https://conda.anaconda.org/conda-forge/noarch/anyio-4.14.2-pyhcf101f3_0.conda - conda: https://conda.anaconda.org/conda-forge/noarch/appnope-0.1.4-pyhd8ed1ab_1.conda @@ -2435,7 +2435,7 @@ environments: indexes: - https://pypi.org/simple packages: - p1: + linux-64-glibc-2-35: - conda: https://conda.anaconda.org/conda-forge/linux-64/_openmp_mutex-4.5-20_gnu.conda - conda: https://conda.anaconda.org/conda-forge/linux-64/argon2-cffi-bindings-25.1.0-py314h5bd0f2a_2.conda - conda: https://conda.anaconda.org/conda-forge/linux-64/bzip2-1.0.8-hda65f42_9.conda @@ -2598,7 +2598,7 @@ environments: - pypi: https://files.pythonhosted.org/packages/fb/76/641ae371508676492379f16e2fa48f4e2c11741bd63c48be4b12a6b09cba/aiosignal-1.4.0-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/fe/3b/8ec5074bcfc450fe84273713b4b0a0dd47c0249358f5d82eb8104ffe2520/multidict-6.7.1-cp314-cp314-manylinux2014_x86_64.manylinux_2_17_x86_64.manylinux_2_28_x86_64.whl - pypi: https://files.pythonhosted.org/packages/fe/d7/a71b12b625c44cfcb95d41397752edccc0a731903081609cffd6720f671d/asteval-1.0.9-py3-none-any.whl - p2: + osx-arm64-macos-14-0: - conda: https://conda.anaconda.org/conda-forge/noarch/_python_abi3_support-1.0-hd8ed1ab_3.conda - conda: https://conda.anaconda.org/conda-forge/noarch/anyio-4.14.2-pyhcf101f3_0.conda - conda: https://conda.anaconda.org/conda-forge/noarch/appnope-0.1.4-pyhd8ed1ab_1.conda diff --git a/src/easyscience/__init__.py b/src/easyscience/__init__.py index ae709682..b5adb47d 100644 --- a/src/easyscience/__init__.py +++ b/src/easyscience/__init__.py @@ -11,6 +11,7 @@ global_object.stack.enabled = False +from . import math # noqa: E402 from .base_classes import ObjBase # noqa: E402 from .fitting import AvailableMinimizers # noqa: E402 from .fitting import Fitter # noqa: E402 @@ -21,12 +22,12 @@ __version__ = version('easyscience') __all__ = [ - __version__, - global_object, - ObjBase, - CollectionBase, - AvailableMinimizers, - Fitter, - DescriptorNumber, - Parameter, + '__version__', + 'global_object', + 'ObjBase', + 'CollectionBase', + 'AvailableMinimizers', + 'Fitter', + 'DescriptorNumber', + 'Parameter', ] diff --git a/src/easyscience/fitting/minimizers/minimizer_bumps.py b/src/easyscience/fitting/minimizers/minimizer_bumps.py index 315d856c..500edce6 100644 --- a/src/easyscience/fitting/minimizers/minimizer_bumps.py +++ b/src/easyscience/fitting/minimizers/minimizer_bumps.py @@ -417,42 +417,45 @@ def mcmc_sample( weights : np.ndarray Flattened weight array. samples : int, default=10000 - Number of raw samples to draw across all chains, before thinning. - A guaranteed minimum, not an exact count: DREAM advances in - blocks of 10 generations (one generation = one draw per chain) - and stops at the first block boundary at or past ``samples``. + Number of raw samples to draw across all chains, before + thinning. A guaranteed minimum, not an exact count: DREAM + advances in blocks of 10 generations (one generation = one + draw per chain) and stops at the first block boundary at or + past ``samples``. burn : int, default=2000 Burn-in generations to discard. BUMPS counts ``burn`` in - generations while ``samples`` counts raw draws, so ``burn=500`` - discards ``500 * n_chains`` raw samples. + generations while ``samples`` counts raw draws, so + ``burn=500`` discards ``500 * n_chains`` raw samples. thin : int, default=10 - Thinning interval — only every ``thin``-th generation is stored. + Thinning interval — only every ``thin``-th generation is + stored. population : int | None, default=None - BUMPS DREAM population count per parameter (number of parallel - chains): BUMPS creates ``ceil(population * n_parameters)`` chains. + BUMPS DREAM population count per parameter (number of + parallel chains): BUMPS creates ``ceil(population * + n_parameters)`` chains. resume_state : MCMCDraw | None, default=None A BUMPS ``MCMCDraw`` state object from a previous - ``mcmc_sample()`` call (e.g. ``PosteriorResults.sampler_state``). - When provided, DREAM **continues** the saved chain instead of - starting cold. The population, parameter count, and parameter - names must match the current model — a ``ValueError`` is raised - otherwise. - - ``samples`` must be the **total** number of raw samples, not an - increment: to extend an existing chain of ``N`` raw samples by - ``M``, pass ``samples=N + M`` (DREAM keeps only the last - ``samples`` draws in its buffer). The `Sampler.extend` helper - computes this for you. - - ``burn`` is forced to 0 on resume: a previously-converged chain is - never re-burned. - - The ``population`` and ``initializer`` parameters - have **no effect** when ``resume_state`` is provided — they - are determined by the saved state. - - Resuming against *different* data is undefined behaviour (the - chain's likelihood changes underneath it). + ``mcmc_sample()`` call (e.g. + ``PosteriorResults.sampler_state``). When provided, DREAM + **continues** the saved chain instead of starting cold. The + population, parameter count, and parameter names must match + the current model — a ``ValueError`` is raised otherwise. + + ``samples`` must be the **total** number of raw samples, not + an increment: to extend an existing chain of ``N`` raw + samples by ``M``, pass ``samples=N + M`` (DREAM keeps only + the last ``samples`` draws in its buffer). The + ``Sampler.extend`` helper computes this for you. + + ``burn`` is forced to 0 on resume: a previously-converged + chain is never re-burned. + + The ``population`` and ``initializer`` parameters have **no + effect** when ``resume_state`` is provided — they are + determined by the saved state. + + Resuming against *different* data is undefined behaviour + (the chain's likelihood changes underneath it). sampler_kwargs : dict | None, default=None Additional keyword arguments forwarded to ``bumps.fitters.fit``. @@ -475,9 +478,9 @@ def mcmc_sample( ------ ValueError If the input shapes or weights are invalid, if - ``progress_callback`` is not callable, or if ``resume_state`` - is incompatible with the current model (parameter count, - names/order, or population mismatch). + ``progress_callback`` is not callable, or if + ``resume_state`` is incompatible with the current model + (parameter count, names/order, or population mismatch). FitError If DREAM sampling was aborted by the user (via ``abort_test``). @@ -608,28 +611,32 @@ def _validate_resume_state( population: int | None, burn: int, ) -> tuple[int, int]: - """Check that ``resume_state`` is compatible with ``problem`` and + """ + Check that ``resume_state`` is compatible with ``problem`` and resolve the population and burn values to use when resuming. Parameters ---------- problem : FitProblem - The freshly built BUMPS ``FitProblem`` for the current model. + The freshly built BUMPS ``FitProblem`` for the current + model. resume_state : MCMCDraw The saved chain state to resume from. population : int | None The caller-supplied population scale factor, or ``None``. burn : int - The caller-supplied burn-in, ignored (with a warning) on resume. + The caller-supplied burn-in, ignored (with a warning) on + resume. Returns ------- tuple[int, int] ``(population, burn)`` to pass to DREAM. The population is returned as a **negative** number, which BUMPS' - ``initpop.generate`` reads as an absolute chain count, exactly - reproducing the saved state's population. ``burn`` is always 0: - a previously converged chain is never re-burned. + ``initpop.generate`` reads as an absolute chain count, + exactly reproducing the saved state's population. ``burn`` + is always 0: a previously converged chain is never + re-burned. Raises ------ diff --git a/src/easyscience/fitting/sampler.py b/src/easyscience/fitting/sampler.py index d3915115..3917a5bc 100644 --- a/src/easyscience/fitting/sampler.py +++ b/src/easyscience/fitting/sampler.py @@ -1,6 +1,8 @@ # SPDX-FileCopyrightText: 2026 EasyScience contributors # SPDX-License-Identifier: BSD-3-Clause -"""Bayesian MCMC sampling — the ``Sampler`` class and persistence helpers.""" +""" +Bayesian MCMC sampling — the ``Sampler`` class and persistence helpers. +""" from __future__ import annotations @@ -38,7 +40,9 @@ def _data_fingerprint( y_list: list, w_list: list, ) -> str | None: - """Return a SHA-256 hex digest of concatenated (x|y|weights), or None.""" + """ + Return a SHA-256 hex digest of concatenated (x|y|weights), or None. + """ try: h = hashlib.sha256() for arr in list(x_list) + list(y_list) + list(w_list): @@ -53,7 +57,8 @@ def _validate_dataset_arrays( data: np.ndarray | list | tuple, allow_none_entries: bool = False, ) -> None: - """Check that ``data`` (an array or list of arrays) holds numeric, + """ + Check that ``data`` (an array or list of arrays) holds numeric, at-least-1-D, non-empty arrays. Structural checks (array vs list, matching dataset counts) are in @@ -101,12 +106,13 @@ def _validate_dataset_arrays( def _copy_data(data): - """Copy an array (or list of arrays) and mark the copies read-only. + """ + Copy an array (or list of arrays) and mark the copies read-only. - The ``Sampler`` binds its data at construction; copying decouples the - bound data from the caller's arrays, and the read-only flag stops - in-place mutation of the copies, so the chain and the ``save()`` - fingerprint always describe the data actually sampled. + The ``Sampler`` binds its data at construction; copying decouples + the bound data from the caller's arrays, and the read-only flag + stops in-place mutation of the copies, so the chain and the + ``save()`` fingerprint always describe the data actually sampled. """ if data is None: return None @@ -118,10 +124,13 @@ def _copy_data(data): def _validate_chain_path(path: str | os.PathLike, skip: int = 0) -> str: - """Validate the persistence arguments shared by the save/load functions. + """ + Validate the persistence arguments shared by the save/load + functions. - Returns ``path`` coerced to ``str``. Raises ``TypeError`` if ``path`` is - not path-like and ``ValueError`` if ``skip`` is not a non-negative int. + Returns ``path`` coerced to ``str``. Raises ``TypeError`` if + ``path`` is not path-like and ``ValueError`` if ``skip`` is not a + non-negative int. """ if not isinstance(path, (str, os.PathLike)): raise TypeError(f'path must be a str or os.PathLike, got {type(path).__name__}.') @@ -131,14 +140,15 @@ def _validate_chain_path(path: str | os.PathLike, skip: int = 0) -> str: def _load_bumps_state(path: str, skip: int = 0) -> MCMCDraw: - """Read a BUMPS chain from disk, working around a BUMPS bug. + """ + Read a BUMPS chain from disk, working around a BUMPS bug. ``bumps.dream.state.load_state`` reads the saved buffers with - ``numpy.loadtxt``, which collapses a single-row file to a 1-D array. A - short chain stores a single CR-weight update row, so ``load_state``'s - subsequent ``stats[:, 0]`` indexing raises ``IndexError: too many indices - for array``. We coerce each buffer read back to 2-D before ``load_state`` - consumes it. + ``numpy.loadtxt``, which collapses a single-row file to a 1-D array. + A short chain stores a single CR-weight update row, so + ``load_state``'s subsequent ``stats[:, 0]`` indexing raises + ``IndexError: too many indices for array``. We coerce each buffer + read back to 2-D before ``load_state`` consumes it. Parameters ---------- @@ -170,28 +180,30 @@ def _loadtxt_2d(file, report: int = 0) -> np.ndarray: def load_chain(path: str | os.PathLike, skip: int = 0) -> tuple[MCMCDraw, list[str] | None, dict]: - """Reload a DREAM chain state saved by ``Sampler.save``. + """ + Reload a DREAM chain state saved by ``Sampler.save``. - This is the standalone reader: unlike ``Sampler.load_state`` it needs no - fitter, model or data, so a saved chain can be inspected or post-processed - on a machine that does not have the model. Parameter names are restored - from the sidecar when available (schema versions 1 and 2), falling back to - the state's labels with the minimizer prefix stripped. + This is the standalone reader: unlike ``Sampler.load_state`` it + needs no fitter, model or data, so a saved chain can be inspected or + post-processed on a machine that does not have the model. Parameter + names are restored from the sidecar when available (schema versions + 1 and 2), falling back to the state's labels with the minimizer + prefix stripped. Parameters ---------- path : str | os.PathLike File path prefix used when saving. skip : int, default=0 - Discard the first ``skip`` saved generations on load, forwarded to - ``bumps.dream.state.load_state``. + Discard the first ``skip`` saved generations on load, forwarded + to ``bumps.dream.state.load_state``. Returns ------- tuple[MCMCDraw, list[str] | None, dict] - The reloaded BUMPS chain state, the parameter names (or ``None`` if - neither sidecar nor labels yielded them), and the raw sidecar dict - (empty if absent/unreadable). + The reloaded BUMPS chain state, the parameter names (or ``None`` + if neither sidecar nor labels yielded them), and the raw sidecar + dict (empty if absent/unreadable). Raises ------ @@ -229,15 +241,18 @@ def load_chain(path: str | os.PathLike, skip: int = 0) -> tuple[MCMCDraw, list[s @dataclass class SamplingResults: - """Structured result of an MCMC sampling run (analogous to ``FitResults``). + """ + Structured result of an MCMC sampling run (analogous to + ``FitResults``). Attributes ---------- draws : np.ndarray - Posterior samples, shape ``(n_samples, n_params)``. For results from - ``Sampler.sample()``/``Sampler.extend()`` this is a *trimmed, - outlier-filtered* view of the chain, not the raw buffer, so it holds - fewer rows than ``samples / thin`` — see the notes on ``Sampler``. + Posterior samples, shape ``(n_samples, n_params)``. For results + from ``Sampler.sample()``/``Sampler.extend()`` this is a + *trimmed, outlier-filtered* view of the chain, not the raw + buffer, so it holds fewer rows than ``samples / thin`` — see the + notes on ``Sampler``. param_names : list[str] Parameter names (one per column of ``draws``). logp : np.ndarray @@ -252,8 +267,10 @@ class SamplingResults: state: MCMCDraw def to_legacy_dict(self) -> dict: - """Return the legacy dict shape produced by the deprecated - ``mcmc_sample()`` APIs.""" + """ + Return the legacy dict shape produced by the deprecated + ``mcmc_sample()`` APIs. + """ return { 'draws': self.draws, 'param_names': self.param_names, @@ -263,58 +280,64 @@ def to_legacy_dict(self) -> dict: class Sampler: - """Bayesian MCMC sampler for one dataset, backed by a Fitter's BUMPS minimizer. - - One ``Sampler`` instance represents one chain over one ``(x, y, weights)`` - dataset. The data is bound at construction; ``sample()`` and ``extend()`` - take no data arguments, so a chain can never be extended against different - data (undefined behaviour in BUMPS). The bound data is a defensive, - read-only copy of the caller's arrays, exposed via the ``x``, ``y`` and - ``weights`` properties — mutating the originals after construction has no - effect on the sampler, and there are deliberately no setters: to sample - different data, create a new ``Sampler``. - - Construct directly with a configured ``Fitter`` (or ``MultiFitter``) whose - minimizer has been switched to ``AvailableMinimizers.Bumps``. **Running a - fit first is not required** — the ``Fitter`` supplies the model and the - minimizer, not a fit result, and sampling from the initial parameter values - works fine. - - It is often worth fitting first anyway. DREAM seeds its whole starting - population inside a tiny ball around the parameters' *current* values - (BUMPS' default ``init='eps'``), so sampling from fitted values starts the - chain in the right region and shortens the burn-in needed to reach the - typical set. From a poor initial guess, expect to burn for longer. - - The sampler is BUMPS/DREAM-specific for now: the BUMPS check in ``_run()`` - is the seam where another backend would plug in. + """ + Bayesian MCMC sampler for one dataset, backed by a Fitter's BUMPS + minimizer. + + One ``Sampler`` instance represents one chain over one ``(x, y, + weights)`` dataset. The data is bound at construction; ``sample()`` + and ``extend()`` take no data arguments, so a chain can never be + extended against different data (undefined behaviour in BUMPS). The + bound data is a defensive, read-only copy of the caller's arrays, + exposed via the ``x``, ``y`` and ``weights`` properties — mutating + the originals after construction has no effect on the sampler, and + there are deliberately no setters: to sample different data, create + a new ``Sampler``. + + Construct directly with a configured ``Fitter`` (or ``MultiFitter``) + whose minimizer has been switched to ``AvailableMinimizers.Bumps``. + **Running a fit first is not required** — the ``Fitter`` supplies + the model and the minimizer, not a fit result, and sampling from the + initial parameter values works fine. + + It is often worth fitting first anyway. DREAM seeds its whole + starting population inside a tiny ball around the parameters' + *current* values (BUMPS' default ``init='eps'``), so sampling from + fitted values starts the chain in the right region and shortens the + burn-in needed to reach the typical set. From a poor initial guess, + expect to burn for longer. + + The sampler is BUMPS/DREAM-specific for now: the BUMPS check in + ``_run()`` is the seam where another backend would plug in. Parameters ---------- - fitter : Fitter - A configured ``Fitter`` (or ``MultiFitter``) whose minimizer has been - switched to ``AvailableMinimizers.Bumps``. + fitter : 'Fitter' + A configured ``Fitter`` (or ``MultiFitter``) whose minimizer has + been switched to ``AvailableMinimizers.Bumps``. x : np.ndarray | list[np.ndarray] - Independent variable array (or list of arrays for ``MultiFitter``). + Independent variable array (or list of arrays for + ``MultiFitter``). y : np.ndarray | list[np.ndarray] - Dependent variable array (or list of arrays for ``MultiFitter``). + Dependent variable array (or list of arrays for + ``MultiFitter``). weights : np.ndarray | list[np.ndarray | None] | None, default=None Weight array (or list of arrays for ``MultiFitter``). vectorized : bool, default=False When ``True``, each x array may be multi-dimensional (e.g. an ``(N, M, 2)`` grid for a 2D model) and is left as-is. sampler_kwargs : dict | None, default=None - Per-instance default keyword arguments forwarded to the BUMPS DREAM - sampler on every run, merged with (and overridden by) per-call - ``sampler_kwargs``. + Per-instance default keyword arguments forwarded to the BUMPS + DREAM sampler on every run, merged with (and overridden by) + per-call ``sampler_kwargs``. Raises ------ TypeError - If ``fitter`` is not Fitter-shaped (no ``minimizer``/``fit_function``), - if any dataset in ``x``/``y``/``weights`` is not a numeric array - (e.g. a string), or ``vectorized``/``sampler_kwargs`` have the wrong - type. + If ``fitter`` is not Fitter-shaped (no + ``minimizer``/``fit_function``), if any dataset in + ``x``/``y``/``weights`` is not a numeric array (e.g. a string), + or ``vectorized``/``sampler_kwargs`` have the wrong type. ValueError If ``x``, ``y`` and ``weights`` do not hold matching structures (all arrays, or lists of the same length), or any dataset is a @@ -322,42 +345,44 @@ class Sampler: Notes ----- - **The retained draws are a trimmed view, not the whole chain.** BUMPS' - DREAM sampler defaults to ``trim=True``: once sampling finishes it runs a - convergence-based burn-point detector over the chain and returns only the - portion after that point. ``state.draw()`` additionally drops chains - flagged as outliers. So ``results.draws`` is usually *smaller* than the - chain, and its length is not deterministic — the detector re-runs from - scratch on every ``sample()`` and ``extend()`` call and may place the burn - point differently each time. Read the count off the array; do not predict - it. - - Nor is the chain itself exactly ``samples / thin`` rows: ``samples`` is a - guaranteed minimum, not an exact count. DREAM advances in blocks of 10 - generations (one generation = one draw per chain) and only checks its - stopping condition between blocks, so the raw chain length is ``samples`` - rounded up to a multiple of ``10 * n_chains``. - - The trimming is only a *view* — nothing is dropped from the buffer. To - read the full untrimmed chain:: + **The retained draws are a trimmed view, not the whole chain.** + BUMPS' DREAM sampler defaults to ``trim=True``: once sampling + finishes it runs a convergence-based burn-point detector over the + chain and returns only the portion after that point. + ``state.draw()`` additionally drops chains flagged as outliers. So + ``results.draws`` is usually *smaller* than the chain, and its + length is not deterministic — the detector re-runs from scratch on + every ``sample()`` and ``extend()`` call and may place the burn + point differently each time. Read the count off the array; do not + predict it. + + Nor is the chain itself exactly ``samples / thin`` rows: ``samples`` + is a guaranteed minimum, not an exact count. DREAM advances in + blocks of 10 generations (one generation = one draw per chain) and + only checks its stopping condition between blocks, so the raw chain + length is ``samples`` rounded up to a multiple of ``10 * n_chains``. + + The trimming is only a *view* — nothing is dropped from the buffer. + To read the full untrimmed chain:: results = sampler.sample(samples=10000, burn=500, thin=2) full = results.state.draw(portion=1.0, outliers=True) full.points # (n_chain_rows, n_params) full.logp - To switch the automatic trimming off entirely, pass BUMPS' own ``trim`` - option straight through ``sampler_kwargs``; ``results.draws`` then holds - the whole chain:: + To switch the automatic trimming off entirely, pass BUMPS' own + ``trim`` option straight through ``sampler_kwargs``; + ``results.draws`` then holds the whole chain:: sampler = Sampler( fitter, x, y, weights=w, sampler_kwargs={'trim': False} ) - Note also that trimming does not survive a ``save()``/``load_state()`` - round-trip: BUMPS does not persist the trim point, so a reloaded chain - comes back untrimmed and ``load_state()`` reports *more* draws than the - ``sample()`` call that created it. The chain itself is identical. + Note also that trimming does not survive a + ``save()``/``load_state()`` round-trip: BUMPS does not persist the + trim point, so a reloaded chain comes back untrimmed and + ``load_state()`` reports *more* draws than the ``sample()`` call + that created it. The chain itself is identical. """ def __init__( @@ -435,7 +460,9 @@ def weights(self) -> np.ndarray | list[np.ndarray | None] | None: @property def state(self) -> MCMCDraw | None: - """Raw BUMPS MCMCDraw state (None before first sample/load_state).""" + """ + Raw BUMPS MCMCDraw state (None before first sample/load_state). + """ return self._state @property @@ -459,7 +486,9 @@ def logp(self) -> np.ndarray | None: return self._results.logp if self._results is not None else None def _fingerprint(self) -> str | None: - """SHA-256 fingerprint of the bound (x, y, weights) data, or None.""" + """ + SHA-256 fingerprint of the bound (x, y, weights) data, or None. + """ x_list = list(self._x) if isinstance(self._x, (list, tuple)) else [self._x] y_list = list(self._y) if isinstance(self._y, (list, tuple)) else [self._y] if self._weights is None: @@ -481,7 +510,8 @@ def _run( progress_callback: Callable[[dict], bool | None] | None, abort_test: Callable[[], bool] | None, ) -> SamplingResults: - """Shared sampling engine for ``sample()`` and ``extend()``. + """ + Shared sampling engine for ``sample()`` and ``extend()``. Argument validation for ``samples``/``burn``/``thin`` lives in ``Bumps.mcmc_sample`` (single source of truth). @@ -545,37 +575,43 @@ def sample( progress_callback: Callable[[dict], bool | None] | None = None, abort_test: Callable[[], bool] | None = None, ) -> SamplingResults: - """Run fresh Bayesian MCMC sampling on the bound data. + """ + Run fresh Bayesian MCMC sampling on the bound data. - Calling ``sample()`` on a sampler that already holds a chain starts a - **fresh** chain — the previous state and results are replaced. Use - ``extend()`` to continue an existing chain. + Calling ``sample()`` on a sampler that already holds a chain + starts a **fresh** chain — the previous state and results are + replaced. Use ``extend()`` to continue an existing chain. Parameters ---------- samples : int, default=10000 - Number of raw samples to draw across all chains, before thinning. - This is a guaranteed minimum, not an exact count: DREAM advances - in blocks of 10 generations (one generation = one draw per chain) - and stops at the first block boundary at or past ``samples``. + Number of raw samples to draw across all chains, before + thinning. This is a guaranteed minimum, not an exact count: + DREAM advances in blocks of 10 generations (one generation = + one draw per chain) and stops at the first block boundary at + or past ``samples``. burn : int, default=2000 - Burn-in generations to discard before collecting samples. Note - BUMPS counts ``burn`` in generations while ``samples`` counts raw - draws, so ``burn=500`` discards ``500 * n_chains`` raw samples. + Burn-in generations to discard before collecting samples. + Note BUMPS counts ``burn`` in generations while ``samples`` + counts raw draws, so ``burn=500`` discards ``500 * + n_chains`` raw samples. thin : int, default=10 - Thinning interval — only every ``thin``-th generation is kept, - which reduces autocorrelation between consecutive draws. + Thinning interval — only every ``thin``-th generation is + kept, which reduces autocorrelation between consecutive + draws. population : int | None, default=None - DREAM population **scale factor** (not an absolute chain count): - BUMPS creates ``ceil(population * n_parameters)`` parallel chains. + DREAM population **scale factor** (not an absolute chain + count): BUMPS creates ``ceil(population * n_parameters)`` + parallel chains. sampler_kwargs : dict | None, default=None Additional keyword arguments forwarded to the BUMPS DREAM sampler (merged over the instance defaults). progress_callback : Callable[[dict], bool | None] | None, default=None - Optional callback invoked at each DREAM generation. The payload - dict includes ``iteration`` and ``sampling: True``. + Optional callback invoked at each DREAM generation. The + payload dict includes ``iteration`` and ``sampling: True``. abort_test : Callable[[], bool] | None, default=None - Optional callable that returns ``True`` to abort sampling early. + Optional callable that returns ``True`` to abort sampling + early. Returns ------- @@ -584,14 +620,16 @@ def sample( Notes ----- - ``results.draws`` is a trimmed, outlier-filtered view of the chain, - so it is usually smaller than the chain and its length is not - predictable from ``samples``/``thin``. See the ``Sampler`` class notes - for how to read the full chain or switch trimming off. + ``results.draws`` is a trimmed, outlier-filtered view of the + chain, so it is usually smaller than the chain and its length is + not predictable from ``samples``/``thin``. See the ``Sampler`` + class notes for how to read the full chain or switch trimming + off. Exceptions propagate from the sampling engine: ``ValueError`` if - ``samples``, ``burn``, or ``thin`` are invalid, and ``RuntimeError`` - if the active minimizer is not a BUMPS instance. + ``samples``, ``burn``, or ``thin`` are invalid, and + ``RuntimeError`` if the active minimizer is not a BUMPS + instance. """ if self._state is not None: global_object.log.getLogger('fitting').warning( @@ -618,38 +656,43 @@ def extend( progress_callback: Callable[[dict], bool | None] | None = None, abort_test: Callable[[], bool] | None = None, ) -> SamplingResults: - """Continue the existing chain with additional samples. + """ + Continue the existing chain with additional samples. DREAM stores draws in a fixed-size ring buffer sized to its - ``samples`` parameter; this method does the ring-buffer arithmetic - for you (``samples = stored_generations * population + - additional_samples``) so no existing draws are dropped from the - buffer, regardless of the thinning interval. Runs with ``burn=0`` — - re-burning a converged chain is usually a mistake, and BUMPS forces - it to 0 on resume in any case. The DREAM population is recovered from - the saved state and cannot be changed on extend. + ``samples`` parameter; this method does the ring-buffer + arithmetic for you (``samples = stored_generations * population + + additional_samples``) so no existing draws are dropped from + the buffer, regardless of the thinning interval. Runs with + ``burn=0`` — re-burning a converged chain is usually a mistake, + and BUMPS forces it to 0 on resume in any case. The DREAM + population is recovered from the saved state and cannot be + changed on extend. Parameters ---------- additional_samples : int, default=5000 - Number of additional DREAM samples to draw, in the same units - as ``samples`` in ``sample()``. This grows the ring buffer by - exactly ``additional_samples``; how many of the new draws become - *visible* in ``results.draws`` depends on the thinning interval - and on BUMPS' automatic trimming (see Notes). + Number of additional DREAM samples to draw, in the same + units as ``samples`` in ``sample()``. This grows the ring + buffer by exactly ``additional_samples``; how many of the + new draws become *visible* in ``results.draws`` depends on + the thinning interval and on BUMPS' automatic trimming (see + Notes). thin : int, default=10 Thinning interval for the retained draws. total_samples : int | None, default=None - Advanced: total retained samples requested from the ring buffer, - **overriding** the ``additional_samples`` arithmetic. With - ``total_samples=N``, only the last N draws are retained. + Advanced: total retained samples requested from the ring + buffer, **overriding** the ``additional_samples`` + arithmetic. With ``total_samples=N``, only the last N draws + are retained. sampler_kwargs : dict | None, default=None Additional keyword arguments forwarded to the BUMPS DREAM sampler (merged over the instance defaults). progress_callback : Callable[[dict], bool | None] | None, default=None Optional callback invoked at each DREAM generation. abort_test : Callable[[], bool] | None, default=None - Optional callable that returns ``True`` to abort sampling early. + Optional callable that returns ``True`` to abort sampling + early. Returns ------- @@ -665,16 +708,17 @@ def extend( Notes ----- ``results.draws`` will **not** grow by exactly - ``additional_samples / thin``. Two things get in the way, neither of - them a re-applied burn-in: BUMPS re-runs its burn-point detector over - the whole extended chain and re-trims the returned view (so the - visible count can even shrink), and DREAM advances in blocks of 10 - generations, so the raw growth is ``additional_samples`` rounded up - to a multiple of ``10 * n_chains`` (i.e. at least - ``additional_samples / thin`` retained rows). To see the chain - itself, compare ``results.state.draw(portion=1.0, outliers=True)`` - before and after, or run with ``sampler_kwargs={'trim': False}``. See - the ``Sampler`` class notes. + ``additional_samples / thin``. Two things get in the way, + neither of them a re-applied burn-in: BUMPS re-runs its + burn-point detector over the whole extended chain and re-trims + the returned view (so the visible count can even shrink), and + DREAM advances in blocks of 10 generations, so the raw growth is + ``additional_samples`` rounded up to a multiple of ``10 * + n_chains`` (i.e. at least ``additional_samples / thin`` retained + rows). To see the chain itself, compare + ``results.state.draw(portion=1.0, outliers=True)`` before and + after, or run with ``sampler_kwargs={'trim': False}``. See the + ``Sampler`` class notes. """ if self._state is None: raise RuntimeError('No chain to extend. Call sample() or load_state() first.') @@ -698,14 +742,15 @@ def extend( ) def save(self, path: str | os.PathLike) -> None: - """Persist the chain state and metadata to disk. + """ + Persist the chain state and metadata to disk. Writes the BUMPS native files (``-chain.mc.gz``, ``-point.mc.gz`` and ``-stats.mc.gz``) plus a ``.params.json`` sidecar with the parameter names, the - easyscience version, and a fingerprint of the bound data (verified - with a warning on ``load_state()``). Use ``load_chain`` to read the - files back without a fitter. + easyscience version, and a fingerprint of the bound data + (verified with a warning on ``load_state()``). Use + ``load_chain`` to read the files back without a fitter. Parameters ---------- @@ -746,13 +791,14 @@ def save(self, path: str | os.PathLike) -> None: json.dump(sidecar, f, indent=2) def load_state(self, path: str | os.PathLike, skip: int = 0) -> SamplingResults: - """Load a previously saved chain into this sampler. + """ + Load a previously saved chain into this sampler. - The sampler must be constructed with the same fitter and data used to - create the chain — ``extend()`` then continues the saved chain. If the - sidecar carries a data fingerprint and it does not match this - sampler's bound data, a warning is logged (extending a chain against - different data is undefined behaviour). + The sampler must be constructed with the same fitter and data + used to create the chain — ``extend()`` then continues the saved + chain. If the sidecar carries a data fingerprint and it does not + match this sampler's bound data, a warning is logged (extending + a chain against different data is undefined behaviour). Populates ``state`` and ``results`` (draws, log-posterior and parameter names) from the saved chain, so summaries and @@ -763,8 +809,8 @@ def load_state(self, path: str | os.PathLike, skip: int = 0) -> SamplingResults: path : str | os.PathLike File path prefix used in ``save()``. skip : int, default=0 - Discard the first ``skip`` saved generations on load. Useful for - trimming additional burn-in without re-sampling. + Discard the first ``skip`` saved generations on load. Useful + for trimming additional burn-in without re-sampling. Returns ------- @@ -780,11 +826,11 @@ def load_state(self, path: str | os.PathLike, skip: int = 0) -> SamplingResults: Notes ----- - BUMPS does not persist its automatic trim point, so a reloaded chain - comes back **untrimmed**: this method reports more draws than the - ``sample()`` call that wrote the file, even though the chain is - identical. Use ``skip`` to discard leading generations explicitly. - See the ``Sampler`` class notes. + BUMPS does not persist its automatic trim point, so a reloaded + chain comes back **untrimmed**: this method reports more draws + than the ``sample()`` call that wrote the file, even though the + chain is identical. Use ``skip`` to discard leading generations + explicitly. See the ``Sampler`` class notes. """ # noqa: DOC502 -- raised in _validate_chain_path via load_chain state, param_names, sidecar = load_chain(path, skip=skip) diff --git a/src/easyscience/math.py b/src/easyscience/math.py new file mode 100644 index 00000000..5a1c7b86 --- /dev/null +++ b/src/easyscience/math.py @@ -0,0 +1,16 @@ +# SPDX-FileCopyrightText: 2026 EasyScience contributors +# SPDX-License-Identifier: BSD-3-Clause + +""" +Mathematical functions for EasyScience ``Parameter`` objects. + +Use them through the package namespace, for example +``es.math.sin(angle)``. Each function takes ``Parameter`` objects and +returns a new ``Parameter`` with the value, the propagated variance and +exact fitting bounds. The implementation lives in +``easyscience.variable.math``; this module only re-exports its public +names. +""" + +from easyscience.variable.math import * # noqa: F403 +from easyscience.variable.math import __all__ # noqa: F401 diff --git a/src/easyscience/variable/__init__.py b/src/easyscience/variable/__init__.py index 121d5cf9..a300efe9 100644 --- a/src/easyscience/variable/__init__.py +++ b/src/easyscience/variable/__init__.py @@ -7,14 +7,17 @@ from .descriptor_bool import DescriptorBool from .descriptor_number import DescriptorNumber from .descriptor_str import DescriptorStr +from .math import * # noqa: F403 +from .math import __all__ as _math_all from .parameter import Parameter __all__ = [ - DescriptorAnyType, - DescriptorArray, - DescriptorBase, - DescriptorBool, - DescriptorNumber, - DescriptorStr, - Parameter, + 'DescriptorAnyType', + 'DescriptorArray', + 'DescriptorBase', + 'DescriptorBool', + 'DescriptorNumber', + 'DescriptorStr', + 'Parameter', + *_math_all, ] diff --git a/src/easyscience/variable/descriptor_number.py b/src/easyscience/variable/descriptor_number.py index 5bcbcb0d..ed8f1cd5 100644 --- a/src/easyscience/variable/descriptor_number.py +++ b/src/easyscience/variable/descriptor_number.py @@ -10,6 +10,7 @@ from typing import Dict from typing import List from typing import Optional +from typing import Self from typing import Union import numpy as np @@ -115,7 +116,7 @@ def __init__( self._convert_unit(self._base_unit()) @classmethod - def from_scipp(cls, name: str, full_value: Variable, **kwargs: Any) -> DescriptorNumber: + def from_scipp(cls, name: str, full_value: Variable, **kwargs: Any) -> Self: """ Create a DescriptorNumber from a scipp constant. @@ -130,8 +131,8 @@ def from_scipp(cls, name: str, full_value: Variable, **kwargs: Any) -> Descripto Returns ------- - DescriptorNumber - DescriptorNumber. + Self + A new instance of the class this is called on. Raises ------ diff --git a/src/easyscience/variable/math.py b/src/easyscience/variable/math.py new file mode 100644 index 00000000..7c77adca --- /dev/null +++ b/src/easyscience/variable/math.py @@ -0,0 +1,1927 @@ +# SPDX-FileCopyrightText: 2026 EasyScience contributors +# SPDX-License-Identifier: BSD-3-Clause + +""" +Mathematical functions for EasyScience ``Parameter`` objects. + +Every public function takes one or more ``Parameter`` objects and +returns a new ``Parameter`` holding + +* the **value**, evaluated by Scipp (which also validates the units), +* the propagated **variance** (first-order "delta method", uncorrelated inputs), and +* the **fitting bounds**, i.e. the exact range the function can take when its + argument stays inside the argument's own bounds. + +Conventions: + +* *Values are validated strictly.* A value outside the function's domain, or a + non-finite result, raises ``ValueError``; such values would break a fit loop. +* *Bounds are permissive.* They never raise: bounds outside the domain are + clipped to it, infinite bounds are accepted, ``tan`` bounds that cross an + asymptote become ``(-inf, inf)``, and bounds that collapse to (nearly) a single + number are padded by a relative 2e-9 around the value. + +Input and result classes: + +Every function accepts a ``DescriptorNumber`` or a ``Parameter`` (a +``Parameter`` is a ``DescriptorNumber``), including user subclasses of +either. The result is a ``Parameter``, with bounds, if any argument is a +``Parameter``, and a ``DescriptorNumber``, without bounds, otherwise. It +is always the base class and never the subclass, as for the arithmetic +operators. A ``DescriptorNumber`` argument counts as the single point at +its value when the bounds of a ``Parameter`` result are derived. + +Using ``es.math`` in a model: + +``es.math`` works on scalar ``Parameter`` objects, so it cannot be +applied over an ``x`` array. Use it to derive scalar quantities, with +their uncertainty, from the parameters of a model or of a fit. Each +result is a new ``Parameter`` that is a *snapshot* of its inputs: it +does not follow later changes. Inside a model, derive it in ``__call__`` +so that it is recomputed on every evaluation of the fit, and not once in +``__init__``: + +```python +class GaussianPeak(ModelBase): + def sigma(self): + two = self._two + return self._fwhm / (2 * es.math.sqrt(2 * es.math.log(two))) + + def __call__(self, x): + sigma = self.sigma().value # derived on every evaluation + shift = (x - self.center.value) ** 2 + return self.amplitude.value * np.exp(-shift / (2 * sigma**2)) +``` + +Two things to keep in mind. Propagated variances assume uncorrelated +inputs, because a ``Parameter`` stores a variance but no covariance; a +function of two strongly correlated fit parameters (for example +``arctan2`` of the slope and intercept of a line) can be off by a large +factor. And a fit reports different uncertainties for different +minimizers, so every derived variance is only as good as the variance it +starts from. + +Implementation: each function is described by a small, immutable +``_MathOp`` specification (which Scipp function to call, which domain it +accepts, how to derive its bounds, ...). A single dispatcher, +``_apply``, does the work for all of them, so adding a function means +adding a specification, not new control flow. + +The public entry point is ``easyscience.math`` (``es.math``), which +re-exports the names in ``__all__``. They are also available from +``easyscience.variable``. +""" + +import math +import warnings +from collections.abc import Callable +from dataclasses import dataclass +from functools import cache +from typing import TypeVar +from typing import overload + +import numpy as np +import scipp as sc + +from easyscience.variable.descriptor_number import DescriptorNumber +from easyscience.variable.parameter import Parameter + +# A closed interval ``(lo, hi)`` of an argument, in the unit the bounds strategy works in. +_Interval = tuple[float, float] +# A bounds strategy maps one interval per argument to the ``(min, max)`` of the result. +_BoundsStrategy = Callable[..., tuple[float, float]] +# A function returns the class of its input (see the module docstring), as Scipp's do. +_D = TypeVar('_D', bound=DescriptorNumber) + +_HALF_PI = np.pi / 2 +_TWO_PI = 2 * np.pi + +# tan: a value closer than this (in rad) to an asymptote is rejected whatever its variance; +# a larger standard deviation widens the margin. +_TAN_ASYMPTOTE_TOLERANCE = 1e-7 + +# sinc: the extrema are the roots of tan(x) = x. Two consecutive roots hold the largest +# maximum and the deepest minimum of an interval and a root just below ``lo`` may be +# filtered out, so four roots (one more than the three strictly needed) are refined. +_SINC_ROOT_COUNT = 4 +# Newton's method converges quadratically from the asymptotic first guess; 8 steps is far +# more than the ~4 needed to reach machine precision. +_SINC_NEWTON_STEPS = 8 +# Beyond this the lobes of sinc are below 1e-15 (they fall off like 1/x), and the lobe index +# would no longer fit numpy's integer range; the endpoints give the bounds to that precision. +_SINC_NEGLIGIBLE_LOBES_FROM = 1e15 + +# tan: a bound within this distance (in units of pi) of an asymptote counts as sitting on it, so +# that 90 deg converted to radians is recognised despite round-off. +_TAN_EDGE_TOLERANCE = 1e-12 + +# ``Parameter`` refuses ``min`` and ``max`` that are equal or within a relative 1e-9 of each +# other. Result bounds that collapse (a saturating or underflowing function, or an extremum +# at which the function is quadratic) are therefore padded around the value by twice that +# relative tolerance, or by an absolute amount if the value is 0. +_COLLAPSE_PAD = 2e-9 +_COLLAPSE_PAD_ZERO = 1e-12 + +__all__ = [ + 'acos', + 'acosh', + 'arccos', + 'arccosh', + 'arcsin', + 'arcsinh', + 'arctan', + 'arctan2', + 'arctanh', + 'asin', + 'asinh', + 'atan', + 'atan2', + 'atanh', + 'cos', + 'cosh', + 'erf', + 'erfc', + 'exp', + 'log', + 'log10', + 'sin', + 'sinc', + 'sinh', + 'sqrt', + 'tan', + 'tanh', +] + + +# ========================================== +# Public functions: standard maths +# ========================================== + + +def exp(x: _D) -> _D: + """ + Calculate the exponential, ``e**x``, of a dimensionless parameter. + + Parameters + ---------- + x : _D + Dimensionless input. + + Returns + ------- + _D + Parameter or DescriptorNumber + ``exp(x)`` with propagated variance and bounds ``(exp(x.min), + exp(x.max))``. The class follows the inputs: a ``Parameter`` + (with bounds) if any argument is one, otherwise a + ``DescriptorNumber`` (without bounds). + + Raises + ------ + NotImplementedError + If ``x`` is not a ``Parameter`` or ``DescriptorNumber``. + UnitError + If ``x`` is not dimensionless. + ValueError + If the result overflows (``x`` greater than about 709.78) or, + when ``x`` has a variance, if the propagated variance overflows. + The variance scales with ``exp(2x)``, so that happens from about + 355 for a unit variance. + + Examples + -------- + + Exponential of a dimensionless parameter, with its bounds: + + >>> import easyscience as es + >>> x = es.Parameter('x', 1.0, min=0.0, max=2.0) + >>> result = es.math.exp(x) + >>> print(f'{result.value:.4f} in [{result.min:.4f}, {result.max:.4f}]') + 2.7183 in [1.0000, 7.3891] + + The variance is propagated (first order): + + >>> x = es.Parameter('x', 0.0, variance=0.01) + >>> result = es.math.exp(x) + >>> print(f'{result.value:.4f}, variance {result.variance:.4f}') + 1.0000, variance 0.0100 + """ # noqa: DOC502, DOC203, DOC105 -- see _apply + return _apply(_EXP, x) + + +def log(x: _D) -> _D: + """ + Calculate the natural logarithm of a dimensionless parameter. + + Parameters + ---------- + x : _D + Dimensionless input. Its value must be strictly positive. + + Returns + ------- + _D + Parameter or DescriptorNumber + ``log(x)`` with propagated variance and bounds. The class + follows the inputs: a ``Parameter`` (with bounds) if any + argument is one, otherwise a ``DescriptorNumber`` (without + bounds). + + Raises + ------ + NotImplementedError + If ``x`` is not a ``Parameter`` or ``DescriptorNumber``. + UnitError + If ``x`` is not dimensionless. + ValueError + If the value of ``x`` is not strictly positive. + + Notes + ----- + Bounds below zero are clipped to the domain, so ``x`` in ``[-1, 5]`` + gives result bounds ``[-inf, log(5)]``. + + Examples + -------- + + Natural logarithm, with its bounds: + + >>> import easyscience as es + >>> x = es.Parameter('x', 10.0, min=1.0, max=100.0) + >>> result = es.math.log(x) + >>> print(f'{result.value:.4f} in [{result.min:.4f}, {result.max:.4f}]') + 2.3026 in [0.0000, 4.6052] + + Bounds below zero are clipped to the domain: + + >>> x = es.Parameter('x', 2.0, min=-1.0, max=5.0) + >>> result = es.math.log(x) + >>> print(f'bounds [{result.min:.4f}, {result.max:.4f}]') + bounds [-inf, 1.6094] + """ # noqa: DOC502, DOC203, DOC105 -- see _apply + return _apply(_LOG, x) + + +def log10(x: _D) -> _D: + """ + Calculate the base-10 logarithm of a dimensionless parameter. + + Parameters + ---------- + x : _D + Dimensionless input. Its value must be strictly positive. + + Returns + ------- + _D + Parameter or DescriptorNumber + ``log10(x)`` with propagated variance and bounds. The class + follows the inputs: a ``Parameter`` (with bounds) if any + argument is one, otherwise a ``DescriptorNumber`` (without + bounds). + + Raises + ------ + NotImplementedError + If ``x`` is not a ``Parameter`` or ``DescriptorNumber``. + UnitError + If ``x`` is not dimensionless. + ValueError + If the value of ``x`` is not strictly positive. + + Notes + ----- + Bounds below zero are clipped to the domain, so ``x`` in ``[0, + 100]`` gives result bounds ``[-inf, 2]``. + + Examples + -------- + + Base-10 logarithm, with its bounds: + + >>> import easyscience as es + >>> x = es.Parameter('x', 100.0, min=10.0, max=1000.0) + >>> result = es.math.log10(x) + >>> print(f'{result.value:.4f} in [{result.min:.4f}, {result.max:.4f}]') + 2.0000 in [1.0000, 3.0000] + + Bounds touching zero are clipped to the domain: + + >>> x = es.Parameter('x', 10.0, min=0.0, max=100.0) + >>> result = es.math.log10(x) + >>> print(f'bounds [{result.min:.4f}, {result.max:.4f}]') + bounds [-inf, 2.0000] + """ # noqa: DOC502, DOC203, DOC105 -- see _apply + return _apply(_LOG10, x) + + +def sqrt(x: _D) -> _D: + """ + Calculate the square root of a parameter. + + Parameters + ---------- + x : _D + Input that is dimensionless, or whose unit is a perfect square + (e.g. ``m**2``). Its value must not be negative. + + Returns + ------- + _D + Parameter or DescriptorNumber + ``sqrt(x)`` with propagated variance and bounds. The unit is the + square root of the input unit (``m**2`` gives ``m``). The class + follows the inputs: a ``Parameter`` (with bounds) if any + argument is one, otherwise a ``DescriptorNumber`` (without + bounds). + + Raises + ------ + NotImplementedError + If ``x`` is not a ``Parameter`` or ``DescriptorNumber``. + UnitError + If the unit of ``x`` has no square root in Scipp (e.g. ``m``). + ValueError + If the value of ``x`` is negative, or if the result is not + finite (for example a non-zero variance at ``x == 0``, where the + derivative is infinite). + + Notes + ----- + Bounds below zero are clipped to the domain. EasyScience + ``Parameter`` objects also support ``x ** 0.5``. + + Examples + -------- + + Square root of a dimensionless parameter, with its bounds: + + >>> import easyscience as es + >>> x = es.Parameter('x', 4.0, min=1.0, max=9.0) + >>> result = es.math.sqrt(x) + >>> print(f'{result.value:.4f} in [{result.min:.4f}, {result.max:.4f}]') + 2.0000 in [1.0000, 3.0000] + + The unit is square-rooted too: + + >>> area = es.Parameter('area', 9.0, unit='m**2') + >>> result = es.math.sqrt(area) + >>> print(result.value, result.unit) + 3.0 m + """ # noqa: DOC502, DOC203, DOC105 -- see _apply + return _apply(_SQRT, x) + + +def erf(x: _D) -> _D: + """ + Calculate the error function of a dimensionless parameter. + + Parameters + ---------- + x : _D + Dimensionless input. + + Returns + ------- + _D + Parameter or DescriptorNumber + ``erf(x)`` with propagated variance and bounds in ``[-1, 1]``. + The class follows the inputs: a ``Parameter`` (with bounds) if + any argument is one, otherwise a ``DescriptorNumber`` (without + bounds). + + Raises + ------ + NotImplementedError + If ``x`` is not a ``Parameter`` or ``DescriptorNumber``. + UnitError + If ``x`` is not dimensionless. + + Notes + ----- + Scipp does not propagate variances through ``erf``. EasyScience does + it with the delta method, using ``d erf / dx = 2 / sqrt(pi) * + exp(-x**2)``. + + Examples + -------- + + Error function, with its bounds: + + >>> import easyscience as es + >>> x = es.Parameter('x', 1.0, min=0.0, max=2.0) + >>> result = es.math.erf(x) + >>> print(f'{result.value:.4f} in [{result.min:.4f}, {result.max:.4f}]') + 0.8427 in [0.0000, 0.9953] + + The variance is propagated, although Scipp cannot do it for erf: + + >>> x = es.Parameter('x', 0.0, variance=0.01) + >>> result = es.math.erf(x) + >>> print(f'{result.value:.4f}, variance {result.variance:.4f}') + 0.0000, variance 0.0127 + """ # noqa: DOC502, DOC203, DOC105 -- see _apply + return _apply(_ERF, x) + + +def erfc(x: _D) -> _D: + """ + Calculate the complementary error function, ``1 - erf(x)``, of a + dimensionless parameter. + + Parameters + ---------- + x : _D + Dimensionless input. + + Returns + ------- + _D + Parameter or DescriptorNumber + ``erfc(x)`` with propagated variance and bounds in ``[0, 2]``. + The class follows the inputs: a ``Parameter`` (with bounds) if + any argument is one, otherwise a ``DescriptorNumber`` (without + bounds). + + Raises + ------ + NotImplementedError + If ``x`` is not a ``Parameter`` or ``DescriptorNumber``. + UnitError + If ``x`` is not dimensionless. + + Notes + ----- + Scipp does not propagate variances through ``erfc``. EasyScience + does it with the delta method, using ``d erfc / dx = -2 / sqrt(pi) * + exp(-x**2)``. The function is decreasing, so the bounds are swapped + relative to the input. + + Examples + -------- + + Complementary error function; it is decreasing, so the bounds swap: + + >>> import easyscience as es + >>> x = es.Parameter('x', 1.0, min=0.0, max=2.0) + >>> result = es.math.erfc(x) + >>> print(f'{result.value:.4f} in [{result.min:.4f}, {result.max:.4f}]') + 0.1573 in [0.0047, 1.0000] + + The variance is propagated, although Scipp cannot do it for erfc: + + >>> x = es.Parameter('x', 0.0, variance=0.01) + >>> result = es.math.erfc(x) + >>> print(f'{result.value:.4f}, variance {result.variance:.4f}') + 1.0000, variance 0.0127 + """ # noqa: DOC502, DOC203, DOC105 -- see _apply + return _apply(_ERFC, x) + + +# ========================================== +# Public functions: trigonometric +# ========================================== + + +def sin(x: _D) -> _D: + """ + Calculate the trigonometric sine of an angle parameter. + + Parameters + ---------- + x : _D + Angle, with unit ``rad`` or ``deg``. + + Returns + ------- + _D + Parameter or DescriptorNumber + Dimensionless ``sin(x)`` with propagated variance and exact + bounds. The class follows the inputs: a ``Parameter`` (with + bounds) if any argument is one, otherwise a ``DescriptorNumber`` + (without bounds). + + Raises + ------ + NotImplementedError + If ``x`` is not a ``Parameter`` or ``DescriptorNumber``. + UnitError + If the unit of ``x`` is not an angle. + + Notes + ----- + The bounds are exact: the endpoints and the interior maxima and + minima (``pi/2 + n*pi``) are evaluated, in radians whatever the + input unit. Infinite bounds, or bounds wider than one period, give + ``[-1, 1]``. + + Examples + -------- + + Angle in radians: + + >>> import easyscience as es + >>> import numpy as np + >>> angle = es.Parameter('angle', np.pi / 2, unit='rad', min=0.0, max=np.pi) + >>> result = es.math.sin(angle) + >>> print(f'{result.value:.4f} in [{result.min:.4f}, {result.max:.4f}]') + 1.0000 in [0.0000, 1.0000] + + Angle in degrees; the bounds are computed in radians internally: + + >>> angle = es.Parameter('angle', 30.0, unit='deg', min=0.0, max=90.0) + >>> result = es.math.sin(angle) + >>> print(f'{result.value:.4f} in [{result.min:.4f}, {result.max:.4f}]') + 0.5000 in [0.0000, 1.0000] + """ # noqa: DOC502, DOC203, DOC105 -- see _apply + return _apply(_SIN, x) + + +def cos(x: _D) -> _D: + """ + Calculate the trigonometric cosine of an angle parameter. + + Parameters + ---------- + x : _D + Angle, with unit ``rad`` or ``deg``. + + Returns + ------- + _D + Parameter or DescriptorNumber + Dimensionless ``cos(x)`` with propagated variance and exact + bounds. The class follows the inputs: a ``Parameter`` (with + bounds) if any argument is one, otherwise a ``DescriptorNumber`` + (without bounds). + + Raises + ------ + NotImplementedError + If ``x`` is not a ``Parameter`` or ``DescriptorNumber``. + UnitError + If the unit of ``x`` is not an angle. + + Notes + ----- + The bounds are exact: the endpoints and the interior maxima and + minima (``n*pi``) are evaluated, in radians whatever the input unit. + Infinite bounds, or bounds wider than one period, give ``[-1, 1]``. + + Examples + -------- + + Angle in radians: + + >>> import easyscience as es + >>> import numpy as np + >>> angle = es.Parameter( + ... 'angle', np.pi / 3, unit='rad', min=0.0, max=np.pi / 2 + ... ) + >>> result = es.math.cos(angle) + >>> print(f'{result.value:.4f} in [{result.min:.4f}, {result.max:.4f}]') + 0.5000 in [0.0000, 1.0000] + + Angle in degrees: + + >>> angle = es.Parameter('angle', 60.0, unit='deg', min=0.0, max=180.0) + >>> result = es.math.cos(angle) + >>> print(f'{result.value:.4f} in [{result.min:.4f}, {result.max:.4f}]') + 0.5000 in [-1.0000, 1.0000] + """ # noqa: DOC502, DOC203, DOC105 -- see _apply + return _apply(_COS, x) + + +def tan(x: _D) -> _D: + """ + Calculate the trigonometric tangent of an angle parameter. + + Parameters + ---------- + x : _D + Angle, with unit ``rad`` or ``deg``. + + Returns + ------- + _D + Parameter or DescriptorNumber + Dimensionless ``tan(x)`` with propagated variance and bounds. + The class follows the inputs: a ``Parameter`` (with bounds) if + any argument is one, otherwise a ``DescriptorNumber`` (without + bounds). + + Raises + ------ + NotImplementedError + If ``x`` is not a ``Parameter`` or ``DescriptorNumber``. + UnitError + If the unit of ``x`` is not an angle. + ValueError + If the value of ``x`` is within 1e-7 rad, or one standard + deviation if that is larger, of an asymptote (``pi/2 + n*pi``), + where the propagated variance would explode. + + Notes + ----- + The *value* is guarded against asymptotes, but the *bounds* are + allowed to cross them: if ``x.min`` and ``x.max`` lie on different + branches of ``tan``, the result bounds are ``(-inf, inf)``. A bound + exactly on an asymptote makes the result bound infinite on that side + only, so ``[0, 90]`` deg gives ``[0, inf)``. + + Examples + -------- + + Angle in radians: + + >>> import easyscience as es + >>> import numpy as np + >>> angle = es.Parameter( + ... 'angle', np.pi / 4, unit='rad', min=0.0, max=np.pi / 3 + ... ) + >>> result = es.math.tan(angle) + >>> print(f'{result.value:.4f} in [{result.min:.4f}, {result.max:.4f}]') + 1.0000 in [0.0000, 1.7321] + + Angle in degrees; bounds that cross an asymptote (90 deg) become + infinite: + + >>> angle = es.Parameter('angle', 45.0, unit='deg', min=0.0, max=120.0) + >>> result = es.math.tan(angle) + >>> print(f'bounds [{result.min:.4f}, {result.max:.4f}]') + bounds [-inf, inf] + """ # noqa: DOC502, DOC203, DOC105 -- see _apply + return _apply(_TAN, x) + + +def arcsin(x: _D, *, unit: str = 'rad') -> _D: + """ + Calculate the inverse sine of a dimensionless parameter. + + Parameters + ---------- + x : _D + Dimensionless input; its value must be in ``[-1, 1]``. + unit : str, default='rad' + Unit of the returned angle, ``'rad'`` or ``'deg'``. The value, + the variance and the bounds are all converted. + + Returns + ------- + _D + Parameter or DescriptorNumber + Angle in ``unit`` (``[-pi/2, pi/2]`` rad), with propagated + variance and bounds. The class follows the inputs: a + ``Parameter`` (with bounds) if any argument is one, otherwise a + ``DescriptorNumber`` (without bounds). + + Raises + ------ + NotImplementedError + If ``x`` is not a ``Parameter`` or ``DescriptorNumber``. + UnitError + If ``x`` is not dimensionless, or ``unit`` is not an angle unit. + ValueError + If the value of ``x`` is outside ``[-1, 1]``, or if the result + is not finite (a non-zero variance at ``x == ±1``, where the + derivative is infinite). + + Notes + ----- + Bounds outside ``[-1, 1]`` are clipped to the domain. + + Examples + -------- + + The result is an angle in radians: + + >>> import easyscience as es + >>> x = es.Parameter('x', 0.5, min=0.0, max=1.0) + >>> result = es.math.arcsin(x) + >>> result.unit + 'rad' + + Ask for degrees with ``unit``; the value, variance and bounds + follow: + + >>> x = es.Parameter('x', 0.5, min=0.0, max=1.0) + >>> result = es.math.arcsin(x, unit='deg') + >>> print( + ... f'{result.value:.4f} {result.unit} in [{result.min:.4f}, {result.max:.4f}]' + ... ) + 30.0000 deg in [0.0000, 90.0000] + """ # noqa: DOC502, DOC203, DOC105 -- see _apply and _angle_in + return _angle_in(_apply(_ARCSIN, x), unit) + + +def arccos(x: _D, *, unit: str = 'rad') -> _D: + """ + Calculate the inverse cosine of a dimensionless parameter. + + Parameters + ---------- + x : _D + Dimensionless input; its value must be in ``[-1, 1]``. + unit : str, default='rad' + Unit of the returned angle, ``'rad'`` or ``'deg'``. The value, + the variance and the bounds are all converted. + + Returns + ------- + _D + Parameter or DescriptorNumber + Angle in ``unit`` (``[0, pi]`` rad), with propagated variance + and bounds. The class follows the inputs: a ``Parameter`` (with + bounds) if any argument is one, otherwise a ``DescriptorNumber`` + (without bounds). + + Raises + ------ + NotImplementedError + If ``x`` is not a ``Parameter`` or ``DescriptorNumber``. + UnitError + If ``x`` is not dimensionless, or ``unit`` is not an angle unit. + ValueError + If the value of ``x`` is outside ``[-1, 1]``, or if the result + is not finite (a non-zero variance at ``x == ±1``, where the + derivative is infinite). + + Notes + ----- + Bounds outside ``[-1, 1]`` are clipped to the domain. The function + is decreasing, so ``x.min`` maps to the result's maximum. + + Examples + -------- + + The result is an angle in radians; the function is decreasing, so + the bounds swap: + + >>> import easyscience as es + >>> x = es.Parameter('x', 0.5, min=0.0, max=1.0) + >>> result = es.math.arccos(x) + >>> print( + ... f'{result.value:.4f} {result.unit} in [{result.min:.4f}, {result.max:.4f}]' + ... ) + 1.0472 rad in [0.0000, 1.5708] + + Ask for degrees with ``unit``; the value, variance and bounds + follow: + + >>> x = es.Parameter('x', 0.5, min=0.0, max=1.0) + >>> result = es.math.arccos(x, unit='deg') + >>> print( + ... f'{result.value:.4f} {result.unit} in [{result.min:.4f}, {result.max:.4f}]' + ... ) + 60.0000 deg in [0.0000, 90.0000] + """ # noqa: DOC502, DOC203, DOC105 -- see _apply and _angle_in + return _angle_in(_apply(_ARCCOS, x), unit) + + +def arctan(x: _D, *, unit: str = 'rad') -> _D: + """ + Calculate the inverse tangent of a dimensionless parameter. + + Parameters + ---------- + x : _D + Dimensionless input. + unit : str, default='rad' + Unit of the returned angle, ``'rad'`` or ``'deg'``. The value, + the variance and the bounds are all converted. + + Returns + ------- + _D + Parameter or DescriptorNumber + Angle in ``unit`` (``[-pi/2, pi/2]`` rad), with propagated + variance and bounds. The class follows the inputs: a + ``Parameter`` (with bounds) if any argument is one, otherwise a + ``DescriptorNumber`` (without bounds). + + Raises + ------ + NotImplementedError + If ``x`` is not a ``Parameter`` or ``DescriptorNumber``. + UnitError + If ``x`` is not dimensionless, or ``unit`` is not an angle unit. + + Examples + -------- + + The result is an angle in radians: + + >>> import easyscience as es + >>> x = es.Parameter('x', 1.0, min=0.0, max=10.0) + >>> result = es.math.arctan(x) + >>> print( + ... f'{result.value:.4f} {result.unit} in [{result.min:.4f}, {result.max:.4f}]' + ... ) + 0.7854 rad in [0.0000, 1.4711] + + Ask for degrees with ``unit``; the value, variance and bounds + follow: + + >>> x = es.Parameter('x', 1.0, min=0.0, max=10.0) + >>> result = es.math.arctan(x, unit='deg') + >>> print( + ... f'{result.value:.4f} {result.unit} in [{result.min:.4f}, {result.max:.4f}]' + ... ) + 45.0000 deg in [0.0000, 84.2894] + """ # noqa: DOC502, DOC203, DOC105 -- see _apply and _angle_in + return _angle_in(_apply(_ARCTAN, x), unit) + + +@overload +def arctan2(y: Parameter, x: DescriptorNumber, *, unit: str = ...) -> Parameter: ... + + +@overload +def arctan2(y: DescriptorNumber, x: Parameter, *, unit: str = ...) -> Parameter: ... + + +@overload +def arctan2(y: DescriptorNumber, x: DescriptorNumber, *, unit: str = ...) -> DescriptorNumber: ... + + +def arctan2(y: DescriptorNumber, x: DescriptorNumber, *, unit: str = 'rad') -> DescriptorNumber: + """ + Calculate the angle of the point ``(x, y)``, using both signs to + find the quadrant. + + Parameters + ---------- + y : DescriptorNumber + Cartesian y-coordinate. + x : DescriptorNumber + Cartesian x-coordinate. Must have the same unit as ``y``. + unit : str, default='rad' + Unit of the returned angle, ``'rad'`` or ``'deg'``. The value, + the variance and the bounds are all converted. + + Returns + ------- + DescriptorNumber + Parameter or DescriptorNumber + Angle in ``unit`` (``[-pi, pi]`` rad), with propagated variance + and bounds. The class follows the inputs: a ``Parameter`` (with + bounds) if any argument is one, otherwise a ``DescriptorNumber`` + (without bounds). + + Raises + ------ + NotImplementedError + If ``y`` or ``x`` is not a ``Parameter`` or + ``DescriptorNumber``. + UnitError + If ``y`` and ``x`` have different units, or ``unit`` is not an + angle unit. + ValueError + If an argument has a variance and the point ``(x, y)`` lies + within the larger of the two standard deviations of the origin, + where the angle is undefined, or if the result is not finite. + + Warns + ----- + UserWarning If all four bounds are finite and the box they span + touches the branch cut. + + Notes + ----- + Scipp does not propagate variances through ``atan2``. EasyScience + does it with the delta method, assuming ``y`` and ``x`` are + uncorrelated. + + The bounds are the exact range of the angle over the box ``[y.min, + y.max] x [x.min, x.max]``. If the box touches the branch cut (``y = + 0`` with ``x < 0``) the angle jumps between ``-pi`` and ``pi``; the + bounds are then ``[-pi, pi]``, and a ``UserWarning`` is raised when + all four bounds are finite, since gradient-based minimisers behave + badly across the jump. + + Examples + -------- + + The angle of the point (x, y) = (1, 1) m, in radians: + + >>> import easyscience as es + >>> y = es.Parameter('y', 1.0, unit='m', min=0.5, max=2.0) + >>> x = es.Parameter('x', 1.0, unit='m', min=1.0, max=3.0) + >>> result = es.math.arctan2(y, x) + >>> print( + ... f'{result.value:.4f} {result.unit} in [{result.min:.4f}, {result.max:.4f}]' + ... ) + 0.7854 rad in [0.1651, 1.1071] + + Ask for degrees with ``unit``; the value, variance and bounds + follow: + + >>> y = es.Parameter('y', 1.0, unit='m', min=0.5, max=2.0) + >>> x = es.Parameter('x', 1.0, unit='m', min=1.0, max=3.0) + >>> result = es.math.arctan2(y, x, unit='deg') + >>> print( + ... f'{result.value:.4f} {result.unit} in [{result.min:.4f}, {result.max:.4f}]' + ... ) + 45.0000 deg in [9.4623, 63.4349] + """ # noqa: DOC502, DOC203, DOC105 -- see _apply and _angle_in + return _angle_in(_apply(_ARCTAN2, y, x), unit) + + +def sinc(x: _D) -> _D: + """ + Calculate the unnormalised sinc function, ``sin(x) / x``, of a + parameter. + + Parameters + ---------- + x : _D + Input, with unit ``rad`` or dimensionless. + + Returns + ------- + _D + Parameter or DescriptorNumber + Dimensionless ``sinc(x)`` (equal to 1 at ``x = 0``) with + propagated variance and exact bounds. The class follows the + inputs: a ``Parameter`` (with bounds) if any argument is one, + otherwise a ``DescriptorNumber`` (without bounds). + + Raises + ------ + NotImplementedError + If ``x`` is not a ``Parameter`` or ``DescriptorNumber``. + UnitError + If the unit of ``x`` is neither ``rad`` nor dimensionless. + + Notes + ----- + This matches Scipp and the physics convention, not the normalised + ``sin(pi x) / (pi x)`` of ``numpy.sinc``. + + The bounds are exact. The extrema of sinc lie at ``x = 0`` and at + the roots of ``tan(x) = x``; the roots inside the interval are found + numerically. The global minimum, about -0.2172, is at ``x`` of about + 4.4934. + + Examples + -------- + + Angle in radians: + + >>> import easyscience as es + >>> import numpy as np + >>> angle = es.Parameter('angle', np.pi / 2, unit='rad', min=0.0, max=10.0) + >>> result = es.math.sinc(angle) + >>> print(f'{result.value:.4f} in [{result.min:.4f}, {result.max:.4f}]') + 0.6366 in [-0.2172, 1.0000] + + sinc does not accept degrees; convert an angle in degrees first: + + >>> angle = es.Parameter('angle', 90.0, unit='deg') + >>> angle.convert_unit('rad') + >>> result = es.math.sinc(angle) + >>> print(f'{result.value:.4f} in [{result.min:.4f}, {result.max:.4f}]') + 0.6366 in [-0.2172, 1.0000] + """ # noqa: DOC502, DOC203, DOC105 -- see _apply + return _apply(_SINC, x) + + +# ========================================== +# Public functions: hyperbolic +# ========================================== + + +def sinh(x: _D) -> _D: + """ + Calculate the hyperbolic sine of a dimensionless parameter. + + Parameters + ---------- + x : _D + Dimensionless input. + + Returns + ------- + _D + Parameter or DescriptorNumber + ``sinh(x)`` with propagated variance and bounds. The class + follows the inputs: a ``Parameter`` (with bounds) if any + argument is one, otherwise a ``DescriptorNumber`` (without + bounds). + + Raises + ------ + NotImplementedError + If ``x`` is not a ``Parameter`` or ``DescriptorNumber``. + UnitError + If ``x`` is not dimensionless. + ValueError + If the result overflows (``|x|`` greater than about 710) or, + when ``x`` has a variance, if the propagated variance overflows + (from ``|x|`` of about 355 for a unit variance). + + Examples + -------- + + Hyperbolic sine, with its bounds: + + >>> import easyscience as es + >>> x = es.Parameter('x', 1.0, min=0.0, max=2.0) + >>> result = es.math.sinh(x) + >>> print(f'{result.value:.4f} in [{result.min:.4f}, {result.max:.4f}]') + 1.1752 in [0.0000, 3.6269] + + The variance is propagated (first order): + + >>> x = es.Parameter('x', 0.0, variance=0.01) + >>> result = es.math.sinh(x) + >>> print(f'{result.value:.4f}, variance {result.variance:.4f}') + 0.0000, variance 0.0100 + """ # noqa: DOC502, DOC203, DOC105 -- see _apply + return _apply(_SINH, x) + + +def cosh(x: _D) -> _D: + """ + Calculate the hyperbolic cosine of a dimensionless parameter. + + Parameters + ---------- + x : _D + Dimensionless input. + + Returns + ------- + _D + Parameter or DescriptorNumber + ``cosh(x)`` with propagated variance and bounds. The class + follows the inputs: a ``Parameter`` (with bounds) if any + argument is one, otherwise a ``DescriptorNumber`` (without + bounds). + + Raises + ------ + NotImplementedError + If ``x`` is not a ``Parameter`` or ``DescriptorNumber``. + UnitError + If ``x`` is not dimensionless. + ValueError + If the result overflows (``|x|`` greater than about 710) or, + when ``x`` has a variance, if the propagated variance overflows + (from ``|x|`` of about 355 for a unit variance). + + Notes + ----- + ``cosh`` is not monotonic: if the bounds contain 0, the lower result + bound is 1. + + Examples + -------- + + Hyperbolic cosine; the minimum, 1, is inside the bounds: + + >>> import easyscience as es + >>> x = es.Parameter('x', 0.5, min=-1.0, max=2.0) + >>> result = es.math.cosh(x) + >>> print(f'{result.value:.4f} in [{result.min:.4f}, {result.max:.4f}]') + 1.1276 in [1.0000, 3.7622] + + Bounds on one side of zero: cosh is monotonic there: + + >>> x = es.Parameter('x', 1.0, min=0.5, max=2.0) + >>> result = es.math.cosh(x) + >>> print(f'{result.value:.4f} in [{result.min:.4f}, {result.max:.4f}]') + 1.5431 in [1.1276, 3.7622] + """ # noqa: DOC502, DOC203, DOC105 -- see _apply + return _apply(_COSH, x) + + +def tanh(x: _D) -> _D: + """ + Calculate the hyperbolic tangent of a dimensionless parameter. + + Parameters + ---------- + x : _D + Dimensionless input. + + Returns + ------- + _D + Parameter or DescriptorNumber + ``tanh(x)`` with propagated variance and bounds in ``[-1, 1]``. + The class follows the inputs: a ``Parameter`` (with bounds) if + any argument is one, otherwise a ``DescriptorNumber`` (without + bounds). + + Raises + ------ + NotImplementedError + If ``x`` is not a ``Parameter`` or ``DescriptorNumber``. + UnitError + If ``x`` is not dimensionless. + + Examples + -------- + + Hyperbolic tangent, with its bounds: + + >>> import easyscience as es + >>> x = es.Parameter('x', 0.5, min=0.0, max=2.0) + >>> result = es.math.tanh(x) + >>> print(f'{result.value:.4f} in [{result.min:.4f}, {result.max:.4f}]') + 0.4621 in [0.0000, 0.9640] + + Without bounds the result is limited to [-1, 1]: + + >>> x = es.Parameter('x', 0.5) + >>> result = es.math.tanh(x) + >>> print(f'{result.value:.4f} in [{result.min:.4f}, {result.max:.4f}]') + 0.4621 in [-1.0000, 1.0000] + """ # noqa: DOC502, DOC203, DOC105 -- see _apply + return _apply(_TANH, x) + + +def arcsinh(x: _D) -> _D: + """ + Calculate the inverse hyperbolic sine of a dimensionless parameter. + + Parameters + ---------- + x : _D + Dimensionless input. + + Returns + ------- + _D + Parameter or DescriptorNumber + ``arcsinh(x)`` with propagated variance and bounds. The class + follows the inputs: a ``Parameter`` (with bounds) if any + argument is one, otherwise a ``DescriptorNumber`` (without + bounds). + + Raises + ------ + NotImplementedError + If ``x`` is not a ``Parameter`` or ``DescriptorNumber``. + UnitError + If ``x`` is not dimensionless. + + Examples + -------- + + Inverse hyperbolic sine, with its bounds: + + >>> import easyscience as es + >>> x = es.Parameter('x', 1.0, min=0.0, max=10.0) + >>> result = es.math.arcsinh(x) + >>> print(f'{result.value:.4f} in [{result.min:.4f}, {result.max:.4f}]') + 0.8814 in [0.0000, 2.9982] + + The variance is propagated (first order): + + >>> x = es.Parameter('x', 0.0, variance=0.01) + >>> result = es.math.arcsinh(x) + >>> print(f'{result.value:.4f}, variance {result.variance:.4f}') + 0.0000, variance 0.0100 + """ # noqa: DOC502, DOC203, DOC105 -- see _apply + return _apply(_ARCSINH, x) + + +def arccosh(x: _D) -> _D: + """ + Calculate the inverse hyperbolic cosine of a dimensionless + parameter. + + Parameters + ---------- + x : _D + Dimensionless input; its value must be at least 1. + + Returns + ------- + _D + Parameter or DescriptorNumber + ``arccosh(x)`` (non-negative) with propagated variance and + bounds. The class follows the inputs: a ``Parameter`` (with + bounds) if any argument is one, otherwise a ``DescriptorNumber`` + (without bounds). + + Raises + ------ + NotImplementedError + If ``x`` is not a ``Parameter`` or ``DescriptorNumber``. + UnitError + If ``x`` is not dimensionless. + ValueError + If the value of ``x`` is below 1, or if the result is not finite + (a non-zero variance at ``x == 1``, where the derivative is + infinite). + + Notes + ----- + Bounds below 1 are clipped to the domain. + + Examples + -------- + + Inverse hyperbolic cosine, with its bounds: + + >>> import easyscience as es + >>> x = es.Parameter('x', 2.0, min=1.0, max=10.0) + >>> result = es.math.arccosh(x) + >>> print(f'{result.value:.4f} in [{result.min:.4f}, {result.max:.4f}]') + 1.3170 in [0.0000, 2.9932] + + Bounds below 1 are clipped to the domain: + + >>> x = es.Parameter('x', 2.0, min=-5.0, max=10.0) + >>> result = es.math.arccosh(x) + >>> print(f'bounds [{result.min:.4f}, {result.max:.4f}]') + bounds [0.0000, 2.9932] + """ # noqa: DOC502, DOC203, DOC105 -- see _apply + return _apply(_ARCCOSH, x) + + +def arctanh(x: _D) -> _D: + """ + Calculate the inverse hyperbolic tangent of a dimensionless + parameter. + + Parameters + ---------- + x : _D + Dimensionless input; its value must be strictly between -1 and + 1. + + Returns + ------- + _D + Parameter or DescriptorNumber + ``arctanh(x)`` with propagated variance and bounds. The class + follows the inputs: a ``Parameter`` (with bounds) if any + argument is one, otherwise a ``DescriptorNumber`` (without + bounds). + + Raises + ------ + NotImplementedError + If ``x`` is not a ``Parameter`` or ``DescriptorNumber``. + UnitError + If ``x`` is not dimensionless. + ValueError + If the value of ``x`` is not strictly inside ``(-1, 1)``. + + Notes + ----- + Bounds beyond ±1 are clipped to the domain, which makes the result + bound infinite. + + Examples + -------- + + Inverse hyperbolic tangent, with its bounds: + + >>> import easyscience as es + >>> x = es.Parameter('x', 0.5, min=0.0, max=0.9) + >>> result = es.math.arctanh(x) + >>> print(f'{result.value:.4f} in [{result.min:.4f}, {result.max:.4f}]') + 0.5493 in [0.0000, 1.4722] + + Bounds beyond ±1 are clipped, which makes the result bounds + infinite: + + >>> x = es.Parameter('x', 0.5, min=-2.0, max=2.0) + >>> result = es.math.arctanh(x) + >>> print(f'bounds [{result.min:.4f}, {result.max:.4f}]') + bounds [-inf, inf] + """ # noqa: DOC502, DOC203, DOC105 -- see _apply + return _apply(_ARCTANH, x) + + +# ========================================== +# Aliases +# ========================================== +# The C / Scipp spellings of the inverse functions. They are the same function objects, so +# they share the docstring, and the parameters they return are named after the canonical +# function (``asin(x)`` returns a parameter called ``arcsin(x)``). + +asin = arcsin +acos = arccos +atan = arctan +atan2 = arctan2 +asinh = arcsinh +acosh = arccosh +atanh = arctanh + + +# ========================================== +# Specification of a maths operation +# ========================================== + + +@dataclass(frozen=True) +class _Domain: + """ + The set of argument values a function accepts: an interval, open or + closed at each end. + """ + + lo: float = -np.inf + hi: float = np.inf + open_lo: bool = False + open_hi: bool = False + + def contains(self, value: float) -> bool: + above = value > self.lo if self.open_lo else value >= self.lo + below = value < self.hi if self.open_hi else value <= self.hi + return above and below + + def clip(self, lo: float, hi: float) -> _Interval: + """ + Restrict ``[lo, hi]`` to the domain (open ends clip to the edge: + the limit). + """ + return max(lo, self.lo), min(hi, self.hi) + + def __str__(self) -> str: + left = '(' if self.open_lo else '[' + right = ')' if self.open_hi else ']' + return f'{left}{self.lo}, {self.hi}{right}' + + +@dataclass(frozen=True) +class _MathOp: + """ + Everything that makes one maths function different from another. + + Attributes + ---------- + name : str + Name used in error messages and in the name of the resulting + ``Parameter``. + sc_func : Callable + The Scipp function (it validates units and propagates variance, + if it can). + bounds : _BoundsStrategy + Maps one ``(lo, hi)`` interval per argument to the ``(min, + max)`` of the result. + domain : _Domain, default=_Domain() + Accepted argument values. Applies to every argument. + native_unit : str | None, default=None + If set, bounds are converted to this unit before ``bounds`` runs + (``'rad'`` for the trigonometric functions, whose numpy versions + assume radians). + jacobian : Callable | None, default=None + Partial derivatives of the function, ``jacobian(*values) -> + (df/dx1, df/dx2, ...)``. Set it for functions whose variance + Scipp cannot propagate; EasyScience then applies the delta + method itself. + guard : Callable[..., None] | None, default=None + Extra check of the arguments' values, ``guard(*params)``, + raising ``ValueError`` where the delta method breaks down (tan + asymptotes, the arctan2 origin). + """ + + name: str + sc_func: Callable + bounds: _BoundsStrategy + domain: _Domain = _Domain() + native_unit: str | None = None + jacobian: Callable | None = None + guard: Callable[..., None] | None = None + + +# ========================================== +# The dispatcher +# ========================================== + + +# NOTE(): Every result is an independent, non-fixed ``Parameter`` that is a snapshot of its inputs +# (as are the results of ``Parameter.__mul__`` and friends): it does not follow later changes and, +# if handed to a ``Fitter``, would be varied as a free parameter of its own. Whether results should +# be ``fixed=True``, or linked to their inputs (compare ``Parameter.from_dependency``), is for the +# lead developer to decide. +# TODO(): Cost. A chain of five es.math / Parameter operations (a Gaussian width derived from a +# FWHM inside a model's __call__) took about 1.3 ms per evaluation, constant over thousands of +# evaluations and without leaking. That is fine for ordinary fits, but it has not been profiled +# where the time goes: profile before using es.math in hot loops (tens of thousands of +# evaluations), and consider a cheaper path if it matters. +def _apply(op: _MathOp, *params: DescriptorNumber) -> DescriptorNumber: + """ + Evaluate ``op`` on ``params``: validate, call Scipp, derive bounds, + repack. + + Scipp validates the units, so a wrong unit raises its ``UnitError`` + before any of the checks below. + + Parameters + ---------- + op : _MathOp + Specification of the function. + *params : DescriptorNumber + The arguments (a ``Parameter`` is a ``DescriptorNumber``). + + Returns + ------- + DescriptorNumber + A new ``Parameter`` named ``'()'`` if any + argument is a ``Parameter``, otherwise a new + ``DescriptorNumber`` with that name. + + Raises + ------ + NotImplementedError + If an argument is not a ``Parameter`` or ``DescriptorNumber``. + ValueError + If an argument's value is outside the domain, or the result is + not finite. + """ + for param in params: + if not isinstance(param, DescriptorNumber): + raise NotImplementedError( + 'EasyScience math functions currently only support `Parameter` and ' + '`DescriptorNumber` objects. Raw numbers and arrays are not yet supported. ' + f'Found: {type(param)}' + ) + + # 1. Scipp first: it owns the unit rules, so a unit error is reported before anything else. + manual_variance = op.jacobian is not None + sc_args = [_to_scipp_scalar(param, strip_variance=manual_variance) for param in params] + sc_result = op.sc_func(*sc_args) + + # 2. Values are validated strictly: a bad value would break a fit loop. + for param in params: + if not op.domain.contains(param.value): + raise ValueError( + f"The value of '{param.name}' ({param.value}) is outside the domain " + f'{op.domain} of {op.name}().' + ) + if op.guard is not None: + op.guard(*params) + + # 3. Variance for functions Scipp cannot handle: the delta method. + if manual_variance: + variance = _delta_method_variance(op, params) + if variance is not None: + sc_result.variance = variance + + label = ', '.join(param.name for param in params) + if not _is_finite(sc_result.value) or not _is_finite(sc_result.variance): + raise ValueError( + f'{op.name}({label}) is not finite: value={sc_result.value}, ' + f'variance={sc_result.variance}. Check the value and variance of the argument(s).' + ) + + # A plain DescriptorNumber has no bounds, so neither has a result made only from those. + if not any(isinstance(param, Parameter) for param in params): + return DescriptorNumber.from_scipp(name=f'{op.name}({label})', full_value=sc_result) + + # 4. Bounds are permissive: clip to the domain, convert, and let the strategy decide. + intervals = [_native_interval(op, param) for param in params] + with np.errstate(all='ignore'): + new_min, new_max = op.bounds(*intervals) + # Guard against round-off at an extremum: the value is always inside its own bounds. + new_min, new_max = min(new_min, sc_result.value), max(new_max, sc_result.value) + # Parameter rejects (nearly) identical bounds, so bounds that collapsed are padded. + new_min, new_max = _pad_if_collapsed(new_min, new_max, sc_result.value) + + return Parameter.from_scipp( + name=f'{op.name}({label})', + full_value=sc_result, + min=new_min, + max=new_max, + ) + + +# TODO(): ``rtol=1e-9`` below mirrors the literal tolerance in the min/max setters of ``Parameter``, +# which is not a named constant. Expose it in parameter.py and import it here, so that the two +# cannot drift apart (if ``Parameter`` tightens its tolerance this padding would silently stop +# preventing the 'identical bounds' ValueError). +def _pad_if_collapsed(lo: float, hi: float, value: float) -> tuple[float, float]: + """ + Widen bounds that ``Parameter`` would reject as identical, + symmetrically around ``value``. + + Parameters + ---------- + lo : float + Lower bound. + hi : float + Upper bound. + value : float + The result value, which lies inside ``[lo, hi]``. + + Returns + ------- + tuple[float, float] + ``(lo, hi)`` unchanged, or padded by ``_COLLAPSE_PAD * |value|`` + (``_COLLAPSE_PAD_ZERO`` if the value is 0) on each side if they + are equal or within ``Parameter``'s tolerance. + """ + if not np.isclose(lo, hi, rtol=1e-9, atol=0.0): + return lo, hi + pad = abs(value) * _COLLAPSE_PAD or _COLLAPSE_PAD_ZERO + return lo - pad, hi + pad + + +def _angle_in(result: DescriptorNumber, unit: str) -> DescriptorNumber: + """ + Express an angle ``result`` (in radians, as Scipp returns it) in + ``unit``. + + The value, the variance and the bounds are converted together by + ``Parameter.convert_unit``. + + Parameters + ---------- + result : DescriptorNumber + An angle in radians. It is converted in place. + unit : str + Target unit, ``'rad'`` or ``'deg'``. + + Returns + ------- + DescriptorNumber + ``result``, in ``unit``. + + Raises + ------ + UnitError + If ``unit`` is not an angle unit. + """ # noqa: DOC502 -- raised by Parameter.convert_unit + if unit != result.unit: + result.convert_unit(unit) + return result + + +def _variance(param: DescriptorNumber) -> float: + """ + Variance of a parameter; a missing (``None``) or zero variance + counts as 0.0. + """ + return float(param.variance) if param.variance else 0.0 + + +def _to_scipp_scalar(param: DescriptorNumber, strip_variance: bool = False) -> sc.Variable: + """ + Convert a Parameter to a Scipp scalar; a zero or missing variance is + dropped. + """ + kwargs = {'value': float(param.value), 'unit': param.unit} + variance = _variance(param) + if not strip_variance and variance > 0: + kwargs['variance'] = variance + return sc.scalar(**kwargs) + + +def _delta_method_variance(op: _MathOp, params: tuple[DescriptorNumber, ...]) -> float | None: + """ + First-order error propagation, ``var(y) = sum_i (dy/dx_i)**2 * + var(x_i)``. + + Only arguments that have a variance contribute (this avoids ``0 * + inf``). Returns ``None`` if no argument has a variance. + """ + variances = [_variance(p) for p in params] + if not any(variances): + return None + with np.errstate(all='ignore'): + partials = op.jacobian(*(np.float64(p.value) for p in params)) + return float(sum(j**2 * v for j, v in zip(partials, variances, strict=True) if v > 0)) + + +def _is_finite(value: float | None) -> bool: + """ + Whether ``value`` is finite; a missing value (no variance) counts as + fine. + """ + return value is None or bool(np.isfinite(value)) + + +@cache +def _native_factor(unit: str, native_unit: str | None) -> float: + """ + Factor converting ``unit`` to ``native_unit`` (1.0 if none is + needed), asked of Scipp once. + """ + if native_unit is None or unit == native_unit: + return 1.0 + return float(sc.scalar(1.0, unit=unit).to(unit=native_unit).value) + + +def _native_interval(op: _MathOp, param: DescriptorNumber) -> _Interval: + """ + The parameter's bounds, clipped to the domain and in the unit the + strategy expects. + """ + if isinstance(param, Parameter): + lo, hi = float(param.min), float(param.max) + else: # a plain DescriptorNumber has no bounds: it is the single point at its value + lo = hi = float(param.value) + lo, hi = op.domain.clip(lo, hi) + factor = _native_factor(str(param.unit), op.native_unit) + return lo * factor, hi * factor + + +# ========================================== +# Bounds strategies +# ========================================== +# Each takes one (lo, hi) interval per argument and returns the exact (min, max) of the +# function over that interval. They run inside ``np.errstate(all='ignore')``, so an edge +# such as log(0) = -inf is returned instead of warning. + + +def _monotonic(func: Callable) -> _BoundsStrategy: + """ + Monotonic function (increasing or decreasing): the extremes sit at + the endpoints. + """ + + def bounds(interval: _Interval) -> tuple[float, float]: + lo, hi = interval + at_lo, at_hi = float(func(lo)), float(func(hi)) + return min(at_lo, at_hi), max(at_lo, at_hi) + + return bounds + + +def _extrema( + func: Callable, critical_points: Callable[[float, float], np.ndarray] +) -> _BoundsStrategy: + """ + Non-monotonic function: the extremes sit at the endpoints or at + critical points. + + ``critical_points(lo, hi)`` returns the points inside ``[lo, hi]`` + where the function may have a local extremum. ``nan`` candidates, + such as ``sin(inf)``, are ignored. + """ + + def bounds(interval: _Interval) -> tuple[float, float]: + lo, hi = interval + candidates = np.concatenate(([lo, hi], critical_points(lo, hi))) + values = np.asarray(func(candidates), dtype=float) + return float(np.nanmin(values)), float(np.nanmax(values)) + + return bounds + + +def _periodic_critical_points(offset: float) -> Callable[[float, float], np.ndarray]: + """ + Critical points ``offset + n*pi`` (n an integer) of sin (``pi/2``) + or cos (``0``). + + A whole period, or an infinite interval, holds a maximum and a + minimum, so those two are returned without enumerating. + """ + + def points(lo: float, hi: float) -> np.ndarray: + if not (np.isfinite(lo) and np.isfinite(hi)) or hi - lo >= _TWO_PI: + return np.array([offset, offset + np.pi]) + n = np.arange(np.ceil((lo - offset) / np.pi), np.floor((hi - offset) / np.pi) + 1) + return offset + n * np.pi + + return points + + +def _cosh_critical_points(lo: float, hi: float) -> np.ndarray: + """The single minimum of cosh, at 0.""" + return np.array([0.0]) if lo <= 0.0 <= hi else np.array([]) + + +def _unnormalized_sinc(x) -> np.ndarray: + """ + sin(x)/x, with sinc(0) = 1 and the limit sinc(±inf) = 0. + + ``np.sinc`` is the *normalised* sin(pi x)/(pi x), so the argument is + scaled by 1/pi. + """ + x = np.asarray(x, dtype=float) + with np.errstate(all='ignore'): + return np.where(np.isinf(x), 0.0, np.sinc(x / np.pi)) + + +def _sinc_positive_extrema(lo: float, hi: float, count: int = _SINC_ROOT_COUNT) -> np.ndarray: + """ + Stationary points of sinc on the positive axis inside ``[lo, hi]``, + nearest to ``lo``. + + They are the positive roots of ``tan(x) = x``, one in each ``(k*pi, + k*pi + pi/2)``. The heights of the lobes fall off like 1/x, so the + largest maximum and the deepest minimum in the interval are among + the first few roots at or after ``lo``. + + Beyond ``_SINC_NEGLIGIBLE_LOBES_FROM`` the lobes are negligible and + none are returned. + + All ``count`` roots are refined at once with Newton's method on + ``g(x) = x cos(x) - sin(x)``, ``g'(x) = -x sin(x)``, with the + vectorised update ``x <- x - g(x) / g'(x)`` (no Python loop over the + roots). + """ + if hi <= 0.0 or lo > _SINC_NEGLIGIBLE_LOBES_FROM: + return np.array([]) + first = max(1, int(np.floor(max(lo, 0.0) / np.pi))) + k = np.arange(first, first + count) + x = (k + 0.5) * np.pi + x = x - 1.0 / x # asymptotic first guess for the root of tan(x) = x + for _ in range(_SINC_NEWTON_STEPS): + x = x - (x * np.cos(x) - np.sin(x)) / (-x * np.sin(x)) + return x[(x >= lo) & (x <= hi)] + + +def _sinc_critical_points(lo: float, hi: float) -> np.ndarray: + """ + Stationary points of sinc: 0 (the global maximum) and the roots of + tan(x) = x on both sides. + """ + points = [_sinc_positive_extrema(lo, hi), -_sinc_positive_extrema(-hi, -lo)] + if lo <= 0.0 <= hi: + points.append(np.array([0.0])) + return np.concatenate(points) + + +def _tan_branch(position: float) -> tuple[float, bool]: + """ + Branch of tan that a point lies on, and whether it sits exactly on + an asymptote. + + Parameters + ---------- + position : float + ``(x - pi/2) / pi``, so that the asymptotes ``pi/2 + n*pi`` are + the integers. + + Returns + ------- + tuple[float, bool] + The branch index (the asymptote's own index if the point sits on + one) and whether the point sits on an asymptote. + """ + nearest = float(np.round(position)) + on_asymptote = abs(position - nearest) <= _TAN_EDGE_TOLERANCE + return (nearest if on_asymptote else float(np.floor(position))), on_asymptote + + +def _tan_bounds(interval: _Interval) -> tuple[float, float]: + """ + Bounds of tan(x) for ``x`` in radians. + + tan is increasing between asymptotes (``pi/2 + n*pi``). An interval + that contains an asymptote in its interior covers all real numbers. + A bound exactly on an asymptote makes the result unbounded on that + side only: ``[0, pi/2]`` gives ``[0, inf)``. + """ + lo, hi = interval + if not (np.isfinite(lo) and np.isfinite(hi)): + return -np.inf, np.inf + branch_lo, lo_on_asymptote = _tan_branch((lo - _HALF_PI) / np.pi) + branch_hi, hi_on_asymptote = _tan_branch((hi - _HALF_PI) / np.pi) + if hi_on_asymptote: + branch_hi -= 1 # an upper bound on an asymptote closes the branch below it + if branch_lo != branch_hi: + return -np.inf, np.inf + lower = -np.inf if lo_on_asymptote else float(np.tan(lo)) + upper = np.inf if hi_on_asymptote else float(np.tan(hi)) + return lower, upper + + +def _std(param: DescriptorNumber) -> float: + """Standard deviation of a parameter; 0 if it has no variance.""" + return math.sqrt(_variance(param)) + + +# NOTE(): Design decision awaiting sign-off. The guards below reject a +# value that lies within one standard deviation of a tan asymptote (or of the arctan2 origin), +# because the delta method is meaningless there. The price is that a legitimate fit result that is +# consistent with the singularity (for example a = 0.1 +- 0.2 and b = 0.1 +- 0.2 in +# arctan2(b, a)) cannot be turned into a derived quantity. Alternatives: only warn, or return a +# very large variance instead of raising. +def _tan_guard(param: DescriptorNumber) -> None: + """ + The value must keep clear of a tan asymptote, where the variance + would explode. + + The delta method is only meaningful if the asymptote lies more than + one standard deviation away, so that distance (never less than + ``_TAN_ASYMPTOTE_TOLERANCE``) is the limit. + """ + factor = _native_factor(str(param.unit), 'rad') + value = float(param.value) * factor + margin = max(_TAN_ASYMPTOTE_TOLERANCE, _std(param) * factor) + distance = abs(value % np.pi - _HALF_PI) + if distance <= margin: + raise ValueError( + f"The value of '{param.name}' ({param.value} {param.unit}) is within " + f'{margin:.3g} rad of a tan asymptote (pi/2 + n*pi); the propagated variance ' + f'would explode.' + ) + + +# NOTE(): The same one-standard-deviation design decision as for _tan_guard applies here. +def _arctan2_guard(y: DescriptorNumber, x: DescriptorNumber) -> None: + """ + The point ``(x, y)`` must keep clear of the origin, where the angle + is undefined. + + Only checked if an argument has a variance: the delta method then + divides by ``x**2 + y**2``. The point must lie more than one + standard deviation from the origin. + """ + margin = max(_std(y), _std(x)) + if margin > 0 and np.hypot(float(y.value), float(x.value)) <= margin: + raise ValueError( + f"The point ('{x.name}', '{y.name}') = ({x.value}, {y.value}) is within one standard " + f'deviation ({margin:.3g}) of the origin, where the angle of arctan2() is undefined ' + f'and its propagated variance would explode.' + ) + + +def _arctan2_bounds(y_interval: _Interval, x_interval: _Interval) -> tuple[float, float]: + """ + Bounds of arctan2(y, x) over the box ``y_interval x x_interval``. + + The angle is a harmonic function, so away from the branch cut its + extremes sit on the corners of the box. If the box touches the cut + (``y = 0`` with ``x < 0``) the angle jumps between -pi and pi and + the range is the full ``[-pi, pi]``. + """ + (y_lo, y_hi), (x_lo, x_hi) = y_interval, x_interval + if y_lo < 0.0 <= y_hi and x_lo < 0.0: + if np.isfinite([y_lo, y_hi, x_lo, x_hi]).all(): + # stacklevel=4: this function -> _apply -> arctan2 -> the user's call. + warnings.warn( + 'The bounds of the arguments touch the arctan2 branch cut (y = 0, x < 0), where ' + 'the angle jumps between -pi and pi. Gradient-based minimisers behave badly ' + 'across such a jump.', + UserWarning, + stacklevel=4, + ) + return -np.pi, np.pi + corners = [np.arctan2(y, x) for y in (y_lo, y_hi) for x in (x_lo, x_hi)] + return float(min(corners)), float(max(corners)) + + +# ========================================== +# Jacobians (for variance Scipp cannot propagate) +# ========================================== + + +# NumPy scalars (not Python floats) are passed in: 0/0 then gives nan under ``np.errstate`` +# and is caught by the finiteness check, instead of raising ``ZeroDivisionError``. + + +def _arctan2_jacobian(y: np.float64, x: np.float64) -> tuple[np.float64, np.float64]: + r2 = x**2 + y**2 + return x / r2, -y / r2 + + +def _erf_jacobian(x: np.float64) -> tuple[np.float64]: + return (2.0 / np.sqrt(np.pi) * np.exp(-(x**2)),) + + +def _erfc_jacobian(x: np.float64) -> tuple[np.float64]: + return (-2.0 / np.sqrt(np.pi) * np.exp(-(x**2)),) + + +# ========================================== +# The specifications +# ========================================== + +_POSITIVE = _Domain(lo=0.0, open_lo=True) +_NON_NEGATIVE = _Domain(lo=0.0) +_AT_LEAST_ONE = _Domain(lo=1.0) +_CLOSED_UNIT_INTERVAL = _Domain(lo=-1.0, hi=1.0) +_OPEN_UNIT_INTERVAL = _Domain(lo=-1.0, hi=1.0, open_lo=True, open_hi=True) + +_EXP = _MathOp('exp', sc.exp, _monotonic(np.exp)) +_LOG = _MathOp('log', sc.log, _monotonic(np.log), domain=_POSITIVE) +_LOG10 = _MathOp('log10', sc.log10, _monotonic(np.log10), domain=_POSITIVE) +_SQRT = _MathOp('sqrt', sc.sqrt, _monotonic(np.sqrt), domain=_NON_NEGATIVE) +_ERF = _MathOp('erf', sc.erf, _monotonic(math.erf), jacobian=_erf_jacobian) +_ERFC = _MathOp('erfc', sc.erfc, _monotonic(math.erfc), jacobian=_erfc_jacobian) + +_SIN = _MathOp( + 'sin', sc.sin, _extrema(np.sin, _periodic_critical_points(_HALF_PI)), native_unit='rad' +) +_COS = _MathOp('cos', sc.cos, _extrema(np.cos, _periodic_critical_points(0.0)), native_unit='rad') +_TAN = _MathOp('tan', sc.tan, _tan_bounds, native_unit='rad', guard=_tan_guard) +_ARCSIN = _MathOp('arcsin', sc.asin, _monotonic(np.arcsin), domain=_CLOSED_UNIT_INTERVAL) +_ARCCOS = _MathOp('arccos', sc.acos, _monotonic(np.arccos), domain=_CLOSED_UNIT_INTERVAL) +_ARCTAN = _MathOp('arctan', sc.atan, _monotonic(np.arctan)) +_ARCTAN2 = _MathOp( + 'arctan2', + lambda y, x: sc.atan2(y=y, x=x), # Scipp's atan2 takes keyword-only arguments + _arctan2_bounds, + jacobian=_arctan2_jacobian, + guard=_arctan2_guard, +) +_SINC = _MathOp('sinc', sc.sinc, _extrema(_unnormalized_sinc, _sinc_critical_points)) + +_SINH = _MathOp('sinh', sc.sinh, _monotonic(np.sinh)) +_COSH = _MathOp('cosh', sc.cosh, _extrema(np.cosh, _cosh_critical_points)) +_TANH = _MathOp('tanh', sc.tanh, _monotonic(np.tanh)) +_ARCSINH = _MathOp('arcsinh', sc.asinh, _monotonic(np.arcsinh)) +_ARCCOSH = _MathOp('arccosh', sc.acosh, _monotonic(np.arccosh), domain=_AT_LEAST_ONE) +_ARCTANH = _MathOp('arctanh', sc.atanh, _monotonic(np.arctanh), domain=_OPEN_UNIT_INTERVAL) diff --git a/tests/integration/fitting/test_math_in_fit.py b/tests/integration/fitting/test_math_in_fit.py new file mode 100644 index 00000000..d2bd11f0 --- /dev/null +++ b/tests/integration/fitting/test_math_in_fit.py @@ -0,0 +1,680 @@ +# SPDX-FileCopyrightText: 2026 EasyScience contributors +# SPDX-License-Identifier: BSD-3-Clause + +""" +``es.math`` applied to the parameters of a fit. + +``es.math`` takes scalar ``Parameter`` objects, so it cannot be applied over an ``x`` array +inside a model. Its use is to derive quantities, with uncertainties, from the parameters a fit +produces. These tests fit small models and check the derived value and variance against closed +forms and against sampling the fitted distribution. + +The models are test-only, like ``AbsSin`` in ``test_fitter.py``: + +* ``StraightLine``: dimensionless parameters, so every dimensionless function can be applied. +* ``AngleLine``: a free parameter in degrees, for the trigonometric functions and ``sinc``. +* ``Quadrature``: ``a cos + b sin`` with uncorrelated ``a`` and ``b``, for amplitude and phase. +* ``GaussianPeak``: parametrised by its FWHM; the width is derived with ``log`` and ``sqrt`` + *inside* ``__call__`` on every evaluation (so a ``Parameter`` is created per call), and an + in-window fraction is derived after the fit with ``erf``. +""" + +import gc +import math +import time + +import numpy as np +import pytest + +import easyscience as es +from easyscience import AvailableMinimizers +from easyscience import Fitter +from easyscience import Parameter +from easyscience import global_object +from easyscience.base_classes import ModelBase + +pytestmark = pytest.mark.fast + +SIGMA = 0.05 # standard deviation of the noise added to the synthetic data +OMEGA = 2.0 * np.pi # angular frequency of the Quadrature model (period 1) + +_erf = np.vectorize(math.erf) +_erfc = np.vectorize(math.erfc) + + +@pytest.fixture(autouse=True) +def clear(): + global_object.map._clear() + + +# ---------------------------------------------------------------------- +# Test-only models +# ---------------------------------------------------------------------- + + +class StraightLine(ModelBase): + """y = slope * x + intercept""" + + def __init__(self, slope: float, intercept: float): + super().__init__() + self._slope = Parameter('slope', slope) + self._intercept = Parameter('intercept', intercept) + + @property + def slope(self) -> Parameter: + return self._slope + + @slope.setter + def slope(self, value: float) -> None: + self._slope.value = value + + @property + def intercept(self) -> Parameter: + return self._intercept + + @intercept.setter + def intercept(self, value: float) -> None: + self._intercept.value = value + + def __call__(self, x: np.ndarray) -> np.ndarray: + return self.slope.value * x + self.intercept.value + + +class AngleLine(ModelBase): + """y = tan(angle) * x + intercept, with the angle a free parameter in degrees.""" + + def __init__(self, angle: float, intercept: float): + super().__init__() + self._angle = Parameter('angle', angle, unit='deg', min=1.0, max=80.0) + self._intercept = Parameter('intercept', intercept) + + @property + def angle(self) -> Parameter: + return self._angle + + @angle.setter + def angle(self, value: float) -> None: + self._angle.value = value + + @property + def intercept(self) -> Parameter: + return self._intercept + + @intercept.setter + def intercept(self, value: float) -> None: + self._intercept.value = value + + def __call__(self, x: np.ndarray) -> np.ndarray: + return np.tan(np.deg2rad(self.angle.value)) * x + self.intercept.value + + +class Quadrature(ModelBase): + """y = a cos(OMEGA x) + b sin(OMEGA x)""" + + def __init__(self, a: float, b: float): + super().__init__() + self._a = Parameter('a', a) + self._b = Parameter('b', b) + + @property + def a(self) -> Parameter: + return self._a + + @a.setter + def a(self, value: float) -> None: + self._a.value = value + + @property + def b(self) -> Parameter: + return self._b + + @b.setter + def b(self, value: float) -> None: + self._b.value = value + + def __call__(self, x: np.ndarray) -> np.ndarray: + return self.a.value * np.cos(OMEGA * x) + self.b.value * np.sin(OMEGA * x) + + +class GaussianPeak(ModelBase): + """ + y = amplitude * exp(-(x - center)**2 / (2 sigma**2)), with sigma derived from the FWHM. + + ``sigma()`` uses ``es.math`` (``log``, ``sqrt``) and so creates new ``Parameter`` objects on + every call; ``__call__`` calls it on every evaluation. + """ + + def __init__(self, amplitude: float, center: float, fwhm: float): + super().__init__() + self._amplitude = Parameter('amplitude', amplitude) + self._center = Parameter('center', center) + self._fwhm = Parameter('fwhm', fwhm, min=0.1, max=10.0) + self._two = Parameter('two', 2.0) # a constant, so that log() receives a Parameter + + @property + def amplitude(self) -> Parameter: + return self._amplitude + + @amplitude.setter + def amplitude(self, value: float) -> None: + self._amplitude.value = value + + @property + def center(self) -> Parameter: + return self._center + + @center.setter + def center(self, value: float) -> None: + self._center.value = value + + @property + def fwhm(self) -> Parameter: + return self._fwhm + + @fwhm.setter + def fwhm(self, value: float) -> None: + self._fwhm.value = value + + def sigma(self) -> Parameter: + """sigma = FWHM / (2 sqrt(2 ln 2)), with value and variance propagated.""" + return self._fwhm / (2 * es.math.sqrt(2 * es.math.log(self._two))) + + def fraction_within(self, window: float) -> Parameter: + """Fraction of the area within +-window of the center: erf(window / (sqrt(2) sigma)).""" + return es.math.erf(window / (self.sigma() * math.sqrt(2.0))) + + def __call__(self, x: np.ndarray) -> np.ndarray: + sigma = self.sigma().value + return self.amplitude.value * np.exp(-((x - self.center.value) ** 2) / (2 * sigma**2)) + + +FWHM_PER_SIGMA = 2.0 * math.sqrt(2.0 * math.log(2.0)) + + +def fit(model: ModelBase, x: np.ndarray, y: np.ndarray) -> None: + """Fit ``model`` to ``y`` with weights = 1/sigma (the repo's convention).""" + result = Fitter(model, model).fit(x, y, weights=np.full_like(x, 1.0 / SIGMA)) + assert result.success + + +def noisy(values: np.ndarray, seed: int) -> np.ndarray: + return values + np.random.default_rng(seed).normal(0.0, SIGMA, values.size) + + +# ---------------------------------------------------------------------- +# Dimensionless functions on the parameters of a fitted line +# ---------------------------------------------------------------------- + +# (function, parameter it is applied to, reference, derivative of the reference). +# slope is about 0.5 (inside (-1, 1)) and intercept about 2.0 (above 1), so every function +# is evaluated well inside its domain. +UNARY_CASES = [ + ('exp', 'slope', np.exp, np.exp), + ('erf', 'slope', _erf, lambda x: 2 / np.sqrt(np.pi) * np.exp(-(x**2))), + ('erfc', 'slope', _erfc, lambda x: -2 / np.sqrt(np.pi) * np.exp(-(x**2))), + ('sinh', 'slope', np.sinh, np.cosh), + ('cosh', 'slope', np.cosh, np.sinh), + ('tanh', 'slope', np.tanh, lambda x: 1 / np.cosh(x) ** 2), + ('arcsinh', 'slope', np.arcsinh, lambda x: 1 / np.sqrt(x**2 + 1)), + ('arcsin', 'slope', np.arcsin, lambda x: 1 / np.sqrt(1 - x**2)), + ('arccos', 'slope', np.arccos, lambda x: -1 / np.sqrt(1 - x**2)), + ('arctan', 'slope', np.arctan, lambda x: 1 / (1 + x**2)), + ('arctanh', 'slope', np.arctanh, lambda x: 1 / (1 - x**2)), + ('log', 'intercept', np.log, lambda x: 1 / x), + ('log10', 'intercept', np.log10, lambda x: 1 / (x * np.log(10))), + ('sqrt', 'intercept', np.sqrt, lambda x: 0.5 / np.sqrt(x)), + ('arccosh', 'intercept', np.arccosh, lambda x: 1 / np.sqrt(x**2 - 1)), +] +UNARY_IDS = [case[0] for case in UNARY_CASES] + + +@pytest.fixture +def fitted_line(): + x = np.linspace(0.0, 5.0, 200) + model = StraightLine(slope=0.1, intercept=1.0) + fit(model, x, noisy(0.5 * x + 2.0, seed=1)) + return model, x + + +class TestDimensionlessFunctionsOnFittedLine: + @pytest.mark.parametrize( + ('name', 'attribute', 'reference', 'derivative'), UNARY_CASES, ids=UNARY_IDS + ) + def test_value_and_variance_follow_the_fitted_parameter( + self, fitted_line, name, attribute, reference, derivative + ): + # WHEN + model, _ = fitted_line + fitted = getattr(model, attribute) + + # THEN + derived = getattr(es.math, name)(fitted) + + # EXPECT + assert fitted.variance > 0 # the fit produced an uncertainty to propagate + assert derived.value == pytest.approx(float(reference(fitted.value)), rel=1e-12) + expected_variance = float(derivative(fitted.value)) ** 2 * fitted.variance + assert derived.variance == pytest.approx(expected_variance, rel=1e-9) + + @pytest.mark.parametrize( + ('name', 'attribute', 'reference', 'derivative'), UNARY_CASES, ids=UNARY_IDS + ) + def test_variance_agrees_with_sampling_the_fitted_distribution( + self, fitted_line, name, attribute, reference, derivative + ): + # WHEN + model, _ = fitted_line + fitted = getattr(model, attribute) + rng = np.random.default_rng(7) + samples = rng.normal(fitted.value, np.sqrt(fitted.variance), 200_000) + + # THEN + derived = getattr(es.math, name)(fitted) + + # EXPECT + # An independent check of the first-order (delta) method: the spread of f(samples). + assert np.std(reference(samples)) == pytest.approx(np.sqrt(derived.variance), rel=0.02) + + @pytest.mark.parametrize(('alias', 'canonical'), [('asin', 'arcsin'), ('atanh', 'arctanh')]) + def test_aliases_give_the_same_result_as_the_canonical_function( + self, fitted_line, alias, canonical + ): + # WHEN + model, _ = fitted_line + + # THEN + via_alias = getattr(es.math, alias)(model.slope) + via_canonical = getattr(es.math, canonical)(model.slope) + + # EXPECT + assert (via_alias.value, via_alias.variance) == ( + via_canonical.value, + via_canonical.variance, + ) + + +# ---------------------------------------------------------------------- +# A fitted angle in degrees +# ---------------------------------------------------------------------- + + +@pytest.fixture +def fitted_angle(): + x = np.linspace(0.0, 5.0, 200) + model = AngleLine(angle=20.0, intercept=1.0) + fit(model, x, noisy(np.tan(np.deg2rad(30.0)) * x + 2.0, seed=2)) + return model + + +class TestFittedAngle: + def test_fit_recovers_the_angle_within_its_uncertainty(self, fitted_angle): + # WHEN THEN + angle = fitted_angle.angle + + # EXPECT + assert abs(angle.value - 30.0) < 5 * np.sqrt(angle.variance) + + @pytest.mark.parametrize( + ('name', 'derivative'), + [ + ('sin', np.cos), + ('cos', lambda t: -np.sin(t)), + ('tan', lambda t: 1 / np.cos(t) ** 2), + ], + ) + def test_trigonometric_value_and_variance_with_the_angle_in_degrees( + self, fitted_angle, name, derivative + ): + # WHEN + angle = fitted_angle.angle + radians = np.deg2rad(angle.value) + + # THEN + derived = getattr(es.math, name)(angle) + + # EXPECT + # d f / d(deg) = f'(rad) * pi / 180, so the variance picks up (pi / 180) ** 2. + expected_variance = float(derivative(radians)) ** 2 * angle.variance * np.deg2rad(1.0) ** 2 + assert derived.unit == 'dimensionless' + assert derived.value == pytest.approx(float(getattr(np, name)(radians)), rel=1e-12) + assert derived.variance == pytest.approx(expected_variance, rel=1e-9) + + @pytest.mark.parametrize( + ('name', 'expected_min', 'expected_max'), + [ + ('sin', np.sin(np.deg2rad(1.0)), np.sin(np.deg2rad(80.0))), + ('cos', np.cos(np.deg2rad(80.0)), np.cos(np.deg2rad(1.0))), + ('tan', np.tan(np.deg2rad(1.0)), np.tan(np.deg2rad(80.0))), + ], + ) + def test_bounds_of_the_fitted_angle_are_mapped_through_the_function( + self, fitted_angle, name, expected_min, expected_max + ): + # WHEN + angle = fitted_angle.angle # bounds [1, 80] deg survive the fit + + # THEN + derived = getattr(es.math, name)(angle) + + # EXPECT + assert derived.min == pytest.approx(expected_min, abs=1e-12) + assert derived.max == pytest.approx(expected_max, abs=1e-12) + + def test_sinc_of_the_fitted_angle_after_converting_to_radians(self, fitted_angle): + # WHEN + angle = fitted_angle.angle + angle.convert_unit('rad') # sinc accepts rad or dimensionless, not deg + + # THEN + derived = es.math.sinc(angle) + + # EXPECT + x = angle.value + expected_variance = ((x * np.cos(x) - np.sin(x)) / x**2) ** 2 * angle.variance + assert derived.value == pytest.approx(np.sin(x) / x, rel=1e-12) + assert derived.variance == pytest.approx(expected_variance, rel=1e-9) + + +# ---------------------------------------------------------------------- +# Amplitude and phase of a fitted a cos + b sin +# ---------------------------------------------------------------------- + + +@pytest.fixture +def fitted_quadrature(): + x = np.arange(0.0, 4.0, 0.01) # four whole periods, endpoint excluded + model = Quadrature(a=1.0, b=0.1) + fit(model, x, noisy(1.5 * np.cos(OMEGA * x) + 0.8 * np.sin(OMEGA * x), seed=3)) + return model, x + + +class TestFittedQuadrature: + def test_premise_a_and_b_are_uncorrelated(self, fitted_quadrature): + # WHEN + _, x = fitted_quadrature + design = np.column_stack([np.cos(OMEGA * x), np.sin(OMEGA * x)]) + + # THEN + covariance = np.linalg.inv(design.T @ design) + correlation = covariance[0, 1] / np.sqrt(covariance[0, 0] * covariance[1, 1]) + + # EXPECT + # Over whole periods cos and sin are orthogonal; the delta method relies on this. + assert abs(correlation) < 1e-9 + + def test_amplitude_from_parameter_arithmetic_and_sqrt(self, fitted_quadrature): + # WHEN + model, _ = fitted_quadrature + a, b = model.a, model.b + + # THEN + amplitude = es.math.sqrt(a**2 + b**2) + + # EXPECT + expected_value = np.hypot(a.value, b.value) + expected_variance = (a.value**2 * a.variance + b.value**2 * b.variance) / expected_value**2 + assert amplitude.value == pytest.approx(expected_value, rel=1e-12) + assert amplitude.variance == pytest.approx(expected_variance, rel=1e-9) + assert abs(amplitude.value - 1.7) < 5 * np.sqrt(amplitude.variance) + + def test_phase_from_arctan2(self, fitted_quadrature): + # WHEN + model, _ = fitted_quadrature + a, b = model.a, model.b + + # THEN + phase = es.math.arctan2(b, a) + + # EXPECT + r2 = a.value**2 + b.value**2 + expected_variance = (a.value**2 * b.variance + b.value**2 * a.variance) / r2**2 + assert phase.unit == 'rad' + assert phase.value == pytest.approx(np.arctan2(b.value, a.value), rel=1e-12) + assert phase.variance == pytest.approx(expected_variance, rel=1e-9) + assert abs(phase.value - np.arctan2(0.8, 1.5)) < 5 * np.sqrt(phase.variance) + + def test_phase_variance_agrees_with_sampling_uncorrelated_parameters(self, fitted_quadrature): + # WHEN + model, _ = fitted_quadrature + a, b = model.a, model.b + rng = np.random.default_rng(11) + sampled_a = rng.normal(a.value, np.sqrt(a.variance), 200_000) + sampled_b = rng.normal(b.value, np.sqrt(b.variance), 200_000) + + # THEN + phase = es.math.arctan2(b, a) + + # EXPECT + assert np.std(np.arctan2(sampled_b, sampled_a)) == pytest.approx( + np.sqrt(phase.variance), rel=0.03 + ) + + def test_phase_in_degrees_scales_value_and_variance(self, fitted_quadrature): + # WHEN + model, _ = fitted_quadrature + a, b = model.a, model.b + + # THEN + in_radians = es.math.arctan2(b, a) + in_degrees = es.math.arctan2(b, a, unit='deg') + + # EXPECT + assert in_degrees.unit == 'deg' + assert in_degrees.value == pytest.approx(np.degrees(in_radians.value), rel=1e-12) + assert in_degrees.variance == pytest.approx( + in_radians.variance * np.degrees(1.0) ** 2, rel=1e-9 + ) + + +# ---------------------------------------------------------------------- +# Deriving values inside __call__ (a Parameter is created per evaluation) +# ---------------------------------------------------------------------- + +# NOTE(): Under DFO the library itself prints "RuntimeWarning: invalid value encountered in sqrt" +# (minimizer_base.py, in the error-matrix calculation). It is not caused by es.math, and it is +# deliberately not filtered here so that it stays visible. +# +# NOTE(): The uncertainty a fit reports depends on the minimizer. For the same data the variance of +# the fitted FWHM was about 6.2e-5 (LMFit), 6.7e-5 (Bumps) and 3.0e-6 (DFO). Whether DFO's error +# estimate is wrong is a question for the lead developer. A derived variance is only as good as the +# variance it starts from, so these tests compare a derived variance only with the variance of the +# same engine's fit, never across engines. +ENGINES = [None, AvailableMinimizers.LMFit, AvailableMinimizers.Bumps, AvailableMinimizers.DFO] + + +def count_parameters() -> int: + """Number of live ``Parameter`` objects (after a garbage collection).""" + gc.collect() + return sum(isinstance(obj, Parameter) for obj in gc.get_objects()) + + +@pytest.fixture +def gaussian_data(): + x = np.linspace(-5.0, 5.0, 200) + truth = GaussianPeak(amplitude=3.0, center=0.4, fwhm=1.8) + return x, noisy(truth(x), seed=5) + + +class TestGaussianPeakWithDerivedWidth: + def test_derived_width_value_and_variance_match_the_closed_form(self): + # WHEN + peak = GaussianPeak(amplitude=1.0, center=0.0, fwhm=1.8) + peak.fwhm.variance = 0.01 + + # THEN + sigma = peak.sigma() + + # EXPECT + # sigma = fwhm / c with c a constant, so var(sigma) = var(fwhm) / c**2. + assert sigma.value == pytest.approx(1.8 / FWHM_PER_SIGMA, rel=1e-12) + assert sigma.variance == pytest.approx(0.01 / FWHM_PER_SIGMA**2, rel=1e-9) + + def test_curve_equals_the_numpy_reference(self): + # WHEN + x = np.linspace(-5.0, 5.0, 50) + peak = GaussianPeak(amplitude=3.0, center=0.4, fwhm=1.8) + + # THEN + y = peak(x) + + # EXPECT + sigma = 1.8 / FWHM_PER_SIGMA + assert y == pytest.approx(3.0 * np.exp(-((x - 0.4) ** 2) / (2 * sigma**2)), rel=1e-12) + + @pytest.mark.parametrize('engine', ENGINES, ids=[str(engine) for engine in ENGINES]) + def test_fit_recovers_the_parameters_with_every_minimizer(self, gaussian_data, engine): + # WHEN + x, y = gaussian_data + peak = GaussianPeak(amplitude=2.0, center=0.0, fwhm=1.0) + fitter = Fitter(peak, peak) + if engine is not None: + try: + fitter.switch_minimizer(engine) + except AttributeError: + pytest.skip(reason=f'{engine} is not installed') + + # THEN + result = fitter.fit(x, y, weights=np.full_like(x, 1.0 / SIGMA)) + + # EXPECT + assert result.success + assert peak.amplitude.value == pytest.approx(3.0, abs=0.05) + assert peak.center.value == pytest.approx(0.4, abs=0.05) + assert peak.fwhm.value == pytest.approx(1.8, abs=0.05) + + def test_in_window_fraction_after_the_fit(self, gaussian_data): + # WHEN + x, y = gaussian_data + peak = GaussianPeak(amplitude=2.0, center=0.0, fwhm=1.0) + fit(peak, x, y) + window = 1.0 + + # THEN + fraction = peak.fraction_within(window) + + # EXPECT + # u = window * c / (sqrt(2) fwhm) depends on the fwhm alone, so the delta method is exact + # to first order: var(f) = (erf'(u) * du/dfwhm)**2 * var(fwhm), with du/dfwhm = -u/fwhm. + fwhm = peak.fwhm + u = window * FWHM_PER_SIGMA / (math.sqrt(2.0) * fwhm.value) + slope = 2 / math.sqrt(math.pi) * math.exp(-(u**2)) * (-u / fwhm.value) + assert fraction.value == pytest.approx(math.erf(u), rel=1e-12) + assert fraction.variance == pytest.approx(slope**2 * fwhm.variance, rel=1e-9) + + def test_in_window_fraction_variance_agrees_with_sampling_the_fwhm(self, gaussian_data): + # WHEN + x, y = gaussian_data + peak = GaussianPeak(amplitude=2.0, center=0.0, fwhm=1.0) + fit(peak, x, y) + rng = np.random.default_rng(13) + fwhm_samples = rng.normal(peak.fwhm.value, np.sqrt(peak.fwhm.variance), 200_000) + + # THEN + fraction = peak.fraction_within(1.0) + + # EXPECT + sampled = _erf(1.0 * FWHM_PER_SIGMA / (math.sqrt(2.0) * fwhm_samples)) + assert np.std(sampled) == pytest.approx(np.sqrt(fraction.variance), rel=0.03) + + +class TestCostOfDerivingInsideCall: + """ + ``GaussianPeak.__call__`` creates several ``Parameter`` objects per evaluation (the + ``log``, ``sqrt`` and arithmetic intermediates). These tests pin that this is affordable and + does not leak. + """ + + def test_parameters_created_per_call_are_garbage_collected(self): + # WHEN + x = np.linspace(-5.0, 5.0, 200) + peak = GaussianPeak(amplitude=3.0, center=0.4, fwhm=1.8) + peak(x) # warm-up + before = count_parameters() + + # THEN + for _ in range(500): + peak(x) + after = count_parameters() + + # EXPECT + assert after == before + + def test_a_single_evaluation_stays_cheap(self): + # WHEN + x = np.linspace(-5.0, 5.0, 200) + peak = GaussianPeak(amplitude=3.0, center=0.4, fwhm=1.8) + peak(x) # warm-up + durations = [] + + # THEN + for _ in range(200): + start = time.perf_counter() + peak(x) + durations.append(time.perf_counter() - start) + + # EXPECT + # About 1 ms on a development machine; the bound is generous for slow CI runners. + assert np.median(durations) < 0.05 + + def test_results_are_snapshots_so_deriving_once_goes_stale(self): + # WHEN + peak = GaussianPeak(amplitude=3.0, center=0.4, fwhm=1.8) + derived_once = peak.sigma() + stale_value = derived_once.value + + # THEN + peak.fwhm = 3.0 + + # EXPECT + # Deriving in __init__ (once) would silently not follow the fit; derive in __call__. + assert derived_once.value == stale_value + assert peak.sigma().value == pytest.approx(3.0 / FWHM_PER_SIGMA, rel=1e-12) + + +# ---------------------------------------------------------------------- +# Known limitation: Parameter variances carry no covariance +# ---------------------------------------------------------------------- + + +class TestCorrelatedParametersKnownLimitation: + """ + A ``Parameter`` stores a variance but no covariance with other parameters. The slope and + intercept of a line are strongly anticorrelated (about -0.87 here), so a function of both + (``arctan2``) is propagated as if they were independent and its variance is wrong. + + ``strict=True``: if covariance support is ever added, this XPASSes and fails, which is the + signal to remove the marker. + """ + + def test_premise_slope_and_intercept_are_strongly_correlated(self, fitted_line): + # WHEN + _, x = fitted_line + design = np.column_stack([x, np.ones_like(x)]) + + # THEN + covariance = np.linalg.inv(design.T @ design) * SIGMA**2 + correlation = covariance[0, 1] / np.sqrt(covariance[0, 0] * covariance[1, 1]) + + # EXPECT + assert correlation < -0.8 + + @pytest.mark.xfail( + strict=True, + raises=AssertionError, + reason='Known limitation: Parameter variances carry no covariance (deferred).', + ) + def test_arctan2_of_correlated_parameters_matches_the_full_covariance(self, fitted_line): + # WHEN + model, x = fitted_line + design = np.column_stack([x, np.ones_like(x)]) + covariance = np.linalg.inv(design.T @ design) * SIGMA**2 # true parameter covariance + slope, intercept = model.slope.value, model.intercept.value + r2 = slope**2 + intercept**2 + jacobian = np.array([-intercept / r2, slope / r2]) # d atan2(intercept, slope) / d(p) + variance_with_covariance = float(jacobian @ covariance @ jacobian) + + # THEN + derived = es.math.arctan2(model.intercept, model.slope) + + # EXPECT + assert derived.variance == pytest.approx(variance_with_covariance, rel=0.2) diff --git a/tests/unit/variable/test_math.py b/tests/unit/variable/test_math.py new file mode 100644 index 00000000..e225c94e --- /dev/null +++ b/tests/unit/variable/test_math.py @@ -0,0 +1,1908 @@ +# SPDX-FileCopyrightText: 2026 EasyScience contributors +# SPDX-License-Identifier: BSD-3-Clause + +""" +Unit tests for ``easyscience.variable.math``. + +Conventions encoded here (agreed with the lead developer): + +* **Values are validated strictly.** A value outside a function's domain, + or a result that is not finite, raises ``ValueError`` (such values break + model-fitting loops). +* **Bounds are permissive.** Bounds never raise: bounds outside a domain are + clipped to it, infinite bounds are accepted and bounds that cross an + asymptote (``tan``) blow up to ``(-inf, inf)``. +* **Bounds are exact**, not sampled: they are compared against a brute-force + dense evaluation of the function over the input interval. +""" + +import importlib +import math +import warnings + +import numpy as np +import pytest +import scipp as sc +from scipp import UnitError + +import easyscience as es +import easyscience.variable as es_variable +from easyscience import DescriptorNumber +from easyscience import Parameter +from easyscience import global_object +from easyscience.variable import DescriptorArray +from easyscience.variable import math as es_math + +INF = np.inf +SINC_GLOBAL_MIN = -0.21723362821122166 # sinc(x) at its first minimum, x ~ 4.4934 + +UNARY_FUNCTIONS = [ + 'exp', + 'log', + 'log10', + 'sqrt', + 'erf', + 'erfc', + 'sin', + 'cos', + 'tan', + 'arcsin', + 'arccos', + 'arctan', + 'sinc', + 'sinh', + 'cosh', + 'tanh', + 'arcsinh', + 'arccosh', + 'arctanh', +] +ALL_FUNCTIONS = [*UNARY_FUNCTIONS, 'arctan2'] + + +# ---------------------------------------------------------------------- +# Helpers / reference implementations +# ---------------------------------------------------------------------- + + +def make(value, unit='', *, min=-INF, max=INF, variance=0.0, name='x') -> Parameter: + """Shorthand for building a ``Parameter`` with the bounds under test.""" + return Parameter(name=name, value=value, unit=unit, variance=variance, min=min, max=max) + + +def fn(name: str): + return getattr(es_math, name) + + +def unnormalized_sinc(x): + """sin(x)/x, the definition Scipp uses (``np.sinc`` is the normalised one).""" + return np.sinc(np.asarray(x) / np.pi) + + +_erf = np.vectorize(math.erf) +_erfc = np.vectorize(math.erfc) + + +def dense_range(np_func, lo, hi, n=200_001): + """Brute-force oracle: the range of ``np_func`` sampled on a dense grid.""" + xs = np.linspace(lo, hi, n) + with np.errstate(all='ignore'): + ys = np_func(xs) + return float(np.min(ys)), float(np.max(ys)) + + +def assert_bounds(result: Parameter, expected_min, expected_max, tol=1e-9): + assert np.isclose(result.min, expected_min, rtol=0, atol=tol), (result.min, expected_min) + assert np.isclose(result.max, expected_max, rtol=0, atol=tol), (result.max, expected_max) + + +# name -> (numpy/math reference, derivative of the reference, value, unit) +VALUE_CASES = { + 'exp': (np.exp, np.exp, 0.3, ''), + 'log': (np.log, lambda x: 1 / x, 2.0, ''), + 'log10': (np.log10, lambda x: 1 / (x * np.log(10)), 2.0, ''), + 'sqrt': (np.sqrt, lambda x: 0.5 / np.sqrt(x), 2.0, ''), + 'erf': (math.erf, lambda x: 2 / np.sqrt(np.pi) * np.exp(-(x**2)), 0.4, ''), + 'erfc': (math.erfc, lambda x: -2 / np.sqrt(np.pi) * np.exp(-(x**2)), 0.4, ''), + 'sin': (np.sin, np.cos, 0.7, 'rad'), + 'cos': (np.cos, lambda x: -np.sin(x), 0.7, 'rad'), + 'tan': (np.tan, lambda x: 1 / np.cos(x) ** 2, 0.7, 'rad'), + 'arcsin': (np.arcsin, lambda x: 1 / np.sqrt(1 - x**2), 0.4, ''), + 'arccos': (np.arccos, lambda x: -1 / np.sqrt(1 - x**2), 0.4, ''), + 'arctan': (np.arctan, lambda x: 1 / (1 + x**2), 0.4, ''), + 'sinc': (unnormalized_sinc, lambda x: (x * np.cos(x) - np.sin(x)) / x**2, 1.3, 'rad'), + 'sinh': (np.sinh, np.cosh, 0.4, ''), + 'cosh': (np.cosh, np.sinh, 0.4, ''), + 'tanh': (np.tanh, lambda x: 1 / np.cosh(x) ** 2, 0.4, ''), + 'arcsinh': (np.arcsinh, lambda x: 1 / np.sqrt(x**2 + 1), 0.4, ''), + 'arccosh': (np.arccosh, lambda x: 1 / np.sqrt(x**2 - 1), 1.7, ''), + 'arctanh': (np.arctanh, lambda x: 1 / (1 - x**2), 0.4, ''), +} + + +@pytest.fixture(autouse=True) +def clear(): + global_object.map._clear() + + +# ---------------------------------------------------------------------- +# Value and variance +# ---------------------------------------------------------------------- + + +class TestValueAndVariance: + @pytest.mark.parametrize('name', UNARY_FUNCTIONS) + def test_value_matches_reference(self, name): + # WHEN + reference, _, value, unit = VALUE_CASES[name] + x = make(value, unit) + + # THEN + result = fn(name)(x) + + # EXPECT + assert isinstance(result, Parameter) + assert result.value == pytest.approx(float(reference(value)), rel=1e-12) + + @pytest.mark.parametrize('name', UNARY_FUNCTIONS) + def test_variance_matches_delta_method(self, name): + # WHEN + # sigma_y^2 = f'(x)^2 * sigma_x^2. This also covers erf/erfc, for which + # Scipp itself has no variance support and EasyScience must supply it. + _, derivative, value, unit = VALUE_CASES[name] + variance = 0.04 + x = make(value, unit, variance=variance) + + # THEN + result = fn(name)(x) + + # EXPECT + expected = float(derivative(value)) ** 2 * variance + assert result.variance == pytest.approx(expected, rel=1e-9) + + @pytest.mark.parametrize('name', ['sin', 'cos', 'tan']) + def test_variance_in_degrees_follows_the_unit_conversion(self, name): + # WHEN + _, derivative, _, _ = VALUE_CASES[name] + angle_in_degrees, variance_in_degrees2 = 30.0, 1.0 + x = make(angle_in_degrees, 'deg', variance=variance_in_degrees2) + + # THEN + result = fn(name)(x) + + # EXPECT + # d f / d(deg) = f'(rad) * pi / 180, so the variance picks up (pi / 180) ** 2. + radians = np.deg2rad(angle_in_degrees) + expected = float(derivative(radians)) ** 2 * variance_in_degrees2 * np.deg2rad(1.0) ** 2 + assert result.variance == pytest.approx(expected, rel=1e-9) + + @pytest.mark.parametrize('name', UNARY_FUNCTIONS) + @pytest.mark.parametrize('variance', [0.0, None]) + def test_no_input_variance_gives_no_output_variance(self, name, variance): + # WHEN + _, _, value, unit = VALUE_CASES[name] + x = make(value, unit, variance=variance) + + # THEN + result = fn(name)(x) + + # EXPECT + assert result.variance is None + + @pytest.mark.parametrize('name', UNARY_FUNCTIONS) + def test_inputs_are_not_modified(self, name): + # WHEN + _, _, value, unit = VALUE_CASES[name] + x = make(value, unit, variance=0.04, min=value - 0.1, max=value + 0.1) + before = (x.value, x.variance, x.min, x.max, x.unit) + + # THEN + fn(name)(x) + + # EXPECT + assert (x.value, x.variance, x.min, x.max, x.unit) == before + + def test_result_name_is_built_from_function_and_argument(self): + # WHEN + angle = make(0.5, 'rad', name='angle') + y = make(1.0, name='y') + x = make(2.0, name='x') + + # THEN + sin_result = es_math.sin(angle) + sinc_result = es_math.sinc(angle) + arctan2_result = es_math.arctan2(y, x) + + # EXPECT + assert sin_result.name == 'sin(angle)' + assert sinc_result.name == 'sinc(angle)' + assert arctan2_result.name == 'arctan2(y, x)' + + +# ---------------------------------------------------------------------- +# Bounds: exact, compared against brute-force evaluation +# ---------------------------------------------------------------------- + +# (function, reference, lo, hi, value, unit) +ORACLE_CASES = [ + # sin: bounds exactly on an extremum (regression: a 100-point grid missed + # sin(pi/2) = 1 and Parameter.from_scipp then raised "value > max") + ('sin', np.sin, 0.0, np.pi, np.pi / 2, 'rad'), + ('sin', np.sin, 0.1, 0.5, 0.3, 'rad'), # monotonic piece only + ('sin', np.sin, -10.0, 10.0, 1.0, 'rad'), # several periods + ('sin', np.sin, 1.0, 5.0, 2.0, 'rad'), # contains a max and a min + ('sin', np.sin, 2.0, 3.0, 2.5, 'rad'), # decreasing piece only + ('cos', np.cos, -1.0, 2.0, 0.5, 'rad'), + ('cos', np.cos, 3.0, 7.0, 4.0, 'rad'), + ('cos', np.cos, 0.2, 0.9, 0.5, 'rad'), + ('sinc', unnormalized_sinc, -10.0, 10.0, 1.0, 'rad'), + ('sinc', unnormalized_sinc, 3.0, 8.0, 5.0, 'rad'), # first minimum and second maximum + ('sinc', unnormalized_sinc, 0.5, 2.0, 1.0, 'rad'), + ('sinc', unnormalized_sinc, -1.0, 20.0, 3.0, ''), # dimensionless is accepted too + ('sinc', unnormalized_sinc, 100.0, 110.0, 105.0, 'rad'), # far-away, tiny lobes + ('cosh', np.cosh, -1.0, 2.0, 0.5, ''), # minimum inside + ('cosh', np.cosh, 1.0, 2.0, 1.5, ''), # monotonic piece only + ('tan', np.tan, -1.0, 1.0, 0.5, 'rad'), + ('exp', np.exp, -2.0, 3.0, 0.5, ''), + ('log', np.log, 0.5, 8.0, 2.0, ''), + ('log10', np.log10, 0.5, 8.0, 2.0, ''), + ('sqrt', np.sqrt, 0.5, 8.0, 2.0, ''), + ('erf', _erf, -1.0, 2.0, 0.5, ''), + ('erfc', _erfc, -1.0, 2.0, 0.5, ''), # decreasing + ('sinh', np.sinh, -1.0, 2.0, 0.5, ''), + ('tanh', np.tanh, -1.0, 2.0, 0.5, ''), + ('arcsinh', np.arcsinh, -1.0, 2.0, 0.5, ''), + ('arccosh', np.arccosh, 1.5, 6.0, 2.0, ''), + ('arctanh', np.arctanh, -0.5, 0.8, 0.1, ''), + ('arcsin', np.arcsin, -0.5, 0.8, 0.1, ''), + ('arccos', np.arccos, -0.5, 0.7, 0.1, ''), # decreasing + ('arctan', np.arctan, -3.0, 5.0, 1.0, ''), +] + + +class TestBounds: + @pytest.mark.parametrize(('name', 'reference', 'lo', 'hi', 'value', 'unit'), ORACLE_CASES) + def test_bounds_match_brute_force_range(self, name, reference, lo, hi, value, unit): + # WHEN + x = make(value, unit, min=lo, max=hi) + + # THEN + result = fn(name)(x) + + # EXPECT + true_min, true_max = dense_range(reference, lo, hi) + # Contains everything the dense sampling found (up to rounding) ... + assert result.min <= true_min + 1e-12 + assert result.max >= true_max - 1e-12 + # ... and is tight: no looser than the sampling resolution. + assert result.min >= true_min - 1e-6 + assert result.max <= true_max + 1e-6 + + @pytest.mark.parametrize( + ('name', 'value', 'unit', 'expected_min', 'expected_max'), + [ + ('sin', 0.5, 'rad', -1.0, 1.0), + ('cos', 0.5, 'rad', -1.0, 1.0), + ('tan', 0.3, 'rad', -INF, INF), + ('sinc', 0.5, 'rad', SINC_GLOBAL_MIN, 1.0), + ('exp', 0.5, '', 0.0, INF), + ('log', 2.0, '', -INF, INF), + ('log10', 2.0, '', -INF, INF), + ('sqrt', 2.0, '', 0.0, INF), + ('erf', 0.5, '', -1.0, 1.0), + ('erfc', 0.5, '', 0.0, 2.0), + ('sinh', 0.5, '', -INF, INF), + ('cosh', 0.5, '', 1.0, INF), + ('tanh', 0.5, '', -1.0, 1.0), + ('arcsinh', 0.5, '', -INF, INF), + ('arccosh', 2.0, '', 0.0, INF), + ('arctanh', 0.3, '', -INF, INF), + ('arcsin', 0.3, '', -np.pi / 2, np.pi / 2), + ('arccos', 0.3, '', 0.0, np.pi), + ('arctan', 0.3, '', -np.pi / 2, np.pi / 2), + ], + ) + def test_default_infinite_bounds_are_accepted( + self, name, value, unit, expected_min, expected_max + ): + # WHEN + # A freshly created Parameter has bounds (-inf, inf); that must just work. + x = make(value, unit) + + # THEN + result = fn(name)(x) + + # EXPECT + assert_bounds(result, expected_min, expected_max) + + def test_sin_with_degrees_bounds_are_computed_in_radians(self): + # WHEN + # Regression: numpy used to be fed 0..90 as if they were radians. + angle = make(45.0, 'deg', min=0.0, max=90.0) + + # THEN + result = es_math.sin(angle) + + # EXPECT + assert result.value == pytest.approx(np.sqrt(0.5)) + assert_bounds(result, 0.0, 1.0) + + def test_sin_bounds_include_an_interior_maximum(self): + # WHEN + # sin peaks at 90 deg, which lies inside [15, 95] deg, so the upper bound is 1.0 and + # NOT sin(95 deg) = 0.9962, the value at the endpoint. + angle = make(45.0, 'deg', min=15.0, max=95.0) + + # THEN + result = es_math.sin(angle) + + # EXPECT + assert result.value == pytest.approx(0.707106781, abs=1e-9) + assert_bounds(result, 0.258819045, 1.0) + + def test_sin_bounds_are_the_endpoint_values_when_no_extremum_is_inside(self): + # WHEN + # sin is increasing on all of [15, 85] deg, so the bounds are just sin(15) and sin(85). + angle = make(45.0, 'deg', min=15.0, max=85.0) + + # THEN + result = es_math.sin(angle) + + # EXPECT + assert result.value == pytest.approx(0.707106781, abs=1e-9) + assert_bounds(result, 0.258819045, 0.996194698) + + def test_cos_bounds_are_computed_in_radians_for_degree_input(self): + # WHEN + angle = make(30.0, 'deg', min=0.0, max=90.0) + + # THEN + result = es_math.cos(angle) + + # EXPECT + assert_bounds(result, 0.0, 1.0) + + def test_sin_with_degrees_and_infinite_bounds(self): + # WHEN + angle = make(45.0, 'deg') + + # THEN + result = es_math.sin(angle) + + # EXPECT + assert_bounds(result, -1.0, 1.0) + + def test_inverse_trig_bounds_are_in_radians(self): + # WHEN + x = make(0.5, min=0.0, max=1.0) + + # THEN + result = es_math.arcsin(x) + + # EXPECT + assert result.unit == 'rad' + assert_bounds(result, 0.0, np.pi / 2) + + def test_sinc_global_minimum_is_found(self): + # WHEN + x = make(1.0, 'rad', min=-20.0, max=20.0) + + # THEN + result = es_math.sinc(x) + + # EXPECT + assert_bounds(result, SINC_GLOBAL_MIN, 1.0) + + def test_sinc_negative_one_sided_infinite_bounds(self): + # WHEN + x = make(-150.0, 'rad', max=-100.0) + + # THEN + result = es_math.sinc(x) + + # EXPECT + true_min, true_max = dense_range(unnormalized_sinc, -1000.0, -100.0) + assert result.min <= true_min + 1e-12 + assert result.max >= true_max - 1e-12 + assert abs(result.min) < 0.01 + assert 0 < result.max < 0.01 + + @pytest.mark.parametrize( + ('value', 'lo', 'hi'), + [(2e20, 1e20, INF), (-2e20, -INF, -1e20), (3e20, 1e20, 5e20)], + ) + def test_sinc_with_huge_bounds_does_not_raise(self, value, lo, hi): + # WHEN + # Regression: int(floor(lo / pi)) overflowed numpy's integer range for lo >~ 5.8e19. + x = make(value, 'rad', min=lo, max=hi) + + # THEN + result = es_math.sinc(x) + + # EXPECT + # The lobes are below 1e-19 out here, so the bounds are essentially [0, 0]. + assert abs(result.min) < 1e-19 + assert abs(result.max) < 1e-19 + + def test_sinc_one_sided_infinite_bounds(self): + # WHEN + # Lobes decay like 1/x; the extremes beyond x=100 are tiny but present. + x = make(150.0, 'rad', min=100.0) + + # THEN + result = es_math.sinc(x) + + # EXPECT + true_min, true_max = dense_range(unnormalized_sinc, 100.0, 1000.0) + assert result.min <= true_min + 1e-12 + assert result.max >= true_max - 1e-12 + assert abs(result.min) < 0.01 + assert 0 < result.max < 0.01 + + +# ---------------------------------------------------------------------- +# Domain: strict on values, permissive on bounds +# ---------------------------------------------------------------------- + + +class TestValueGuard: + @pytest.mark.parametrize( + ('name', 'value', 'unit'), + [ + ('log', 0.0, ''), + ('log', -1.0, ''), + ('log10', 0.0, ''), + ('log10', -2.0, ''), + ('sqrt', -1.0, ''), + ('arcsin', 1.5, ''), + ('arcsin', -1.5, ''), + ('arccos', 1.5, ''), + ('arccos', -1.5, ''), + ('arccosh', 0.5, ''), + ('arctanh', 1.0, ''), + ('arctanh', -1.0, ''), + ('arctanh', 1.5, ''), + ], + ) + def test_value_outside_domain_raises(self, name, value, unit): + # WHEN THEN + x = make(value, unit) + + # EXPECT + with pytest.raises(ValueError, match='outside the domain'): + fn(name)(x) + + @pytest.mark.parametrize('name', ['exp', 'sinh', 'cosh']) + def test_overflowing_result_raises(self, name): + # WHEN THEN + x = make(1000.0) + + # EXPECT + with pytest.raises(ValueError, match='not finite'): + fn(name)(x) + + @pytest.mark.parametrize('name', ['exp', 'sinh', 'cosh']) + def test_overflowing_variance_raises_although_the_value_is_finite(self, name): + # WHEN THEN + # The value is finite at x=400, but the variance scales with exp(2x) and overflows. + x = make(400.0, variance=1.0) + with np.errstate(all='ignore'): + assert np.isfinite(getattr(np, name)(400.0)) + + # EXPECT + with pytest.raises(ValueError, match='not finite'): + fn(name)(x) + + @pytest.mark.parametrize('name', ['sqrt', 'arcsin', 'arccos', 'arccosh']) + def test_infinite_variance_at_domain_edge_raises(self, name): + # WHEN THEN + # d/dx is infinite at the edge, so a non-zero input variance gives a + # non-finite output variance. + edge = {'sqrt': 0.0, 'arcsin': 1.0, 'arccos': 1.0, 'arccosh': 1.0}[name] + x = make(edge, variance=0.01) + + # EXPECT + with pytest.raises(ValueError, match='not finite'): + fn(name)(x) + + @pytest.mark.parametrize( + ('name', 'value', 'expected'), + [ + ('sqrt', 0.0, 0.0), + ('arcsin', 1.0, np.pi / 2), + ('arcsin', -1.0, -np.pi / 2), + ('arccos', 1.0, 0.0), + ('arccos', -1.0, np.pi), + ('arccosh', 1.0, 0.0), + ], + ) + def test_domain_edge_value_is_allowed_without_variance(self, name, value, expected): + # WHEN + x = make(value) + + # THEN + result = fn(name)(x) + + # EXPECT + assert result.value == pytest.approx(expected) + + def test_error_message_names_the_parameter(self): + # WHEN THEN + x = make(-1.0, name='my_param') + + # EXPECT + with pytest.raises(ValueError, match='my_param'): + es_math.log(x) + + +class TestPermissiveBounds: + """Bounds that leave the domain are clipped to it, never an error or NaN.""" + + @pytest.fixture(autouse=True) + def no_warnings(self): + with warnings.catch_warnings(): + warnings.simplefilter('error') + yield + + def test_log_bounds_below_zero(self): + # WHEN + x = make(2.0, min=-1.0, max=5.0) + + # THEN + result = es_math.log(x) + + # EXPECT + assert_bounds(result, -INF, np.log(5.0)) + + def test_log10_bounds_touching_zero(self): + # WHEN + x = make(10.0, min=0.0, max=100.0) + + # THEN + result = es_math.log10(x) + + # EXPECT + assert_bounds(result, -INF, 2.0) + + def test_sqrt_bounds_below_zero(self): + # WHEN + x = make(4.0, min=-4.0, max=9.0) + + # THEN + result = es_math.sqrt(x) + + # EXPECT + assert_bounds(result, 0.0, 3.0) + + def test_arcsin_bounds_beyond_one(self): + # WHEN + x = make(0.2, min=-3.0, max=3.0) + + # THEN + result = es_math.arcsin(x) + + # EXPECT + assert_bounds(result, -np.pi / 2, np.pi / 2) + + def test_arccos_bounds_beyond_one(self): + # WHEN + x = make(0.2, min=-3.0, max=3.0) + + # THEN + result = es_math.arccos(x) + + # EXPECT + assert_bounds(result, 0.0, np.pi) + + def test_arccosh_bounds_below_one(self): + # WHEN + x = make(3.0, min=-5.0, max=10.0) + + # THEN + result = es_math.arccosh(x) + + # EXPECT + assert_bounds(result, 0.0, np.arccosh(10.0)) + + def test_arctanh_bounds_beyond_one(self): + # WHEN + x = make(0.3, min=-5.0, max=5.0) + + # THEN + result = es_math.arctanh(x) + + # EXPECT + assert_bounds(result, -INF, INF) + + @pytest.mark.parametrize('name', ['log', 'log10', 'sqrt', 'arcsin', 'arccos', 'arctanh']) + def test_bounds_are_never_nan(self, name): + # WHEN + x = make(0.5, min=-10.0, max=10.0) + + # THEN + result = fn(name)(x) + + # EXPECT + assert not np.isnan(result.min) + assert not np.isnan(result.max) + + +# ---------------------------------------------------------------------- +# tan +# ---------------------------------------------------------------------- + + +class TestTan: + def test_bounds_without_asymptote(self): + # WHEN + angle = make(0.3, 'rad', min=-1.0, max=1.0) + + # THEN + result = es_math.tan(angle) + + # EXPECT + assert_bounds(result, np.tan(-1.0), np.tan(1.0)) + + def test_bounds_crossing_asymptote_blow_up(self): + # WHEN + angle = make(0.3, 'rad', min=0.0, max=2.0) + + # THEN + result = es_math.tan(angle) + + # EXPECT + assert_bounds(result, -INF, INF) + + def test_bounds_crossing_asymptote_in_degrees(self): + # WHEN + angle = make(30.0, 'deg', min=0.0, max=100.0) + + # THEN + result = es_math.tan(angle) + + # EXPECT + assert_bounds(result, -INF, INF) + + def test_bounds_in_degrees_without_asymptote(self): + # WHEN + angle = make(30.0, 'deg', min=0.0, max=60.0) + + # THEN + result = es_math.tan(angle) + + # EXPECT + assert_bounds(result, 0.0, np.tan(np.pi / 3)) + + def test_upper_bound_exactly_on_an_asymptote_is_infinite_on_that_side_only(self): + # WHEN + # tan rises from tan(0) = 0 to +inf as the angle approaches 90 deg. + angle = make(45.0, 'deg', min=0.0, max=90.0) + + # THEN + result = es_math.tan(angle) + + # EXPECT + assert_bounds(result, 0.0, INF) + + def test_upper_bound_exactly_on_an_asymptote_in_radians(self): + # WHEN + angle = make(0.5, 'rad', min=0.0, max=np.pi / 2) + + # THEN + result = es_math.tan(angle) + + # EXPECT + assert_bounds(result, 0.0, INF) + + def test_lower_bound_exactly_on_an_asymptote_is_infinite_on_that_side_only(self): + # WHEN + # Just after 90 deg tan is -inf and rises to tan(120 deg) = -sqrt(3). + angle = make(100.0, 'deg', min=90.0, max=120.0) + + # THEN + result = es_math.tan(angle) + + # EXPECT + assert_bounds(result, -INF, -np.sqrt(3.0)) + + def test_bounds_on_two_asymptotes_cover_all_real_numbers(self): + # WHEN + angle = make(0.0, 'deg', min=-90.0, max=90.0) + + # THEN + result = es_math.tan(angle) + + # EXPECT + assert_bounds(result, -INF, INF) + + @pytest.mark.parametrize( + ('value', 'unit'), [(np.pi / 2, 'rad'), (90.0, 'deg'), (-90.0, 'deg')] + ) + def test_value_on_asymptote_raises(self, value, unit): + # WHEN THEN + angle = make(value, unit) + + # EXPECT + with pytest.raises(ValueError, match='asymptote'): + es_math.tan(angle) + + def test_value_within_one_sigma_of_asymptote_raises(self): + # WHEN THEN + angle = make(np.pi / 2 - 1e-6, 'rad', variance=1e-4) + + # EXPECT + with pytest.raises(ValueError, match='asymptote'): + es_math.tan(angle) + + def test_value_within_one_sigma_of_asymptote_in_degrees_raises(self): + # WHEN THEN + angle = make(89.9, 'deg', variance=1.0) + + # EXPECT + with pytest.raises(ValueError, match='asymptote'): + es_math.tan(angle) + + @pytest.mark.parametrize('value', [np.pi / 2 - 5e-8, np.pi / 2 + 5e-8]) + def test_value_inside_the_numerical_floor_raises(self, value): + # WHEN THEN + # With no variance the margin is the 1e-7 rad floor. + angle = make(value, 'rad') + + # EXPECT + with pytest.raises(ValueError, match='asymptote'): + es_math.tan(angle) + + @pytest.mark.parametrize('value', [np.pi / 2 - 2e-7, np.pi / 2 + 2e-7]) + def test_value_just_outside_the_numerical_floor_is_allowed(self, value): + # WHEN + angle = make(value, 'rad') + + # THEN + result = es_math.tan(angle) + + # EXPECT + assert np.isfinite(result.value) + + def test_value_within_one_sigma_of_asymptote_raises_when_sigma_exceeds_the_floor(self): + # WHEN THEN + # variance=1e-4 -> sigma = 0.01 rad, far above the 1e-7 rad floor; 0.009 rad away. + angle = make(np.pi / 2 - 0.009, 'rad', variance=1e-4) + + # EXPECT + with pytest.raises(ValueError, match='asymptote'): + es_math.tan(angle) + + def test_value_just_beyond_one_sigma_of_asymptote_is_allowed(self): + # WHEN + # sigma = 0.01 rad; the asymptote is 0.011 rad away. + angle = make(np.pi / 2 - 0.011, 'rad', variance=1e-4) + + # THEN + result = es_math.tan(angle) + + # EXPECT + assert np.isfinite(result.variance) + + @pytest.mark.parametrize( + ('value', 'unit'), [(3 * np.pi / 2, 'rad'), (270.0, 'deg'), (-3 * np.pi / 2, 'rad')] + ) + def test_asymptotes_other_than_plus_minus_pi_over_two_raise(self, value, unit): + # WHEN THEN + angle = make(value, unit) + + # EXPECT + with pytest.raises(ValueError, match='asymptote'): + es_math.tan(angle) + + @pytest.mark.parametrize('variance', [0.0, None]) + def test_missing_variance_is_guarded_like_zero_variance(self, variance): + # WHEN THEN + angle = make(np.pi / 2, 'rad', variance=variance) + + # EXPECT + with pytest.raises(ValueError, match='asymptote'): + es_math.tan(angle) + + def test_value_far_from_asymptote_in_sigma_is_allowed(self): + # WHEN + angle = make(1.0, 'rad', variance=1e-4) + + # THEN + result = es_math.tan(angle) + + # EXPECT + assert np.isfinite(result.variance) + + def test_value_near_but_not_on_asymptote_is_allowed(self): + # WHEN + angle = make(np.pi / 2 - 1e-3, 'rad') + + # THEN + result = es_math.tan(angle) + + # EXPECT + assert result.value == pytest.approx(np.tan(np.pi / 2 - 1e-3)) + + +# ---------------------------------------------------------------------- +# arctan2 (binary; Scipp has no variance support for it) +# ---------------------------------------------------------------------- + + +class TestArctan2: + def test_value_is_the_angle_in_radians_and_unit_is_rad(self): + # WHEN + y = make(1.0, 'm', name='y') + x = make(2.0, 'm', name='x') + + # THEN + result = es_math.arctan2(y, x) + + # EXPECT + assert result.value == pytest.approx(np.arctan2(1.0, 2.0)) + assert result.unit == 'rad' + + @pytest.mark.parametrize(('y', 'x'), [(1.0, 2.0), (1.0, -2.0), (-1.0, -2.0), (-1.0, 2.0)]) + def test_value_is_in_the_correct_quadrant(self, y, x): + # WHEN + y_param = make(y) + x_param = make(x) + + # THEN + result = es_math.arctan2(y_param, x_param) + + # EXPECT + assert result.value == pytest.approx(np.arctan2(y, x)) + + def test_variance_matches_delta_method(self): + # WHEN + y, x, var_y, var_x = 1.0, 2.0, 0.04, 0.09 + y_param = make(y, variance=var_y) + x_param = make(x, variance=var_x) + + # THEN + result = es_math.arctan2(y_param, x_param) + + # EXPECT + expected = (x**2 * var_y + y**2 * var_x) / (x**2 + y**2) ** 2 + assert result.variance == pytest.approx(expected, rel=1e-12) + + def test_variance_with_only_one_argument_uncertain(self): + # WHEN + y = make(1.0, variance=0.04) + x = make(2.0) + + # THEN + result = es_math.arctan2(y, x) + + # EXPECT + assert result.variance == pytest.approx(2.0**2 * 0.04 / (1.0 + 4.0) ** 2, rel=1e-12) + + def test_no_variance_when_neither_argument_has_one(self): + # WHEN + y = make(1.0) + x = make(2.0) + + # THEN + result = es_math.arctan2(y, x) + + # EXPECT + assert result.variance is None + + def test_variance_at_origin_raises(self): + # WHEN THEN + y = make(0.0, variance=0.01) + x = make(0.0, variance=0.01) + + # EXPECT + with pytest.raises(ValueError, match='origin'): + es_math.arctan2(y, x) + + @pytest.mark.parametrize('var_y', [0.0, None]) + @pytest.mark.parametrize('var_x', [0.0, None]) + def test_zero_and_missing_variance_are_equivalent(self, var_y, var_x): + # WHEN + y = make(1.0, variance=var_y) + x = make(2.0, variance=var_x) + + # THEN + result = es_math.arctan2(y, x) + + # EXPECT + assert result.variance is None + + def test_variance_with_only_x_uncertain(self): + # WHEN + y = make(1.0) + x = make(2.0, variance=0.09) + + # THEN + result = es_math.arctan2(y, x) + + # EXPECT + assert result.variance == pytest.approx(1.0**2 * 0.09 / (4.0 + 1.0) ** 2, rel=1e-12) + + def test_variance_with_one_argument_none_and_the_other_uncertain(self): + # WHEN + y = make(1.0, variance=None) + x = make(2.0, variance=0.09) + + # THEN + result = es_math.arctan2(y, x) + + # EXPECT + assert result.variance == pytest.approx(0.09 / 25.0, rel=1e-12) + + def test_only_x_uncertain_and_close_to_origin_raises(self): + # WHEN THEN + y = make(0.0) + x = make(0.01, variance=0.01) + + # EXPECT + with pytest.raises(ValueError, match='origin'): + es_math.arctan2(y, x) + + def test_point_exactly_one_sigma_from_origin_raises(self): + # WHEN THEN + # sigma = sqrt(0.25) = 0.5 exactly, and the point is exactly 0.5 from the origin. + y = make(0.0) + x = make(0.5, variance=0.25) + + # EXPECT + with pytest.raises(ValueError, match='origin'): + es_math.arctan2(y, x) + + def test_point_just_outside_one_sigma_from_origin_is_allowed(self): + # WHEN + y = make(0.0) + x = make(0.5 + 1e-6, variance=0.25) + + # THEN + result = es_math.arctan2(y, x) + + # EXPECT + assert np.isfinite(result.variance) + + def test_origin_without_variance_is_allowed(self): + # WHEN + y = make(0.0) + x = make(0.0) + + # THEN + result = es_math.arctan2(y, x) + + # EXPECT + assert result.value == 0.0 + + def test_point_within_one_sigma_of_origin_raises(self): + # WHEN THEN + y = make(0.01, variance=0.01) + x = make(0.0) + + # EXPECT + with pytest.raises(ValueError, match='origin'): + es_math.arctan2(y, x) + + def test_point_far_from_origin_in_sigma_is_allowed(self): + # WHEN + y = make(1.0, variance=0.01) + x = make(0.0) + + # THEN + result = es_math.arctan2(y, x) + + # EXPECT + assert np.isfinite(result.variance) + + def test_mismatched_units_raise(self): + # WHEN THEN + y = make(1.0, 'm') + x = make(2.0, 's') + + # EXPECT + with pytest.raises(UnitError): + es_math.arctan2(y, x) + + @pytest.mark.parametrize( + ('y_bounds', 'x_bounds', 'y', 'x'), + [ + ((0.5, 2.0), (1.0, 3.0), 1.0, 2.0), # first quadrant + ((-2.0, -0.5), (1.0, 3.0), -1.0, 2.0), # fourth quadrant + ((0.5, 2.0), (-3.0, -1.0), 1.0, -2.0), # second quadrant + ((-2.0, -0.5), (-3.0, -1.0), -1.0, -2.0), # third quadrant + ((0.0, 1.0), (0.0, 2.0), 0.5, 1.0), # origin on a corner + ((0.0, 1.0), (-1.0, 2.0), 0.5, 1.0), # y >= 0, origin on an edge + ((-1.0, 1.0), (0.0, 2.0), 0.5, 1.0), # x >= 0, origin on an edge + ], + ) + def test_bounds_without_branch_cut_match_brute_force(self, y_bounds, x_bounds, y, x): + # WHEN + y_param = make(y, min=y_bounds[0], max=y_bounds[1]) + x_param = make(x, min=x_bounds[0], max=x_bounds[1]) + + # THEN + with warnings.catch_warnings(): + warnings.simplefilter('error') + result = es_math.arctan2(y_param, x_param) + + # EXPECT + ys, xs = np.meshgrid( + np.linspace(*y_bounds, 801), np.linspace(*x_bounds, 801), indexing='ij' + ) + angles = np.arctan2(ys, xs) + assert result.min <= angles.min() + 1e-12 + assert result.max >= angles.max() - 1e-12 + assert result.min >= angles.min() - 1e-2 + assert result.max <= angles.max() + 1e-2 + + def test_bounds_crossing_branch_cut_warn_and_span_full_range(self): + # WHEN + # y crosses 0 while x < 0: the angle jumps between -pi and +pi. + y = make(0.3, min=-1.0, max=1.0) + x = make(-1.0, min=-2.0, max=-0.5) + + # THEN + with pytest.warns(UserWarning, match='branch cut'): + result = es_math.arctan2(y, x) + + # EXPECT + assert_bounds(result, -np.pi, np.pi) + + def test_bounds_containing_origin_span_full_range(self): + # WHEN + y = make(0.3, min=-1.0, max=1.0) + x = make(0.3, min=-1.0, max=1.0) + + # THEN + with pytest.warns(UserWarning): + result = es_math.arctan2(y, x) + + # EXPECT + assert_bounds(result, -np.pi, np.pi) + + def test_one_sided_infinite_bounds_span_full_range_without_warning(self): + # WHEN + # y crosses 0 while x extends to -inf: the warning is only for fully finite bounds. + y = make(0.3, min=-1.0, max=1.0) + x = make(-2.0, min=-INF, max=-1.0) + + # THEN + with warnings.catch_warnings(): + warnings.simplefilter('error') + result = es_math.arctan2(y, x) + + # EXPECT + assert_bounds(result, -np.pi, np.pi) + + def test_default_infinite_bounds_span_full_range_without_warning(self): + # WHEN + # Unbounded parameters are the default; warning on every call would be noise. + y = make(1.0) + x = make(-2.0) + + # THEN + with warnings.catch_warnings(): + warnings.simplefilter('error') + result = es_math.arctan2(y, x) + + # EXPECT + assert_bounds(result, -np.pi, np.pi) + + def test_value_is_always_inside_bounds(self): + # WHEN + y = make(1.0, min=0.5, max=2.0) + x = make(2.0, min=1.0, max=3.0) + + # THEN + result = es_math.arctan2(y, x) + + # EXPECT + assert result.min <= result.value <= result.max + + +# ---------------------------------------------------------------------- +# Units +# ---------------------------------------------------------------------- + + +class TestUnits: + @pytest.mark.parametrize('name', ['sin', 'cos', 'tan']) + @pytest.mark.parametrize('unit', ['m', '']) + def test_trig_requires_an_angle_unit(self, name, unit): + # WHEN THEN + x = make(0.5, unit) + + # EXPECT + with pytest.raises(UnitError): + fn(name)(x) + + @pytest.mark.parametrize('name', ['sin', 'cos']) + def test_trig_result_is_dimensionless(self, name): + # WHEN + radians = make(0.5, 'rad') + degrees = make(30.0, 'deg') + + # THEN + result_rad = fn(name)(radians) + result_deg = fn(name)(degrees) + + # EXPECT + assert result_rad.unit == 'dimensionless' + assert result_deg.unit == 'dimensionless' + + @pytest.mark.parametrize('name', ['arcsin', 'arccos', 'arctan']) + def test_inverse_trig_returns_radians(self, name): + # WHEN + x = make(0.5) + + # THEN + result = fn(name)(x) + + # EXPECT + assert result.unit == 'rad' + + @pytest.mark.parametrize('name', ['arcsin', 'arccos', 'arctan']) + def test_inverse_trig_requires_dimensionless_input(self, name): + # WHEN THEN + x = make(0.5, 'rad') + + # EXPECT + with pytest.raises(UnitError): + fn(name)(x) + + @pytest.mark.parametrize( + 'name', + ['exp', 'log', 'log10', 'erf', 'erfc', 'sinh', 'cosh', 'tanh', 'arcsinh', 'arccosh'], + ) + def test_dimensionless_only_functions_reject_units(self, name): + # WHEN THEN + x = make(1.5, 'm') + + # EXPECT + with pytest.raises(UnitError): + fn(name)(x) + + def test_arctanh_rejects_units(self): + # WHEN THEN + x = make(0.5, 'm') + + # EXPECT + with pytest.raises(UnitError): + es_math.arctanh(x) + + def test_sinc_accepts_radians_and_dimensionless(self): + # WHEN + radians = make(0.5, 'rad') + dimensionless = make(0.5, '') + + # THEN + result_rad = es_math.sinc(radians) + result_dimensionless = es_math.sinc(dimensionless) + + # EXPECT + assert result_rad.value == pytest.approx(np.sin(0.5) / 0.5) + assert result_dimensionless.value == pytest.approx(np.sin(0.5) / 0.5) + + def test_sinc_rejects_degrees(self): + # WHEN THEN + angle = make(30.0, 'deg') + + # EXPECT + with pytest.raises(UnitError): + es_math.sinc(angle) + + def test_sqrt_of_a_square_unit(self): + # WHEN + area = make(4.0, 'm**2', min=1.0, max=9.0) + + # THEN + result = es_math.sqrt(area) + + # EXPECT + assert result.unit == 'm' + assert result.value == pytest.approx(2.0) + assert_bounds(result, 1.0, 3.0) + + def test_unit_error_wins_over_value_error(self): + # WHEN THEN + # Scipp validates units first, so the more fundamental error is reported. + x = make(-1.0, 'm') + + # EXPECT + with pytest.raises(UnitError): + es_math.log(x) + + +# ---------------------------------------------------------------------- +# Input types +# ---------------------------------------------------------------------- + + +class TestInputTypes: + @pytest.mark.parametrize('name', UNARY_FUNCTIONS) + def test_plain_numbers_are_not_supported(self, name): + # WHEN THEN + x = 0.5 + + # EXPECT + with pytest.raises(NotImplementedError): + fn(name)(x) + + def test_scipp_variable_is_not_supported(self): + # WHEN THEN + x = sc.scalar(0.5, unit='rad') + + # EXPECT + with pytest.raises(NotImplementedError): + es_math.sin(x) + + @pytest.mark.parametrize(('y', 'x'), [(1.0, 2.0), (1.0, make(2.0)), (make(1.0), 2.0)]) + def test_arctan2_requires_descriptor_numbers_or_parameters(self, y, x): + # WHEN THEN + # (y, x) contain at least one plain number + + # EXPECT + with pytest.raises(NotImplementedError): + es_math.arctan2(y, x) + + +# ---------------------------------------------------------------------- +# Invariants over random inputs +# ---------------------------------------------------------------------- + +# (value range, unit). Bounds are drawn around the value and may leave the domain. +RANDOM_RANGES = { + 'exp': ((-3.0, 3.0), ''), + 'log': ((0.1, 10.0), ''), + 'log10': ((0.1, 10.0), ''), + 'sqrt': ((0.1, 10.0), ''), + 'erf': ((-2.0, 2.0), ''), + 'erfc': ((-2.0, 2.0), ''), + 'sin': ((-10.0, 10.0), 'rad'), + 'cos': ((-10.0, 10.0), 'rad'), + 'tan': ((-1.4, 1.4), 'rad'), + 'arcsin': ((-0.9, 0.9), ''), + 'arccos': ((-0.9, 0.9), ''), + 'arctan': ((-5.0, 5.0), ''), + 'sinc': ((-10.0, 10.0), 'rad'), + 'sinh': ((-3.0, 3.0), ''), + 'cosh': ((-3.0, 3.0), ''), + 'tanh': ((-3.0, 3.0), ''), + 'arcsinh': ((-3.0, 3.0), ''), + 'arccosh': ((1.1, 10.0), ''), + 'arctanh': ((-0.9, 0.9), ''), +} + + +# Array-capable references where VALUE_CASES holds a scalar-only one (math.erf). +ARRAY_REFERENCE = {'erf': _erf, 'erfc': _erfc} + + +class TestInvariants: + @pytest.mark.parametrize('name', UNARY_FUNCTIONS) + def test_value_lies_within_result_bounds_and_bounds_are_not_nan(self, name): + # WHEN + (low, high), unit = RANDOM_RANGES[name] + rng = np.random.default_rng(seed=1234) + for _ in range(50): + value = rng.uniform(low, high) + lo = value - rng.uniform(0.01, 3.0) + hi = value + rng.uniform(0.01, 3.0) + x = make(value, unit, min=lo, max=hi, variance=0.001) + + # THEN + with warnings.catch_warnings(): + warnings.simplefilter('error') + result = fn(name)(x) + + # EXPECT + assert not np.isnan(result.min) and not np.isnan(result.max) + assert result.min <= result.value <= result.max + + @pytest.mark.parametrize('name', UNARY_FUNCTIONS) + def test_random_bounds_match_brute_force_range(self, name): + # WHEN + # Unlike the invariant above (which _apply enforces by widening the bounds), this + # compares against an independent dense evaluation, so a wrong strategy is caught. + (low, high), unit = RANDOM_RANGES[name] + reference = ARRAY_REFERENCE.get(name, VALUE_CASES[name][0]) + rng = np.random.default_rng(seed=2024) + for _ in range(25): + value = rng.uniform(low, high) + lo, hi = rng.uniform(low, value), rng.uniform(value, high) + if hi - lo < 0.05: + continue + x = make(value, unit, min=lo, max=hi) + + # THEN + result = fn(name)(x) + + # EXPECT + true_min, true_max = dense_range(reference, lo, hi, n=20_001) + assert result.min <= true_min + 1e-9 + assert result.max >= true_max - 1e-9 + assert result.min >= true_min - 1e-3 + assert result.max <= true_max + 1e-3 + + def test_arctan2_value_lies_within_result_bounds(self): + # WHEN + rng = np.random.default_rng(seed=99) + for _ in range(50): + y, x = rng.uniform(-3.0, 3.0, size=2) + y_lo, y_hi = y - rng.uniform(0.01, 2.0), y + rng.uniform(0.01, 2.0) + x_lo, x_hi = x - rng.uniform(0.01, 2.0), x + rng.uniform(0.01, 2.0) + y_param = make(y, min=y_lo, max=y_hi) + x_param = make(x, min=x_lo, max=x_hi) + + # THEN + with warnings.catch_warnings(): + warnings.simplefilter('ignore', UserWarning) + result = es_math.arctan2(y_param, x_param) + + # EXPECT + assert result.min <= result.value <= result.max + + +# ---------------------------------------------------------------------- +# Public API +# ---------------------------------------------------------------------- + + +class TestPublicApi: + @pytest.mark.parametrize('name', ALL_FUNCTIONS) + def test_function_is_exported_from_easyscience_variable(self, name): + # WHEN THEN + exported = getattr(es_variable, name) + + # EXPECT + assert callable(exported) + assert exported is getattr(es_math, name) + + @pytest.mark.parametrize('name', ALL_FUNCTIONS) + def test_function_has_a_docstring(self, name): + # WHEN + function = fn(name) + + # THEN + docstring = function.__doc__ + + # EXPECT + assert docstring and docstring.strip() + + def test_variable_all_contains_only_strings(self): + # WHEN + exported_names = es_variable.__all__ + + # THEN + all_are_strings = all(isinstance(entry, str) for entry in exported_names) + + # EXPECT + assert all_are_strings + for name in ALL_FUNCTIONS: + assert name in exported_names + + +# ---------------------------------------------------------------------- +# Output unit of the inverse trigonometric functions +# ---------------------------------------------------------------------- + +INVERSE_TRIG = ['arcsin', 'arccos', 'arctan'] + + +class TestInverseOutputUnit: + @pytest.mark.parametrize('name', INVERSE_TRIG) + def test_default_unit_is_radians(self, name): + # WHEN + x = make(0.4, min=0.2, max=0.8) + + # THEN + result = fn(name)(x) + + # EXPECT + assert result.unit == 'rad' + + @pytest.mark.parametrize('name', INVERSE_TRIG) + def test_explicit_radians_equals_the_default(self, name): + # WHEN + x = make(0.4, min=0.2, max=0.8, variance=0.04) + + # THEN + default = fn(name)(x) + explicit = fn(name)(x, unit='rad') + + # EXPECT + assert explicit.unit == 'rad' + assert (explicit.value, explicit.variance, explicit.min, explicit.max) == ( + default.value, + default.variance, + default.min, + default.max, + ) + + @pytest.mark.parametrize('name', INVERSE_TRIG) + def test_degrees_convert_value_and_bounds(self, name): + # WHEN + x = make(0.4, min=0.2, max=0.8) + in_radians = fn(name)(x) + + # THEN + in_degrees = fn(name)(x, unit='deg') + + # EXPECT + assert in_degrees.unit == 'deg' + assert in_degrees.value == pytest.approx(np.degrees(in_radians.value), rel=1e-12) + assert in_degrees.min == pytest.approx(np.degrees(in_radians.min), rel=1e-12) + assert in_degrees.max == pytest.approx(np.degrees(in_radians.max), rel=1e-12) + + @pytest.mark.parametrize('name', INVERSE_TRIG) + def test_degrees_convert_variance_with_the_squared_factor(self, name): + # WHEN + _, derivative, value, _ = VALUE_CASES[name] + variance = 0.04 + x = make(value, variance=variance) + + # THEN + result = fn(name)(x, unit='deg') + + # EXPECT + variance_in_rad2 = float(derivative(value)) ** 2 * variance + assert result.variance == pytest.approx(variance_in_rad2 * np.degrees(1.0) ** 2, rel=1e-9) + + def test_arcsin_example_in_degrees(self): + # WHEN + x = make(0.5, min=0.0, max=1.0) + + # THEN + result = es_math.arcsin(x, unit='deg') + + # EXPECT + assert result.value == pytest.approx(30.0, abs=1e-9) + assert_bounds(result, 0.0, 90.0) + + def test_arctan2_in_degrees(self): + # WHEN + y = make(1.0, 'm', min=0.5, max=2.0, name='y') + x = make(1.0, 'm', min=1.0, max=3.0, name='x') + + # THEN + result = es_math.arctan2(y, x, unit='deg') + + # EXPECT + assert result.unit == 'deg' + assert result.value == pytest.approx(45.0, abs=1e-9) + assert result.min == pytest.approx(np.degrees(np.arctan2(0.5, 3.0)), rel=1e-12) + assert result.max == pytest.approx(np.degrees(np.arctan2(2.0, 1.0)), rel=1e-12) + + def test_arctan2_variance_in_degrees(self): + # WHEN + y, x, var_y, var_x = 1.0, 2.0, 0.04, 0.09 + y_param = make(y, 'm', variance=var_y) + x_param = make(x, 'm', variance=var_x) + + # THEN + result = es_math.arctan2(y_param, x_param, unit='deg') + + # EXPECT + variance_in_rad2 = (x**2 * var_y + y**2 * var_x) / (x**2 + y**2) ** 2 + assert result.variance == pytest.approx(variance_in_rad2 * np.degrees(1.0) ** 2, rel=1e-9) + + @pytest.mark.parametrize('name', INVERSE_TRIG) + def test_a_unit_that_is_not_an_angle_raises(self, name): + # WHEN THEN + x = make(0.4) + + # EXPECT + with pytest.raises(UnitError): + fn(name)(x, unit='m') + + def test_arctan2_with_a_unit_that_is_not_an_angle_raises(self): + # WHEN THEN + y = make(1.0, 'm') + x = make(2.0, 'm') + + # EXPECT + with pytest.raises(UnitError): + es_math.arctan2(y, x, unit='m') + + @pytest.mark.parametrize('name', INVERSE_TRIG) + def test_unit_is_keyword_only(self, name): + # WHEN THEN + x = make(0.4) + + # EXPECT + with pytest.raises(TypeError): + fn(name)(x, 'deg') + + +# ---------------------------------------------------------------------- +# Aliases +# ---------------------------------------------------------------------- + +ALIASES = { + 'asin': 'arcsin', + 'acos': 'arccos', + 'atan': 'arctan', + 'atan2': 'arctan2', + 'asinh': 'arcsinh', + 'acosh': 'arccosh', + 'atanh': 'arctanh', +} + + +class TestAliases: + @pytest.mark.parametrize(('alias', 'canonical'), ALIASES.items()) + def test_alias_is_the_same_function_object(self, alias, canonical): + # WHEN THEN + alias_function = getattr(es_math, alias) + + # EXPECT + assert alias_function is getattr(es_math, canonical) + + @pytest.mark.parametrize(('alias', 'canonical'), ALIASES.items()) + def test_alias_is_exported_from_easyscience_variable(self, alias, canonical): + # WHEN THEN + exported = getattr(es_variable, alias) + + # EXPECT + assert exported is getattr(es_variable, canonical) + assert alias in es_variable.__all__ + + def test_result_is_named_after_the_canonical_function(self): + # WHEN + x = make(0.4, name='x') + + # THEN + result = es_variable.asin(x) + + # EXPECT + assert result.name == 'arcsin(x)' + + def test_alias_accepts_the_unit_keyword(self): + # WHEN + x = make(0.5, min=0.0, max=1.0) + + # THEN + result = es_variable.asin(x, unit='deg') + + # EXPECT + assert result.unit == 'deg' + assert result.value == pytest.approx(30.0, abs=1e-9) + + def test_atan2_alias_keeps_the_y_then_x_argument_order(self): + # WHEN + y = make(1.0, name='y') + x = make(-1.0, name='x') + + # THEN + result = es_variable.atan2(y, x) + + # EXPECT + assert result.value == pytest.approx(3 * np.pi / 4) + + +# ---------------------------------------------------------------------- +# Inheritance: Parameter and DescriptorNumber +# ---------------------------------------------------------------------- + + +class MyParameter(Parameter): + """A user subclass whose constructor needs an extra argument.""" + + def __init__(self, name, value, tag, **kwargs): + super().__init__(name, value, **kwargs) + self.tag = tag + + +class MyDescriptor(DescriptorNumber): + """A user subclass whose constructor needs an extra argument.""" + + def __init__(self, name, value, tag, **kwargs): + super().__init__(name, value, **kwargs) + self.tag = tag + + +def descriptor(value, unit='', *, variance=0.0, name='x') -> DescriptorNumber: + """Shorthand for a plain ``DescriptorNumber`` (no bounds).""" + return DescriptorNumber(name=name, value=value, unit=unit, variance=variance) + + +class TestInheritance: + @pytest.mark.parametrize('name', UNARY_FUNCTIONS) + def test_descriptor_number_input_returns_a_descriptor_number(self, name): + # WHEN + _, _, value, unit = VALUE_CASES[name] + plain = descriptor(value, unit, variance=0.04) + parameter = make(value, unit, variance=0.04) + + # THEN + result = fn(name)(plain) + reference = fn(name)(parameter) + + # EXPECT + assert type(result) is DescriptorNumber + assert not hasattr(result, 'min') + assert result.value == pytest.approx(reference.value, rel=1e-12) + assert result.variance == pytest.approx(reference.variance, rel=1e-12) + + @pytest.mark.parametrize('name', UNARY_FUNCTIONS) + def test_parameter_input_still_returns_a_parameter(self, name): + # WHEN + _, _, value, unit = VALUE_CASES[name] + x = make(value, unit) + + # THEN + result = fn(name)(x) + + # EXPECT + assert type(result) is Parameter + + def test_mixed_arctan2_returns_a_parameter_with_point_bounds(self): + # WHEN + y = make(1.0, min=0.5, max=2.0, name='y') + x = descriptor(2.0, name='x') + point_y = descriptor(1.0, name='y') + range_x = make(2.0, min=1.0, max=3.0, name='x') + + # THEN + result = es_math.arctan2(y, x) + swapped = es_math.arctan2(point_y, range_x) + + # EXPECT + assert type(result) is Parameter + assert_bounds(result, np.arctan2(0.5, 2.0), np.arctan2(2.0, 2.0)) + assert type(swapped) is Parameter + assert_bounds(swapped, np.arctan2(1.0, 3.0), np.arctan2(1.0, 1.0)) + + def test_arctan2_of_two_descriptor_numbers_returns_a_descriptor_number(self): + # WHEN + y = descriptor(1.0, name='y') + x = descriptor(1.0, name='x') + + # THEN + result = es_math.arctan2(y, x) + + # EXPECT + assert type(result) is DescriptorNumber + assert result.value == pytest.approx(np.pi / 4) + + def test_user_subclasses_are_accepted_and_return_base_classes(self): + # WHEN + sub_parameter = MyParameter('p', 0.5, 'tag', unit='rad', min=0.0, max=1.0) + sub_descriptor = MyDescriptor('d', 0.5, 'tag', unit='rad') + + # THEN + from_parameter = es_math.sin(sub_parameter) + from_descriptor = es_math.sin(sub_descriptor) + + # EXPECT + assert type(from_parameter) is Parameter + assert type(from_descriptor) is DescriptorNumber + assert from_parameter.value == pytest.approx(np.sin(0.5)) + assert from_descriptor.value == pytest.approx(np.sin(0.5)) + + def test_unit_keyword_on_a_descriptor_number(self): + # WHEN + x = descriptor(0.5, name='x') + + # THEN + result = es_math.arcsin(x, unit='deg') + + # EXPECT + assert type(result) is DescriptorNumber + assert result.unit == 'deg' + assert result.value == pytest.approx(30.0, abs=1e-9) + + def test_variance_in_degrees_for_a_descriptor_number(self): + # WHEN + angle = descriptor(30.0, 'deg', variance=1.0) + + # THEN + result = es_math.sin(angle) + + # EXPECT + expected = np.cos(np.deg2rad(30.0)) ** 2 * np.deg2rad(1.0) ** 2 + assert result.variance == pytest.approx(expected, rel=1e-9) + + @pytest.mark.parametrize( + ('name', 'value', 'unit', 'message'), + [ + ('tan', np.pi / 2, 'rad', 'asymptote'), + ('log', -1.0, '', 'outside the domain'), + ('exp', 1000.0, '', 'not finite'), + ], + ) + def test_guards_apply_to_a_descriptor_number(self, name, value, unit, message): + # WHEN THEN + x = descriptor(value, unit) + + # EXPECT + with pytest.raises(ValueError, match=message): + fn(name)(x) + + @pytest.mark.parametrize( + 'unsupported', + [ + 0.5, + np.array([0.5]), + sc.scalar(0.5, unit='rad'), + DescriptorArray(name='arr', value=[0.5], unit='rad', dimensions=['x']), + ], + ids=['float', 'ndarray', 'scipp-variable', 'descriptor-array'], + ) + def test_other_types_still_raise_not_implemented(self, unsupported): + # WHEN THEN + argument = unsupported + + # EXPECT + with pytest.raises(NotImplementedError, match='DescriptorNumber'): + es_math.sin(argument) + + +# ---------------------------------------------------------------------- +# The es.math namespace +# ---------------------------------------------------------------------- + + +class TestEsMathNamespace: + def test_es_math_is_the_easyscience_math_module(self): + # WHEN THEN + namespace = es.math + + # EXPECT + assert namespace.__name__ == 'easyscience.math' + + def test_es_math_is_not_the_standard_library_math(self): + # WHEN + import math as stdlib_math + + # THEN + namespace = es.math + + # EXPECT + assert stdlib_math is not namespace + assert stdlib_math.sin(0.0) == 0.0 + + def test_es_math_all_equals_the_implementation_all(self): + # WHEN THEN + public_names = es.math.__all__ + + # EXPECT + assert public_names == es_math.__all__ + assert len(public_names) == 27 + assert all(isinstance(name, str) for name in public_names) + + @pytest.mark.parametrize('name', es_math.__all__) + def test_es_math_exposes_the_same_function_objects(self, name): + # WHEN THEN + exposed = getattr(es.math, name) + + # EXPECT + assert exposed is getattr(es_math, name) + + def test_es_math_exposes_nothing_beyond_the_public_names(self): + # WHEN THEN + public_attributes = {name for name in dir(es.math) if not name.startswith('_')} + + # EXPECT + assert public_attributes == set(es_math.__all__) + + def test_math_is_not_exported_by_a_star_import_of_the_package(self): + # WHEN THEN + top_level_all = es.__all__ + + # EXPECT + # `from easyscience import *` would otherwise replace the standard-library `math`. + assert 'math' not in top_level_all + + def test_variable_still_exports_every_math_name(self): + # WHEN THEN + math_names = set(es_math.__all__) + + # EXPECT + assert math_names <= set(es_variable.__all__) + + def test_from_import_works(self): + # WHEN THEN + module = importlib.import_module('easyscience.math') + + # EXPECT + assert module.sin is es_math.sin + assert module.asin is es_math.arcsin + + +# ---------------------------------------------------------------------- +# Collapsed result bounds +# ---------------------------------------------------------------------- + +COLLAPSED_BOUNDS_CASES = [ + # (function, value, unit, min, max, why the exact result bounds collapse) + ('erf', 15.0, '', 10.0, 20.0, 'erf(10) == erf(20) == 1.0'), + ('erf', 5.5, '', 5.0, 6.0, 'erf(5) and erf(6) differ by 1.5e-12 (< rtol 1e-9)'), + ('tanh', 20.0, '', 15.0, 25.0, 'tanh saturates at 1.0'), + ('exp', -900.0, '', -1000.0, -800.0, 'both bounds underflow to 0.0'), + ('erfc', 35.0, '', 30.0, 40.0, 'both bounds underflow to 0.0'), + ('cos', 0.0, 'rad', -1e-5, 1e-5, 'cos is quadratic at its maximum: [1 - 5e-11, 1]'), + ('cosh', 0.0, '', -1e-5, 1e-5, 'cosh is quadratic at its minimum: [1, 1 + 5e-11]'), + ('sin', np.pi / 2, 'rad', np.pi / 2 - 1e-5, np.pi / 2 + 1e-5, 'quadratic at the maximum'), + ('arctan', 1e11, '', 1e10, 1e12, 'arctan saturates at pi/2'), +] +COLLAPSED_IDS = [f'{case[0]}[{case[3]:g},{case[4]:g}]' for case in COLLAPSED_BOUNDS_CASES] + + +class TestCollapsedBounds: + """ + ``Parameter`` rejects ``min`` and ``max`` that are equal or within ``rtol=1e-9``. Valid + inputs whose function saturates, underflows or is quadratic at an extremum would make the + exact result bounds collapse, so such bounds are padded around the value (by twice that + relative tolerance, or by an absolute 1e-12 if the value is 0). + """ + + @pytest.mark.parametrize( + ('name', 'value', 'unit', 'lo', 'hi', 'why'), COLLAPSED_BOUNDS_CASES, ids=COLLAPSED_IDS + ) + def test_collapsed_bounds_do_not_raise_and_contain_the_value( + self, name, value, unit, lo, hi, why + ): + # WHEN + x = make(value, unit, min=lo, max=hi) + + # THEN + result = fn(name)(x) + + # EXPECT + assert not np.isnan(result.min) and not np.isnan(result.max) + assert result.min <= result.value <= result.max + assert result.min < result.max + + @pytest.mark.parametrize( + ('name', 'value', 'unit', 'lo', 'hi', 'why'), COLLAPSED_BOUNDS_CASES, ids=COLLAPSED_IDS + ) + def test_padded_bounds_contain_the_true_range(self, name, value, unit, lo, hi, why): + # WHEN + x = make(value, unit, min=lo, max=hi) + reference = ARRAY_REFERENCE.get(name, VALUE_CASES[name][0]) + + # THEN + result = fn(name)(x) + + # EXPECT + true_min, true_max = dense_range(reference, lo, hi, n=20_001) + assert result.min <= true_min + assert result.max >= true_max + + @pytest.mark.parametrize( + ('name', 'value', 'unit', 'lo', 'hi', 'why'), COLLAPSED_BOUNDS_CASES, ids=COLLAPSED_IDS + ) + def test_padding_is_tiny(self, name, value, unit, lo, hi, why): + # WHEN + x = make(value, unit, min=lo, max=hi) + + # THEN + result = fn(name)(x) + + # EXPECT + # Twice the 1e-9 relative tolerance on each side, i.e. a width of about 4e-9 * |value| + # (an absolute 2e-12 around a value of 0). + assert result.max - result.min <= 1e-8 * max(abs(result.value), 1e-3) + + def test_zero_value_gets_an_absolute_pad(self): + # WHEN + # exp(-900) underflows to exactly 0.0, so a relative pad would again be zero. + x = make(-900.0, min=-1000.0, max=-800.0) + + # THEN + result = es_math.exp(x) + + # EXPECT + assert result.value == 0.0 + assert_bounds(result, -1e-12, 1e-12, tol=1e-15) + + def test_bounds_that_do_not_collapse_are_left_exact(self): + # WHEN + x = make(0.5, min=0.0, max=1.0) + + # THEN + result = es_math.erf(x) + + # EXPECT + assert_bounds(result, 0.0, math.erf(1.0), tol=1e-15)