diff --git a/solarwindpy/fitfunctions/CONTRIBUTING.md b/solarwindpy/fitfunctions/CONTRIBUTING.md new file mode 100644 index 00000000..8e84c717 --- /dev/null +++ b/solarwindpy/fitfunctions/CONTRIBUTING.md @@ -0,0 +1,374 @@ +# Contributing to fitfunctions + +This document defines the standards, conventions, and quality requirements for contributing +to the `solarwindpy.fitfunctions` module. It is standalone and will be integrated into +unified project documentation once all submodules have contribution standards. + +## 1. Overview + +The `fitfunctions` module provides a framework for fitting mathematical models to data +using `scipy.optimize.curve_fit`. Each fit function is a class that inherits from +`FitFunction` and implements three required abstract properties. + +**Key files:** +- `core.py` - Base `FitFunction` class and exceptions +- `hinge.py`, `lines.py`, `gaussians.py`, etc. - Concrete implementations +- `tests/fitfunctions/` - Test suite + +## 2. Development Workflow (TDD) + +Follow Test-Driven Development with separate commits for tests and implementation: + +``` +1. Requirements → What does this function model? What are the parameters? +2. Test Writing → Commit: test(fitfunctions): add tests for +3. Implementation → Commit: feat(fitfunctions): add +4. Verification → All tests pass, including existing tests +``` + +**Commit order matters:** Tests are committed before implementation. This documents the +expected behavior and ensures tests are not written to pass existing code. + +## 3. FitFunction Class Requirements + +### 3.1 Required Abstract Properties + +Every `FitFunction` subclass MUST implement these three properties: + +| Property | Returns | Purpose | +|----------|---------|---------| +| `function` | callable | The mathematical function `f(x, *params)` to fit | +| `p0` | list | Initial parameter guesses (data-driven) | +| `TeX_function` | str | LaTeX representation for plotting | + +**Minimal implementation:** + +```python +from .core import FitFunction + +class MyFunction(FitFunction): + r"""One-line description. + + Extended description with math: + + .. math:: + + f(x) = m \cdot x + b + + Parameters + ---------- + xobs : array-like + Independent variable observations. + yobs : array-like + Dependent variable observations. + **kwargs + Additional arguments passed to :class:`FitFunction`. + """ + + @property + def function(self): + def my_func(x, m, b): + return m * x + b + return my_func + + @property + def p0(self) -> list: + assert self.sufficient_data + x = self.observations.used.x + y = self.observations.used.y + # Data-driven estimation (see §3.2) + m = (y[-1] - y[0]) / (x[-1] - x[0]) + b = y[0] - m * x[0] + return [m, b] + + @property + def TeX_function(self) -> str: + return r"f(x) = m \cdot x + b" +``` + +### 3.2 p0 Estimation (Data-Driven) + +Initial parameter guesses MUST be data-driven. Hardcoded domain values are prohibited. + +**REQUIRED pattern:** + +```python +@property +def p0(self) -> list: + assert self.sufficient_data + x = self.observations.used.x + y = self.observations.used.y + + # Data-driven estimation examples: + x0 = (x.max() + x.min()) / 2 # Midpoint for transitions + y0 = np.median(y[x > x0]) # Baseline from data + m = np.polyfit(x[:10], y[:10], 1)[0] # Slope from segment + A = y.max() - y.min() # Amplitude from range + + return [x0, y0, m, A] +``` + +**PROHIBITED (hardcoded values):** + +```python +# BAD: Domain-specific hardcoded values +x0 = 425 # Solar wind speed +m1 = 0.0163 # Kasper 2007 value +``` + +**Why data-driven?** +- Works on arbitrary datasets, not just solar wind +- Enables reuse across scientific domains +- Reduces "magic number" bugs + +### 3.3 Optional Overrides + +**Custom `__init__` with guess parameters:** + +```python +def __init__( + self, + xobs, + yobs, + guess_x0: float | None = None, # Optional user hint + **kwargs, +): + self._guess_x0 = guess_x0 + super().__init__(xobs, yobs, **kwargs) +``` + +**Derived properties:** + +```python +@property +def xs(self) -> float: + """Saturation x-coordinate (derived from fitted params).""" + return self.popt["x1"] + self.popt["yh"] / self.popt["m1"] +``` + +### 3.4 Code Conventions + +| Convention | Standard | Example | +|------------|----------|---------| +| Class names | PascalCase | `HingeSaturation`, `GaussianPlusHeavySide` | +| Method names | snake_case | `make_fit()`, `build_plotter()` | +| Property names | snake_case | `popt`, `TeX_function` | +| Docstrings | NumPy style with `r"""` | See example above | +| Type hints | Selective (params with defaults) | `guess_x0: float = None` | +| Imports | Relative, future annotations | `from .core import FitFunction` | +| LaTeX | Raw strings | `r"$\chi^2_\nu$"` | + +**Import template:** + +```python +r"""Module docstring.""" + +from __future__ import annotations + +import numpy as np + +from .core import FitFunction +``` + +## 4. Test Requirements + +### 4.1 Test Categories (E1-E7) + +Every FitFunction MUST have tests in categories E1-E5. E6-E7 are recommended where applicable. + +| Category | Purpose | Tolerance | Required? | +|----------|---------|-----------|-----------| +| E1. Function Evaluation | Verify exact f(x) values | `rtol=1e-10` | YES | +| E2. Parameter Recovery (Clean) | Fit recovers known params | `rel_error < 2%` | YES | +| E3. Parameter Recovery (Noisy) | Statistical precision | `deviation < 2σ` | YES | +| E4. Initial Parameter (p0) | p0 enables convergence | `isfinite(popt)` | YES | +| E5. Edge Cases | Error handling | `raises Exception` | YES | +| E6. Derived Properties | Internal consistency | `rtol=1e-6` | If applicable | +| E7. Behavioral | Continuity, transitions | `rtol=0.1` | If applicable | + +### 4.2 Fixture Pattern + +All fixtures MUST return `(x, y, w, true_params)`: + +```python +@pytest.fixture +def clean_gaussian_data(): + """Clean Gaussian data with known parameters.""" + true_params = {"mu": 5.0, "sigma": 1.0, "A": 10.0} + x = np.linspace(0, 10, 200) + y = gaussian(x, **true_params) + w = np.ones_like(x) + return x, y, w, true_params + + +@pytest.fixture +def noisy_gaussian_data(): + """Noisy Gaussian data with known parameters.""" + rng = np.random.default_rng(42) # Deterministic seed + true_params = {"mu": 5.0, "sigma": 1.0, "A": 10.0} + noise_std = 0.5 # 5% of amplitude + + x = np.linspace(0, 10, 200) + y_true = gaussian(x, **true_params) + y = y_true + rng.normal(0, noise_std, len(x)) + w = np.ones_like(x) / noise_std + + return x, y, w, true_params +``` + +**Conventions:** +- Random seed: `np.random.default_rng(42)` for reproducibility +- Noise level: 3-5% of signal amplitude +- Weights: `w = np.ones_like(x) / noise_std` for noisy data + +### 4.3 Assertion Patterns + +**REQUIRED: Use `np.testing.assert_allclose` with `err_msg`:** + +```python +np.testing.assert_allclose( + result, expected, rtol=0.02, + err_msg=f"param: fitted={result:.4f}, expected={expected:.4f}" +) +``` + +**Tolerance Reference:** + +| Test Type | Tolerance | Justification | +|-----------|-----------|---------------| +| Exact math (E1) | `rtol=1e-10` | Floating point precision | +| Clean fitting (E2) | `rel_error < 0.02` | curve_fit convergence | +| Noisy fitting (E3) | `deviation < 2*sigma` | 95% confidence interval | +| Derived quantities (E6) | `rtol=1e-6` | Computed from fitted params | +| Behavioral (E7) | `rtol=0.1` | Approximate behavior | + +### 4.4 Test Parameterization (REQUIRED for multi-case tests) + +Use `@pytest.mark.parametrize` to avoid code duplication: + +**Pattern 1: Multiple parameter sets** + +```python +@pytest.mark.parametrize( + "true_params", + [ + {"mu": 5.0, "sigma": 1.0, "A": 10.0}, # Standard case + {"mu": 0.0, "sigma": 0.5, "A": 5.0}, # Edge: mu at origin + {"mu": 10.0, "sigma": 2.0, "A": 1.0}, # Edge: small amplitude + ], + ids=["standard", "mu_at_origin", "small_amplitude"], +) +def test_gaussian_recovers_parameters(true_params): + """Test parameter recovery across configurations.""" + # Single test logic, multiple cases +``` + +**Pattern 2: Multiple classes** + +```python +@pytest.mark.parametrize("cls", [Line, LineXintercept]) +def test_make_fit_success(cls, simple_linear_data): + """All line classes should fit successfully.""" + x, y, w, _ = simple_linear_data + fit = cls(x, y) + fit.make_fit() + assert np.isfinite(fit.popt['m']) +``` + +**Best practices:** +- Use `ids=` for readable test names +- Use dict for multiple parameters +- Document each case with comments +- Include edge cases + +## 5. Test Patterns (DO THIS) + +| Pattern | Purpose | Example | +|---------|---------|---------| +| Analytic expected values | Verify math is correct | `expected = m * x + b` | +| Known parameter recovery | Verify fitting works | Generate data, fit, compare | +| Statistical bounds | Handle noise properly | `assert deviation < 2 * sigma` | +| Boundary conditions | Verify edge behavior | Test at x=x0 for step functions | +| Derived property consistency | Verify internal math | `assert m2 == (y2-y1)/(x2-x1)` | +| Documented tolerances | Explain precision | `rtol=0.02 # curve_fit convergence` | +| Error messages with context | Enable debugging | `err_msg=f"fitted={x:.3f}"` | +| Deterministic random seeds | Reproducibility | `rng = np.random.default_rng(42)` | + +## 6. Test Anti-Patterns (DO NOT DO THIS) + +| Anti-Pattern | Why It's Bad | Good Alternative | +|--------------|--------------|------------------| +| `assert fit.popt is not None` | Proves nothing about correctness | `assert_allclose(fit.popt['m'], 2.0, rtol=0.02)` | +| `assert isinstance(fit.popt, dict)` | Verifies structure, not behavior | Verify actual parameter values | +| `assert len(fit.popt) == 3` | Trivial, no math validation | Verify each parameter value | +| `rtol=0.1 # works` | Unexplained, arbitrary | `rtol=0.02 # curve_fit convergence` | +| `assert a == b` (no message) | Hard to debug failures | `assert a == b, f"got {a}, expected {b}"` | +| `rtol=1e-15` for noisy data | Flaky tests | `rtol=0.02` for fitting | +| Only test clean data | Misses real-world behavior | Include noisy data with 2σ bounds | +| `np.random.normal(...)` | Non-reproducible failures | `rng.normal(...)` with fixed seed | + +## 7. Non-Trivial Test Criteria + +**Definition:** A test is **non-trivial** if it would FAIL on a plausible incorrect implementation. + +Every test MUST satisfy ALL of these criteria: + +| Criterion | Requirement | Anti-Example | +|-----------|-------------|--------------| +| Numeric assertion | Uses `assert_allclose` with explicit tolerance | `assert popt is not None` | +| Known expected value | Expected value is analytically computed | `assert result` (truthy) | +| Justified tolerance | rtol/atol documented with reasoning | `rtol=0.1 # seems to work` | +| Failure diagnostic | Error message shows actual vs expected | Bare `AssertionError` | +| Mathematical meaning | Tests a model property, not structure | `assert len(popt) == 4` | +| Would fail if broken | A plausible bug would cause failure | Test that always passes | + +**Example non-trivial test:** + +```python +def test_line_evaluates_correctly(): + """Line: f(x) = m*x + b should give exact values. + + Non-trivial because: + - Tests specific numeric values, not just "runs" + - Would fail if formula were m*x - b (sign error) + - Tolerance is 1e-10 (floating point, not fitting) + """ + m, b = 2.0, 1.0 + x = np.array([0.0, 1.0, 2.0]) + expected = np.array([1.0, 3.0, 5.0]) # Analytically computed + + fit = Line(x, expected) + result = fit.function(x, m, b) + + np.testing.assert_allclose( + result, expected, rtol=1e-10, + err_msg=f"Line(x, m={m}, b={b}) should equal m*x + b" + ) +``` + +## 8. Quality Checklist + +Before submitting a PR, verify: + +- [ ] All 3 abstract properties implemented (`function`, `p0`, `TeX_function`) +- [ ] p0 is data-driven (no hardcoded domain values) +- [ ] Tests cover categories E1-E5 minimum +- [ ] All tests are non-trivial (pass criteria in §7) +- [ ] Docstrings complete with `.. math::` blocks +- [ ] `make_fit()` converges on clean data +- [ ] `make_fit()` within 2σ on noisy data +- [ ] Class exported in `__init__.py` +- [ ] No regressions (all existing tests pass) + +## 9. Complete Examples + +For complete implementation examples, see: + +- **Implementation:** `hinge.py:HingeSaturation` (lines 20-180) +- **Tests:** `tests/fitfunctions/test_hinge.py:TestHingeSaturation` + +For test organization patterns: + +- **Parameterization:** `tests/fitfunctions/test_composite.py` (lines 359-365) +- **Fixtures:** `tests/fitfunctions/test_hinge.py` (fixture definitions) +- **Categories E1-E7:** `tests/fitfunctions/test_heaviside.py` (section headers) diff --git a/solarwindpy/fitfunctions/__init__.py b/solarwindpy/fitfunctions/__init__.py index f0e0a5a8..6311a5c1 100644 --- a/solarwindpy/fitfunctions/__init__.py +++ b/solarwindpy/fitfunctions/__init__.py @@ -7,8 +7,10 @@ from . import power_laws from . import moyal -# from . import hinge +from . import hinge +from . import heaviside from . import trend_fits +from . import composite FitFunction = core.FitFunction Gaussian = gaussians.Gaussian @@ -16,8 +18,17 @@ Line = lines.Line PowerLaw = power_laws.PowerLaw Moyal = moyal.Moyal -# Hinge = hinge.Hinge +HingeSaturation = hinge.HingeSaturation +TwoLine = hinge.TwoLine +Saturation = hinge.Saturation +HingeMin = hinge.HingeMin +HingeMax = hinge.HingeMax +HingeAtPoint = hinge.HingeAtPoint +HeavySide = heaviside.HeavySide TrendFit = trend_fits.TrendFit +GaussianPlusHeavySide = composite.GaussianPlusHeavySide +GaussianTimesHeavySide = composite.GaussianTimesHeavySide +GaussianTimesHeavySidePlusHeavySide = composite.GaussianTimesHeavySidePlusHeavySide # Exception classes for better error handling FitFunctionError = core.FitFunctionError diff --git a/solarwindpy/fitfunctions/composite.py b/solarwindpy/fitfunctions/composite.py new file mode 100644 index 00000000..1dc3d8cd --- /dev/null +++ b/solarwindpy/fitfunctions/composite.py @@ -0,0 +1,559 @@ +r"""Composite fit functions combining Gaussians with Heaviside step functions. + +This module provides fit functions that combine Gaussian distributions with +Heaviside step functions for modeling distributions with sharp transitions +or truncations. + +Classes +------- +GaussianPlusHeavySide + Gaussian peak with additive step function offset. +GaussianTimesHeavySide + Gaussian truncated at a threshold (zero below x0). +GaussianTimesHeavySidePlusHeavySide + Gaussian for x >= x0 with constant plateau for x < x0. +""" + +from __future__ import annotations + +import numpy as np + +from .core import FitFunction + + +class GaussianPlusHeavySide(FitFunction): + r"""Gaussian plus Heaviside step function for offset peak modeling. + + Models a Gaussian peak with a step function that adds an offset + below a threshold: + + .. math:: + + f(x) = A \cdot e^{-\frac{1}{2}\left(\frac{x-\mu}{\sigma}\right)^2} + + y_1 \cdot H(x_0 - x) + y_0 + + where :math:`H(z)` is the Heaviside step function. + + Parameters + ---------- + xobs : array-like + Independent variable observations. + yobs : array-like + Dependent variable observations. + **kwargs + Additional arguments passed to :class:`FitFunction`. + + Attributes + ---------- + x0 : float + Transition x-coordinate (fitted parameter). + y0 : float + Constant offset applied everywhere (fitted parameter). + y1 : float + Additional offset for x < x0 (fitted parameter). + mu : float + Gaussian mean (fitted parameter). + sigma : float + Gaussian standard deviation (fitted parameter). + A : float + Gaussian amplitude (fitted parameter). + + Examples + -------- + >>> import numpy as np + >>> from solarwindpy.fitfunctions import GaussianPlusHeavySide + >>> x = np.linspace(0, 10, 100) + >>> # Gaussian peak at mu=5 with step down at x0=2 + >>> y = 4*np.exp(-0.5*((x-5)/1)**2) + 3*np.heaviside(2-x, 0.5) + 1 + >>> fit = GaussianPlusHeavySide(x, y) + >>> fit.make_fit() + >>> print(f"mu={fit.popt['mu']:.2f}, x0={fit.popt['x0']:.2f}") + mu=5.00, x0=2.00 + """ + + def __init__(self, xobs, yobs, **kwargs): + super().__init__(xobs, yobs, **kwargs) + + @property + def function(self): + r"""The Gaussian plus Heaviside function. + + Returns + ------- + callable + Function with signature ``f(x, x0, y0, y1, mu, sigma, A)``. + """ + + def gaussian_heavy_side(x, x0, y0, y1, mu, sigma, A): + r"""Evaluate Gaussian plus Heaviside model. + + Parameters + ---------- + x : array-like + Independent variable values. + x0 : float + Transition x-coordinate. + y0 : float + Constant offset everywhere. + y1 : float + Additional offset for x < x0. + mu : float + Gaussian mean. + sigma : float + Gaussian standard deviation. + A : float + Gaussian amplitude. + + Returns + ------- + numpy.ndarray + Model values at x. + """ + + def gaussian(x, mu, sigma, A): + arg = -0.5 * (((x - mu) / sigma) ** 2.0) + return A * np.exp(arg) + + def heavy_side(x, x0, y0, y1): + out = (y1 * np.heaviside(x0 - x, 0.5)) + y0 + return out + + out = gaussian(x, mu, sigma, A) + heavy_side(x, x0, y0, y1) + return out + + return gaussian_heavy_side + + @property + def p0(self) -> list: + r"""Calculate initial parameter guess. + + The initial guess is derived from: + - Weighted mean and std for Gaussian parameters + - x0 estimated as 0.75 * mean + - y0 = 0, y1 = 0.8 * peak + - Gaussian parameters recalculated for x > x0 + + Returns + ------- + list + Initial guesses as [x0, y0, y1, mu, sigma, A]. + + Raises + ------ + AssertionError + If insufficient data for estimation. + + Notes + ----- + # TODO: Convert to data-driven p0 estimation (see GH issue #XX) + """ + assert self.sufficient_data + + x, y = self.observations.used.x, self.observations.used.y + mean = (x * y).sum() / y.sum() + std = np.sqrt(((x - mean) ** 2.0 * y).sum() / y.sum()) + + try: + peak = y.max() + except ValueError as e: + chk = ( + r"zero-size array to reduction operation maximum " + "which has no identity" + ) + if str(e).startswith(chk): + msg = ( + "There is no maximum of a zero-size array. " + "Please check input data." + ) + raise ValueError(msg) + raise + + x0 = 0.75 * mean + y1 = 0.8 * peak + + tk = x > x0 + x, y = x[tk], y[tk] + mean = (x * y).sum() / y.sum() + std = np.sqrt(((x - mean) ** 2.0 * y).sum() / y.sum()) + + p0 = [x0, 0, y1, mean, std, peak] + return p0 + + @property + def TeX_function(self) -> str: + r"""LaTeX representation of the model. + + Returns + ------- + str + Multi-line LaTeX string describing the function. + """ + tex = "\n".join( + [ + r"f(x)=A \cdot e^{-\frac{1}{2} (\frac{x-\mu}{\sigma})^2} +", + r"\left(y_1 \cdot H(x_0 - x) + y_0\right)", + ] + ) + return tex + + +class GaussianTimesHeavySide(FitFunction): + r"""Gaussian multiplied by Heaviside for truncated distribution modeling. + + Models a Gaussian that is truncated (zeroed) below a threshold: + + .. math:: + + f(x) = A \cdot e^{-\frac{1}{2}\left(\frac{x-\mu}{\sigma}\right)^2} + \cdot H(x - x_0) + + where :math:`H(z)` is the Heaviside step function with :math:`H(0) = 1`. + + Parameters + ---------- + xobs : array-like + Independent variable observations. + yobs : array-like + Dependent variable observations. + guess_x0 : float, optional + Initial guess for the transition x-coordinate. If None, must be + provided for fitting to work properly. + **kwargs + Additional arguments passed to :class:`FitFunction`. + + Attributes + ---------- + x0 : float + Transition x-coordinate (fitted parameter). + mu : float + Gaussian mean (fitted parameter). + sigma : float + Gaussian standard deviation (fitted parameter). + A : float + Gaussian amplitude (fitted parameter). + + Examples + -------- + >>> import numpy as np + >>> from solarwindpy.fitfunctions import GaussianTimesHeavySide + >>> x = np.linspace(0, 10, 100) + >>> # Truncated Gaussian: zero for x < 3 + >>> y = 4*np.exp(-0.5*((x-5)/1)**2) * np.heaviside(x-3, 1.0) + >>> fit = GaussianTimesHeavySide(x, y, guess_x0=3.0) + >>> fit.make_fit() + >>> print(f"x0={fit.popt['x0']:.2f}, mu={fit.popt['mu']:.2f}") + x0=3.00, mu=5.00 + """ + + def __init__(self, xobs, yobs, guess_x0=None, **kwargs): + if guess_x0 is None: + raise ValueError( + "guess_x0 is required for GaussianTimesHeavySide. " + "Provide an initial estimate for the transition x-coordinate." + ) + super().__init__(xobs, yobs, **kwargs) + self._guess_x0 = guess_x0 + + @property + def guess_x0(self) -> float: + r"""Initial guess for transition x-coordinate used in p0 calculation.""" + return self._guess_x0 + + @property + def function(self): + r"""The Gaussian times Heaviside function. + + Returns + ------- + callable + Function with signature ``f(x, x0, mu, sigma, A)``. + """ + + def gaussian_heavy_side(x, x0, mu, sigma, A): + r"""Evaluate Gaussian times Heaviside model. + + Parameters + ---------- + x : array-like + Independent variable values. + x0 : float + Transition x-coordinate. + mu : float + Gaussian mean. + sigma : float + Gaussian standard deviation. + A : float + Gaussian amplitude. + + Returns + ------- + numpy.ndarray + Model values at x. + """ + + def gaussian(x, mu, sigma, A): + arg = -0.5 * (((x - mu) / sigma) ** 2.0) + return A * np.exp(arg) + + out = gaussian(x, mu, sigma, A) * np.heaviside(x - x0, 1.0) + return out + + return gaussian_heavy_side + + @property + def p0(self) -> list: + r"""Calculate initial parameter guess. + + The initial guess is derived from: + - ``guess_x0`` for x0 + - Weighted mean and std for Gaussian parameters (using x > x0) + - Peak amplitude from data + + Returns + ------- + list + Initial guesses as [x0, mu, sigma, A]. + + Raises + ------ + AssertionError + If insufficient data for estimation. + + Notes + ----- + # TODO: Convert to data-driven p0 estimation (see GH issue #XX) + """ + assert self.sufficient_data + + x, y = self.observations.used.x, self.observations.used.y + x0 = self.guess_x0 + tk = x > x0 + + x, y = x[tk], y[tk] + mean = (x * y).sum() / y.sum() + std = np.sqrt(((x - mean) ** 2.0 * y).sum() / y.sum()) + + try: + peak = y.max() + except ValueError as e: + chk = ( + r"zero-size array to reduction operation maximum " + "which has no identity" + ) + if str(e).startswith(chk): + msg = ( + "There is no maximum of a zero-size array. " + "Please check input data." + ) + raise ValueError(msg) + raise + + p0 = [x0, mean, std, peak] + return p0 + + @property + def TeX_function(self) -> str: + r"""LaTeX representation of the model. + + Returns + ------- + str + LaTeX string describing the function. + """ + tex = ( + r"f(x)=A \cdot e^{-\frac{1}{2} (\frac{x-\mu}{\sigma})^2} \times H(x - x_0)" + ) + return tex + + +class GaussianTimesHeavySidePlusHeavySide(FitFunction): + r"""Gaussian times Heaviside plus Heaviside for plateau-to-peak modeling. + + Models a distribution that transitions from a constant plateau to a + Gaussian peak: + + .. math:: + + f(x) = A \cdot e^{-\frac{1}{2}\left(\frac{x-\mu}{\sigma}\right)^2} + \cdot H(x - x_0) + y_1 \cdot H(x_0 - x) + + where :math:`H(z)` is the Heaviside step function. + + This gives: + - Constant :math:`y_1` for :math:`x < x_0` + - Gaussian for :math:`x > x_0` + - Transition at :math:`x = x_0` + + Parameters + ---------- + xobs : array-like + Independent variable observations. + yobs : array-like + Dependent variable observations. + guess_x0 : float, optional + Initial guess for the transition x-coordinate. If None, must be + provided for fitting to work properly. + **kwargs + Additional arguments passed to :class:`FitFunction`. + + Attributes + ---------- + x0 : float + Transition x-coordinate (fitted parameter). + y1 : float + Constant plateau value for x < x0 (fitted parameter). + mu : float + Gaussian mean (fitted parameter). + sigma : float + Gaussian standard deviation (fitted parameter). + A : float + Gaussian amplitude (fitted parameter). + + Examples + -------- + >>> import numpy as np + >>> from solarwindpy.fitfunctions import GaussianTimesHeavySidePlusHeavySide + >>> x = np.linspace(0, 10, 100) + >>> # Plateau at y=2 for x<3, then Gaussian peak + >>> y = 4*np.exp(-0.5*((x-5)/1)**2)*np.heaviside(x-3, 0.5) + 2*np.heaviside(3-x, 0.5) + >>> fit = GaussianTimesHeavySidePlusHeavySide(x, y, guess_x0=3.0) + >>> fit.make_fit() + >>> print(f"x0={fit.popt['x0']:.2f}, y1={fit.popt['y1']:.2f}") + x0=3.00, y1=2.00 + """ + + def __init__(self, xobs, yobs, guess_x0=None, **kwargs): + if guess_x0 is None: + raise ValueError( + "guess_x0 is required for GaussianTimesHeavySidePlusHeavySide. " + "Provide an initial estimate for the transition x-coordinate." + ) + super().__init__(xobs, yobs, **kwargs) + self._guess_x0 = guess_x0 + + @property + def guess_x0(self) -> float: + r"""Initial guess for transition x-coordinate used in p0 calculation.""" + return self._guess_x0 + + @property + def function(self): + r"""The Gaussian times Heaviside plus Heaviside function. + + Returns + ------- + callable + Function with signature ``f(x, x0, y1, mu, sigma, A)``. + """ + + def gaussian_heavy_side(x, x0, y1, mu, sigma, A): + r"""Evaluate Gaussian times Heaviside plus Heaviside model. + + Parameters + ---------- + x : array-like + Independent variable values. + x0 : float + Transition x-coordinate. + y1 : float + Constant plateau value for x < x0. + mu : float + Gaussian mean. + sigma : float + Gaussian standard deviation. + A : float + Gaussian amplitude. + + Returns + ------- + numpy.ndarray + Model values at x. + """ + + def gaussian(x, mu, sigma, A): + arg = -0.5 * (((x - mu) / sigma) ** 2.0) + return A * np.exp(arg) + + def heavy_side(x, x0, y1): + out = y1 * np.heaviside(x0 - x, 1.0) + return out + + out = gaussian(x, mu, sigma, A) * np.heaviside(x - x0, 1.0) + heavy_side( + x, x0, y1 + ) + return out + + return gaussian_heavy_side + + @property + def p0(self) -> list: + r"""Calculate initial parameter guess. + + The initial guess is derived from: + - ``guess_x0`` for x0 + - Mean of y values for x < x0 as y1 estimate + - Weighted mean and std for Gaussian parameters (using x > x0) + - Peak amplitude from data + + Returns + ------- + list + Initial guesses as [x0, y1, mu, sigma, A]. + + Raises + ------ + AssertionError + If insufficient data for estimation. + + Notes + ----- + # TODO: Convert to data-driven p0 estimation (see GH issue #XX) + """ + assert self.sufficient_data + + x, y = self.observations.used.x, self.observations.used.y + x0 = self.guess_x0 + tk = x > x0 + + y1 = y[~tk].mean() + if np.isnan(y1): + y1 = 0 + + x, y = x[tk], y[tk] + mean = (x * y).sum() / y.sum() + std = np.sqrt(((x - mean) ** 2.0 * y).sum() / y.sum()) + + try: + peak = y.max() + except ValueError as e: + chk = ( + r"zero-size array to reduction operation maximum " + "which has no identity" + ) + if str(e).startswith(chk): + msg = ( + "There is no maximum of a zero-size array. " + "Please check input data." + ) + raise ValueError(msg) + raise + + p0 = [x0, y1, mean, std, peak] + return p0 + + @property + def TeX_function(self) -> str: + r"""LaTeX representation of the model. + + Returns + ------- + str + Multi-line LaTeX string describing the function. + """ + tex = "\n".join( + [ + r"f(x)=A \cdot e^{-\frac{1}{2} (\frac{x-\mu}{\sigma})^2} \times H(x - x_0) + ", + r"y_1 \cdot H(x_0 - x)", + ] + ) + return tex diff --git a/solarwindpy/fitfunctions/core.py b/solarwindpy/fitfunctions/core.py index 64cae010..cb5b2009 100644 --- a/solarwindpy/fitfunctions/core.py +++ b/solarwindpy/fitfunctions/core.py @@ -810,7 +810,7 @@ def make_fit(self, return_exception=False, **kwargs): try: res, p0 = self._run_least_squares(**kwargs) - except (RuntimeError, ValueError) as e: + except (RuntimeError, ValueError, FitFailedError) as e: # print("fitting failed", flush=True) # raise if return_exception: diff --git a/solarwindpy/fitfunctions/heaviside.py b/solarwindpy/fitfunctions/heaviside.py new file mode 100644 index 00000000..8655ebeb --- /dev/null +++ b/solarwindpy/fitfunctions/heaviside.py @@ -0,0 +1,184 @@ +r"""Heaviside step function fit. + +This module provides a fit function for a Heaviside (step) function, +commonly used for modeling abrupt transitions in data. +""" + +from __future__ import annotations + +import numpy as np + +from .core import FitFunction + + +class HeavySide(FitFunction): + r"""Heaviside step function for modeling abrupt transitions. + + The model is a step function with transition at x0: + + .. math:: + + f(x) = y_1 \cdot H(x_0 - x, \tfrac{1}{2}(y_0 + y_1)) + y_0 + + where H is the Heaviside step function. The behavior is: + + - For x < x0: :math:`f(x) = y_1 + y_0` + - For x > x0: :math:`f(x) = y_0` + - For x == x0: :math:`f(x) = y_1 \cdot \tfrac{1}{2}(y_0 + y_1) + y_0` + + Parameters + ---------- + xobs : array-like + Independent variable observations. + yobs : array-like + Dependent variable observations. + guess_x0 : float, optional + Initial guess for transition x-coordinate. If not provided, + estimated as the midpoint of the x data range. + guess_y0 : float, optional + Initial guess for baseline level (value for x > x0). If not + provided, estimated from data above the transition. + guess_y1 : float, optional + Initial guess for step height. If not provided, estimated + from data below and above the transition. + **kwargs + Additional arguments passed to :class:`FitFunction`. + + Attributes + ---------- + x0 : float + Step transition x-coordinate (fitted parameter). + y0 : float + Baseline level for x > x0 (fitted parameter). + y1 : float + Step height (fitted parameter). The value for x < x0 is y0 + y1. + + Examples + -------- + >>> import numpy as np + >>> from solarwindpy.fitfunctions import HeavySide + >>> x = np.linspace(0, 10, 100) + >>> y = np.where(x < 5, 5, 2) # Step down at x=5 + >>> fit = HeavySide(x, y, guess_x0=5, guess_y0=2, guess_y1=3) + >>> fit.make_fit() + >>> print(f"Step at x={fit.popt['x0']:.2f}") + Step at x=5.00 + + Notes + ----- + The initial parameter estimation (p0) uses heuristics based on + the data distribution. For best results with noisy or complex + data, providing manual guesses or passing p0 directly to + make_fit() is recommended. + """ + + def __init__( + self, + xobs, + yobs, + guess_x0: float | None = None, + guess_y0: float | None = None, + guess_y1: float | None = None, + **kwargs, + ): + self._guess_x0 = guess_x0 + self._guess_y0 = guess_y0 + self._guess_y1 = guess_y1 + super().__init__(xobs, yobs, **kwargs) + + @property + def function(self): + r"""The Heaviside step function. + + Returns + ------- + callable + Function with signature ``f(x, x0, y0, y1)``. + """ + + def heavy_side(x, x0, y0, y1): + r"""Evaluate Heaviside step function. + + Parameters + ---------- + x : array-like + Independent variable values. + x0 : float + Step transition x-coordinate. + y0 : float + Baseline level (value for x > x0). + y1 : float + Step height (y_left - y0 = y1, so y_left = y0 + y1). + + Returns + ------- + numpy.ndarray + Model values at x. + """ + out = y1 * np.heaviside(x0 - x, 0.5 * (y0 + y1)) + y0 + return out + + return heavy_side + + @property + def p0(self) -> list: + r"""Calculate initial parameter guess. + + The initial guess is derived from: + - User-provided guesses if available + - Otherwise, heuristic estimates from the data + + Returns + ------- + list + Initial guesses as [x0, y0, y1]. + + Raises + ------ + AssertionError + If insufficient data for estimation. + """ + assert self.sufficient_data + + x = self.observations.used.x + y = self.observations.used.y + + # Use guesses if provided, otherwise estimate from data + if self._guess_x0 is not None: + x0 = self._guess_x0 + else: + # Estimate x0 as midpoint of data range + x0 = (x.max() + x.min()) / 2 + + if self._guess_y0 is not None and self._guess_y1 is not None: + y0 = self._guess_y0 + y1 = self._guess_y1 + else: + # Estimate y0 and y1 from data above and below x0 + y0 = np.median(y[x > x0]) # Value after step + y1 = np.median(y[x < x0]) - y0 # Step height (y_left - y0) + if np.isnan(y0): + y0 = y.mean() + if np.isnan(y1): + y1 = 0.0 + + p0 = [x0, y0, y1] + return p0 + + @property + def TeX_function(self) -> str: + r"""LaTeX representation of the model. + + Returns + ------- + str + Multi-line LaTeX string describing the Heaviside function. + """ + tex = "\n".join( + [ + r"f(x) = y_1 \cdot H(x_0 - x) + y_0", + r"x < x_0: f(x) = y_0 + y_1", + r"x > x_0: f(x) = y_0", + ] + ) + return tex diff --git a/solarwindpy/fitfunctions/hinge.py b/solarwindpy/fitfunctions/hinge.py new file mode 100644 index 00000000..a1aa0819 --- /dev/null +++ b/solarwindpy/fitfunctions/hinge.py @@ -0,0 +1,1297 @@ +r"""Hinge (piecewise linear) fit functions. + +This module provides fit functions for piecewise linear models with a +hinge point, commonly used for modeling saturation behavior. +""" + +from __future__ import annotations + +from collections import namedtuple + +import numpy as np + +from .core import FitFunction + + +# Named tuple for x-intercepts used by HingeAtPoint +XIntercepts = namedtuple("XIntercepts", "x1,x2") + + +class HingeSaturation(FitFunction): + r"""Piecewise linear function with hinge point for saturation modeling. + + The model consists of two linear segments joined at a hinge point (xh, yh): + + - Rising region (x < xh): :math:`f(x) = m_1 (x - x_1)` + - Plateau region (x >= xh): :math:`f(x) = m_2 (x - x_2)` + + where the slopes and intercepts are related by continuity at the hinge: + + - :math:`m_1 = y_h / (x_h - x_1)` + - :math:`x_2 = x_h - y_h / m_2` + + Parameters + ---------- + xobs : array-like + Independent variable observations. + yobs : array-like + Dependent variable observations. + guess_xh : float, optional + Initial guess for hinge x-coordinate. Default is 326. + guess_yh : float, optional + Initial guess for hinge y-coordinate. Default is 0.5. + **kwargs + Additional arguments passed to :class:`FitFunction`. + + Attributes + ---------- + xh : float + Hinge x-coordinate (fitted parameter). + yh : float + Hinge y-coordinate (fitted parameter). + x1 : float + x-intercept of rising line (fitted parameter). + m2 : float + Slope of plateau region (fitted parameter). m2=0 gives constant saturation. + + Examples + -------- + >>> import numpy as np + >>> from solarwindpy.fitfunctions import HingeSaturation + >>> x = np.linspace(0, 15, 100) + >>> y = np.where(x < 5, 2*x, 10) # Saturation at y=10 for x>=5 + >>> fit = HingeSaturation(x, y, guess_xh=5, guess_yh=10) + >>> fit.make_fit() + >>> print(f"Hinge at ({fit.popt['xh']:.2f}, {fit.popt['yh']:.2f})") + Hinge at (5.00, 10.00) + """ + + def __init__( + self, + xobs, + yobs, + guess_xh: float = 326, + guess_yh: float = 0.5, + **kwargs, + ): + super().__init__(xobs, yobs, **kwargs) + self._saturation_guess = (guess_xh, guess_yh) + + @property + def saturation_guess(self) -> tuple[float, float]: + r"""Guess for saturation transition (xh, yh) used in p0 calculation.""" + return self._saturation_guess + + @property + def function(self): + r"""The hinge saturation function. + + Returns + ------- + callable + Function with signature ``f(x, xh, yh, x1, m2)``. + """ + + def hinge_saturation(x, xh, yh, x1, m2): + r"""Evaluate hinge saturation model. + + Parameters + ---------- + x : array-like + Independent variable values. + xh : float + Hinge x-coordinate. + yh : float + Hinge y-coordinate. + x1 : float + x-intercept of rising line. + m2 : float + Slope of plateau region. + + Returns + ------- + numpy.ndarray + Model values at x. + """ + m1 = yh / (xh - x1) + x2 = xh - (yh / m2) if abs(m2) > 1e-15 else np.inf + + y1 = m1 * (x - x1) + y2 = m2 * (x - x2) if abs(m2) > 1e-15 else yh * np.ones_like(x) + + out = np.minimum(y1, y2) + return out + + return hinge_saturation + + @property + def p0(self) -> list: + r"""Calculate initial parameter guess. + + The initial guess is derived from: + - ``saturation_guess`` for (xh, yh) + - Estimated x1 from linear fit to rising region, or data minimum + - Median slope in plateau region for m2 + + Returns + ------- + list + Initial guesses as [xh, yh, x1, m2]. + + Raises + ------ + AssertionError + If insufficient data for estimation. + """ + assert self.sufficient_data + + xh, yh = self.saturation_guess + + x = self.observations.used.x + y = self.observations.used.y + + # Estimate x1 from data in rising region + # m1 = yh / (xh - x1), so x1 = xh - yh/m1 + rising_mask = x < xh + if rising_mask.sum() >= 2: + x_rising = x[rising_mask] + y_rising = y[rising_mask] + # Simple linear regression to estimate slope m1 + m1_est = np.polyfit(x_rising, y_rising, 1)[0] + if abs(m1_est) > 1e-10: + x1 = xh - yh / m1_est + else: + x1 = x.min() + else: + # Fall back to minimum x value + x1 = x.min() + + # Estimate m2 from slope in plateau region + plateau_mask = x >= xh + if plateau_mask.sum() >= 2: + m2 = np.median(np.ediff1d(y[plateau_mask]) / np.ediff1d(x[plateau_mask])) + else: + m2 = 0.0 + + p0 = [xh, yh, x1, m2] + return p0 + + @property + def TeX_function(self) -> str: + r"""LaTeX representation of the model. + + Returns + ------- + str + Multi-line LaTeX string describing the piecewise function. + """ + tex = "\n".join( + [ + r"f(x)=\min(y_1, \, y_2)", + r"y_i = m_i(x-x_i)", + r"m_1 = \frac{y_h}{x_h - x_1}", + r"x_2 = x_h - \frac{y_h}{m_2}", + ] + ) + return tex + + +class TwoLine(FitFunction): + r"""Piecewise linear function with two intersecting lines using minimum. + + The model consists of two linear segments: + + .. math:: + + f(x) = \min(y_1, y_2) + + where: + + - :math:`y_1 = m_1 (x - x_1)` + - :math:`y_2 = m_2 (x - x_2)` + + The lines intersect at the saturation point :math:`(x_s, s)` where: + + - :math:`x_s = \frac{m_1 x_1 - m_2 x_2}{m_1 - m_2}` + - :math:`s = m_1 (x_s - x_1) = m_2 (x_s - x_2)` + + Parameters + ---------- + xobs : array-like + Independent variable observations. + yobs : array-like + Dependent variable observations. + guess_xs : float, optional + Initial guess for saturation x-coordinate. Default is 425.0. + **kwargs + Additional arguments passed to :class:`FitFunction`. + + Attributes + ---------- + x1 : float + x-intercept of first line (fitted parameter). + x2 : float + x-intercept of second line (fitted parameter). + m1 : float + Slope of first line (fitted parameter). + m2 : float + Slope of second line (fitted parameter). + xs : float + x-coordinate of intersection point (derived property). + s : float + y-coordinate of intersection point (derived property). + theta : float + Angle between the two lines in radians (derived property). + + Examples + -------- + >>> import numpy as np + >>> from solarwindpy.fitfunctions import TwoLine + >>> x = np.linspace(0, 15, 100) + >>> y = np.minimum(2*(x-0), -1*(x-15)) # Two lines intersecting at (5, 10) + >>> fit = TwoLine(x, y, guess_xs=5.0) + >>> fit.make_fit() + >>> print(f"Intersection at ({fit.xs:.2f}, {fit.s:.2f})") + Intersection at (5.00, 10.00) + """ + + def __init__( + self, + xobs, + yobs, + guess_xs: float = 425.0, + **kwargs, + ): + super().__init__(xobs, yobs, **kwargs) + self._guess_xs = guess_xs + + @property + def guess_xs(self) -> float: + r"""Initial guess for saturation x-coordinate used in p0 calculation.""" + return self._guess_xs + + @property + def function(self): + r"""The two-line minimum function. + + Returns + ------- + callable + Function with signature ``f(x, x1, x2, m1, m2)``. + """ + + def twoline(x, x1, x2, m1, m2): + r"""Evaluate two-line minimum model. + + Parameters + ---------- + x : array-like + Independent variable values. + x1 : float + x-intercept of first line. + x2 : float + x-intercept of second line. + m1 : float + Slope of first line. + m2 : float + Slope of second line. + + Returns + ------- + numpy.ndarray + Model values at x. + """ + l1 = m1 * (x - x1) + l2 = m2 * (x - x2) + out = np.minimum(l1, l2) + return out + + return twoline + + @property + def xs(self) -> float: + r"""x-coordinate of the intersection (saturation) point. + + Calculated as: + + .. math:: + + x_s = \frac{m_1 x_1 - m_2 x_2}{m_1 - m_2} + """ + popt = self.popt + x1 = popt["x1"] + x2 = popt["x2"] + m1 = popt["m1"] + m2 = popt["m2"] + n = (m1 * x1) - (m2 * x2) + d = m1 - m2 + return n / d + + @property + def s(self) -> float: + r"""y-coordinate of the intersection (saturation) point. + + Calculated as: + + .. math:: + + s = m_1 (x_s - x_1) + """ + popt = self.popt + x1 = popt["x1"] + m1 = popt["m1"] + xs = self.xs + return m1 * (xs - x1) + + @property + def theta(self) -> float: + r"""Angle between the two lines in radians. + + Calculated as: + + .. math:: + + \theta = \arctan(m_1) - \arctan(m_2) + """ + m1 = self.popt["m1"] + m2 = self.popt["m2"] + return np.arctan(m1) - np.arctan(m2) + + @property + def p0(self) -> list: + r"""Calculate initial parameter guess. + + The initial guess is derived from the data by estimating slopes + and intercepts in regions separated by ``guess_xs``. + + Returns + ------- + list + Initial guesses as [x1, x2, m1, m2]. + + Raises + ------ + AssertionError + If insufficient data for estimation. + + Notes + ----- + # TODO: Convert to data-driven p0 estimation (see GH issue #XX) + Uses hardcoded xs=425 as the default separation point, which is + appropriate for solar wind speed analysis. + """ + assert self.sufficient_data + + def estimate_line(x, y, tk, xs): + x = x[tk] + y = y[tk] + + m_set = np.ediff1d(y) / np.ediff1d(x) + m = np.nanmedian(m_set) + x0 = np.nanmedian(x - (y / m)) + return x0, m + + x = self.observations.used.x + y = self.observations.used.y + + # TODO: Convert to data-driven p0 estimation (see GH issue #XX) + xs = self._guess_xs + tk = x <= xs + + x1, m1 = estimate_line(x, y, tk, xs) + x2, m2 = estimate_line(x, y, ~tk, xs) + + p0 = [x1, x2, m1, m2] + return p0 + + @property + def TeX_function(self) -> str: + r"""LaTeX representation of the model. + + Returns + ------- + str + Multi-line LaTeX string describing the piecewise function. + """ + tex = "\n".join( + [ + r"f(x) \, =\min\left(y_1, \, y_2\right)", + r"y_i = m_i(x - x_i)", + ] + ) + return tex + + +class Saturation(FitFunction): + r"""Piecewise linear function reparameterized for saturation analysis. + + This is an alternative parameterization of :class:`TwoLine` where the + saturation point coordinates and the angle between lines are used + directly as parameters: + + .. math:: + + f(x) = \min(y_1, y_2) + + Parameters are :math:`(x_1, x_s, s, \theta)` where: + + - :math:`x_1`: x-intercept of the rising line + - :math:`x_s`: x-coordinate of saturation point + - :math:`s`: y-coordinate of saturation point + - :math:`\theta`: angle between the two lines (radians) + + The slopes are derived as: + + - :math:`m_1 = s / (x_s - x_1)` + - :math:`m_2 = \tan(\arctan(m_1) - \theta)` + + Parameters + ---------- + xobs : array-like + Independent variable observations. + yobs : array-like + Dependent variable observations. + guess_xs : float, optional + Initial guess for saturation x-coordinate. Default is 425.0. + guess_s : float, optional + Initial guess for saturation y-value. Default is 0.5. + **kwargs + Additional arguments passed to :class:`FitFunction`. + + Attributes + ---------- + x1 : float + x-intercept of rising line (fitted parameter). + xs : float + x-coordinate of saturation point (fitted parameter). + s : float + y-coordinate of saturation point (fitted parameter). + theta : float + Angle between lines (fitted parameter, radians). + m1 : float + Slope of rising line (derived property). + m2 : float + Slope of plateau line (derived property). + x2 : float + x-intercept of plateau line (derived property). + + Examples + -------- + >>> import numpy as np + >>> from solarwindpy.fitfunctions import Saturation + >>> x = np.linspace(0, 15, 100) + >>> y = np.minimum(2*(x-0), -1*(x-15)) + >>> fit = Saturation(x, y, guess_xs=5.0, guess_s=10.0) + >>> fit.make_fit() + >>> print(f"Saturation at ({fit.popt['xs']:.2f}, {fit.popt['s']:.2f})") + Saturation at (5.00, 10.00) + """ + + def __init__( + self, + xobs, + yobs, + guess_xs: float = 425.0, + guess_s: float = 0.5, + **kwargs, + ): + super().__init__(xobs, yobs, **kwargs) + self._saturation_guess = (guess_xs, guess_s) + + @property + def saturation_guess(self) -> tuple[float, float]: + r"""Guess for saturation transition (xs, s) used in p0 calculation.""" + return self._saturation_guess + + @property + def function(self): + r"""The saturation function. + + Returns + ------- + callable + Function with signature ``f(x, x1, xs, s, theta)``. + """ + + def saturation(x, x1, xs, s, theta): + r"""Evaluate saturation model. + + Parameters + ---------- + x : array-like + Independent variable values. + x1 : float + x-intercept of rising line. + xs : float + x-coordinate of saturation point. + s : float + y-coordinate of saturation point. + theta : float + Angle between lines in radians. + + Returns + ------- + numpy.ndarray + Model values at x. + """ + m1 = s / (xs - x1) + m2 = np.tan(np.arctan(m1) - theta) + x2 = xs - (s / m2) + + l1 = m1 * (x - x1) + l2 = m2 * (x - x2) + out = np.minimum(l1, l2) + return out + + return saturation + + @property + def m1(self) -> float: + r"""Slope of the rising line. + + Calculated as: + + .. math:: + + m_1 = \frac{s}{x_s - x_1} + """ + popt = self.popt + s = popt["s"] + xs = popt["xs"] + x1 = popt["x1"] + return s / (xs - x1) + + @property + def m2(self) -> float: + r"""Slope of the plateau line. + + Calculated as: + + .. math:: + + m_2 = \tan(\arctan(m_1) - \theta) + """ + popt = self.popt + theta = popt["theta"] + m1 = self.m1 + return np.tan(np.arctan(m1) - theta) + + @property + def x2(self) -> float: + r"""x-intercept of the plateau line. + + Calculated as: + + .. math:: + + x_2 = x_s - \frac{s}{m_2} + """ + popt = self.popt + s = popt["s"] + xs = popt["xs"] + m2 = self.m2 + return xs - (s / m2) + + @property + def p0(self) -> list: + r"""Calculate initial parameter guess. + + The initial guess is derived from the data by estimating slopes + and intercepts in regions separated by the saturation guess. + + Returns + ------- + list + Initial guesses as [x1, xs, s, theta]. + + Raises + ------ + AssertionError + If insufficient data for estimation. + + Notes + ----- + # TODO: Convert to data-driven p0 estimation (see GH issue #XX) + Uses hardcoded xs=425 as the default separation point, which is + appropriate for solar wind speed analysis. + """ + assert self.sufficient_data + + def estimate_line(x, y, tk, xs): + x = x[tk] + y = y[tk] + + m_set = np.ediff1d(y) / np.ediff1d(x) + m = np.nanmedian(m_set) + x0 = np.nanmedian(x - (y / m)) + s = m * (xs - x0) + return x0, m, s + + x = self.observations.used.x + y = self.observations.used.y + + # TODO: Convert to data-driven p0 estimation (see GH issue #XX) + xs, _ = self.saturation_guess + tk = x <= xs + + x1, m1, s1 = estimate_line(x, y, tk, xs) + x2, m2, s2 = estimate_line(x, y, ~tk, xs) + + s = np.nanmedian([s1, s2]) + theta = np.arctan((m1 - m2) / (1 + m1 * m2)) + + p0 = [x1, xs, s, theta] + return p0 + + @property + def TeX_function(self) -> str: + r"""LaTeX representation of the model. + + Returns + ------- + str + Multi-line LaTeX string describing the piecewise function. + """ + tex = "\n".join( + [ + r"f \, \left(x, x_1, x_s, s, \theta\right)=\min\left(y_1, \, y_2\right)", + r"y_i = m_i(x - x_i)", + r"x_s = \frac{m_2 x_2 - m_1 x_1}{m_2 - m_1}", + r"s = m_i (x_s - x_i)", + r"\theta = \arctan\left(\frac{m_1 - m_2}{1 + m_1 m_2}\right)", + ] + ) + return tex + + +class HingeMin(FitFunction): + r"""Piecewise linear function with hinge point using minimum. + + The model consists of two linear segments joined at a hinge point: + + .. math:: + + f(x) = \min(y_1, y_2) + + where: + + - :math:`y_1 = m_1 (x - x_1)` + - :math:`y_2 = m_2 (x - x_2)` + + Both lines pass through the hinge point :math:`(h, y_h)` where + :math:`y_h = m_1 (h - x_1)`. The second slope is constrained by: + + .. math:: + + m_2 = m_1 \frac{h - x_1}{h - x_2} + + Parameters + ---------- + xobs : array-like + Independent variable observations. + yobs : array-like + Dependent variable observations. + guess_h : float, optional + Initial guess for hinge x-coordinate. Default is 400.0. + **kwargs + Additional arguments passed to :class:`FitFunction`. + + Attributes + ---------- + m1 : float + Slope of first line (fitted parameter). + x1 : float + x-intercept of first line (fitted parameter). + x2 : float + x-intercept of second line (fitted parameter). + h : float + x-coordinate of hinge point (fitted parameter). + m2 : float + Slope of second line (derived property). + theta : float + Angle between the two lines in radians (derived property). + + Examples + -------- + >>> import numpy as np + >>> from solarwindpy.fitfunctions import HingeMin + >>> x = np.linspace(0, 15, 100) + >>> y = np.minimum(2*(x-0), -2*(x-10)) # Two lines meeting at (5, 10) + >>> fit = HingeMin(x, y, guess_h=5.0) + >>> fit.make_fit() + >>> print(f"Hinge at x={fit.popt['h']:.2f}") + Hinge at x=5.00 + """ + + def __init__( + self, + xobs, + yobs, + guess_h: float = 400.0, + **kwargs, + ): + super().__init__(xobs, yobs, **kwargs) + self._guess_h = guess_h + + @property + def guess_h(self) -> float: + r"""Initial guess for hinge x-coordinate used in p0 calculation.""" + return self._guess_h + + @property + def function(self): + r"""The hinge minimum function. + + Returns + ------- + callable + Function with signature ``f(x, m1, x1, x2, h)``. + """ + + def hinge(x, m1, x1, x2, h): + r"""Evaluate hinge minimum model. + + Parameters + ---------- + x : array-like + Independent variable values. + m1 : float + Slope of first line. + x1 : float + x-intercept of first line. + x2 : float + x-intercept of second line. + h : float + x-coordinate of hinge point. + + Returns + ------- + numpy.ndarray + Model values at x. + """ + m2 = m1 * (h - x1) / (h - x2) + l1 = m1 * (x - x1) + l2 = m2 * (x - x2) + out = np.minimum(l1, l2) + return out + + return hinge + + @property + def m2(self) -> float: + r"""Slope of the second line. + + Derived from the constraint that both lines pass through the hinge: + + .. math:: + + m_2 = m_1 \frac{h - x_1}{h - x_2} + """ + popt = self.popt + h = popt["h"] + m1 = popt["m1"] + x1 = popt["x1"] + x2 = popt["x2"] + return m1 * (h - x1) / (h - x2) + + @property + def theta(self) -> float: + r"""Angle between the two lines in radians. + + Calculated using arctan2 for proper quadrant handling: + + .. math:: + + \theta = \arctan2(m_1 - m_2, 1 + m_1 m_2) + """ + m1 = self.popt["m1"] + m2 = self.m2 + top = m1 - m2 + bottom = 1 + (m1 * m2) + return np.arctan2(top, bottom) + + @property + def p0(self) -> list: + r"""Calculate initial parameter guess. + + The initial guess estimates slopes and intercepts from the data + in regions separated by the hinge guess. + + Returns + ------- + list + Initial guesses as [m1, x1, x2, h]. + + Raises + ------ + AssertionError + If insufficient data for estimation. + + Notes + ----- + # TODO: Convert to data-driven p0 estimation (see GH issue #XX) + Default guess_h=400 is appropriate for solar wind speed analysis. + """ + assert self.sufficient_data + + x = self.observations.used.x + y = self.observations.used.y + h = self._guess_h + + # Estimate m1 and x1 from region below hinge + tk_below = x < h + if tk_below.sum() >= 2: + m1_set = np.ediff1d(y[tk_below]) / np.ediff1d(x[tk_below]) + m1 = np.nanmedian(m1_set) + x1 = np.nanmedian(x[tk_below] - (y[tk_below] / m1)) + else: + # Fall back to simple estimate + m1 = (y.max() - y.min()) / (x.max() - x.min()) + x1 = x.min() + + # Estimate m2 and x2 from plateau region + tk_above = x >= h + if tk_above.sum() >= 2: + m2 = np.median(np.ediff1d(y[tk_above]) / np.ediff1d(x[tk_above])) + x2 = np.median(x[tk_above] - (y[tk_above] / m2)) + else: + m2 = 0.0 + x2 = h + + p0 = [m1, x1, x2, h] + return p0 + + @property + def TeX_function(self) -> str: + r"""LaTeX representation of the model. + + Returns + ------- + str + Multi-line LaTeX string describing the piecewise function. + """ + tex = "\n".join( + [ + r"f(x)=\min(m_1(x-x_1), \, m_2(x-x_2))", + r"m_2 = m_1 \frac{h - x_1}{h - x_2}", + ] + ) + return tex + + +class HingeMax(FitFunction): + r"""Piecewise linear function with hinge point using maximum. + + The model consists of two linear segments joined at a hinge point: + + .. math:: + + f(x) = \max(y_1, y_2) + + where: + + - :math:`y_1 = m_1 (x - x_1)` + - :math:`y_2 = m_2 (x - x_2)` + + Both lines pass through the hinge point :math:`(h, y_h)` where + :math:`y_h = m_1 (h - x_1)`. The second slope is constrained by: + + .. math:: + + m_2 = m_1 \frac{h - x_1}{h - x_2} + + This is the same as :class:`HingeMin` but uses ``np.maximum`` instead + of ``np.minimum``, suitable for V-shaped patterns opening upward. + + Parameters + ---------- + xobs : array-like + Independent variable observations. + yobs : array-like + Dependent variable observations. + guess_h : float, optional + Initial guess for hinge x-coordinate. Default is 400.0. + **kwargs + Additional arguments passed to :class:`FitFunction`. + + Attributes + ---------- + m1 : float + Slope of first line (fitted parameter). + x1 : float + x-intercept of first line (fitted parameter). + x2 : float + x-intercept of second line (fitted parameter). + h : float + x-coordinate of hinge point (fitted parameter). + m2 : float + Slope of second line (derived property). + theta : float + Angle between the two lines in radians (derived property). + + Examples + -------- + >>> import numpy as np + >>> from solarwindpy.fitfunctions import HingeMax + >>> x = np.linspace(0, 15, 100) + >>> y = np.maximum(-2*(x-0), 2*(x-10)) # V-shape with vertex at (5, -10) + >>> fit = HingeMax(x, y, guess_h=5.0) + >>> fit.make_fit() + >>> print(f"Hinge at x={fit.popt['h']:.2f}") + Hinge at x=5.00 + """ + + def __init__( + self, + xobs, + yobs, + guess_h: float = 400.0, + **kwargs, + ): + super().__init__(xobs, yobs, **kwargs) + self._guess_h = guess_h + + @property + def guess_h(self) -> float: + r"""Initial guess for hinge x-coordinate used in p0 calculation.""" + return self._guess_h + + @property + def function(self): + r"""The hinge maximum function. + + Returns + ------- + callable + Function with signature ``f(x, m1, x1, x2, h)``. + """ + + def hinge(x, m1, x1, x2, h): + r"""Evaluate hinge maximum model. + + Parameters + ---------- + x : array-like + Independent variable values. + m1 : float + Slope of first line. + x1 : float + x-intercept of first line. + x2 : float + x-intercept of second line. + h : float + x-coordinate of hinge point. + + Returns + ------- + numpy.ndarray + Model values at x. + """ + m2 = m1 * (h - x1) / (h - x2) + l1 = m1 * (x - x1) + l2 = m2 * (x - x2) + out = np.maximum(l1, l2) + return out + + return hinge + + @property + def m2(self) -> float: + r"""Slope of the second line. + + Derived from the constraint that both lines pass through the hinge: + + .. math:: + + m_2 = m_1 \frac{h - x_1}{h - x_2} + """ + popt = self.popt + h = popt["h"] + m1 = popt["m1"] + x1 = popt["x1"] + x2 = popt["x2"] + return m1 * (h - x1) / (h - x2) + + @property + def theta(self) -> float: + r"""Angle between the two lines in radians. + + Calculated using arctan2 for proper quadrant handling: + + .. math:: + + \theta = \arctan2(m_1 - m_2, 1 + m_1 m_2) + """ + m1 = self.popt["m1"] + m2 = self.m2 + top = m1 - m2 + bottom = 1 + (m1 * m2) + return np.arctan2(top, bottom) + + @property + def p0(self) -> list: + r"""Calculate initial parameter guess. + + The initial guess estimates slopes and intercepts from the data + in regions separated by the hinge guess. + + Returns + ------- + list + Initial guesses as [m1, x1, x2, h]. + + Raises + ------ + AssertionError + If insufficient data for estimation. + + Notes + ----- + # TODO: Convert to data-driven p0 estimation (see GH issue #XX) + Default guess_h=400 is appropriate for solar wind speed analysis. + """ + assert self.sufficient_data + + x = self.observations.used.x + y = self.observations.used.y + h = self._guess_h + + # Estimate m1 and x1 from region below hinge + tk_below = x < h + if tk_below.sum() >= 2: + m1_set = np.ediff1d(y[tk_below]) / np.ediff1d(x[tk_below]) + m1 = np.nanmedian(m1_set) + x1 = np.nanmedian(x[tk_below] - (y[tk_below] / m1)) + else: + # Fall back to simple estimate + m1 = (y.max() - y.min()) / (x.max() - x.min()) + x1 = x.min() + + # Estimate m2 and x2 from plateau region + tk_above = x >= h + if tk_above.sum() >= 2: + m2 = np.median(np.ediff1d(y[tk_above]) / np.ediff1d(x[tk_above])) + x2 = np.median(x[tk_above] - (y[tk_above] / m2)) + else: + m2 = 0.0 + x2 = h + + p0 = [m1, x1, x2, h] + return p0 + + @property + def TeX_function(self) -> str: + r"""LaTeX representation of the model. + + Returns + ------- + str + Multi-line LaTeX string describing the piecewise function. + """ + tex = "\n".join( + [ + r"f(x)=\max(m_1(x-x_1), \, m_2(x-x_2))", + r"m_2 = m_1 \frac{h - x_1}{h - x_2}", + ] + ) + return tex + + +class HingeAtPoint(FitFunction): + r"""Piecewise linear function passing through a specified hinge point. + + The model consists of two linear segments that both pass through the + hinge point :math:`(x_h, y_h)`: + + .. math:: + + f(x) = \min(y_1, y_2) + + where: + + - :math:`y_1 = m_1 (x - x_1)` with :math:`x_1 = x_h - y_h / m_1` + - :math:`y_2 = m_2 (x - x_2)` with :math:`x_2 = x_h - y_h / m_2` + + Parameters + ---------- + xobs : array-like + Independent variable observations. + yobs : array-like + Dependent variable observations. + guess_xh : float, optional + Initial guess for hinge x-coordinate. Default is 400.0. + guess_yh : float, optional + Initial guess for hinge y-coordinate. Default is 0.5. + **kwargs + Additional arguments passed to :class:`FitFunction`. + + Attributes + ---------- + xh : float + x-coordinate of hinge point (fitted parameter). + yh : float + y-coordinate of hinge point (fitted parameter). + m1 : float + Slope of first line (fitted parameter). + m2 : float + Slope of second line (fitted parameter). + x_intercepts : XIntercepts + Named tuple with x1 and x2 attributes (derived property). + + Examples + -------- + >>> import numpy as np + >>> from solarwindpy.fitfunctions import HingeAtPoint + >>> x = np.linspace(0, 15, 100) + >>> y = np.minimum(2*(x-0), -1*(x-15)) # Hinge at (5, 10) + >>> fit = HingeAtPoint(x, y, guess_xh=5.0, guess_yh=10.0) + >>> fit.make_fit() + >>> print(f"Hinge at ({fit.popt['xh']:.2f}, {fit.popt['yh']:.2f})") + Hinge at (5.00, 10.00) + """ + + def __init__( + self, + xobs, + yobs, + guess_xh: float = 400.0, + guess_yh: float = 0.5, + **kwargs, + ): + super().__init__(xobs, yobs, **kwargs) + self._hinge_guess = (guess_xh, guess_yh) + + @property + def hinge_guess(self) -> tuple[float, float]: + r"""Guess for hinge point (xh, yh) used in p0 calculation.""" + return self._hinge_guess + + @property + def function(self): + r"""The hinge-at-point function. + + Returns + ------- + callable + Function with signature ``f(x, xh, yh, m1, m2)``. + """ + + def hinge_at_point(x, xh, yh, m1, m2): + r"""Evaluate hinge-at-point model. + + Parameters + ---------- + x : array-like + Independent variable values. + xh : float + x-coordinate of hinge point. + yh : float + y-coordinate of hinge point. + m1 : float + Slope of first line. + m2 : float + Slope of second line. + + Returns + ------- + numpy.ndarray + Model values at x. + """ + x1 = xh - (yh / m1) + x2 = xh - (yh / m2) + + y1 = m1 * (x - x1) + y2 = m2 * (x - x2) + + out = np.minimum(y1, y2) + return out + + return hinge_at_point + + @property + def x_intercepts(self) -> XIntercepts: + r"""x-intercepts of the two lines. + + Returns a named tuple with: + + - x1 = xh - yh / m1 + - x2 = xh - yh / m2 + """ + popt = self.popt + xh = popt["xh"] + yh = popt["yh"] + m1 = popt["m1"] + m2 = popt["m2"] + x1 = xh - (yh / m1) + x2 = xh - (yh / m2) + return XIntercepts(x1, x2) + + @property + def p0(self) -> list: + r"""Calculate initial parameter guess. + + The initial guess uses the hinge_guess for (xh, yh) and estimates + slopes from the data in regions separated by the hinge guess. + + Returns + ------- + list + Initial guesses as [xh, yh, m1, m2]. + + Raises + ------ + AssertionError + If insufficient data for estimation. + + Notes + ----- + # TODO: Convert to data-driven p0 estimation (see GH issue #XX) + Default guess_xh=400 and guess_yh=0.5 are appropriate for solar + wind speed analysis. + """ + assert self.sufficient_data + + xh, yh_guess = self._hinge_guess + + x = self.observations.used.x + y = self.observations.used.y + + # Estimate yh from data near the hinge point + yh = y[np.argmin(np.abs(x - xh))] + + # Estimate m1 from region below hinge + tk_below = x < xh + if tk_below.sum() >= 2: + m1_set = np.ediff1d(y[tk_below]) / np.ediff1d(x[tk_below]) + m1 = np.nanmedian(m1_set) + else: + # Fall back to simple estimate + m1 = yh / (xh - x.min()) if xh > x.min() else 1.0 + + # Estimate m2 from region above hinge + tk_above = x >= xh + if tk_above.sum() >= 2: + m2 = np.median(np.ediff1d(y[tk_above]) / np.ediff1d(x[tk_above])) + else: + m2 = 0.0 + + p0 = [xh, yh, m1, m2] + return p0 + + @property + def TeX_function(self) -> str: + r"""LaTeX representation of the model. + + Returns + ------- + str + Multi-line LaTeX string describing the piecewise function. + """ + tex = "\n".join( + [ + r"f(x)=\min(y_1, \, y_2)", + r"y_i = m_i(x-x_i)", + r"x_i = x_h - \frac{y_h}{m_i}", + ] + ) + return tex diff --git a/tests/fitfunctions/test_composite.py b/tests/fitfunctions/test_composite.py new file mode 100644 index 00000000..240fb14c --- /dev/null +++ b/tests/fitfunctions/test_composite.py @@ -0,0 +1,1350 @@ +"""Tests for Gaussian-Heaviside composite fit functions. + +This module tests three composite functions that combine Gaussian and Heaviside +step functions: + +1. GaussianPlusHeavySide: + f(x) = Gaussian(x, mu, sigma, A) + y1 * H(x0-x) + y0 + - Gaussian peak plus a step function that adds y1 for x < x0 + +2. GaussianTimesHeavySide: + f(x) = Gaussian(x, mu, sigma, A) * H(x-x0) + - Gaussian truncated at x0; zero for x < x0 + +3. GaussianTimesHeavySidePlusHeavySide: + f(x) = Gaussian(x, mu, sigma, A) * H(x-x0) + y1 * H(x0-x) + - Gaussian for x >= x0, constant y1 for x < x0 + +Where: +- Gaussian(x, mu, sigma, A) = A * exp(-0.5 * ((x-mu)/sigma)^2) +- H(z) is the Heaviside step function (H(z) = 0 for z < 0, 0.5 at z=0, 1 for z > 0) + +Mathematical derivations for expected values: + +For a standard Gaussian at x=mu: + Gaussian(mu, mu, sigma, A) = A * exp(0) = A + +For Heaviside transitions: + H(x0 - x) = 1 when x < x0, 0.5 when x = x0, 0 when x > x0 + H(x - x0) = 0 when x < x0, 0.5 when x = x0, 1 when x > x0 +""" + +import inspect + +import numpy as np +import pytest + +from solarwindpy.fitfunctions.composite import ( + GaussianPlusHeavySide, + GaussianTimesHeavySide, + GaussianTimesHeavySidePlusHeavySide, +) +from solarwindpy.fitfunctions.core import InsufficientDataError + + +# ============================================================================= +# Helper Functions +# ============================================================================= + + +def gaussian(x, mu, sigma, A): + """Standard Gaussian function for test reference calculations.""" + return A * np.exp(-0.5 * ((x - mu) / sigma) ** 2) + + +# ============================================================================= +# GaussianPlusHeavySide Fixtures +# ============================================================================= + + +@pytest.fixture +def gph_clean_data(): + """Clean data for GaussianPlusHeavySide. + + Parameters: x0=2.0, y0=1.0, y1=3.0, mu=5.0, sigma=1.0, A=4.0 + + Function behavior: + - For x < x0: f(x) = Gaussian(x) + y1 + y0 = Gaussian(x) + 4.0 + - For x = x0: f(x) = Gaussian(x) + 0.5*y1 + y0 = Gaussian(x) + 2.5 + - For x > x0: f(x) = Gaussian(x) + y0 = Gaussian(x) + 1.0 + """ + true_params = {"x0": 2.0, "y0": 1.0, "y1": 3.0, "mu": 5.0, "sigma": 1.0, "A": 4.0} + x = np.linspace(0, 10, 200) + + # Build y: Gaussian + y1*H(x0-x) + y0 + gauss = gaussian(x, true_params["mu"], true_params["sigma"], true_params["A"]) + heaviside_term = true_params["y1"] * np.heaviside(true_params["x0"] - x, 0.5) + y = gauss + heaviside_term + true_params["y0"] + + w = np.ones_like(x) + return x, y, w, true_params + + +@pytest.fixture +def gph_noisy_data(): + """Noisy data for GaussianPlusHeavySide with 3% noise. + + Parameters: x0=2.0, y0=1.0, y1=3.0, mu=5.0, sigma=1.0, A=4.0 + Noise std = 0.15 (approximately 3% of peak amplitude A+y0+y1) + """ + rng = np.random.default_rng(42) + true_params = {"x0": 2.0, "y0": 1.0, "y1": 3.0, "mu": 5.0, "sigma": 1.0, "A": 4.0} + x = np.linspace(0, 10, 200) + + gauss = gaussian(x, true_params["mu"], true_params["sigma"], true_params["A"]) + heaviside_term = true_params["y1"] * np.heaviside(true_params["x0"] - x, 0.5) + y_true = gauss + heaviside_term + true_params["y0"] + + noise_std = 0.15 + y = y_true + rng.normal(0, noise_std, len(x)) + w = np.ones_like(x) / noise_std + return x, y, w, true_params + + +# ============================================================================= +# GaussianTimesHeavySide Fixtures +# ============================================================================= + + +@pytest.fixture +def gth_clean_data(): + """Clean data for GaussianTimesHeavySide. + + Parameters: x0=3.0, mu=5.0, sigma=1.0, A=4.0 + + Function behavior: + - For x < x0: f(x) = 0 (Heaviside is 0) + - For x = x0: f(x) = Gaussian(x0) * 1.0 = Gaussian(x0) + - For x > x0: f(x) = Gaussian(x) (Heaviside is 1) + """ + true_params = {"x0": 3.0, "mu": 5.0, "sigma": 1.0, "A": 4.0} + x = np.linspace(0, 10, 200) + + # Build y: Gaussian * H(x-x0) + gauss = gaussian(x, true_params["mu"], true_params["sigma"], true_params["A"]) + y = gauss * np.heaviside(x - true_params["x0"], 1.0) + + w = np.ones_like(x) + return x, y, w, true_params + + +@pytest.fixture +def gth_noisy_data(): + """Noisy data for GaussianTimesHeavySide with 3% noise. + + Parameters: x0=3.0, mu=5.0, sigma=1.0, A=4.0 + Noise std = 0.12 (approximately 3% of peak amplitude) + """ + rng = np.random.default_rng(42) + true_params = {"x0": 3.0, "mu": 5.0, "sigma": 1.0, "A": 4.0} + x = np.linspace(0, 10, 200) + + gauss = gaussian(x, true_params["mu"], true_params["sigma"], true_params["A"]) + y_true = gauss * np.heaviside(x - true_params["x0"], 1.0) + + noise_std = 0.12 + y = y_true + rng.normal(0, noise_std, len(x)) + # Ensure y >= 0 for x < x0 (the function should be ~0 there) + y[x < true_params["x0"]] = np.abs(y[x < true_params["x0"]]) + w = np.ones_like(x) / noise_std + return x, y, w, true_params + + +# ============================================================================= +# GaussianTimesHeavySidePlusHeavySide Fixtures +# ============================================================================= + + +@pytest.fixture +def gthph_clean_data(): + """Clean data for GaussianTimesHeavySidePlusHeavySide. + + Parameters: x0=3.0, y1=2.0, mu=5.0, sigma=1.0, A=4.0 + + Function behavior: + - For x < x0: f(x) = y1 (constant plateau) + - For x = x0: f(x) = Gaussian(x0) + y1 (both H(0) = 1.0) + - For x > x0: f(x) = Gaussian(x) (pure Gaussian) + """ + true_params = {"x0": 3.0, "y1": 2.0, "mu": 5.0, "sigma": 1.0, "A": 4.0} + x = np.linspace(0, 10, 200) + + # Build y: Gaussian * H(x-x0) + y1 * H(x0-x) with H(0) = 1.0 + gauss = gaussian(x, true_params["mu"], true_params["sigma"], true_params["A"]) + y = gauss * np.heaviside(x - true_params["x0"], 1.0) + true_params[ + "y1" + ] * np.heaviside(true_params["x0"] - x, 1.0) + + w = np.ones_like(x) + return x, y, w, true_params + + +@pytest.fixture +def gthph_noisy_data(): + """Noisy data for GaussianTimesHeavySidePlusHeavySide with 3% noise. + + Parameters: x0=3.0, y1=2.0, mu=5.0, sigma=1.0, A=4.0 + Noise std = 0.12 (approximately 3% of peak amplitude) + """ + rng = np.random.default_rng(42) + true_params = {"x0": 3.0, "y1": 2.0, "mu": 5.0, "sigma": 1.0, "A": 4.0} + x = np.linspace(0, 10, 200) + + gauss = gaussian(x, true_params["mu"], true_params["sigma"], true_params["A"]) + y_true = gauss * np.heaviside(x - true_params["x0"], 1.0) + true_params[ + "y1" + ] * np.heaviside(true_params["x0"] - x, 1.0) + + noise_std = 0.12 + y = y_true + rng.normal(0, noise_std, len(x)) + w = np.ones_like(x) / noise_std + return x, y, w, true_params + + +# ============================================================================= +# GaussianPlusHeavySide Tests +# ============================================================================= + + +class TestGaussianPlusHeavySide: + """Tests for GaussianPlusHeavySide fit function. + + GaussianPlusHeavySide models: + f(x) = Gaussian(x, mu, sigma, A) + y1 * H(x0-x) + y0 + + This produces a Gaussian peak with: + - A constant offset y0 everywhere + - An additional step of height y1 for x < x0 + """ + + # ------------------------------------------------------------------------- + # E1. Function Evaluation Tests (Exact Values) + # ------------------------------------------------------------------------- + + def test_func_evaluates_below_x0_correctly(self): + """For x < x0: f(x) = Gaussian(x) + y1 + y0. + + At x=0 with params x0=2, y0=1, y1=3, mu=5, sigma=1, A=4: + - Gaussian(0) = 4 * exp(-0.5 * ((0-5)/1)^2) = 4 * exp(-12.5) ~ 1.48e-5 + - H(2-0) = H(2) = 1 + - f(0) = Gaussian(0) + 3*1 + 1 ~ 4.0 (Gaussian contribution negligible) + """ + x0, y0, y1, mu, sigma, A = 2.0, 1.0, 3.0, 5.0, 1.0, 4.0 + + x_test = np.array([0.0, 1.0]) # Both below x0=2 + gauss_vals = gaussian(x_test, mu, sigma, A) + expected = gauss_vals + y1 * 1.0 + y0 # H(x0-x)=1 for x x0: f(x) = Gaussian(x) + y0 (Heaviside term is 0). + + At x=5 (Gaussian peak) with params x0=2, y0=1, y1=3, mu=5, sigma=1, A=4: + - Gaussian(5) = 4 * exp(0) = 4.0 + - H(2-5) = H(-3) = 0 + - f(5) = 4.0 + 0 + 1.0 = 5.0 + """ + x0, y0, y1, mu, sigma, A = 2.0, 1.0, 3.0, 5.0, 1.0, 4.0 + + x_test = np.array([3.0, 5.0, 7.0]) # All above x0=2 + gauss_vals = gaussian(x_test, mu, sigma, A) + expected = gauss_vals + y0 # H(x0-x)=0 for x>x0 + + x_dummy = np.array([0.0, 10.0]) + y_dummy = np.array([4.0, 1.0]) + obj = GaussianPlusHeavySide(x_dummy, y_dummy) + result = obj.function(x_test, x0, y0, y1, mu, sigma, A) + + np.testing.assert_allclose( + result, + expected, + rtol=1e-10, + err_msg="Above x0: f(x) should equal Gaussian(x) + y0", + ) + + def test_func_evaluates_at_x0_correctly(self): + """At x = x0: f(x) = Gaussian(x0) + 0.5*y1 + y0. + + At x=2 with params x0=2, y0=1, y1=3, mu=5, sigma=1, A=4: + - Gaussian(2) = 4 * exp(-0.5 * ((2-5)/1)^2) = 4 * exp(-4.5) ~ 0.0446 + - H(0) = 0.5 + - f(2) = Gaussian(2) + 3*0.5 + 1 = Gaussian(2) + 2.5 + """ + x0, y0, y1, mu, sigma, A = 2.0, 1.0, 3.0, 5.0, 1.0, 4.0 + + x_test = np.array([x0]) + gauss_val = gaussian(x_test, mu, sigma, A) + expected = gauss_val + 0.5 * y1 + y0 + + x_dummy = np.array([0.0, 10.0]) + y_dummy = np.array([4.0, 1.0]) + obj = GaussianPlusHeavySide(x_dummy, y_dummy) + result = obj.function(x_test, x0, y0, y1, mu, sigma, A) + + np.testing.assert_allclose( + result, + expected, + rtol=1e-10, + err_msg="At x0: f(x) should equal Gaussian(x0) + 0.5*y1 + y0", + ) + + def test_func_evaluates_at_gaussian_peak(self): + """At x = mu (Gaussian peak): f(mu) = A + y0 (assuming mu > x0). + + At x=5 with params x0=2, y0=1, y1=3, mu=5, sigma=1, A=4: + - Gaussian(5) = A = 4.0 + - H(2-5) = 0 + - f(5) = 4.0 + 0 + 1.0 = 5.0 + """ + x0, y0, y1, mu, sigma, A = 2.0, 1.0, 3.0, 5.0, 1.0, 4.0 + + x_test = np.array([mu]) + expected = np.array([A + y0]) # mu > x0, so H term is 0 + + x_dummy = np.array([0.0, 10.0]) + y_dummy = np.array([4.0, 1.0]) + obj = GaussianPlusHeavySide(x_dummy, y_dummy) + result = obj.function(x_test, x0, y0, y1, mu, sigma, A) + + np.testing.assert_allclose( + result, + expected, + rtol=1e-10, + err_msg="At Gaussian peak (mu > x0): f(mu) = A + y0", + ) + + # ------------------------------------------------------------------------- + # E2. Parameter Recovery Tests (Clean Data) + # ------------------------------------------------------------------------- + + def test_fit_recovers_gaussian_parameters_from_clean_data(self, gph_clean_data): + """Fitting noise-free data should recover Gaussian parameters within 2%. + + Note: The step function parameters (x0, y0, y1) are inherently difficult + to constrain due to the zero gradient of Heaviside functions. Only the + Gaussian parameters (mu, sigma, A) are expected to be recovered precisely. + """ + x, y, w, true_params = gph_clean_data + + obj = GaussianPlusHeavySide(x, y) + obj.make_fit() + + # Only test Gaussian parameters which can be well-constrained + for param in ["mu", "sigma", "A"]: + true_val = true_params[param] + fitted_val = obj.popt[param] + if abs(true_val) < 0.1: + assert abs(fitted_val - true_val) < 0.05, ( + f"{param}: fitted={fitted_val:.4f}, true={true_val:.4f}, " + f"absolute error exceeds 0.05" + ) + else: + rel_error = abs(fitted_val - true_val) / abs(true_val) + assert rel_error < 0.02, ( + f"{param}: fitted={fitted_val:.4f}, true={true_val:.4f}, " + f"rel_error={rel_error:.2%} exceeds 2% tolerance" + ) + + @pytest.mark.parametrize( + "true_params", + [ + {"x0": 2.0, "y0": 1.0, "y1": 3.0, "mu": 5.0, "sigma": 1.0, "A": 4.0}, + {"x0": 4.0, "y0": 0.5, "y1": 2.0, "mu": 6.0, "sigma": 0.8, "A": 3.0}, + {"x0": 1.0, "y0": 2.0, "y1": 1.0, "mu": 4.0, "sigma": 1.5, "A": 5.0}, + ], + ) + def test_fit_recovers_gaussian_parameters_for_various_combinations( + self, true_params + ): + """Gaussian parameters should be recovered for diverse parameter combinations. + + Note: Only Gaussian parameters (mu, sigma, A) are tested since step function + parameters are inherently difficult to constrain. Tolerance is 10% due to + parameter coupling effects in this complex model. + """ + x = np.linspace(0, 10, 200) + gauss = gaussian(x, true_params["mu"], true_params["sigma"], true_params["A"]) + heaviside_term = true_params["y1"] * np.heaviside(true_params["x0"] - x, 0.5) + y = gauss + heaviside_term + true_params["y0"] + + obj = GaussianPlusHeavySide(x, y) + obj.make_fit() + + # Only test Gaussian parameters which can be well-constrained + for param in ["mu", "sigma", "A"]: + true_val = true_params[param] + fitted_val = obj.popt[param] + if abs(true_val) < 0.1: + assert ( + abs(fitted_val - true_val) < 0.1 + ), f"{param}: fitted={fitted_val:.4f}, true={true_val:.4f}" + else: + rel_error = abs(fitted_val - true_val) / abs(true_val) + assert rel_error < 0.10, ( + f"{param}: fitted={fitted_val:.4f}, true={true_val:.4f}, " + f"rel_error={rel_error:.2%}" + ) + + # ------------------------------------------------------------------------- + # E3. Noisy Data Tests (Precision Bounds) + # ------------------------------------------------------------------------- + + def test_fit_with_noise_recovers_gaussian_parameters_within_2sigma( + self, gph_noisy_data + ): + """Gaussian parameters should be within 2sigma of true values (95% confidence). + + Note: Step function parameters (x0, y0, y1) are excluded because Heaviside + functions have zero gradient almost everywhere, making their uncertainties + unreliable for this type of test. + """ + x, y, w, true_params = gph_noisy_data + + obj = GaussianPlusHeavySide(x, y) + obj.make_fit() + + # Only test Gaussian parameters which have well-defined uncertainties + for param in ["mu", "sigma", "A"]: + true_val = true_params[param] + fitted_val = obj.popt[param] + sigma = obj.psigma[param] + deviation = abs(fitted_val - true_val) + + assert deviation < 2 * sigma, ( + f"{param}: |fitted({fitted_val:.4f}) - true({true_val:.4f})| = " + f"{deviation:.4f} exceeds 2sigma = {2*sigma:.4f}" + ) + + def test_fit_uncertainty_scales_with_noise(self): + """Higher noise should produce larger parameter uncertainties.""" + rng = np.random.default_rng(42) + true_params = { + "x0": 2.0, + "y0": 1.0, + "y1": 3.0, + "mu": 5.0, + "sigma": 1.0, + "A": 4.0, + } + x = np.linspace(0, 10, 200) + gauss = gaussian(x, true_params["mu"], true_params["sigma"], true_params["A"]) + heaviside_term = true_params["y1"] * np.heaviside(true_params["x0"] - x, 0.5) + y_true = gauss + heaviside_term + true_params["y0"] + + # Low noise fit + y_low = y_true + rng.normal(0, 0.05, len(x)) + fit_low = GaussianPlusHeavySide(x, y_low) + fit_low.make_fit() + + # High noise fit + rng2 = np.random.default_rng(43) + y_high = y_true + rng2.normal(0, 0.25, len(x)) + fit_high = GaussianPlusHeavySide(x, y_high) + fit_high.make_fit() + + # High noise should have larger uncertainties for most parameters + for param in ["mu", "sigma", "A"]: + assert fit_high.psigma[param] > fit_low.psigma[param], ( + f"{param}: high_noise_sigma={fit_high.psigma[param]:.4f} should be " + f"> low_noise_sigma={fit_low.psigma[param]:.4f}" + ) + + # ------------------------------------------------------------------------- + # E4. Initial Parameter Estimation Tests + # ------------------------------------------------------------------------- + + def test_p0_returns_list_with_correct_length(self, gph_clean_data): + """p0 should return a list with 6 elements.""" + x, y, w, true_params = gph_clean_data + obj = GaussianPlusHeavySide(x, y) + + p0 = obj.p0 + assert isinstance(p0, list), f"p0 should be a list, got {type(p0)}" + assert len(p0) == 6, f"p0 should have 6 elements, got {len(p0)}" + + def test_p0_enables_successful_convergence(self, gph_noisy_data): + """Fit should converge when initialized with estimated p0.""" + x, y, w, true_params = gph_noisy_data + + obj = GaussianPlusHeavySide(x, y) + obj.make_fit() + + assert all( + np.isfinite(v) for v in obj.popt.values() + ), f"Fit did not converge: popt={obj.popt}" + + # ------------------------------------------------------------------------- + # E5. Edge Case and Error Handling Tests + # ------------------------------------------------------------------------- + + def test_insufficient_data_raises_error(self): + """Fitting with insufficient data should raise InsufficientDataError.""" + x = np.array([1.0, 2.0, 3.0]) # Only 3 points for 6 parameters + y = np.array([2.0, 4.0, 3.0]) + + obj = GaussianPlusHeavySide(x, y) + + with pytest.raises(InsufficientDataError): + obj.make_fit() + + def test_function_signature(self): + """Test that function has correct parameter signature.""" + x = np.array([0.0, 5.0, 10.0]) + y = np.array([4.0, 5.0, 1.0]) + obj = GaussianPlusHeavySide(x, y) + + sig = inspect.signature(obj.function) + params = tuple(sig.parameters.keys()) + + assert params == ( + "x", + "x0", + "y0", + "y1", + "mu", + "sigma", + "A", + ), f"Function should have signature (x, x0, y0, y1, mu, sigma, A), got {params}" + + def test_callable_interface(self, gph_clean_data): + """Test that fitted object is callable and returns correct shape.""" + x, y, w, true_params = gph_clean_data + + obj = GaussianPlusHeavySide(x, y) + obj.make_fit() + + x_test = np.array([1.0, 5.0, 10.0]) + y_pred = obj(x_test) + + assert ( + y_pred.shape == x_test.shape + ), f"Predicted shape {y_pred.shape} should match input shape {x_test.shape}" + assert np.all(np.isfinite(y_pred)), "All predicted values should be finite" + + def test_popt_has_correct_keys(self, gph_clean_data): + """Test that popt contains expected parameter names.""" + x, y, w, true_params = gph_clean_data + + obj = GaussianPlusHeavySide(x, y) + obj.make_fit() + + expected_keys = {"x0", "y0", "y1", "mu", "sigma", "A"} + actual_keys = set(obj.popt.keys()) + + assert ( + actual_keys == expected_keys + ), f"popt keys should be {expected_keys}, got {actual_keys}" + + def test_psigma_has_same_keys_as_popt(self, gph_noisy_data): + """Test that psigma has same keys as popt.""" + x, y, w, true_params = gph_noisy_data + + obj = GaussianPlusHeavySide(x, y) + obj.make_fit() + + assert set(obj.psigma.keys()) == set(obj.popt.keys()), ( + f"psigma keys {set(obj.psigma.keys())} should match " + f"popt keys {set(obj.popt.keys())}" + ) + + def test_psigma_values_are_nonnegative(self, gph_noisy_data): + """Test that parameter uncertainties are non-negative. + + Note: x0 uncertainty can be zero or very small due to the step function + nature (zero gradient almost everywhere), so we only check for non-negative + values. Gaussian parameters (mu, sigma, A) should have positive uncertainties. + """ + x, y, w, true_params = gph_noisy_data + + obj = GaussianPlusHeavySide(x, y) + obj.make_fit() + + for param, sigma in obj.psigma.items(): + assert sigma >= 0, f"psigma['{param}'] = {sigma} should be non-negative" + + # Gaussian parameters should have positive uncertainties + for param in ["mu", "sigma", "A"]: + assert obj.psigma[param] > 0, f"psigma['{param}'] should be positive" + + +# ============================================================================= +# GaussianTimesHeavySide Tests +# ============================================================================= + + +class TestGaussianTimesHeavySide: + """Tests for GaussianTimesHeavySide fit function. + + GaussianTimesHeavySide models: + f(x) = Gaussian(x, mu, sigma, A) * H(x-x0) + + This produces a truncated Gaussian that is: + - Zero for x < x0 + - Gaussian(x) for x > x0 + - Gaussian(x0) at x = x0 (since H(0) = 1.0 in this implementation) + """ + + # ------------------------------------------------------------------------- + # E1. Function Evaluation Tests (Exact Values) + # ------------------------------------------------------------------------- + + def test_func_evaluates_below_x0_as_zero(self): + """For x < x0: f(x) = 0 (Heaviside is 0). + + With x0=3, the function should be exactly 0 for all x < 3. + """ + x0, mu, sigma, A = 3.0, 5.0, 1.0, 4.0 + + x_test = np.array([0.0, 1.0, 2.0, 2.99]) # All below x0=3 + expected = np.zeros_like(x_test) + + x_dummy = np.array([0.0, 10.0]) + y_dummy = np.array([0.0, 4.0]) + obj = GaussianTimesHeavySide(x_dummy, y_dummy, guess_x0=x0) + result = obj.function(x_test, x0, mu, sigma, A) + + np.testing.assert_allclose( + result, + expected, + rtol=1e-10, + err_msg="Below x0: f(x) should be exactly 0", + ) + + def test_func_evaluates_above_x0_as_gaussian(self): + """For x > x0: f(x) = Gaussian(x). + + With x0=3, mu=5, sigma=1, A=4: + - f(5) = 4 * exp(0) = 4.0 (Gaussian peak) + - f(4) = 4 * exp(-0.5) ~ 2.426 + - f(6) = 4 * exp(-0.5) ~ 2.426 + """ + x0, mu, sigma, A = 3.0, 5.0, 1.0, 4.0 + + x_test = np.array([4.0, 5.0, 6.0, 8.0]) # All above x0=3 + expected = gaussian(x_test, mu, sigma, A) + + x_dummy = np.array([0.0, 10.0]) + y_dummy = np.array([0.0, 4.0]) + obj = GaussianTimesHeavySide(x_dummy, y_dummy, guess_x0=x0) + result = obj.function(x_test, x0, mu, sigma, A) + + np.testing.assert_allclose( + result, + expected, + rtol=1e-10, + err_msg="Above x0: f(x) should equal Gaussian(x)", + ) + + def test_func_evaluates_at_x0_correctly(self): + """At x = x0: f(x) = Gaussian(x0) * 1.0 = Gaussian(x0). + + This implementation uses H(0) = 1.0, so at the transition point + the function equals the full Gaussian value. + + With x0=3, mu=5, sigma=1, A=4: + - Gaussian(3) = 4 * exp(-0.5 * ((3-5)/1)^2) = 4 * exp(-2) ~ 0.541 + """ + x0, mu, sigma, A = 3.0, 5.0, 1.0, 4.0 + + x_test = np.array([x0]) + expected = gaussian(x_test, mu, sigma, A) # H(0) = 1.0 + + x_dummy = np.array([0.0, 10.0]) + y_dummy = np.array([0.0, 4.0]) + obj = GaussianTimesHeavySide(x_dummy, y_dummy, guess_x0=x0) + result = obj.function(x_test, x0, mu, sigma, A) + + np.testing.assert_allclose( + result, + expected, + rtol=1e-10, + err_msg="At x0: f(x0) should equal Gaussian(x0) since H(0)=1.0", + ) + + def test_func_evaluates_at_gaussian_peak(self): + """At x = mu: f(mu) = A (assuming mu > x0). + + With x0=3, mu=5, A=4: + - Gaussian(5) = 4 * exp(0) = 4.0 + - H(5-3) = H(2) = 1.0 + - f(5) = 4.0 * 1.0 = 4.0 + """ + x0, mu, sigma, A = 3.0, 5.0, 1.0, 4.0 + + x_test = np.array([mu]) + expected = np.array([A]) # Full Gaussian amplitude + + x_dummy = np.array([0.0, 10.0]) + y_dummy = np.array([0.0, 4.0]) + obj = GaussianTimesHeavySide(x_dummy, y_dummy, guess_x0=x0) + result = obj.function(x_test, x0, mu, sigma, A) + + np.testing.assert_allclose( + result, + expected, + rtol=1e-10, + err_msg="At Gaussian peak (mu > x0): f(mu) = A", + ) + + # ------------------------------------------------------------------------- + # E2. Parameter Recovery Tests (Clean Data) + # ------------------------------------------------------------------------- + + def test_fit_recovers_exact_parameters_from_clean_data(self, gth_clean_data): + """Fitting noise-free data should recover parameters within 2%.""" + x, y, w, true_params = gth_clean_data + + obj = GaussianTimesHeavySide(x, y, guess_x0=true_params["x0"]) + obj.make_fit() + + for param, true_val in true_params.items(): + fitted_val = obj.popt[param] + if abs(true_val) < 0.1: + assert abs(fitted_val - true_val) < 0.05, ( + f"{param}: fitted={fitted_val:.4f}, true={true_val:.4f}, " + f"absolute error exceeds 0.05" + ) + else: + rel_error = abs(fitted_val - true_val) / abs(true_val) + assert rel_error < 0.02, ( + f"{param}: fitted={fitted_val:.4f}, true={true_val:.4f}, " + f"rel_error={rel_error:.2%} exceeds 2% tolerance" + ) + + @pytest.mark.parametrize( + "true_params", + [ + {"x0": 3.0, "mu": 5.0, "sigma": 1.0, "A": 4.0}, + {"x0": 2.0, "mu": 6.0, "sigma": 0.8, "A": 3.0}, + {"x0": 4.0, "mu": 7.0, "sigma": 1.5, "A": 5.0}, + ], + ) + def test_fit_recovers_various_parameter_combinations(self, true_params): + """Fitting should work for diverse parameter combinations.""" + x = np.linspace(0, 12, 250) + gauss = gaussian(x, true_params["mu"], true_params["sigma"], true_params["A"]) + y = gauss * np.heaviside(x - true_params["x0"], 1.0) + + obj = GaussianTimesHeavySide(x, y, guess_x0=true_params["x0"]) + obj.make_fit() + + for param, true_val in true_params.items(): + fitted_val = obj.popt[param] + if abs(true_val) < 0.1: + assert ( + abs(fitted_val - true_val) < 0.1 + ), f"{param}: fitted={fitted_val:.4f}, true={true_val:.4f}" + else: + rel_error = abs(fitted_val - true_val) / abs(true_val) + assert rel_error < 0.05, ( + f"{param}: fitted={fitted_val:.4f}, true={true_val:.4f}, " + f"rel_error={rel_error:.2%}" + ) + + # ------------------------------------------------------------------------- + # E3. Noisy Data Tests (Precision Bounds) + # ------------------------------------------------------------------------- + + def test_fit_with_noise_recovers_gaussian_parameters_within_2sigma( + self, gth_noisy_data + ): + """Gaussian parameters should be within 2sigma of true values (95% confidence). + + Note: x0 is excluded because Heaviside functions have zero gradient almost + everywhere, making its uncertainty unreliable for this type of test. + """ + x, y, w, true_params = gth_noisy_data + + obj = GaussianTimesHeavySide(x, y, guess_x0=true_params["x0"]) + obj.make_fit() + + # Only test Gaussian parameters which have well-defined uncertainties + for param in ["mu", "sigma", "A"]: + true_val = true_params[param] + fitted_val = obj.popt[param] + sigma = obj.psigma[param] + deviation = abs(fitted_val - true_val) + + assert deviation < 2 * sigma, ( + f"{param}: |fitted({fitted_val:.4f}) - true({true_val:.4f})| = " + f"{deviation:.4f} exceeds 2sigma = {2*sigma:.4f}" + ) + + def test_fit_uncertainty_scales_with_noise(self): + """Higher noise should produce larger parameter uncertainties.""" + rng = np.random.default_rng(42) + true_params = {"x0": 3.0, "mu": 5.0, "sigma": 1.0, "A": 4.0} + x = np.linspace(0, 10, 200) + gauss = gaussian(x, true_params["mu"], true_params["sigma"], true_params["A"]) + y_true = gauss * np.heaviside(x - true_params["x0"], 1.0) + + # Low noise fit + y_low = y_true + rng.normal(0, 0.05, len(x)) + y_low[x < true_params["x0"]] = np.abs(y_low[x < true_params["x0"]]) + fit_low = GaussianTimesHeavySide(x, y_low, guess_x0=true_params["x0"]) + fit_low.make_fit() + + # High noise fit + rng2 = np.random.default_rng(43) + y_high = y_true + rng2.normal(0, 0.25, len(x)) + y_high[x < true_params["x0"]] = np.abs(y_high[x < true_params["x0"]]) + fit_high = GaussianTimesHeavySide(x, y_high, guess_x0=true_params["x0"]) + fit_high.make_fit() + + # High noise should have larger uncertainties for most parameters + for param in ["mu", "sigma", "A"]: + assert fit_high.psigma[param] > fit_low.psigma[param], ( + f"{param}: high_noise_sigma={fit_high.psigma[param]:.4f} should be " + f"> low_noise_sigma={fit_low.psigma[param]:.4f}" + ) + + # ------------------------------------------------------------------------- + # E4. Initial Parameter Estimation Tests + # ------------------------------------------------------------------------- + + def test_p0_returns_list_with_correct_length(self, gth_clean_data): + """p0 should return a list with 4 elements.""" + x, y, w, true_params = gth_clean_data + obj = GaussianTimesHeavySide(x, y, guess_x0=true_params["x0"]) + + p0 = obj.p0 + assert isinstance(p0, list), f"p0 should be a list, got {type(p0)}" + assert len(p0) == 4, f"p0 should have 4 elements, got {len(p0)}" + + def test_p0_enables_successful_convergence(self, gth_noisy_data): + """Fit should converge when initialized with estimated p0.""" + x, y, w, true_params = gth_noisy_data + + obj = GaussianTimesHeavySide(x, y, guess_x0=true_params["x0"]) + obj.make_fit() + + assert all( + np.isfinite(v) for v in obj.popt.values() + ), f"Fit did not converge: popt={obj.popt}" + + def test_guess_x0_is_required(self): + """Test that guess_x0 parameter is required for initialization.""" + x = np.array([0.0, 5.0, 10.0]) + y = np.array([0.0, 4.0, 1.0]) + + # This should either raise an error or require guess_x0 + # The exact behavior depends on implementation + with pytest.raises((TypeError, ValueError)): + obj = GaussianTimesHeavySide(x, y) # Missing guess_x0 + + # ------------------------------------------------------------------------- + # E5. Edge Case and Error Handling Tests + # ------------------------------------------------------------------------- + + def test_insufficient_data_raises_error(self): + """Fitting with insufficient data should raise InsufficientDataError.""" + x = np.array([1.0, 2.0, 3.0]) # Only 3 points for 4 parameters + y = np.array([0.0, 0.0, 1.0]) + + obj = GaussianTimesHeavySide(x, y, guess_x0=2.0) + + with pytest.raises(InsufficientDataError): + obj.make_fit() + + def test_function_signature(self): + """Test that function has correct parameter signature.""" + x = np.array([0.0, 5.0, 10.0]) + y = np.array([0.0, 4.0, 1.0]) + obj = GaussianTimesHeavySide(x, y, guess_x0=3.0) + + sig = inspect.signature(obj.function) + params = tuple(sig.parameters.keys()) + + assert params == ( + "x", + "x0", + "mu", + "sigma", + "A", + ), f"Function should have signature (x, x0, mu, sigma, A), got {params}" + + def test_callable_interface(self, gth_clean_data): + """Test that fitted object is callable and returns correct shape.""" + x, y, w, true_params = gth_clean_data + + obj = GaussianTimesHeavySide(x, y, guess_x0=true_params["x0"]) + obj.make_fit() + + x_test = np.array([1.0, 5.0, 10.0]) + y_pred = obj(x_test) + + assert ( + y_pred.shape == x_test.shape + ), f"Predicted shape {y_pred.shape} should match input shape {x_test.shape}" + assert np.all(np.isfinite(y_pred)), "All predicted values should be finite" + + def test_popt_has_correct_keys(self, gth_clean_data): + """Test that popt contains expected parameter names.""" + x, y, w, true_params = gth_clean_data + + obj = GaussianTimesHeavySide(x, y, guess_x0=true_params["x0"]) + obj.make_fit() + + expected_keys = {"x0", "mu", "sigma", "A"} + actual_keys = set(obj.popt.keys()) + + assert ( + actual_keys == expected_keys + ), f"popt keys should be {expected_keys}, got {actual_keys}" + + def test_psigma_values_are_nonnegative(self, gth_noisy_data): + """Test that parameter uncertainties are non-negative. + + Note: x0 uncertainty can be zero or very small due to the step function + nature (zero gradient almost everywhere), so we only check for non-negative + values. Gaussian parameters (mu, sigma, A) should have positive uncertainties. + """ + x, y, w, true_params = gth_noisy_data + + obj = GaussianTimesHeavySide(x, y, guess_x0=true_params["x0"]) + obj.make_fit() + + for param, sigma in obj.psigma.items(): + assert sigma >= 0, f"psigma['{param}'] = {sigma} should be non-negative" + + # Gaussian parameters should have positive uncertainties + for param in ["mu", "sigma", "A"]: + assert obj.psigma[param] > 0, f"psigma['{param}'] should be positive" + + +# ============================================================================= +# GaussianTimesHeavySidePlusHeavySide Tests +# ============================================================================= + + +class TestGaussianTimesHeavySidePlusHeavySide: + """Tests for GaussianTimesHeavySidePlusHeavySide fit function. + + GaussianTimesHeavySidePlusHeavySide models: + f(x) = Gaussian(x, mu, sigma, A) * H(x-x0) + y1 * H(x0-x) + + This produces: + - Constant y1 for x < x0 + - Gaussian(x) for x > x0 + - Transition at x = x0 (depends on Heaviside convention) + """ + + # ------------------------------------------------------------------------- + # E1. Function Evaluation Tests (Exact Values) + # ------------------------------------------------------------------------- + + def test_func_evaluates_below_x0_as_constant(self): + """For x < x0: f(x) = y1 (constant plateau). + + With x0=3, y1=2: + - H(3-x) = 1 for x < 3 + - H(x-3) = 0 for x < 3 + - f(x) = 0 + y1*1 = 2 + """ + x0, y1, mu, sigma, A = 3.0, 2.0, 5.0, 1.0, 4.0 + + x_test = np.array([0.0, 1.0, 2.0, 2.99]) # All below x0=3 + expected = np.full_like(x_test, y1) + + x_dummy = np.array([0.0, 10.0]) + y_dummy = np.array([2.0, 4.0]) + obj = GaussianTimesHeavySidePlusHeavySide(x_dummy, y_dummy, guess_x0=x0) + result = obj.function(x_test, x0, y1, mu, sigma, A) + + np.testing.assert_allclose( + result, + expected, + rtol=1e-10, + err_msg="Below x0: f(x) should be constant y1", + ) + + def test_func_evaluates_above_x0_as_gaussian(self): + """For x > x0: f(x) = Gaussian(x) (H(x0-x) = 0). + + With x0=3, mu=5, sigma=1, A=4: + - H(3-x) = 0 for x > 3 + - H(x-3) = 1 for x > 3 + - f(x) = Gaussian(x) * 1 + y1 * 0 = Gaussian(x) + """ + x0, y1, mu, sigma, A = 3.0, 2.0, 5.0, 1.0, 4.0 + + x_test = np.array([4.0, 5.0, 6.0, 8.0]) # All above x0=3 + expected = gaussian(x_test, mu, sigma, A) + + x_dummy = np.array([0.0, 10.0]) + y_dummy = np.array([2.0, 4.0]) + obj = GaussianTimesHeavySidePlusHeavySide(x_dummy, y_dummy, guess_x0=x0) + result = obj.function(x_test, x0, y1, mu, sigma, A) + + np.testing.assert_allclose( + result, + expected, + rtol=1e-10, + err_msg="Above x0: f(x) should equal Gaussian(x)", + ) + + def test_func_evaluates_at_x0_correctly(self): + """At x = x0: f(x) = Gaussian(x0)*1.0 + y1*1.0. + + This implementation uses H(0) = 1.0 for both Heaviside terms, so at the + transition point both contribute fully. + With x0=3, y1=2, mu=5, sigma=1, A=4: + - Gaussian(3) = 4 * exp(-2) ~ 0.541 + - f(3) = 0.541*1.0 + 2*1.0 = 2.541 + """ + x0, y1, mu, sigma, A = 3.0, 2.0, 5.0, 1.0, 4.0 + + x_test = np.array([x0]) + gauss_val = gaussian(x_test, mu, sigma, A) + expected = gauss_val * 1.0 + y1 * 1.0 # H(0) = 1.0 for both terms + + x_dummy = np.array([0.0, 10.0]) + y_dummy = np.array([2.0, 4.0]) + obj = GaussianTimesHeavySidePlusHeavySide(x_dummy, y_dummy, guess_x0=x0) + result = obj.function(x_test, x0, y1, mu, sigma, A) + + np.testing.assert_allclose( + result, + expected, + rtol=1e-10, + err_msg="At x0: f(x0) should equal Gaussian(x0)*1.0 + y1*1.0", + ) + + def test_func_evaluates_at_gaussian_peak(self): + """At x = mu: f(mu) = A (assuming mu > x0). + + With x0=3, mu=5, A=4: + - Gaussian(5) = 4 * exp(0) = 4.0 + - H(5-3) = H(2) = 1.0 + - H(3-5) = H(-2) = 0.0 + - f(5) = 4.0 * 1.0 + y1 * 0.0 = 4.0 + """ + x0, y1, mu, sigma, A = 3.0, 2.0, 5.0, 1.0, 4.0 + + x_test = np.array([mu]) + expected = np.array([A]) # Full Gaussian amplitude + + x_dummy = np.array([0.0, 10.0]) + y_dummy = np.array([2.0, 4.0]) + obj = GaussianTimesHeavySidePlusHeavySide(x_dummy, y_dummy, guess_x0=x0) + result = obj.function(x_test, x0, y1, mu, sigma, A) + + np.testing.assert_allclose( + result, + expected, + rtol=1e-10, + err_msg="At Gaussian peak (mu > x0): f(mu) = A", + ) + + # ------------------------------------------------------------------------- + # E2. Parameter Recovery Tests (Clean Data) + # ------------------------------------------------------------------------- + + def test_fit_recovers_exact_parameters_from_clean_data(self, gthph_clean_data): + """Fitting noise-free data should recover parameters within 2%.""" + x, y, w, true_params = gthph_clean_data + + obj = GaussianTimesHeavySidePlusHeavySide(x, y, guess_x0=true_params["x0"]) + obj.make_fit() + + for param, true_val in true_params.items(): + fitted_val = obj.popt[param] + if abs(true_val) < 0.1: + assert abs(fitted_val - true_val) < 0.05, ( + f"{param}: fitted={fitted_val:.4f}, true={true_val:.4f}, " + f"absolute error exceeds 0.05" + ) + else: + rel_error = abs(fitted_val - true_val) / abs(true_val) + assert rel_error < 0.02, ( + f"{param}: fitted={fitted_val:.4f}, true={true_val:.4f}, " + f"rel_error={rel_error:.2%} exceeds 2% tolerance" + ) + + @pytest.mark.parametrize( + "true_params", + [ + {"x0": 3.0, "y1": 2.0, "mu": 5.0, "sigma": 1.0, "A": 4.0}, + {"x0": 2.0, "y1": 1.5, "mu": 6.0, "sigma": 0.8, "A": 3.0}, + {"x0": 4.0, "y1": 3.0, "mu": 7.0, "sigma": 1.5, "A": 5.0}, + ], + ) + def test_fit_recovers_various_parameter_combinations(self, true_params): + """Fitting should work for diverse parameter combinations.""" + x = np.linspace(0, 12, 250) + gauss = gaussian(x, true_params["mu"], true_params["sigma"], true_params["A"]) + y = gauss * np.heaviside(x - true_params["x0"], 0.5) + true_params[ + "y1" + ] * np.heaviside(true_params["x0"] - x, 0.5) + + obj = GaussianTimesHeavySidePlusHeavySide(x, y, guess_x0=true_params["x0"]) + obj.make_fit() + + for param, true_val in true_params.items(): + fitted_val = obj.popt[param] + if abs(true_val) < 0.1: + assert ( + abs(fitted_val - true_val) < 0.1 + ), f"{param}: fitted={fitted_val:.4f}, true={true_val:.4f}" + else: + rel_error = abs(fitted_val - true_val) / abs(true_val) + assert rel_error < 0.05, ( + f"{param}: fitted={fitted_val:.4f}, true={true_val:.4f}, " + f"rel_error={rel_error:.2%}" + ) + + # ------------------------------------------------------------------------- + # E3. Noisy Data Tests (Precision Bounds) + # ------------------------------------------------------------------------- + + def test_fit_with_noise_recovers_gaussian_parameters_within_2sigma( + self, gthph_noisy_data + ): + """Gaussian parameters should be within 2sigma of true values (95% confidence). + + Note: x0 and y1 are excluded because Heaviside functions have zero gradient + almost everywhere, making their uncertainties unreliable for this type of test. + """ + x, y, w, true_params = gthph_noisy_data + + obj = GaussianTimesHeavySidePlusHeavySide(x, y, guess_x0=true_params["x0"]) + obj.make_fit() + + # Only test Gaussian parameters which have well-defined uncertainties + for param in ["mu", "sigma", "A"]: + true_val = true_params[param] + fitted_val = obj.popt[param] + sigma = obj.psigma[param] + deviation = abs(fitted_val - true_val) + + assert deviation < 2 * sigma, ( + f"{param}: |fitted({fitted_val:.4f}) - true({true_val:.4f})| = " + f"{deviation:.4f} exceeds 2sigma = {2*sigma:.4f}" + ) + + def test_fit_uncertainty_scales_with_noise(self): + """Higher noise should produce larger parameter uncertainties.""" + rng = np.random.default_rng(42) + true_params = {"x0": 3.0, "y1": 2.0, "mu": 5.0, "sigma": 1.0, "A": 4.0} + x = np.linspace(0, 10, 200) + gauss = gaussian(x, true_params["mu"], true_params["sigma"], true_params["A"]) + y_true = gauss * np.heaviside(x - true_params["x0"], 0.5) + true_params[ + "y1" + ] * np.heaviside(true_params["x0"] - x, 0.5) + + # Low noise fit + y_low = y_true + rng.normal(0, 0.05, len(x)) + fit_low = GaussianTimesHeavySidePlusHeavySide( + x, y_low, guess_x0=true_params["x0"] + ) + fit_low.make_fit() + + # High noise fit + rng2 = np.random.default_rng(43) + y_high = y_true + rng2.normal(0, 0.25, len(x)) + fit_high = GaussianTimesHeavySidePlusHeavySide( + x, y_high, guess_x0=true_params["x0"] + ) + fit_high.make_fit() + + # High noise should have larger uncertainties for most parameters + for param in ["mu", "sigma", "A"]: + assert fit_high.psigma[param] > fit_low.psigma[param], ( + f"{param}: high_noise_sigma={fit_high.psigma[param]:.4f} should be " + f"> low_noise_sigma={fit_low.psigma[param]:.4f}" + ) + + # ------------------------------------------------------------------------- + # E4. Initial Parameter Estimation Tests + # ------------------------------------------------------------------------- + + def test_p0_returns_list_with_correct_length(self, gthph_clean_data): + """p0 should return a list with 5 elements.""" + x, y, w, true_params = gthph_clean_data + obj = GaussianTimesHeavySidePlusHeavySide(x, y, guess_x0=true_params["x0"]) + + p0 = obj.p0 + assert isinstance(p0, list), f"p0 should be a list, got {type(p0)}" + assert len(p0) == 5, f"p0 should have 5 elements, got {len(p0)}" + + def test_p0_enables_successful_convergence(self, gthph_noisy_data): + """Fit should converge when initialized with estimated p0.""" + x, y, w, true_params = gthph_noisy_data + + obj = GaussianTimesHeavySidePlusHeavySide(x, y, guess_x0=true_params["x0"]) + obj.make_fit() + + assert all( + np.isfinite(v) for v in obj.popt.values() + ), f"Fit did not converge: popt={obj.popt}" + + def test_guess_x0_is_required(self): + """Test that guess_x0 parameter is required for initialization.""" + x = np.array([0.0, 5.0, 10.0]) + y = np.array([2.0, 4.0, 1.0]) + + # This should either raise an error or require guess_x0 + with pytest.raises((TypeError, ValueError)): + obj = GaussianTimesHeavySidePlusHeavySide(x, y) # Missing guess_x0 + + # ------------------------------------------------------------------------- + # E5. Edge Case and Error Handling Tests + # ------------------------------------------------------------------------- + + def test_insufficient_data_raises_error(self): + """Fitting with insufficient data should raise InsufficientDataError.""" + x = np.array([1.0, 2.0, 3.0, 4.0]) # Only 4 points for 5 parameters + y = np.array([2.0, 2.0, 2.0, 3.0]) + + obj = GaussianTimesHeavySidePlusHeavySide(x, y, guess_x0=2.5) + + with pytest.raises(InsufficientDataError): + obj.make_fit() + + def test_function_signature(self): + """Test that function has correct parameter signature.""" + x = np.array([0.0, 5.0, 10.0]) + y = np.array([2.0, 4.0, 1.0]) + obj = GaussianTimesHeavySidePlusHeavySide(x, y, guess_x0=3.0) + + sig = inspect.signature(obj.function) + params = tuple(sig.parameters.keys()) + + assert params == ( + "x", + "x0", + "y1", + "mu", + "sigma", + "A", + ), f"Function should have signature (x, x0, y1, mu, sigma, A), got {params}" + + def test_callable_interface(self, gthph_clean_data): + """Test that fitted object is callable and returns correct shape.""" + x, y, w, true_params = gthph_clean_data + + obj = GaussianTimesHeavySidePlusHeavySide(x, y, guess_x0=true_params["x0"]) + obj.make_fit() + + x_test = np.array([1.0, 5.0, 10.0]) + y_pred = obj(x_test) + + assert ( + y_pred.shape == x_test.shape + ), f"Predicted shape {y_pred.shape} should match input shape {x_test.shape}" + assert np.all(np.isfinite(y_pred)), "All predicted values should be finite" + + def test_popt_has_correct_keys(self, gthph_clean_data): + """Test that popt contains expected parameter names.""" + x, y, w, true_params = gthph_clean_data + + obj = GaussianTimesHeavySidePlusHeavySide(x, y, guess_x0=true_params["x0"]) + obj.make_fit() + + expected_keys = {"x0", "y1", "mu", "sigma", "A"} + actual_keys = set(obj.popt.keys()) + + assert ( + actual_keys == expected_keys + ), f"popt keys should be {expected_keys}, got {actual_keys}" + + def test_psigma_values_are_nonnegative(self, gthph_noisy_data): + """Test that parameter uncertainties are non-negative. + + Note: x0 uncertainty can be zero or very small due to the step function + nature (zero gradient almost everywhere), so we only check for non-negative + values. Gaussian parameters (mu, sigma, A) should have positive uncertainties. + """ + x, y, w, true_params = gthph_noisy_data + + obj = GaussianTimesHeavySidePlusHeavySide(x, y, guess_x0=true_params["x0"]) + obj.make_fit() + + for param, sigma in obj.psigma.items(): + assert sigma >= 0, f"psigma['{param}'] = {sigma} should be non-negative" + + # Gaussian parameters should have positive uncertainties + for param in ["mu", "sigma", "A"]: + assert obj.psigma[param] > 0, f"psigma['{param}'] should be positive" + + # ------------------------------------------------------------------------- + # E6. Behavior Consistency Tests + # ------------------------------------------------------------------------- + + def test_transition_continuity(self, gthph_clean_data): + """Test that the function shows expected behavior at transition. + + The function transitions from constant y1 (for x < x0) to Gaussian (for x > x0). + At x = x0, both Heaviside functions contribute 0.5. + """ + x, y, w, true_params = gthph_clean_data + + obj = GaussianTimesHeavySidePlusHeavySide(x, y, guess_x0=true_params["x0"]) + obj.make_fit() + + x0 = obj.popt["x0"] + y1 = obj.popt["y1"] + mu = obj.popt["mu"] + sigma = obj.popt["sigma"] + A = obj.popt["A"] + + # Value just below x0 + x_below = np.array([x0 - 0.01]) + y_below = obj(x_below)[0] + + # Value just above x0 + x_above = np.array([x0 + 0.01]) + y_above = obj(x_above)[0] + + # Below x0 should be close to y1 + np.testing.assert_allclose( + y_below, + y1, + rtol=0.1, + err_msg=f"Value just below x0 ({y_below:.4f}) should be close to y1 ({y1:.4f})", + ) + + # Above x0 should be close to Gaussian(x0) + gauss_at_x0 = gaussian(np.array([x0 + 0.01]), mu, sigma, A)[0] + np.testing.assert_allclose( + y_above, + gauss_at_x0, + rtol=0.1, + err_msg=f"Value just above x0 ({y_above:.4f}) should be close to Gaussian ({gauss_at_x0:.4f})", + ) + + def test_plateau_region_is_constant(self, gthph_clean_data): + """Test that the region x < x0 is a constant plateau at y1.""" + x, y, w, true_params = gthph_clean_data + + obj = GaussianTimesHeavySidePlusHeavySide(x, y, guess_x0=true_params["x0"]) + obj.make_fit() + + x0 = obj.popt["x0"] + y1 = obj.popt["y1"] + + # Test multiple points below x0 + x_plateau = np.array([0.5, 1.0, 1.5, 2.0, 2.5]) + x_plateau = x_plateau[x_plateau < x0] # Ensure all below x0 + + if len(x_plateau) > 0: + y_plateau = obj(x_plateau) + expected = np.full_like(y_plateau, y1) + + np.testing.assert_allclose( + y_plateau, + expected, + rtol=1e-6, + err_msg="Plateau region (x < x0) should be constant at y1", + ) diff --git a/tests/fitfunctions/test_heaviside.py b/tests/fitfunctions/test_heaviside.py new file mode 100644 index 00000000..c3885326 --- /dev/null +++ b/tests/fitfunctions/test_heaviside.py @@ -0,0 +1,653 @@ +"""Tests for HeavySide fit function. + +HeavySide models a step function (Heaviside step function) with parameters: +- x0: transition point (step location) +- y0: baseline level (value for x > x0) +- y1: step height (added to y0 for x < x0) + +The function is defined as: + f(x) = y1 * heaviside(x0 - x, 0.5*(y0+y1)) + y0 + +Behavior: +- For x < x0: heaviside(x0-x) = 1, so f(x) = y1 + y0 +- For x > x0: heaviside(x0-x) = 0, so f(x) = y0 +- For x == x0: heaviside(0, 0.5*(y0+y1)) = 0.5*(y0+y1), so f(x) = y1*0.5*(y0+y1) + y0 + +Note: The p0 property provides heuristic-based initial guesses, but for +best results you may want to provide manual p0 via make_fit(p0=[x0, y0, y1]). +""" + +import inspect + +import numpy as np +import pytest + +from solarwindpy.fitfunctions.heaviside import HeavySide +from solarwindpy.fitfunctions.core import InsufficientDataError + + +# ============================================================================= +# Fixtures +# ============================================================================= + + +@pytest.fixture +def clean_step_data(): + """Perfect step function data (no noise). + + Parameters: x0=5.0, y0=2.0, y1=3.0 + - For x < 5: f(x) = 3 + 2 = 5 + - For x > 5: f(x) = 2 + """ + true_params = {"x0": 5.0, "y0": 2.0, "y1": 3.0} + x = np.linspace(0, 10, 201) # Odd number to avoid x=5 exactly except at midpoint + # Build y using the heaviside formula + y = ( + true_params["y1"] + * np.heaviside( + true_params["x0"] - x, 0.5 * (true_params["y0"] + true_params["y1"]) + ) + + true_params["y0"] + ) + w = np.ones_like(x) + return x, y, w, true_params + + +@pytest.fixture +def noisy_step_data(): + """Step function data with 5% Gaussian noise. + + Parameters: x0=5.0, y0=2.0, y1=3.0 + Noise std = 0.25 (5% of step height + baseline) + """ + rng = np.random.default_rng(42) + true_params = {"x0": 5.0, "y0": 2.0, "y1": 3.0} + x = np.linspace(0, 10, 200) + y_true = ( + true_params["y1"] + * np.heaviside( + true_params["x0"] - x, 0.5 * (true_params["y0"] + true_params["y1"]) + ) + + true_params["y0"] + ) + noise_std = 0.25 + y = y_true + rng.normal(0, noise_std, len(x)) + w = np.ones_like(x) / noise_std + return x, y, w, true_params + + +@pytest.fixture +def negative_step_data(): + """Step function with negative y1 (step down instead of step up). + + Parameters: x0=5.0, y0=8.0, y1=-3.0 + - For x < 5: f(x) = -3 + 8 = 5 + - For x > 5: f(x) = 8 + """ + true_params = {"x0": 5.0, "y0": 8.0, "y1": -3.0} + x = np.linspace(0, 10, 201) + y = ( + true_params["y1"] + * np.heaviside( + true_params["x0"] - x, 0.5 * (true_params["y0"] + true_params["y1"]) + ) + + true_params["y0"] + ) + w = np.ones_like(x) + return x, y, w, true_params + + +@pytest.fixture +def zero_baseline_data(): + """Step function with y0=0 (baseline at zero). + + Parameters: x0=3.0, y0=0.0, y1=4.0 + - For x < 3: f(x) = 4 + 0 = 4 + - For x > 3: f(x) = 0 + """ + true_params = {"x0": 3.0, "y0": 0.0, "y1": 4.0} + x = np.linspace(0, 10, 201) + y = ( + true_params["y1"] + * np.heaviside( + true_params["x0"] - x, 0.5 * (true_params["y0"] + true_params["y1"]) + ) + + true_params["y0"] + ) + w = np.ones_like(x) + return x, y, w, true_params + + +# ============================================================================= +# E1. Function Evaluation Tests (Exact Values) +# ============================================================================= + + +def test_func_evaluates_below_step_correctly(): + """For x < x0: f(x) = y1 + y0. + + With x0=5, y0=2, y1=3: f(x<5) = 3 + 2 = 5. + """ + x0, y0, y1 = 5.0, 2.0, 3.0 + + # Test specific points below the step transition + x_test = np.array([0.0, 1.0, 2.5, 4.0, 4.99]) + expected = np.full_like(x_test, y0 + y1) # All should equal 5.0 + + # Create minimal instance to access function + x_dummy = np.array([0.0, 10.0]) + y_dummy = np.array([5.0, 2.0]) + obj = HeavySide(x_dummy, y_dummy) + result = obj.function(x_test, x0, y0, y1) + + np.testing.assert_allclose( + result, + expected, + rtol=1e-10, + err_msg="Below step (x < x0): f(x) should equal y0 + y1", + ) + + +def test_func_evaluates_above_step_correctly(): + """For x > x0: f(x) = y0. + + With x0=5, y0=2, y1=3: f(x>5) = 2. + """ + x0, y0, y1 = 5.0, 2.0, 3.0 + + # Test points above the step transition + x_test = np.array([5.01, 6.0, 7.5, 10.0, 100.0]) + expected = np.full_like(x_test, y0) # All should equal 2.0 + + x_dummy = np.array([0.0, 10.0]) + y_dummy = np.array([5.0, 2.0]) + obj = HeavySide(x_dummy, y_dummy) + result = obj.function(x_test, x0, y0, y1) + + np.testing.assert_allclose( + result, + expected, + rtol=1e-10, + err_msg="Above step (x > x0): f(x) should equal y0", + ) + + +def test_func_evaluates_at_step_transition(): + """At x == x0: f(x) = y1 * 0.5*(y0+y1) + y0. + + With x0=5, y0=2, y1=3: + f(5) = 3 * 0.5*(2+3) + 2 = 3 * 2.5 + 2 = 7.5 + 2 = 9.5 + + Note: This is unusual behavior for a step function. The typical + midpoint would be 0.5*(y0 + y0+y1) = y0 + 0.5*y1 = 3.5. + """ + x0, y0, y1 = 5.0, 2.0, 3.0 + + x_test = np.array([5.0]) + # f(x0) = y1 * heaviside(0, 0.5*(y0+y1)) + y0 + # = y1 * 0.5*(y0+y1) + y0 + expected_at_transition = y1 * 0.5 * (y0 + y1) + y0 + expected = np.array([expected_at_transition]) + + x_dummy = np.array([0.0, 10.0]) + y_dummy = np.array([5.0, 2.0]) + obj = HeavySide(x_dummy, y_dummy) + result = obj.function(x_test, x0, y0, y1) + + np.testing.assert_allclose( + result, + expected, + rtol=1e-10, + err_msg=f"At step (x == x0): f(x) should equal {expected_at_transition:.2f}", + ) + + +def test_func_with_negative_step_height(): + """Step function with negative y1 (step down from left to right). + + With x0=5, y0=8, y1=-3: + - For x < 5: f(x) = -3 + 8 = 5 + - For x > 5: f(x) = 8 + """ + x0, y0, y1 = 5.0, 8.0, -3.0 + + x_test = np.array([2.0, 4.9, 5.1, 8.0]) + expected = np.array([y0 + y1, y0 + y1, y0, y0]) # [5, 5, 8, 8] + + x_dummy = np.array([0.0, 10.0]) + y_dummy = np.array([5.0, 8.0]) + obj = HeavySide(x_dummy, y_dummy) + result = obj.function(x_test, x0, y0, y1) + + np.testing.assert_allclose( + result, + expected, + rtol=1e-10, + err_msg="Negative y1: step down should give y0+y1 below, y0 above", + ) + + +def test_func_with_zero_baseline(): + """Step function with y0=0. + + With x0=3, y0=0, y1=4: + - For x < 3: f(x) = 4 + 0 = 4 + - For x > 3: f(x) = 0 + """ + x0, y0, y1 = 3.0, 0.0, 4.0 + + x_test = np.array([1.0, 2.9, 3.1, 5.0]) + expected = np.array([y0 + y1, y0 + y1, y0, y0]) # [4, 4, 0, 0] + + x_dummy = np.array([0.0, 10.0]) + y_dummy = np.array([4.0, 0.0]) + obj = HeavySide(x_dummy, y_dummy) + result = obj.function(x_test, x0, y0, y1) + + np.testing.assert_allclose( + result, + expected, + rtol=1e-10, + err_msg="Zero baseline: should give y1 below, 0 above", + ) + + +# ============================================================================= +# E2. Parameter Recovery Tests (Clean Data with Manual p0) +# ============================================================================= + + +def test_fit_recovers_exact_parameters_from_clean_data(clean_step_data): + """Fitting noise-free data with manual p0 should recover parameters within 1%.""" + x, y, w, true_params = clean_step_data + + obj = HeavySide(x, y) + # Must provide p0 manually since automatic p0 raises NotImplementedError + p0 = [true_params["x0"], true_params["y0"], true_params["y1"]] + obj.make_fit(p0=p0) + + for param, true_val in true_params.items(): + fitted_val = obj.popt[param] + # Use absolute tolerance for values near zero + if abs(true_val) < 0.1: + assert abs(fitted_val - true_val) < 0.05, ( + f"{param}: fitted={fitted_val:.4f}, true={true_val:.4f}, " + f"absolute error exceeds 0.05" + ) + else: + rel_error = abs(fitted_val - true_val) / abs(true_val) + assert rel_error < 0.01, ( + f"{param}: fitted={fitted_val:.4f}, true={true_val:.4f}, " + f"rel_error={rel_error:.2%} exceeds 1% tolerance" + ) + + +def test_fit_recovers_negative_step_parameters(negative_step_data): + """Fitting clean data with negative y1 should recover parameters within 2%.""" + x, y, w, true_params = negative_step_data + + obj = HeavySide(x, y) + p0 = [true_params["x0"], true_params["y0"], true_params["y1"]] + obj.make_fit(p0=p0) + + for param, true_val in true_params.items(): + fitted_val = obj.popt[param] + if abs(true_val) < 0.1: + assert ( + abs(fitted_val - true_val) < 0.05 + ), f"{param}: fitted={fitted_val:.4f}, true={true_val:.4f}" + else: + rel_error = abs(fitted_val - true_val) / abs(true_val) + assert rel_error < 0.02, ( + f"{param}: fitted={fitted_val:.4f}, true={true_val:.4f}, " + f"rel_error={rel_error:.2%}" + ) + + +def test_fit_recovers_zero_baseline_parameters(zero_baseline_data): + """Fitting clean data with y0=0 should recover parameters within 2%.""" + x, y, w, true_params = zero_baseline_data + + obj = HeavySide(x, y) + p0 = [true_params["x0"], true_params["y0"], true_params["y1"]] + obj.make_fit(p0=p0) + + for param, true_val in true_params.items(): + fitted_val = obj.popt[param] + if abs(true_val) < 0.1: + # For y0=0, check absolute tolerance + assert ( + abs(fitted_val - true_val) < 0.05 + ), f"{param}: fitted={fitted_val:.4f}, true={true_val:.4f}" + else: + rel_error = abs(fitted_val - true_val) / abs(true_val) + assert rel_error < 0.02, ( + f"{param}: fitted={fitted_val:.4f}, true={true_val:.4f}, " + f"rel_error={rel_error:.2%}" + ) + + +@pytest.mark.parametrize( + "true_params", + [ + {"x0": 5.0, "y0": 2.0, "y1": 3.0}, # Standard step up + {"x0": 3.0, "y0": 10.0, "y1": -5.0}, # Step down + {"x0": 7.0, "y0": 0.0, "y1": 6.0}, # Zero baseline + {"x0": 2.0, "y0": 1.0, "y1": 0.5}, # Small step + ], +) +def test_fit_recovers_various_parameter_combinations(true_params): + """Fitting should work for diverse parameter combinations.""" + x = np.linspace(0, 10, 200) + + # Build y from parameters using the heaviside formula + y = ( + true_params["y1"] + * np.heaviside( + true_params["x0"] - x, 0.5 * (true_params["y0"] + true_params["y1"]) + ) + + true_params["y0"] + ) + + obj = HeavySide(x, y) + p0 = [true_params["x0"], true_params["y0"], true_params["y1"]] + obj.make_fit(p0=p0) + + for param, true_val in true_params.items(): + fitted_val = obj.popt[param] + if abs(true_val) < 0.1: + assert ( + abs(fitted_val - true_val) < 0.05 + ), f"{param}: fitted={fitted_val:.4f}, true={true_val:.4f}" + else: + rel_error = abs(fitted_val - true_val) / abs(true_val) + assert rel_error < 0.02, ( + f"{param}: fitted={fitted_val:.4f}, true={true_val:.4f}, " + f"rel_error={rel_error:.2%}" + ) + + +# ============================================================================= +# E3. Noisy Data Tests (Precision Bounds) +# ============================================================================= + + +def test_fit_with_noise_recovers_parameters_within_tolerance(noisy_step_data): + """Fitted parameters should be close to true values. + + Note: For step functions, the x0 parameter uncertainty can be zero or + very small because the step location is essentially a discrete choice. + We check y0 and y1 against their uncertainties, but use absolute + tolerance for x0. + """ + x, y, w, true_params = noisy_step_data + + obj = HeavySide(x, y) + p0 = [true_params["x0"], true_params["y0"], true_params["y1"]] + obj.make_fit(p0=p0) + + # Check y0 and y1 against their uncertainties (if non-zero) + for param in ["y0", "y1"]: + true_val = true_params[param] + fitted_val = obj.popt[param] + sigma = obj.psigma[param] + deviation = abs(fitted_val - true_val) + + if sigma > 0: + # 2sigma gives 95% confidence + assert deviation < 2 * sigma, ( + f"{param}: |fitted({fitted_val:.4f}) - true({true_val:.4f})| = " + f"{deviation:.4f} exceeds 2sigma = {2*sigma:.4f}" + ) + else: + # If sigma is 0, check absolute tolerance + assert deviation < 0.5, ( + f"{param}: |fitted({fitted_val:.4f}) - true({true_val:.4f})| = " + f"{deviation:.4f} exceeds absolute tolerance 0.5" + ) + + # For x0, check absolute tolerance (step location) + x0_deviation = abs(obj.popt["x0"] - true_params["x0"]) + assert x0_deviation < 0.5, ( + f"x0: |fitted({obj.popt['x0']:.4f}) - true({true_params['x0']:.4f})| = " + f"{x0_deviation:.4f} exceeds tolerance 0.5" + ) + + +def test_fit_uncertainty_scales_with_noise(): + """Higher noise should produce larger parameter uncertainties.""" + rng = np.random.default_rng(42) + true_params = {"x0": 5.0, "y0": 2.0, "y1": 3.0} + x = np.linspace(0, 10, 200) + y_true = ( + true_params["y1"] + * np.heaviside( + true_params["x0"] - x, 0.5 * (true_params["y0"] + true_params["y1"]) + ) + + true_params["y0"] + ) + + p0 = [true_params["x0"], true_params["y0"], true_params["y1"]] + + # Low noise fit + y_low = y_true + rng.normal(0, 0.1, len(x)) + fit_low = HeavySide(x, y_low) + fit_low.make_fit(p0=p0) + + # High noise fit (different seed for independence) + rng2 = np.random.default_rng(43) + y_high = y_true + rng2.normal(0, 0.5, len(x)) + fit_high = HeavySide(x, y_high) + fit_high.make_fit(p0=p0) + + # High noise should have larger uncertainties for at least some parameters + # (y0 and y1 are the primary parameters affected by noise away from transition) + for param in ["y0", "y1"]: + assert fit_high.psigma[param] > fit_low.psigma[param], ( + f"{param}: high_noise_sigma={fit_high.psigma[param]:.4f} should be " + f"> low_noise_sigma={fit_low.psigma[param]:.4f}" + ) + + +# ============================================================================= +# E4. Initial Parameter Estimation Tests +# ============================================================================= + + +def test_p0_returns_list_with_correct_length(clean_step_data): + """p0 should return a list with 3 elements.""" + x, y, w, true_params = clean_step_data + obj = HeavySide(x, y) + + p0 = obj.p0 + assert isinstance(p0, list), f"p0 should be a list, got {type(p0)}" + assert len(p0) == 3, f"p0 should have 3 elements (x0, y0, y1), got {len(p0)}" + + +def test_p0_provides_reasonable_initial_guesses(clean_step_data): + """p0 should provide reasonable heuristic-based initial guesses. + + For clean step data with x0=5, y0=2, y1=3: + - x0 guess should be near midpoint of x range + - y0 guess should be near minimum y value (2) + - y1 guess should be positive (step height estimate) + + Note: The y1 estimate may not be accurate because the HeavySide function + has an unusual value at the transition point x0 (not the simple midpoint). + The heuristic uses max(y) - min(y), which can be inflated by the + transition value y1*0.5*(y0+y1) + y0. + """ + x, y, w, true_params = clean_step_data + obj = HeavySide(x, y) + + p0 = obj.p0 + + # x0 guess should be reasonable (within data range) + assert ( + min(x) <= p0[0] <= max(x) + ), f"x0 guess {p0[0]} should be within data range [{min(x)}, {max(x)}]" + + # y0 guess should be close to minimum y (baseline) + np.testing.assert_allclose( + p0[1], + true_params["y0"], + atol=0.5, + err_msg=f"y0 guess {p0[1]} should be near true y0={true_params['y0']}", + ) + + # y1 guess should be positive and finite (allows fitting to converge) + assert p0[2] > 0, f"y1 guess {p0[2]} should be positive" + assert np.isfinite(p0[2]), f"y1 guess {p0[2]} should be finite" + + +# ============================================================================= +# E5. Derived Quantity Tests (Internal Consistency) +# ============================================================================= + + +def test_step_discontinuity_magnitude(clean_step_data): + """Verify that the step magnitude equals y1. + + The difference between values just below and just above x0 should be y1. + """ + x, y, w, true_params = clean_step_data + + obj = HeavySide(x, y) + p0 = [true_params["x0"], true_params["y0"], true_params["y1"]] + obj.make_fit(p0=p0) + + x0 = obj.popt["x0"] + y1_expected = obj.popt["y1"] + + # Evaluate just below and above the step + epsilon = 0.001 + x_below = np.array([x0 - epsilon]) + x_above = np.array([x0 + epsilon]) + + y_below = obj(x_below)[0] + y_above = obj(x_above)[0] + + step_magnitude = y_below - y_above + + np.testing.assert_allclose( + step_magnitude, + y1_expected, + rtol=1e-3, + err_msg=f"Step magnitude {step_magnitude:.4f} should equal y1={y1_expected:.4f}", + ) + + +# ============================================================================= +# E6. Edge Case and Error Handling Tests +# ============================================================================= + + +def test_insufficient_data_raises_error(): + """Fitting with insufficient data should raise InsufficientDataError.""" + x = np.array([1.0, 2.0]) # Only 2 points for 3 parameters + y = np.array([5.0, 5.0]) + + obj = HeavySide(x, y) + + with pytest.raises(InsufficientDataError): + obj.make_fit(p0=[5.0, 2.0, 3.0]) + + +def test_function_signature(): + """Test that function has correct parameter signature.""" + x = np.array([0.0, 5.0, 10.0]) + y = np.array([5.0, 3.5, 2.0]) + obj = HeavySide(x, y) + + sig = inspect.signature(obj.function) + params = tuple(sig.parameters.keys()) + + assert params == ( + "x", + "x0", + "y0", + "y1", + ), f"Function should have signature (x, x0, y0, y1), got {params}" + + +def test_tex_function_property(): + """Test that TeX_function returns a string. + + Note: The current implementation's TeX_function appears to be copied + from another class and may not accurately represent the Heaviside function. + """ + x = np.array([0.0, 5.0, 10.0]) + y = np.array([5.0, 3.5, 2.0]) + obj = HeavySide(x, y) + + tex = obj.TeX_function + assert isinstance(tex, str), f"TeX_function should return str, got {type(tex)}" + # The TeX string exists even if it's not specific to Heaviside + assert len(tex) > 0, "TeX_function should return non-empty string" + + +def test_callable_interface(clean_step_data): + """Test that fitted object is callable and returns correct shape.""" + x, y, w, true_params = clean_step_data + + obj = HeavySide(x, y) + p0 = [true_params["x0"], true_params["y0"], true_params["y1"]] + obj.make_fit(p0=p0) + + # Test callable interface + x_test = np.array([1.0, 5.0, 9.0]) + y_pred = obj(x_test) + + assert ( + y_pred.shape == x_test.shape + ), f"Predicted shape {y_pred.shape} should match input shape {x_test.shape}" + assert np.all(np.isfinite(y_pred)), "All predicted values should be finite" + + +def test_popt_has_correct_keys(clean_step_data): + """Test that popt contains expected parameter names.""" + x, y, w, true_params = clean_step_data + + obj = HeavySide(x, y) + p0 = [true_params["x0"], true_params["y0"], true_params["y1"]] + obj.make_fit(p0=p0) + + expected_keys = {"x0", "y0", "y1"} + actual_keys = set(obj.popt.keys()) + + assert ( + actual_keys == expected_keys + ), f"popt keys should be {expected_keys}, got {actual_keys}" + + +def test_psigma_has_same_keys_as_popt(noisy_step_data): + """Test that psigma has same keys as popt.""" + x, y, w, true_params = noisy_step_data + + obj = HeavySide(x, y) + p0 = [true_params["x0"], true_params["y0"], true_params["y1"]] + obj.make_fit(p0=p0) + + assert set(obj.psigma.keys()) == set(obj.popt.keys()), ( + f"psigma keys {set(obj.psigma.keys())} should match " + f"popt keys {set(obj.popt.keys())}" + ) + + +def test_psigma_values_are_nonnegative(noisy_step_data): + """Test that all parameter uncertainties are non-negative. + + Note: For step functions, the x0 parameter uncertainty can be zero + because the step location is essentially a discrete choice that the + optimizer converges to exactly. The y0 and y1 uncertainties should + typically be positive when there is noise in the data. + """ + x, y, w, true_params = noisy_step_data + + obj = HeavySide(x, y) + p0 = [true_params["x0"], true_params["y0"], true_params["y1"]] + obj.make_fit(p0=p0) + + for param, sigma in obj.psigma.items(): + assert sigma >= 0, f"psigma['{param}'] = {sigma} should be non-negative" + assert np.isfinite(sigma), f"psigma['{param}'] = {sigma} should be finite" diff --git a/tests/fitfunctions/test_hinge.py b/tests/fitfunctions/test_hinge.py new file mode 100644 index 00000000..7f1c0a34 --- /dev/null +++ b/tests/fitfunctions/test_hinge.py @@ -0,0 +1,2505 @@ +"""Tests for HingeSaturation fit function. + +HingeSaturation models a piecewise linear function with a hinge point (xh, yh): +- Rising region (x < xh): f(x) = m1 * (x - x1) where m1 = yh / (xh - x1) +- Plateau region (x >= xh): f(x) = m2 * (x - x2) where x2 = xh - yh / m2 + +Parameters: +- xh: x-coordinate of hinge point +- yh: y-coordinate of hinge point +- x1: x-intercept of rising line +- m2: slope of plateau (m2=0 gives constant saturation) +""" + +import inspect + +import numpy as np +import pytest + +from solarwindpy.fitfunctions.hinge import HingeSaturation +from solarwindpy.fitfunctions.core import InsufficientDataError + + +# ============================================================================= +# Fixtures +# ============================================================================= + + +@pytest.fixture +def clean_saturation_data(): + """Perfect hinge saturation data (no noise) with m2=0 (flat plateau). + + Parameters: xh=5.0, yh=10.0, x1=0.0, m2=0.0 + This gives m1 = 10/(5-0) = 2.0, so rising region is f(x) = 2x. + """ + true_params = {"xh": 5.0, "yh": 10.0, "x1": 0.0, "m2": 0.0} + x = np.linspace(0.1, 15, 200) + # Build y piecewise to avoid numerical issues + m1 = true_params["yh"] / (true_params["xh"] - true_params["x1"]) + y = np.where( + x < true_params["xh"], + m1 * (x - true_params["x1"]), + true_params["yh"], # m2=0 means flat plateau at yh + ) + w = np.ones_like(x) + return x, y, w, true_params + + +@pytest.fixture +def clean_sloped_plateau_data(): + """Perfect hinge data with non-zero m2 (sloped plateau). + + Parameters: xh=5.0, yh=10.0, x1=0.0, m2=0.5 + Rising: f(x) = 2*(x-0) = 2x + Plateau: f(x) = 0.5*(x - x2) where x2 = 5 - 10/0.5 = -15 + """ + true_params = {"xh": 5.0, "yh": 10.0, "x1": 0.0, "m2": 0.5} + x = np.linspace(0.1, 15, 200) + m1 = true_params["yh"] / (true_params["xh"] - true_params["x1"]) + x2 = true_params["xh"] - true_params["yh"] / true_params["m2"] + y = np.where( + x < true_params["xh"], + m1 * (x - true_params["x1"]), + true_params["m2"] * (x - x2), + ) + w = np.ones_like(x) + return x, y, w, true_params + + +@pytest.fixture +def noisy_saturation_data(): + """Hinge saturation data with 5% Gaussian noise. + + Parameters: xh=5.0, yh=10.0, x1=0.0, m2=0.0 + Noise std = 0.5 (5% of yh) + """ + rng = np.random.default_rng(42) + true_params = {"xh": 5.0, "yh": 10.0, "x1": 0.0, "m2": 0.0} + x = np.linspace(0.5, 15, 200) + m1 = true_params["yh"] / (true_params["xh"] - true_params["x1"]) + y_true = np.where( + x < true_params["xh"], + m1 * (x - true_params["x1"]), + true_params["yh"], + ) + noise_std = 0.5 + y = y_true + rng.normal(0, noise_std, len(x)) + w = np.ones_like(x) / noise_std + return x, y, w, true_params + + +@pytest.fixture +def offset_x1_data(): + """Hinge data with non-zero x1 (offset x-intercept). + + Parameters: xh=5.0, yh=10.0, x1=1.0, m2=0.0 + m1 = 10/(5-1) = 2.5, so rising region is f(x) = 2.5*(x-1) + """ + true_params = {"xh": 5.0, "yh": 10.0, "x1": 1.0, "m2": 0.0} + x = np.linspace(1.1, 15, 200) + m1 = true_params["yh"] / (true_params["xh"] - true_params["x1"]) + y = np.where( + x < true_params["xh"], + m1 * (x - true_params["x1"]), + true_params["yh"], + ) + w = np.ones_like(x) + return x, y, w, true_params + + +# ============================================================================= +# E1. Function Evaluation Tests (Exact Values) +# ============================================================================= + + +def test_func_evaluates_rising_region_correctly(): + """Before hinge: f(x) = m1*(x-x1) where m1 = yh/(xh-x1).""" + # Parameters: xh=5, yh=10, x1=0, m2=0 → m1 = 10/(5-0) = 2.0 + xh, yh, x1, m2 = 5.0, 10.0, 0.0, 0.0 + + # Test specific points in rising region (x < xh) + x_test = np.array([0.0, 1.0, 2.5, 4.0, 5.0]) # includes hinge point + # m1 = 2.0, so f(x) = 2*(x-0) = 2x + expected = np.array([0.0, 2.0, 5.0, 8.0, 10.0]) + + # Create minimal instance to access function + x_dummy = np.array([0.0, 10.0]) + y_dummy = np.array([0.0, 10.0]) + obj = HingeSaturation(x_dummy, y_dummy, guess_xh=xh, guess_yh=yh) + result = obj.function(x_test, xh, yh, x1, m2) + + np.testing.assert_allclose( + result, + expected, + rtol=1e-10, + err_msg="Rising region should follow f(x) = m1*(x-x1)", + ) + + +def test_func_evaluates_saturated_region_correctly(): + """After hinge with m2=0: f(x) = yh (constant plateau).""" + xh, yh, x1, m2 = 5.0, 10.0, 0.0, 0.0 + + # Test points beyond hinge (x > xh) + x_test = np.array([5.5, 6.0, 10.0, 100.0]) + expected = np.array([10.0, 10.0, 10.0, 10.0]) # saturated at yh + + x_dummy = np.array([0.0, 10.0]) + y_dummy = np.array([0.0, 10.0]) + obj = HingeSaturation(x_dummy, y_dummy, guess_xh=xh, guess_yh=yh) + result = obj.function(x_test, xh, yh, x1, m2) + + np.testing.assert_allclose( + result, + expected, + rtol=1e-10, + err_msg="Saturated region with m2=0 should be constant at yh", + ) + + +def test_func_evaluates_sloped_plateau_correctly(): + """After hinge with m2≠0: f(x) = m2*(x-x2) where x2 = xh - yh/m2.""" + xh, yh, x1, m2 = 5.0, 10.0, 0.0, 0.5 + # x2 = 5 - 10/0.5 = 5 - 20 = -15 + # After hinge: f(x) = 0.5*(x - (-15)) = 0.5*(x + 15) + + x_test = np.array([6.0, 8.0, 10.0]) + # f(6) = 0.5*(6+15) = 10.5 + # f(8) = 0.5*(8+15) = 11.5 + # f(10) = 0.5*(10+15) = 12.5 + expected = np.array([10.5, 11.5, 12.5]) + + x_dummy = np.array([0.0, 10.0]) + y_dummy = np.array([0.0, 10.0]) + obj = HingeSaturation(x_dummy, y_dummy, guess_xh=xh, guess_yh=yh) + result = obj.function(x_test, xh, yh, x1, m2) + + np.testing.assert_allclose( + result, + expected, + rtol=1e-10, + err_msg="Sloped plateau should follow f(x) = m2*(x-x2)", + ) + + +# ============================================================================= +# E2. Parameter Recovery Tests (Clean Data) +# ============================================================================= + + +def test_fit_recovers_exact_parameters_from_clean_data(clean_saturation_data): + """Fitting noise-free data should recover parameters within 1%.""" + x, y, w, true_params = clean_saturation_data + + obj = HingeSaturation(x, y, guess_xh=true_params["xh"], guess_yh=true_params["yh"]) + obj.make_fit() + + for param, true_val in true_params.items(): + fitted_val = obj.popt[param] + # Use absolute tolerance for values near zero + if abs(true_val) < 0.1: + assert abs(fitted_val - true_val) < 0.05, ( + f"{param}: fitted={fitted_val:.4f}, true={true_val:.4f}, " + f"absolute error exceeds 0.05" + ) + else: + rel_error = abs(fitted_val - true_val) / abs(true_val) + assert rel_error < 0.01, ( + f"{param}: fitted={fitted_val:.4f}, true={true_val:.4f}, " + f"rel_error={rel_error:.2%} exceeds 1% tolerance" + ) + + +def test_fit_recovers_sloped_plateau_parameters(clean_sloped_plateau_data): + """Fitting clean data with m2≠0 should recover parameters within 2%.""" + x, y, w, true_params = clean_sloped_plateau_data + + obj = HingeSaturation(x, y, guess_xh=true_params["xh"], guess_yh=true_params["yh"]) + obj.make_fit() + + for param, true_val in true_params.items(): + fitted_val = obj.popt[param] + if abs(true_val) < 0.1: + assert ( + abs(fitted_val - true_val) < 0.05 + ), f"{param}: fitted={fitted_val:.4f}, true={true_val:.4f}" + else: + rel_error = abs(fitted_val - true_val) / abs(true_val) + assert rel_error < 0.02, ( + f"{param}: fitted={fitted_val:.4f}, true={true_val:.4f}, " + f"rel_error={rel_error:.2%}" + ) + + +def test_fit_recovers_offset_x1_parameters(offset_x1_data): + """Fitting clean data with x1≠0 should recover parameters within 2%.""" + x, y, w, true_params = offset_x1_data + + obj = HingeSaturation(x, y, guess_xh=true_params["xh"], guess_yh=true_params["yh"]) + obj.make_fit() + + for param, true_val in true_params.items(): + fitted_val = obj.popt[param] + if abs(true_val) < 0.1: + assert ( + abs(fitted_val - true_val) < 0.05 + ), f"{param}: fitted={fitted_val:.4f}, true={true_val:.4f}" + else: + rel_error = abs(fitted_val - true_val) / abs(true_val) + assert rel_error < 0.02, ( + f"{param}: fitted={fitted_val:.4f}, true={true_val:.4f}, " + f"rel_error={rel_error:.2%}" + ) + + +@pytest.mark.parametrize( + "true_params", + [ + {"xh": 5.0, "yh": 10.0, "x1": 0.0, "m2": 0.0}, # Classic saturation + {"xh": 3.0, "yh": 6.0, "x1": 1.0, "m2": 0.0}, # Offset x-intercept + {"xh": 8.0, "yh": 4.0, "x1": 2.0, "m2": 0.3}, # Sloped plateau + {"xh": 5.0, "yh": 10.0, "x1": 0.0, "m2": -0.2}, # Declining plateau + ], +) +def test_fit_recovers_various_parameter_combinations(true_params): + """Fitting should work for diverse parameter combinations.""" + x_start = true_params["x1"] + 0.1 + x_end = true_params["xh"] + 10 + x = np.linspace(x_start, x_end, 200) + + # Build y from parameters + m1 = true_params["yh"] / (true_params["xh"] - true_params["x1"]) + if abs(true_params["m2"]) > 1e-10: + x2 = true_params["xh"] - true_params["yh"] / true_params["m2"] + y = np.where( + x < true_params["xh"], + m1 * (x - true_params["x1"]), + true_params["m2"] * (x - x2), + ) + else: + y = np.where( + x < true_params["xh"], + m1 * (x - true_params["x1"]), + true_params["yh"], + ) + + obj = HingeSaturation(x, y, guess_xh=true_params["xh"], guess_yh=true_params["yh"]) + obj.make_fit() + + for param, true_val in true_params.items(): + fitted_val = obj.popt[param] + if abs(true_val) < 0.1: + assert ( + abs(fitted_val - true_val) < 0.05 + ), f"{param}: fitted={fitted_val:.4f}, true={true_val:.4f}" + else: + rel_error = abs(fitted_val - true_val) / abs(true_val) + assert rel_error < 0.02, ( + f"{param}: fitted={fitted_val:.4f}, true={true_val:.4f}, " + f"rel_error={rel_error:.2%}" + ) + + +# ============================================================================= +# E3. Noisy Data Tests (Precision Bounds) +# ============================================================================= + + +def test_fit_with_noise_recovers_parameters_within_2sigma(noisy_saturation_data): + """Fitted parameters should be within 2σ of true values (95% confidence). + + With 4 parameters, testing each at 1σ (68%) gives only (0.68)^4 ≈ 21% joint + probability of all passing. Using 2σ (95%) gives (0.95)^4 ≈ 81% joint + probability, which is robust for automated testing. + + For well-behaved Gaussian noise, we expect deviations < 2σ with high confidence. + """ + x, y, w, true_params = noisy_saturation_data + + obj = HingeSaturation(x, y, guess_xh=true_params["xh"], guess_yh=true_params["yh"]) + obj.make_fit() + + for param, true_val in true_params.items(): + fitted_val = obj.popt[param] + sigma = obj.psigma[param] + deviation = abs(fitted_val - true_val) + + # 2σ gives 95% confidence per parameter, ~81% joint for 4 parameters + assert deviation < 2 * sigma, ( + f"{param}: |fitted({fitted_val:.4f}) - true({true_val:.4f})| = " + f"{deviation:.4f} exceeds 2σ = {2*sigma:.4f}" + ) + + +def test_fit_uncertainty_scales_with_noise(): + """Higher noise should produce larger parameter uncertainties.""" + rng = np.random.default_rng(42) + true_params = {"xh": 5.0, "yh": 10.0, "x1": 0.0, "m2": 0.0} + x = np.linspace(0.5, 15, 200) + m1 = true_params["yh"] / (true_params["xh"] - true_params["x1"]) + y_true = np.where( + x < true_params["xh"], + m1 * (x - true_params["x1"]), + true_params["yh"], + ) + + # Low noise fit + y_low = y_true + rng.normal(0, 0.2, len(x)) + fit_low = HingeSaturation( + x, y_low, guess_xh=true_params["xh"], guess_yh=true_params["yh"] + ) + fit_low.make_fit() + + # High noise fit (different seed for independence) + rng2 = np.random.default_rng(43) + y_high = y_true + rng2.normal(0, 1.0, len(x)) + fit_high = HingeSaturation( + x, y_high, guess_xh=true_params["xh"], guess_yh=true_params["yh"] + ) + fit_high.make_fit() + + # High noise should have larger uncertainties for most parameters + # (xh and yh are the primary parameters affected by noise) + for param in ["xh", "yh"]: + assert fit_high.psigma[param] > fit_low.psigma[param], ( + f"{param}: high_noise_sigma={fit_high.psigma[param]:.4f} should be " + f"> low_noise_sigma={fit_low.psigma[param]:.4f}" + ) + + +# ============================================================================= +# E4. Initial Parameter Estimation Tests +# ============================================================================= + + +def test_p0_returns_list_with_correct_length(clean_saturation_data): + """p0 should return a list with 4 elements.""" + x, y, w, true_params = clean_saturation_data + obj = HingeSaturation(x, y, guess_xh=true_params["xh"], guess_yh=true_params["yh"]) + + p0 = obj.p0 + assert isinstance(p0, list), f"p0 should be a list, got {type(p0)}" + assert len(p0) == 4, f"p0 should have 4 elements (xh, yh, x1, m2), got {len(p0)}" + + +def test_p0_enables_successful_convergence(noisy_saturation_data): + """Fit should converge when initialized with estimated p0.""" + x, y, w, true_params = noisy_saturation_data + + obj = HingeSaturation(x, y, guess_xh=true_params["xh"], guess_yh=true_params["yh"]) + obj.make_fit() + + # Fit should have converged (popt should exist and be finite) + assert all( + np.isfinite(v) for v in obj.popt.values() + ), f"Fit did not converge: popt={obj.popt}" + + +# ============================================================================= +# E5. Derived Quantity Tests (Internal Consistency) +# ============================================================================= + + +def test_fitted_m1_is_consistent_with_xh_yh_x1(offset_x1_data): + """Verify m1 = yh / (xh - x1) relationship holds for fitted params.""" + x, y, w, true_params = offset_x1_data + # True m1 = 10 / (5 - 1) = 2.5 + + obj = HingeSaturation(x, y, guess_xh=true_params["xh"], guess_yh=true_params["yh"]) + obj.make_fit() + + xh, yh, x1 = obj.popt["xh"], obj.popt["yh"], obj.popt["x1"] + m1_from_params = yh / (xh - x1) + + # Verify by evaluating function in rising region + x_rising = np.array([2.0, 3.0]) # Well before hinge + y_rising = obj(x_rising) + + # y(x) = m1 * (x - x1) → m1 = y / (x - x1) + m1_from_values = y_rising / (x_rising - x1) + + np.testing.assert_allclose( + m1_from_values, + m1_from_params, + rtol=1e-6, + err_msg="m1 derived from function values should match m1 from parameters", + ) + + +def test_hinge_point_is_continuous(clean_sloped_plateau_data): + """Function value at xh should equal yh (continuity at hinge).""" + x, y, w, true_params = clean_sloped_plateau_data + + obj = HingeSaturation(x, y, guess_xh=true_params["xh"], guess_yh=true_params["yh"]) + obj.make_fit() + + # Evaluate exactly at fitted xh + xh = obj.popt["xh"] + yh_expected = obj.popt["yh"] + y_at_hinge = obj(np.array([xh]))[0] + + np.testing.assert_allclose( + y_at_hinge, + yh_expected, + rtol=1e-6, + err_msg=f"f(xh={xh:.3f}) = {y_at_hinge:.3f} should equal yh={yh_expected:.3f}", + ) + + +# ============================================================================= +# E6. Edge Case and Error Handling Tests +# ============================================================================= + + +def test_weighted_fit_respects_weights(): + """Weighted fitting should correctly use sigma to weight observations. + + In FitFunction, weights are interpreted as uncertainties (sigma). Points + with larger sigma have MORE uncertainty and thus LESS influence on the fit. + + Test strategy: Apply non-uniform sigma values and verify the fit converges + correctly. HingeSaturation's min() function makes it inherently robust to + plateau outliers, so we test that weighting works by verifying accurate + parameter recovery with realistic sigma values. + """ + rng = np.random.default_rng(42) + true_params = {"xh": 5.0, "yh": 10.0, "x1": 0.0, "m2": 0.0} + + x = np.linspace(0.5, 15, 100) + m1 = true_params["yh"] / (true_params["xh"] - true_params["x1"]) + y_true = np.where( + x < true_params["xh"], + m1 * (x - true_params["x1"]), + true_params["yh"], + ) + + # Add heteroscedastic noise: larger noise in rising region, smaller in plateau + sigma_true = np.where(x < true_params["xh"], 0.5, 0.1) + noise = rng.normal(0, 1, len(x)) * sigma_true + y = y_true + noise + + # Fit with correct sigma values + fit_weighted = HingeSaturation( + x, y, weights=sigma_true, guess_xh=true_params["xh"], guess_yh=true_params["yh"] + ) + fit_weighted.make_fit() + + # Verify fit converged and parameters are accurate + assert all( + np.isfinite(v) for v in fit_weighted.popt.values() + ), f"Weighted fit did not converge: popt={fit_weighted.popt}" + + # With proper weighting, should recover true parameters within 5% + for param, true_val in true_params.items(): + fitted_val = fit_weighted.popt[param] + if abs(true_val) < 0.1: + # Absolute tolerance for values near zero + assert ( + abs(fitted_val - true_val) < 0.1 + ), f"{param}: fitted={fitted_val:.4f}, true={true_val:.4f}" + else: + rel_error = abs(fitted_val - true_val) / abs(true_val) + assert rel_error < 0.05, ( + f"{param}: fitted={fitted_val:.4f}, true={true_val:.4f}, " + f"rel_error={rel_error:.2%} exceeds 5%" + ) + + +def test_insufficient_data_raises_error(): + """Fitting with insufficient data should raise InsufficientDataError.""" + x = np.array([1.0, 2.0]) # Only 2 points for 4 parameters + y = np.array([2.0, 4.0]) + + obj = HingeSaturation(x, y, guess_xh=5.0, guess_yh=10.0) + + with pytest.raises(InsufficientDataError): + obj.make_fit() + + +def test_function_signature(): + """Test that function has correct parameter signature.""" + x = np.array([0.0, 5.0, 10.0]) + y = np.array([0.0, 10.0, 10.0]) + obj = HingeSaturation(x, y, guess_xh=5.0, guess_yh=10.0) + + sig = inspect.signature(obj.function) + params = tuple(sig.parameters.keys()) + + assert params == ( + "x", + "xh", + "yh", + "x1", + "m2", + ), f"Function should have signature (x, xh, yh, x1, m2), got {params}" + + +def test_tex_function_property(): + """Test that TeX_function returns expected LaTeX string.""" + x = np.array([0.0, 5.0, 10.0]) + y = np.array([0.0, 10.0, 10.0]) + obj = HingeSaturation(x, y, guess_xh=5.0, guess_yh=10.0) + + tex = obj.TeX_function + assert isinstance(tex, str), f"TeX_function should return str, got {type(tex)}" + assert r"\min" in tex, "TeX_function should contain \\min" + assert "y_1" in tex or "y_i" in tex, "TeX_function should reference y components" + + +def test_callable_interface(clean_saturation_data): + """Test that fitted object is callable and returns correct shape.""" + x, y, w, true_params = clean_saturation_data + + obj = HingeSaturation(x, y, guess_xh=true_params["xh"], guess_yh=true_params["yh"]) + obj.make_fit() + + # Test callable interface + x_test = np.array([1.0, 5.0, 10.0]) + y_pred = obj(x_test) + + assert ( + y_pred.shape == x_test.shape + ), f"Predicted shape {y_pred.shape} should match input shape {x_test.shape}" + assert np.all(np.isfinite(y_pred)), "All predicted values should be finite" + + +def test_popt_has_correct_keys(clean_saturation_data): + """Test that popt contains expected parameter names.""" + x, y, w, true_params = clean_saturation_data + + obj = HingeSaturation(x, y, guess_xh=true_params["xh"], guess_yh=true_params["yh"]) + obj.make_fit() + + expected_keys = {"xh", "yh", "x1", "m2"} + actual_keys = set(obj.popt.keys()) + + assert ( + actual_keys == expected_keys + ), f"popt keys should be {expected_keys}, got {actual_keys}" + + +def test_psigma_has_same_keys_as_popt(noisy_saturation_data): + """Test that psigma has same keys as popt.""" + x, y, w, true_params = noisy_saturation_data + + obj = HingeSaturation(x, y, guess_xh=true_params["xh"], guess_yh=true_params["yh"]) + obj.make_fit() + + assert set(obj.psigma.keys()) == set(obj.popt.keys()), ( + f"psigma keys {set(obj.psigma.keys())} should match " + f"popt keys {set(obj.popt.keys())}" + ) + + +def test_psigma_values_are_positive(noisy_saturation_data): + """Test that all parameter uncertainties are positive.""" + x, y, w, true_params = noisy_saturation_data + + obj = HingeSaturation(x, y, guess_xh=true_params["xh"], guess_yh=true_params["yh"]) + obj.make_fit() + + for param, sigma in obj.psigma.items(): + assert sigma > 0, f"psigma['{param}'] = {sigma} should be positive" + + +# ============================================================================= +# ============================================================================= +# +# TESTS FOR NEW HINGE FIT FUNCTION CLASSES +# +# The following sections test five new FitFunction subclasses: +# - TwoLine: Two intersecting lines with np.minimum +# - Saturation: Reparameterized TwoLine with xs, s, theta +# - HingeMin: Hinge with specified intersection point (minimum) +# - HingeMax: Hinge with specified intersection point (maximum) +# - HingeAtPoint: Hinge through a specified (xh, yh) point +# +# ============================================================================= +# ============================================================================= + + +from solarwindpy.fitfunctions.hinge import ( + TwoLine, + Saturation, + HingeMin, + HingeMax, + HingeAtPoint, +) + + +# ============================================================================= +# TwoLine Tests +# ============================================================================= +# TwoLine: f(x) = np.minimum(m1*(x-x1), m2*(x-x2)) +# Parameters: x1, x2, m1, m2 +# Derived: xs = (m1*x1 - m2*x2)/(m1 - m2), s = m1*(xs - x1), theta +# ============================================================================= + + +# ----------------------------------------------------------------------------- +# TwoLine Fixtures +# ----------------------------------------------------------------------------- + + +@pytest.fixture +def clean_twoline_data(): + """Perfect TwoLine data (no noise) with lines intersecting at (5, 10). + + Parameters: x1=0, x2=15, m1=2, m2=-1 + Line1: y = 2*(x - 0) = 2x, passes through (0, 0) and (5, 10) + Line2: y = -1*(x - 15) = -x + 15, passes through (15, 0) and (5, 10) + Intersection: xs = (2*0 - (-1)*15)/(2 - (-1)) = 15/3 = 5 + s = 2*(5 - 0) = 10 + """ + true_params = {"x1": 0.0, "x2": 15.0, "m1": 2.0, "m2": -1.0} + x = np.linspace(-2, 20, 300) + y1 = true_params["m1"] * (x - true_params["x1"]) + y2 = true_params["m2"] * (x - true_params["x2"]) + y = np.minimum(y1, y2) + w = np.ones_like(x) + return x, y, w, true_params + + +@pytest.fixture +def noisy_twoline_data(): + """TwoLine data with 5% Gaussian noise. + + Parameters: x1=0, x2=15, m1=2, m2=-1 + Noise std = 0.5 (5% of max value ~10) + """ + rng = np.random.default_rng(42) + true_params = {"x1": 0.0, "x2": 15.0, "m1": 2.0, "m2": -1.0} + x = np.linspace(-2, 20, 300) + y1 = true_params["m1"] * (x - true_params["x1"]) + y2 = true_params["m2"] * (x - true_params["x2"]) + y_true = np.minimum(y1, y2) + noise_std = 0.5 + y = y_true + rng.normal(0, noise_std, len(x)) + w = np.ones_like(x) / noise_std + return x, y, w, true_params + + +@pytest.fixture +def twoline_parallel_slopes_data(): + """TwoLine where slopes are same sign but different magnitude. + + Parameters: x1=0, x2=10, m1=3, m2=0.5 + Lines intersect where: 3*(x-0) = 0.5*(x-10) + 3x = 0.5x - 5 => 2.5x = -5 => x = -2 + y = 3*(-2) = -6 + """ + true_params = {"x1": 0.0, "x2": 10.0, "m1": 3.0, "m2": 0.5} + x = np.linspace(-5, 15, 300) + y1 = true_params["m1"] * (x - true_params["x1"]) + y2 = true_params["m2"] * (x - true_params["x2"]) + y = np.minimum(y1, y2) + w = np.ones_like(x) + return x, y, w, true_params + + +# ----------------------------------------------------------------------------- +# TwoLine E1. Function Evaluation Tests +# ----------------------------------------------------------------------------- + + +def test_twoline_func_evaluates_line1_region_correctly(): + """TwoLine should follow line1 where m1*(x-x1) < m2*(x-x2). + + With x1=0, x2=15, m1=2, m2=-1: + Line1: y = 2x, Line2: y = -x + 15 + Intersection at x=5. For x<5, line1 < line2. + """ + x1, x2, m1, m2 = 0.0, 15.0, 2.0, -1.0 + + x_test = np.array([0.0, 1.0, 2.0, 3.0, 4.0]) + # y = 2x for these points + expected = np.array([0.0, 2.0, 4.0, 6.0, 8.0]) + + x_dummy = np.linspace(0, 15, 50) + y_dummy = np.minimum(m1 * x_dummy, -m2 * (x_dummy - x2)) + obj = TwoLine(x_dummy, y_dummy, guess_xs=5.0) + result = obj.function(x_test, x1, x2, m1, m2) + + np.testing.assert_allclose( + result, + expected, + rtol=1e-10, + err_msg="TwoLine should follow line1 in left region", + ) + + +def test_twoline_func_evaluates_line2_region_correctly(): + """TwoLine should follow line2 where m2*(x-x2) < m1*(x-x1). + + With x1=0, x2=15, m1=2, m2=-1: + For x>5, line2 = -x + 15 < line1 = 2x. + """ + x1, x2, m1, m2 = 0.0, 15.0, 2.0, -1.0 + + x_test = np.array([6.0, 8.0, 10.0, 12.0, 15.0]) + # y = -x + 15 for these points + expected = np.array([9.0, 7.0, 5.0, 3.0, 0.0]) + + x_dummy = np.linspace(0, 15, 50) + y_dummy = np.minimum(m1 * x_dummy, -m2 * (x_dummy - x2)) + obj = TwoLine(x_dummy, y_dummy, guess_xs=5.0) + result = obj.function(x_test, x1, x2, m1, m2) + + np.testing.assert_allclose( + result, + expected, + rtol=1e-10, + err_msg="TwoLine should follow line2 in right region", + ) + + +def test_twoline_func_evaluates_intersection_correctly(): + """TwoLine should equal both lines at intersection point. + + Intersection at x=5: y = 2*5 = 10 = -5 + 15 = 10 + """ + x1, x2, m1, m2 = 0.0, 15.0, 2.0, -1.0 + xs = 5.0 # Intersection x + + x_test = np.array([xs]) + expected = np.array([10.0]) + + x_dummy = np.linspace(0, 15, 50) + y_dummy = np.minimum(m1 * x_dummy, -m2 * (x_dummy - x2)) + obj = TwoLine(x_dummy, y_dummy, guess_xs=xs) + result = obj.function(x_test, x1, x2, m1, m2) + + np.testing.assert_allclose( + result, + expected, + rtol=1e-10, + err_msg="TwoLine should be continuous at intersection", + ) + + +# ----------------------------------------------------------------------------- +# TwoLine E2. Parameter Recovery Tests (Clean Data) +# ----------------------------------------------------------------------------- + + +def test_twoline_fit_recovers_exact_parameters_from_clean_data(clean_twoline_data): + """Fitting noise-free TwoLine data should recover parameters within 2%.""" + x, y, w, true_params = clean_twoline_data + + obj = TwoLine(x, y, guess_xs=5.0) + obj.make_fit() + + for param, true_val in true_params.items(): + fitted_val = obj.popt[param] + if abs(true_val) < 0.1: + assert abs(fitted_val - true_val) < 0.05, ( + f"{param}: fitted={fitted_val:.4f}, true={true_val:.4f}, " + f"absolute error exceeds 0.05" + ) + else: + rel_error = abs(fitted_val - true_val) / abs(true_val) + assert rel_error < 0.02, ( + f"{param}: fitted={fitted_val:.4f}, true={true_val:.4f}, " + f"rel_error={rel_error:.2%} exceeds 2% tolerance" + ) + + +def test_twoline_fit_recovers_parallel_slopes_parameters(twoline_parallel_slopes_data): + """Fitting TwoLine with same-sign slopes should recover parameters within 2%.""" + x, y, w, true_params = twoline_parallel_slopes_data + + obj = TwoLine(x, y, guess_xs=-2.0) + obj.make_fit() + + for param, true_val in true_params.items(): + fitted_val = obj.popt[param] + if abs(true_val) < 0.1: + assert ( + abs(fitted_val - true_val) < 0.1 + ), f"{param}: fitted={fitted_val:.4f}, true={true_val:.4f}" + else: + rel_error = abs(fitted_val - true_val) / abs(true_val) + assert rel_error < 0.02, ( + f"{param}: fitted={fitted_val:.4f}, true={true_val:.4f}, " + f"rel_error={rel_error:.2%}" + ) + + +# ----------------------------------------------------------------------------- +# TwoLine E3. Noisy Data Tests +# ----------------------------------------------------------------------------- + + +def test_twoline_fit_with_noise_recovers_parameters_within_2sigma(noisy_twoline_data): + """Fitted TwoLine parameters should be within 2sigma of true values. + + With 4 parameters at 2sigma (95%), joint probability is (0.95)^4 = 81%. + """ + x, y, w, true_params = noisy_twoline_data + + obj = TwoLine(x, y, guess_xs=5.0) + obj.make_fit() + + for param, true_val in true_params.items(): + fitted_val = obj.popt[param] + sigma = obj.psigma[param] + deviation = abs(fitted_val - true_val) + + assert deviation < 2 * sigma, ( + f"{param}: |fitted({fitted_val:.4f}) - true({true_val:.4f})| = " + f"{deviation:.4f} exceeds 2sigma = {2*sigma:.4f}" + ) + + +# ----------------------------------------------------------------------------- +# TwoLine E4. Derived Property Tests +# ----------------------------------------------------------------------------- + + +def test_twoline_xs_property_is_consistent(clean_twoline_data): + """Verify xs = (m1*x1 - m2*x2)/(m1 - m2) for fitted params. + + True params: x1=0, x2=15, m1=2, m2=-1 + xs = (2*0 - (-1)*15)/(2 - (-1)) = 15/3 = 5 + """ + x, y, w, true_params = clean_twoline_data + + obj = TwoLine(x, y, guess_xs=5.0) + obj.make_fit() + + m1 = obj.popt["m1"] + m2 = obj.popt["m2"] + x1 = obj.popt["x1"] + x2 = obj.popt["x2"] + + xs_expected = (m1 * x1 - m2 * x2) / (m1 - m2) + + np.testing.assert_allclose( + obj.xs, + xs_expected, + rtol=1e-6, + err_msg="TwoLine.xs should equal (m1*x1 - m2*x2)/(m1 - m2)", + ) + + +def test_twoline_s_property_is_consistent(clean_twoline_data): + """Verify s = m1*(xs - x1) for fitted params. + + True params: xs=5, x1=0, m1=2 + s = 2*(5 - 0) = 10 + """ + x, y, w, true_params = clean_twoline_data + + obj = TwoLine(x, y, guess_xs=5.0) + obj.make_fit() + + m1 = obj.popt["m1"] + x1 = obj.popt["x1"] + xs = obj.xs + + s_expected = m1 * (xs - x1) + + np.testing.assert_allclose( + obj.s, + s_expected, + rtol=1e-6, + err_msg="TwoLine.s should equal m1*(xs - x1)", + ) + + +def test_twoline_theta_property_is_positive_for_converging_lines(clean_twoline_data): + """Theta (angle between lines) should be positive for converging lines. + + Lines y=2x and y=-x+15 form an angle. theta = arctan(m1) - arctan(m2). + theta = arctan(2) - arctan(-1) = 1.107 - (-0.785) = 1.892 rad ~ 108 deg + """ + x, y, w, true_params = clean_twoline_data + + obj = TwoLine(x, y, guess_xs=5.0) + obj.make_fit() + + m1 = true_params["m1"] + m2 = true_params["m2"] + theta_expected = np.arctan(m1) - np.arctan(m2) + + assert ( + obj.theta > 0 + ), f"theta={obj.theta:.4f} should be positive for converging lines" + np.testing.assert_allclose( + obj.theta, + theta_expected, + rtol=0.02, + err_msg=f"TwoLine.theta should be arctan(m1)-arctan(m2)={theta_expected:.4f}", + ) + + +# ----------------------------------------------------------------------------- +# TwoLine E5. Edge Cases and Interface Tests +# ----------------------------------------------------------------------------- + + +def test_twoline_function_signature(): + """Test that TwoLine function has correct parameter signature.""" + x = np.array([0.0, 5.0, 10.0]) + y = np.array([0.0, 10.0, 5.0]) + obj = TwoLine(x, y, guess_xs=5.0) + + sig = inspect.signature(obj.function) + params = tuple(sig.parameters.keys()) + + assert params == ( + "x", + "x1", + "x2", + "m1", + "m2", + ), f"Function should have signature (x, x1, x2, m1, m2), got {params}" + + +def test_twoline_popt_has_correct_keys(clean_twoline_data): + """Test that TwoLine popt contains expected parameter names.""" + x, y, w, true_params = clean_twoline_data + + obj = TwoLine(x, y, guess_xs=5.0) + obj.make_fit() + + expected_keys = {"x1", "x2", "m1", "m2"} + actual_keys = set(obj.popt.keys()) + + assert ( + actual_keys == expected_keys + ), f"popt keys should be {expected_keys}, got {actual_keys}" + + +def test_twoline_callable_interface(clean_twoline_data): + """Test that fitted TwoLine object is callable and returns correct shape.""" + x, y, w, true_params = clean_twoline_data + + obj = TwoLine(x, y, guess_xs=5.0) + obj.make_fit() + + x_test = np.array([1.0, 5.0, 10.0]) + y_pred = obj(x_test) + + assert ( + y_pred.shape == x_test.shape + ), f"Predicted shape {y_pred.shape} should match input shape {x_test.shape}" + assert np.all(np.isfinite(y_pred)), "All predicted values should be finite" + + +def test_twoline_insufficient_data_raises_error(): + """Fitting TwoLine with insufficient data should raise InsufficientDataError.""" + x = np.array([1.0, 2.0, 3.0]) # Only 3 points for 4 parameters + y = np.array([2.0, 4.0, 6.0]) + + obj = TwoLine(x, y, guess_xs=5.0) + + with pytest.raises(InsufficientDataError): + obj.make_fit() + + +# ============================================================================= +# Saturation Tests +# ============================================================================= +# Saturation: Reparameterized TwoLine with parameters (x1, xs, s, theta) +# function: np.minimum(l1, l2) where m1 = s/(xs-x1), m2 = tan(arctan(m1) - theta) +# Derived: m1, m2, x2 +# ============================================================================= + + +# ----------------------------------------------------------------------------- +# Saturation Fixtures +# ----------------------------------------------------------------------------- + + +@pytest.fixture +def clean_saturation_twoline_data(): + """Perfect Saturation data with known x1, xs, s, theta. + + Parameters: x1=0, xs=5, s=10, theta=pi/3 (60 degrees) + m1 = s/(xs-x1) = 10/(5-0) = 2 + m2 = tan(arctan(2) - pi/3) = tan(1.107 - 1.047) = tan(0.060) = 0.060 + Line1: y = 2*(x - 0) = 2x + Line2: y = m2*(x - x2) where x2 = xs - s/m2 + """ + true_params = {"x1": 0.0, "xs": 5.0, "s": 10.0, "theta": np.pi / 3} + m1 = true_params["s"] / (true_params["xs"] - true_params["x1"]) + m2 = np.tan(np.arctan(m1) - true_params["theta"]) + x2 = true_params["xs"] - true_params["s"] / m2 + + x = np.linspace(-2, 20, 300) + y1 = m1 * (x - true_params["x1"]) + y2 = m2 * (x - x2) + y = np.minimum(y1, y2) + w = np.ones_like(x) + return x, y, w, true_params, {"m1": m1, "m2": m2, "x2": x2} + + +@pytest.fixture +def noisy_saturation_twoline_data(): + """Saturation data with 5% Gaussian noise. + + Parameters: x1=0, xs=5, s=10, theta=pi/4 (45 degrees) + """ + rng = np.random.default_rng(42) + true_params = {"x1": 0.0, "xs": 5.0, "s": 10.0, "theta": np.pi / 4} + m1 = true_params["s"] / (true_params["xs"] - true_params["x1"]) + m2 = np.tan(np.arctan(m1) - true_params["theta"]) + x2 = true_params["xs"] - true_params["s"] / m2 + + x = np.linspace(-2, 20, 300) + y1 = m1 * (x - true_params["x1"]) + y2 = m2 * (x - x2) + y_true = np.minimum(y1, y2) + noise_std = 0.5 + y = y_true + rng.normal(0, noise_std, len(x)) + w = np.ones_like(x) / noise_std + return x, y, w, true_params + + +@pytest.fixture +def saturation_small_theta_data(): + """Saturation with small theta (nearly parallel lines after hinge). + + Parameters: x1=0, xs=5, s=10, theta=0.1 (about 5.7 degrees) + """ + true_params = {"x1": 0.0, "xs": 5.0, "s": 10.0, "theta": 0.1} + m1 = true_params["s"] / (true_params["xs"] - true_params["x1"]) + m2 = np.tan(np.arctan(m1) - true_params["theta"]) + x2 = true_params["xs"] - true_params["s"] / m2 + + x = np.linspace(-2, 25, 300) + y1 = m1 * (x - true_params["x1"]) + y2 = m2 * (x - x2) + y = np.minimum(y1, y2) + w = np.ones_like(x) + return x, y, w, true_params + + +# ----------------------------------------------------------------------------- +# Saturation E1. Function Evaluation Tests +# ----------------------------------------------------------------------------- + + +def test_saturation_func_evaluates_rising_region_correctly(): + """Saturation should follow line1 (rising) before saturation point. + + With x1=0, xs=5, s=10: m1 = 10/5 = 2 + Before xs=5: f(x) = 2*(x - 0) = 2x + """ + x1, xs, s, theta = 0.0, 5.0, 10.0, np.pi / 4 + m1 = s / (xs - x1) + + x_test = np.array([0.0, 1.0, 2.0, 3.0, 4.0]) + expected = m1 * (x_test - x1) # [0, 2, 4, 6, 8] + + x_dummy = np.linspace(0, 15, 50) + y_dummy = 2 * x_dummy + obj = Saturation(x_dummy, y_dummy, guess_xs=xs, guess_s=s) + result = obj.function(x_test, x1, xs, s, theta) + + np.testing.assert_allclose( + result, + expected, + rtol=1e-10, + err_msg="Saturation should follow rising line before xs", + ) + + +def test_saturation_func_passes_through_saturation_point(): + """Saturation function should pass through (xs, s). + + At x=xs: both lines should equal s. + """ + x1, xs, s, theta = 0.0, 5.0, 10.0, np.pi / 4 + + x_test = np.array([xs]) + expected = np.array([s]) + + x_dummy = np.linspace(0, 15, 50) + y_dummy = 2 * x_dummy + obj = Saturation(x_dummy, y_dummy, guess_xs=xs, guess_s=s) + result = obj.function(x_test, x1, xs, s, theta) + + np.testing.assert_allclose( + result, + expected, + rtol=1e-10, + err_msg="Saturation function should pass through (xs, s)", + ) + + +def test_saturation_func_theta_controls_plateau_slope(): + """Different theta values should produce different plateau slopes. + + theta=0: m2 = m1 (plateau continues at same slope) + theta>0: m2 < m1 (plateau is less steep) + """ + x1, xs, s = 0.0, 5.0, 10.0 + m1 = s / (xs - x1) # = 2 + + theta_small = 0.1 + theta_large = np.pi / 3 + + m2_small = np.tan(np.arctan(m1) - theta_small) + m2_large = np.tan(np.arctan(m1) - theta_large) + + # Larger theta should give smaller m2 (less steep plateau) + assert ( + m2_small > m2_large + ), f"m2_small={m2_small:.4f} should be > m2_large={m2_large:.4f}" + + +# ----------------------------------------------------------------------------- +# Saturation E2. Parameter Recovery Tests (Clean Data) +# ----------------------------------------------------------------------------- + + +def test_saturation_fit_recovers_exact_parameters_from_clean_data( + clean_saturation_twoline_data, +): + """Fitting noise-free Saturation data should recover parameters within 2%.""" + x, y, w, true_params, derived = clean_saturation_twoline_data + + obj = Saturation(x, y, guess_xs=true_params["xs"], guess_s=true_params["s"]) + obj.make_fit() + + for param, true_val in true_params.items(): + fitted_val = obj.popt[param] + if abs(true_val) < 0.1: + assert ( + abs(fitted_val - true_val) < 0.05 + ), f"{param}: fitted={fitted_val:.4f}, true={true_val:.4f}" + else: + rel_error = abs(fitted_val - true_val) / abs(true_val) + assert rel_error < 0.02, ( + f"{param}: fitted={fitted_val:.4f}, true={true_val:.4f}, " + f"rel_error={rel_error:.2%} exceeds 2% tolerance" + ) + + +def test_saturation_fit_recovers_small_theta_parameters(saturation_small_theta_data): + """Fitting Saturation with small theta should recover parameters within 5%.""" + x, y, w, true_params = saturation_small_theta_data + + obj = Saturation(x, y, guess_xs=true_params["xs"], guess_s=true_params["s"]) + obj.make_fit() + + # Small theta is harder to fit precisely, use 5% tolerance + for param, true_val in true_params.items(): + fitted_val = obj.popt[param] + if abs(true_val) < 0.1: + assert ( + abs(fitted_val - true_val) < 0.1 + ), f"{param}: fitted={fitted_val:.4f}, true={true_val:.4f}" + else: + rel_error = abs(fitted_val - true_val) / abs(true_val) + assert rel_error < 0.05, ( + f"{param}: fitted={fitted_val:.4f}, true={true_val:.4f}, " + f"rel_error={rel_error:.2%}" + ) + + +# ----------------------------------------------------------------------------- +# Saturation E3. Noisy Data Tests +# ----------------------------------------------------------------------------- + + +def test_saturation_fit_with_noise_recovers_parameters_within_2sigma( + noisy_saturation_twoline_data, +): + """Fitted Saturation parameters should be within 2sigma of true values.""" + x, y, w, true_params = noisy_saturation_twoline_data + + obj = Saturation(x, y, guess_xs=true_params["xs"], guess_s=true_params["s"]) + obj.make_fit() + + for param, true_val in true_params.items(): + fitted_val = obj.popt[param] + sigma = obj.psigma[param] + deviation = abs(fitted_val - true_val) + + assert deviation < 2 * sigma, ( + f"{param}: |fitted({fitted_val:.4f}) - true({true_val:.4f})| = " + f"{deviation:.4f} exceeds 2sigma = {2*sigma:.4f}" + ) + + +# ----------------------------------------------------------------------------- +# Saturation E4. Derived Property Tests +# ----------------------------------------------------------------------------- + + +def test_saturation_m1_property_is_consistent(clean_saturation_twoline_data): + """Verify m1 = s/(xs - x1) for fitted params.""" + x, y, w, true_params, derived = clean_saturation_twoline_data + + obj = Saturation(x, y, guess_xs=true_params["xs"], guess_s=true_params["s"]) + obj.make_fit() + + s = obj.popt["s"] + xs = obj.popt["xs"] + x1 = obj.popt["x1"] + + m1_expected = s / (xs - x1) + + np.testing.assert_allclose( + obj.m1, + m1_expected, + rtol=1e-6, + err_msg="Saturation.m1 should equal s/(xs - x1)", + ) + + +def test_saturation_m2_property_is_consistent(clean_saturation_twoline_data): + """Verify m2 = tan(arctan(m1) - theta) for fitted params.""" + x, y, w, true_params, derived = clean_saturation_twoline_data + + obj = Saturation(x, y, guess_xs=true_params["xs"], guess_s=true_params["s"]) + obj.make_fit() + + m1 = obj.m1 + theta = obj.popt["theta"] + + m2_expected = np.tan(np.arctan(m1) - theta) + + np.testing.assert_allclose( + obj.m2, + m2_expected, + rtol=1e-6, + err_msg="Saturation.m2 should equal tan(arctan(m1) - theta)", + ) + + +def test_saturation_x2_property_is_consistent(clean_saturation_twoline_data): + """Verify x2 = xs - s/m2 for fitted params.""" + x, y, w, true_params, derived = clean_saturation_twoline_data + + obj = Saturation(x, y, guess_xs=true_params["xs"], guess_s=true_params["s"]) + obj.make_fit() + + xs = obj.popt["xs"] + s = obj.popt["s"] + m2 = obj.m2 + + x2_expected = xs - s / m2 + + np.testing.assert_allclose( + obj.x2, + x2_expected, + rtol=1e-6, + err_msg="Saturation.x2 should equal xs - s/m2", + ) + + +# ----------------------------------------------------------------------------- +# Saturation E5. Edge Cases and Interface Tests +# ----------------------------------------------------------------------------- + + +def test_saturation_function_signature(): + """Test that Saturation function has correct parameter signature.""" + x = np.array([0.0, 5.0, 10.0]) + y = np.array([0.0, 10.0, 12.0]) + obj = Saturation(x, y, guess_xs=5.0, guess_s=10.0) + + sig = inspect.signature(obj.function) + params = tuple(sig.parameters.keys()) + + assert params == ( + "x", + "x1", + "xs", + "s", + "theta", + ), f"Function should have signature (x, x1, xs, s, theta), got {params}" + + +def test_saturation_popt_has_correct_keys(clean_saturation_twoline_data): + """Test that Saturation popt contains expected parameter names.""" + x, y, w, true_params, derived = clean_saturation_twoline_data + + obj = Saturation(x, y, guess_xs=true_params["xs"], guess_s=true_params["s"]) + obj.make_fit() + + expected_keys = {"x1", "xs", "s", "theta"} + actual_keys = set(obj.popt.keys()) + + assert ( + actual_keys == expected_keys + ), f"popt keys should be {expected_keys}, got {actual_keys}" + + +def test_saturation_callable_interface(clean_saturation_twoline_data): + """Test that fitted Saturation object is callable and returns correct shape.""" + x, y, w, true_params, derived = clean_saturation_twoline_data + + obj = Saturation(x, y, guess_xs=true_params["xs"], guess_s=true_params["s"]) + obj.make_fit() + + x_test = np.array([1.0, 5.0, 10.0]) + y_pred = obj(x_test) + + assert ( + y_pred.shape == x_test.shape + ), f"Predicted shape {y_pred.shape} should match input shape {x_test.shape}" + assert np.all(np.isfinite(y_pred)), "All predicted values should be finite" + + +def test_saturation_insufficient_data_raises_error(): + """Fitting Saturation with insufficient data should raise InsufficientDataError.""" + x = np.array([1.0, 2.0, 3.0]) + y = np.array([2.0, 4.0, 6.0]) + + obj = Saturation(x, y, guess_xs=5.0, guess_s=10.0) + + with pytest.raises(InsufficientDataError): + obj.make_fit() + + +# ============================================================================= +# HingeMin Tests +# ============================================================================= +# HingeMin: f(x) = np.minimum(l1, l2) where l1 = m1*(x-x1), l2 = m2*(x-x2) +# Parameters: m1, x1, x2, h (hinge x-coordinate) +# Constraint: both lines pass through (h, m1*(h-x1)) +# m2 = m1*(h-x1)/(h-x2) +# Derived: m2, theta +# ============================================================================= + + +# ----------------------------------------------------------------------------- +# HingeMin Fixtures +# ----------------------------------------------------------------------------- + + +@pytest.fixture +def clean_hingemin_data(): + """Perfect HingeMin data with lines meeting at hinge point. + + Parameters: m1=2, x1=0, x2=10, h=5 + At hinge h=5: yh = m1*(h-x1) = 2*(5-0) = 10 + m2 = m1*(h-x1)/(h-x2) = 2*(5-0)/(5-10) = 10/(-5) = -2 + Line1: y = 2*(x - 0) = 2x + Line2: y = -2*(x - 10) = -2x + 20 + Intersection: 2x = -2x + 20 => 4x = 20 => x = 5, y = 10 + """ + true_params = {"m1": 2.0, "x1": 0.0, "x2": 10.0, "h": 5.0} + yh = true_params["m1"] * (true_params["h"] - true_params["x1"]) + m2 = ( + true_params["m1"] + * (true_params["h"] - true_params["x1"]) + / (true_params["h"] - true_params["x2"]) + ) + + x = np.linspace(-2, 15, 300) + y1 = true_params["m1"] * (x - true_params["x1"]) + y2 = m2 * (x - true_params["x2"]) + y = np.minimum(y1, y2) + w = np.ones_like(x) + return x, y, w, true_params, {"m2": m2, "yh": yh} + + +@pytest.fixture +def noisy_hingemin_data(): + """HingeMin data with 5% Gaussian noise. + + Parameters: m1=2, x1=0, x2=10, h=5 + """ + rng = np.random.default_rng(42) + true_params = {"m1": 2.0, "x1": 0.0, "x2": 10.0, "h": 5.0} + m2 = ( + true_params["m1"] + * (true_params["h"] - true_params["x1"]) + / (true_params["h"] - true_params["x2"]) + ) + + x = np.linspace(-2, 15, 300) + y1 = true_params["m1"] * (x - true_params["x1"]) + y2 = m2 * (x - true_params["x2"]) + y_true = np.minimum(y1, y2) + noise_std = 0.5 + y = y_true + rng.normal(0, noise_std, len(x)) + w = np.ones_like(x) / noise_std + return x, y, w, true_params + + +@pytest.fixture +def hingemin_positive_slopes_data(): + """HingeMin where both slopes are positive but different. + + Parameters: m1=3, x1=0, x2=-5, h=5 + yh = 3*(5-0) = 15 + m2 = 3*(5-0)/(5-(-5)) = 15/10 = 1.5 + Line1: y = 3x + Line2: y = 1.5*(x + 5) = 1.5x + 7.5 + """ + true_params = {"m1": 3.0, "x1": 0.0, "x2": -5.0, "h": 5.0} + m2 = ( + true_params["m1"] + * (true_params["h"] - true_params["x1"]) + / (true_params["h"] - true_params["x2"]) + ) + + x = np.linspace(-3, 15, 300) + y1 = true_params["m1"] * (x - true_params["x1"]) + y2 = m2 * (x - true_params["x2"]) + y = np.minimum(y1, y2) + w = np.ones_like(x) + return x, y, w, true_params + + +# ----------------------------------------------------------------------------- +# HingeMin E1. Function Evaluation Tests +# ----------------------------------------------------------------------------- + + +def test_hingemin_func_evaluates_line1_region_correctly(): + """HingeMin should follow line1 where m1*(x-x1) < m2*(x-x2). + + With m1=2, x1=0, x2=10, h=5: m2=-2 + For x<5: 2x < -2x+20, so line1 dominates minimum. + """ + m1, x1, x2, h = 2.0, 0.0, 10.0, 5.0 + m2 = m1 * (h - x1) / (h - x2) # = -2 + + x_test = np.array([0.0, 1.0, 2.0, 3.0, 4.0]) + expected = m1 * (x_test - x1) # [0, 2, 4, 6, 8] + + x_dummy = np.linspace(0, 15, 50) + y_dummy = np.minimum(m1 * x_dummy, m2 * (x_dummy - x2)) + obj = HingeMin(x_dummy, y_dummy, guess_h=h) + result = obj.function(x_test, m1, x1, x2, h) + + np.testing.assert_allclose( + result, + expected, + rtol=1e-10, + err_msg="HingeMin should follow line1 in left region", + ) + + +def test_hingemin_func_evaluates_hinge_point_correctly(): + """HingeMin should pass through hinge point (h, yh). + + At h=5: yh = m1*(h-x1) = 2*5 = 10 + """ + m1, x1, x2, h = 2.0, 0.0, 10.0, 5.0 + + x_test = np.array([h]) + expected = np.array([m1 * (h - x1)]) # [10] + + x_dummy = np.linspace(0, 15, 50) + m2 = m1 * (h - x1) / (h - x2) + y_dummy = np.minimum(m1 * x_dummy, m2 * (x_dummy - x2)) + obj = HingeMin(x_dummy, y_dummy, guess_h=h) + result = obj.function(x_test, m1, x1, x2, h) + + np.testing.assert_allclose( + result, + expected, + rtol=1e-10, + err_msg="HingeMin should pass through hinge point", + ) + + +def test_hingemin_func_evaluates_line2_region_correctly(): + """HingeMin should follow line2 where m2*(x-x2) < m1*(x-x1). + + For x>5 with our parameters: -2x+20 < 2x + """ + m1, x1, x2, h = 2.0, 0.0, 10.0, 5.0 + m2 = m1 * (h - x1) / (h - x2) # = -2 + + x_test = np.array([6.0, 7.0, 8.0, 9.0, 10.0]) + expected = m2 * (x_test - x2) # [-2*(x-10)] = [8, 6, 4, 2, 0] + + x_dummy = np.linspace(0, 15, 50) + y_dummy = np.minimum(m1 * x_dummy, m2 * (x_dummy - x2)) + obj = HingeMin(x_dummy, y_dummy, guess_h=h) + result = obj.function(x_test, m1, x1, x2, h) + + np.testing.assert_allclose( + result, + expected, + rtol=1e-10, + err_msg="HingeMin should follow line2 in right region", + ) + + +# ----------------------------------------------------------------------------- +# HingeMin E2. Parameter Recovery Tests (Clean Data) +# ----------------------------------------------------------------------------- + + +def test_hingemin_fit_recovers_exact_parameters_from_clean_data(clean_hingemin_data): + """Fitting noise-free HingeMin data should recover parameters within 2%.""" + x, y, w, true_params, derived = clean_hingemin_data + + obj = HingeMin(x, y, guess_h=true_params["h"]) + obj.make_fit() + + for param, true_val in true_params.items(): + fitted_val = obj.popt[param] + if abs(true_val) < 0.1: + assert ( + abs(fitted_val - true_val) < 0.05 + ), f"{param}: fitted={fitted_val:.4f}, true={true_val:.4f}" + else: + rel_error = abs(fitted_val - true_val) / abs(true_val) + assert rel_error < 0.02, ( + f"{param}: fitted={fitted_val:.4f}, true={true_val:.4f}, " + f"rel_error={rel_error:.2%} exceeds 2% tolerance" + ) + + +def test_hingemin_fit_recovers_positive_slopes_parameters( + hingemin_positive_slopes_data, +): + """Fitting HingeMin with positive slopes should recover parameters within 2%.""" + x, y, w, true_params = hingemin_positive_slopes_data + + obj = HingeMin(x, y, guess_h=true_params["h"]) + obj.make_fit() + + for param, true_val in true_params.items(): + fitted_val = obj.popt[param] + if abs(true_val) < 0.1: + assert ( + abs(fitted_val - true_val) < 0.1 + ), f"{param}: fitted={fitted_val:.4f}, true={true_val:.4f}" + else: + rel_error = abs(fitted_val - true_val) / abs(true_val) + assert rel_error < 0.02, ( + f"{param}: fitted={fitted_val:.4f}, true={true_val:.4f}, " + f"rel_error={rel_error:.2%}" + ) + + +# ----------------------------------------------------------------------------- +# HingeMin E3. Noisy Data Tests +# ----------------------------------------------------------------------------- + + +def test_hingemin_fit_with_noise_recovers_parameters_within_2sigma(noisy_hingemin_data): + """Fitted HingeMin parameters should be within 2sigma of true values.""" + x, y, w, true_params = noisy_hingemin_data + + obj = HingeMin(x, y, guess_h=true_params["h"]) + obj.make_fit() + + for param, true_val in true_params.items(): + fitted_val = obj.popt[param] + sigma = obj.psigma[param] + deviation = abs(fitted_val - true_val) + + assert deviation < 2 * sigma, ( + f"{param}: |fitted({fitted_val:.4f}) - true({true_val:.4f})| = " + f"{deviation:.4f} exceeds 2sigma = {2*sigma:.4f}" + ) + + +# ----------------------------------------------------------------------------- +# HingeMin E4. Derived Property Tests +# ----------------------------------------------------------------------------- + + +def test_hingemin_m2_property_is_consistent(clean_hingemin_data): + """Verify m2 = m1*(h-x1)/(h-x2) for fitted params. + + True: m2 = 2*(5-0)/(5-10) = 10/(-5) = -2 + """ + x, y, w, true_params, derived = clean_hingemin_data + + obj = HingeMin(x, y, guess_h=true_params["h"]) + obj.make_fit() + + m1 = obj.popt["m1"] + x1 = obj.popt["x1"] + x2 = obj.popt["x2"] + h = obj.popt["h"] + + m2_expected = m1 * (h - x1) / (h - x2) + + np.testing.assert_allclose( + obj.m2, + m2_expected, + rtol=1e-6, + err_msg="HingeMin.m2 should equal m1*(h-x1)/(h-x2)", + ) + + +def test_hingemin_theta_property_is_consistent(clean_hingemin_data): + """Verify theta = arctan(m1) - arctan(m2) for fitted params.""" + x, y, w, true_params, derived = clean_hingemin_data + + obj = HingeMin(x, y, guess_h=true_params["h"]) + obj.make_fit() + + m1 = obj.popt["m1"] + m2 = obj.m2 + theta_expected = np.arctan(m1) - np.arctan(m2) + + np.testing.assert_allclose( + obj.theta, + theta_expected, + rtol=1e-6, + err_msg="HingeMin.theta should equal arctan(m1) - arctan(m2)", + ) + + +def test_hingemin_lines_intersect_at_hinge(clean_hingemin_data): + """Verify both lines pass through (h, yh) where yh = m1*(h-x1).""" + x, y, w, true_params, derived = clean_hingemin_data + + obj = HingeMin(x, y, guess_h=true_params["h"]) + obj.make_fit() + + h = obj.popt["h"] + m1 = obj.popt["m1"] + x1 = obj.popt["x1"] + x2 = obj.popt["x2"] + m2 = obj.m2 + + yh_from_line1 = m1 * (h - x1) + yh_from_line2 = m2 * (h - x2) + + np.testing.assert_allclose( + yh_from_line1, + yh_from_line2, + rtol=1e-6, + err_msg="Both lines should pass through hinge point", + ) + + +# ----------------------------------------------------------------------------- +# HingeMin E5. Edge Cases and Interface Tests +# ----------------------------------------------------------------------------- + + +def test_hingemin_function_signature(): + """Test that HingeMin function has correct parameter signature.""" + x = np.array([0.0, 5.0, 10.0]) + y = np.array([0.0, 10.0, 0.0]) + obj = HingeMin(x, y, guess_h=5.0) + + sig = inspect.signature(obj.function) + params = tuple(sig.parameters.keys()) + + assert params == ( + "x", + "m1", + "x1", + "x2", + "h", + ), f"Function should have signature (x, m1, x1, x2, h), got {params}" + + +def test_hingemin_popt_has_correct_keys(clean_hingemin_data): + """Test that HingeMin popt contains expected parameter names.""" + x, y, w, true_params, derived = clean_hingemin_data + + obj = HingeMin(x, y, guess_h=true_params["h"]) + obj.make_fit() + + expected_keys = {"m1", "x1", "x2", "h"} + actual_keys = set(obj.popt.keys()) + + assert ( + actual_keys == expected_keys + ), f"popt keys should be {expected_keys}, got {actual_keys}" + + +def test_hingemin_callable_interface(clean_hingemin_data): + """Test that fitted HingeMin object is callable and returns correct shape.""" + x, y, w, true_params, derived = clean_hingemin_data + + obj = HingeMin(x, y, guess_h=true_params["h"]) + obj.make_fit() + + x_test = np.array([1.0, 5.0, 9.0]) + y_pred = obj(x_test) + + assert ( + y_pred.shape == x_test.shape + ), f"Predicted shape {y_pred.shape} should match input shape {x_test.shape}" + assert np.all(np.isfinite(y_pred)), "All predicted values should be finite" + + +def test_hingemin_insufficient_data_raises_error(): + """Fitting HingeMin with insufficient data should raise InsufficientDataError.""" + x = np.array([1.0, 2.0, 3.0]) + y = np.array([2.0, 4.0, 3.0]) + + obj = HingeMin(x, y, guess_h=2.0) + + with pytest.raises(InsufficientDataError): + obj.make_fit() + + +# ============================================================================= +# HingeMax Tests +# ============================================================================= +# HingeMax: f(x) = np.maximum(l1, l2) where l1 = m1*(x-x1), l2 = m2*(x-x2) +# Parameters: m1, x1, x2, h (hinge x-coordinate) +# Constraint: both lines pass through (h, m1*(h-x1)) +# m2 = m1*(h-x1)/(h-x2) +# Derived: m2, theta +# ============================================================================= + + +# ----------------------------------------------------------------------------- +# HingeMax Fixtures +# ----------------------------------------------------------------------------- + + +@pytest.fixture +def clean_hingemax_data(): + """Perfect HingeMax data with lines meeting at hinge point. + + Parameters: m1=-2, x1=0, x2=10, h=5 + At hinge h=5: yh = m1*(h-x1) = -2*(5-0) = -10 + m2 = m1*(h-x1)/(h-x2) = -2*(5-0)/(5-10) = -10/(-5) = 2 + Line1: y = -2*(x - 0) = -2x + Line2: y = 2*(x - 10) = 2x - 20 + max(-2x, 2x-20) forms a V-shape opening upward with vertex at (5, -10) + """ + true_params = {"m1": -2.0, "x1": 0.0, "x2": 10.0, "h": 5.0} + yh = true_params["m1"] * (true_params["h"] - true_params["x1"]) + m2 = ( + true_params["m1"] + * (true_params["h"] - true_params["x1"]) + / (true_params["h"] - true_params["x2"]) + ) + + x = np.linspace(-2, 15, 300) + y1 = true_params["m1"] * (x - true_params["x1"]) + y2 = m2 * (x - true_params["x2"]) + y = np.maximum(y1, y2) + w = np.ones_like(x) + return x, y, w, true_params, {"m2": m2, "yh": yh} + + +@pytest.fixture +def noisy_hingemax_data(): + """HingeMax data with 5% Gaussian noise. + + Parameters: m1=-2, x1=0, x2=10, h=5 + """ + rng = np.random.default_rng(42) + true_params = {"m1": -2.0, "x1": 0.0, "x2": 10.0, "h": 5.0} + m2 = ( + true_params["m1"] + * (true_params["h"] - true_params["x1"]) + / (true_params["h"] - true_params["x2"]) + ) + + x = np.linspace(-2, 15, 300) + y1 = true_params["m1"] * (x - true_params["x1"]) + y2 = m2 * (x - true_params["x2"]) + y_true = np.maximum(y1, y2) + noise_std = 0.5 + y = y_true + rng.normal(0, noise_std, len(x)) + w = np.ones_like(x) / noise_std + return x, y, w, true_params + + +@pytest.fixture +def hingemax_negative_slopes_data(): + """HingeMax where both slopes are negative. + + Parameters: m1=-3, x1=0, x2=-5, h=5 + yh = -3*(5-0) = -15 + m2 = -3*(5-0)/(5-(-5)) = -15/10 = -1.5 + Line1: y = -3x + Line2: y = -1.5*(x + 5) = -1.5x - 7.5 + """ + true_params = {"m1": -3.0, "x1": 0.0, "x2": -5.0, "h": 5.0} + m2 = ( + true_params["m1"] + * (true_params["h"] - true_params["x1"]) + / (true_params["h"] - true_params["x2"]) + ) + + x = np.linspace(-3, 15, 300) + y1 = true_params["m1"] * (x - true_params["x1"]) + y2 = m2 * (x - true_params["x2"]) + y = np.maximum(y1, y2) + w = np.ones_like(x) + return x, y, w, true_params + + +# ----------------------------------------------------------------------------- +# HingeMax E1. Function Evaluation Tests +# ----------------------------------------------------------------------------- + + +def test_hingemax_func_evaluates_line1_region_correctly(): + """HingeMax should follow line1 where m1*(x-x1) > m2*(x-x2). + + With m1=-2, x1=0, x2=10, h=5: m2=2 + For x<5: -2x > 2x-20, so line1 dominates maximum. + """ + m1, x1, x2, h = -2.0, 0.0, 10.0, 5.0 + m2 = m1 * (h - x1) / (h - x2) # = 2 + + x_test = np.array([0.0, 1.0, 2.0, 3.0, 4.0]) + expected = m1 * (x_test - x1) # [0, -2, -4, -6, -8] + + x_dummy = np.linspace(0, 15, 50) + y_dummy = np.maximum(m1 * x_dummy, m2 * (x_dummy - x2)) + obj = HingeMax(x_dummy, y_dummy, guess_h=h) + result = obj.function(x_test, m1, x1, x2, h) + + np.testing.assert_allclose( + result, + expected, + rtol=1e-10, + err_msg="HingeMax should follow line1 in left region", + ) + + +def test_hingemax_func_evaluates_hinge_point_correctly(): + """HingeMax should pass through hinge point (h, yh). + + At h=5: yh = m1*(h-x1) = -2*5 = -10 + """ + m1, x1, x2, h = -2.0, 0.0, 10.0, 5.0 + + x_test = np.array([h]) + expected = np.array([m1 * (h - x1)]) # [-10] + + x_dummy = np.linspace(0, 15, 50) + m2 = m1 * (h - x1) / (h - x2) + y_dummy = np.maximum(m1 * x_dummy, m2 * (x_dummy - x2)) + obj = HingeMax(x_dummy, y_dummy, guess_h=h) + result = obj.function(x_test, m1, x1, x2, h) + + np.testing.assert_allclose( + result, + expected, + rtol=1e-10, + err_msg="HingeMax should pass through hinge point", + ) + + +def test_hingemax_func_evaluates_line2_region_correctly(): + """HingeMax should follow line2 where m2*(x-x2) > m1*(x-x1). + + For x>5 with our parameters: 2x-20 > -2x + """ + m1, x1, x2, h = -2.0, 0.0, 10.0, 5.0 + m2 = m1 * (h - x1) / (h - x2) # = 2 + + x_test = np.array([6.0, 7.0, 8.0, 9.0, 10.0]) + expected = m2 * (x_test - x2) # [2*(x-10)] = [-8, -6, -4, -2, 0] + + x_dummy = np.linspace(0, 15, 50) + y_dummy = np.maximum(m1 * x_dummy, m2 * (x_dummy - x2)) + obj = HingeMax(x_dummy, y_dummy, guess_h=h) + result = obj.function(x_test, m1, x1, x2, h) + + np.testing.assert_allclose( + result, + expected, + rtol=1e-10, + err_msg="HingeMax should follow line2 in right region", + ) + + +# ----------------------------------------------------------------------------- +# HingeMax E2. Parameter Recovery Tests (Clean Data) +# ----------------------------------------------------------------------------- + + +def test_hingemax_fit_recovers_exact_parameters_from_clean_data(clean_hingemax_data): + """Fitting noise-free HingeMax data should recover parameters within 2%.""" + x, y, w, true_params, derived = clean_hingemax_data + + obj = HingeMax(x, y, guess_h=true_params["h"]) + obj.make_fit() + + for param, true_val in true_params.items(): + fitted_val = obj.popt[param] + if abs(true_val) < 0.1: + assert ( + abs(fitted_val - true_val) < 0.05 + ), f"{param}: fitted={fitted_val:.4f}, true={true_val:.4f}" + else: + rel_error = abs(fitted_val - true_val) / abs(true_val) + assert rel_error < 0.02, ( + f"{param}: fitted={fitted_val:.4f}, true={true_val:.4f}, " + f"rel_error={rel_error:.2%} exceeds 2% tolerance" + ) + + +def test_hingemax_fit_recovers_negative_slopes_parameters( + hingemax_negative_slopes_data, +): + """Fitting HingeMax with negative slopes should recover parameters within 2%.""" + x, y, w, true_params = hingemax_negative_slopes_data + + obj = HingeMax(x, y, guess_h=true_params["h"]) + obj.make_fit() + + for param, true_val in true_params.items(): + fitted_val = obj.popt[param] + if abs(true_val) < 0.1: + assert ( + abs(fitted_val - true_val) < 0.1 + ), f"{param}: fitted={fitted_val:.4f}, true={true_val:.4f}" + else: + rel_error = abs(fitted_val - true_val) / abs(true_val) + assert rel_error < 0.02, ( + f"{param}: fitted={fitted_val:.4f}, true={true_val:.4f}, " + f"rel_error={rel_error:.2%}" + ) + + +# ----------------------------------------------------------------------------- +# HingeMax E3. Noisy Data Tests +# ----------------------------------------------------------------------------- + + +def test_hingemax_fit_with_noise_recovers_parameters_within_2sigma(noisy_hingemax_data): + """Fitted HingeMax parameters should be within 2sigma of true values.""" + x, y, w, true_params = noisy_hingemax_data + + obj = HingeMax(x, y, guess_h=true_params["h"]) + obj.make_fit() + + for param, true_val in true_params.items(): + fitted_val = obj.popt[param] + sigma = obj.psigma[param] + deviation = abs(fitted_val - true_val) + + assert deviation < 2 * sigma, ( + f"{param}: |fitted({fitted_val:.4f}) - true({true_val:.4f})| = " + f"{deviation:.4f} exceeds 2sigma = {2*sigma:.4f}" + ) + + +# ----------------------------------------------------------------------------- +# HingeMax E4. Derived Property Tests +# ----------------------------------------------------------------------------- + + +def test_hingemax_m2_property_is_consistent(clean_hingemax_data): + """Verify m2 = m1*(h-x1)/(h-x2) for fitted params. + + True: m2 = -2*(5-0)/(5-10) = -10/(-5) = 2 + """ + x, y, w, true_params, derived = clean_hingemax_data + + obj = HingeMax(x, y, guess_h=true_params["h"]) + obj.make_fit() + + m1 = obj.popt["m1"] + x1 = obj.popt["x1"] + x2 = obj.popt["x2"] + h = obj.popt["h"] + + m2_expected = m1 * (h - x1) / (h - x2) + + np.testing.assert_allclose( + obj.m2, + m2_expected, + rtol=1e-6, + err_msg="HingeMax.m2 should equal m1*(h-x1)/(h-x2)", + ) + + +def test_hingemax_theta_property_is_consistent(clean_hingemax_data): + """Verify theta = arctan(m1) - arctan(m2) for fitted params.""" + x, y, w, true_params, derived = clean_hingemax_data + + obj = HingeMax(x, y, guess_h=true_params["h"]) + obj.make_fit() + + m1 = obj.popt["m1"] + m2 = obj.m2 + theta_expected = np.arctan(m1) - np.arctan(m2) + + np.testing.assert_allclose( + obj.theta, + theta_expected, + rtol=1e-6, + err_msg="HingeMax.theta should equal arctan(m1) - arctan(m2)", + ) + + +def test_hingemax_lines_intersect_at_hinge(clean_hingemax_data): + """Verify both lines pass through (h, yh) where yh = m1*(h-x1).""" + x, y, w, true_params, derived = clean_hingemax_data + + obj = HingeMax(x, y, guess_h=true_params["h"]) + obj.make_fit() + + h = obj.popt["h"] + m1 = obj.popt["m1"] + x1 = obj.popt["x1"] + x2 = obj.popt["x2"] + m2 = obj.m2 + + yh_from_line1 = m1 * (h - x1) + yh_from_line2 = m2 * (h - x2) + + np.testing.assert_allclose( + yh_from_line1, + yh_from_line2, + rtol=1e-6, + err_msg="Both lines should pass through hinge point", + ) + + +# ----------------------------------------------------------------------------- +# HingeMax E5. Edge Cases and Interface Tests +# ----------------------------------------------------------------------------- + + +def test_hingemax_function_signature(): + """Test that HingeMax function has correct parameter signature.""" + x = np.array([0.0, 5.0, 10.0]) + y = np.array([0.0, -10.0, 0.0]) + obj = HingeMax(x, y, guess_h=5.0) + + sig = inspect.signature(obj.function) + params = tuple(sig.parameters.keys()) + + assert params == ( + "x", + "m1", + "x1", + "x2", + "h", + ), f"Function should have signature (x, m1, x1, x2, h), got {params}" + + +def test_hingemax_popt_has_correct_keys(clean_hingemax_data): + """Test that HingeMax popt contains expected parameter names.""" + x, y, w, true_params, derived = clean_hingemax_data + + obj = HingeMax(x, y, guess_h=true_params["h"]) + obj.make_fit() + + expected_keys = {"m1", "x1", "x2", "h"} + actual_keys = set(obj.popt.keys()) + + assert ( + actual_keys == expected_keys + ), f"popt keys should be {expected_keys}, got {actual_keys}" + + +def test_hingemax_callable_interface(clean_hingemax_data): + """Test that fitted HingeMax object is callable and returns correct shape.""" + x, y, w, true_params, derived = clean_hingemax_data + + obj = HingeMax(x, y, guess_h=true_params["h"]) + obj.make_fit() + + x_test = np.array([1.0, 5.0, 9.0]) + y_pred = obj(x_test) + + assert ( + y_pred.shape == x_test.shape + ), f"Predicted shape {y_pred.shape} should match input shape {x_test.shape}" + assert np.all(np.isfinite(y_pred)), "All predicted values should be finite" + + +def test_hingemax_insufficient_data_raises_error(): + """Fitting HingeMax with insufficient data should raise InsufficientDataError.""" + x = np.array([1.0, 2.0, 3.0]) + y = np.array([-2.0, -4.0, -3.0]) + + obj = HingeMax(x, y, guess_h=2.0) + + with pytest.raises(InsufficientDataError): + obj.make_fit() + + +# ============================================================================= +# HingeAtPoint Tests +# ============================================================================= +# HingeAtPoint: f(x) = np.minimum(y1, y2) where lines pass through (xh, yh) +# Parameters: xh, yh, m1, m2 +# x1 = xh - yh/m1, x2 = xh - yh/m2 +# Derived: x_intercepts (namedtuple with x1, x2) +# ============================================================================= + + +# ----------------------------------------------------------------------------- +# HingeAtPoint Fixtures +# ----------------------------------------------------------------------------- + + +@pytest.fixture +def clean_hingeatpoint_data(): + """Perfect HingeAtPoint data with lines meeting at (xh, yh). + + Parameters: xh=5, yh=10, m1=2, m2=-1 + x1 = 5 - 10/2 = 0 + x2 = 5 - 10/(-1) = 15 + Line1: y = 2*(x - 0) = 2x + Line2: y = -1*(x - 15) = -x + 15 + """ + true_params = {"xh": 5.0, "yh": 10.0, "m1": 2.0, "m2": -1.0} + x1 = true_params["xh"] - true_params["yh"] / true_params["m1"] + x2 = true_params["xh"] - true_params["yh"] / true_params["m2"] + + x = np.linspace(-2, 20, 300) + y1 = true_params["m1"] * (x - x1) + y2 = true_params["m2"] * (x - x2) + y = np.minimum(y1, y2) + w = np.ones_like(x) + return x, y, w, true_params, {"x1": x1, "x2": x2} + + +@pytest.fixture +def noisy_hingeatpoint_data(): + """HingeAtPoint data with 5% Gaussian noise. + + Parameters: xh=5, yh=10, m1=2, m2=-1 + """ + rng = np.random.default_rng(42) + true_params = {"xh": 5.0, "yh": 10.0, "m1": 2.0, "m2": -1.0} + x1 = true_params["xh"] - true_params["yh"] / true_params["m1"] + x2 = true_params["xh"] - true_params["yh"] / true_params["m2"] + + x = np.linspace(-2, 20, 300) + y1 = true_params["m1"] * (x - x1) + y2 = true_params["m2"] * (x - x2) + y_true = np.minimum(y1, y2) + noise_std = 0.5 + y = y_true + rng.normal(0, noise_std, len(x)) + w = np.ones_like(x) / noise_std + return x, y, w, true_params + + +@pytest.fixture +def hingeatpoint_positive_slopes_data(): + """HingeAtPoint where both slopes are positive. + + Parameters: xh=5, yh=10, m1=3, m2=0.5 + x1 = 5 - 10/3 = 1.667 + x2 = 5 - 10/0.5 = -15 + """ + true_params = {"xh": 5.0, "yh": 10.0, "m1": 3.0, "m2": 0.5} + x1 = true_params["xh"] - true_params["yh"] / true_params["m1"] + x2 = true_params["xh"] - true_params["yh"] / true_params["m2"] + + x = np.linspace(-5, 15, 300) + y1 = true_params["m1"] * (x - x1) + y2 = true_params["m2"] * (x - x2) + y = np.minimum(y1, y2) + w = np.ones_like(x) + return x, y, w, true_params + + +# ----------------------------------------------------------------------------- +# HingeAtPoint E1. Function Evaluation Tests +# ----------------------------------------------------------------------------- + + +def test_hingeatpoint_func_evaluates_line1_region_correctly(): + """HingeAtPoint should follow line1 where m1*(x-x1) < m2*(x-x2). + + With xh=5, yh=10, m1=2, m2=-1: x1=0, x2=15 + For x<5: 2x < -x+15, so line1 dominates minimum. + """ + xh, yh, m1, m2 = 5.0, 10.0, 2.0, -1.0 + x1 = xh - yh / m1 # = 0 + x2 = xh - yh / m2 # = 15 + + x_test = np.array([0.0, 1.0, 2.0, 3.0, 4.0]) + expected = m1 * (x_test - x1) # [0, 2, 4, 6, 8] + + x_dummy = np.linspace(0, 15, 50) + y_dummy = np.minimum(m1 * x_dummy, m2 * (x_dummy - x2)) + obj = HingeAtPoint(x_dummy, y_dummy, guess_xh=xh, guess_yh=yh) + result = obj.function(x_test, xh, yh, m1, m2) + + np.testing.assert_allclose( + result, + expected, + rtol=1e-10, + err_msg="HingeAtPoint should follow line1 in left region", + ) + + +def test_hingeatpoint_func_passes_through_hinge_point(): + """HingeAtPoint should pass through (xh, yh). + + At x=xh: f(xh) = yh + """ + xh, yh, m1, m2 = 5.0, 10.0, 2.0, -1.0 + + x_test = np.array([xh]) + expected = np.array([yh]) + + x_dummy = np.linspace(0, 15, 50) + y_dummy = np.minimum(m1 * x_dummy, m2 * (x_dummy - 15)) + obj = HingeAtPoint(x_dummy, y_dummy, guess_xh=xh, guess_yh=yh) + result = obj.function(x_test, xh, yh, m1, m2) + + np.testing.assert_allclose( + result, + expected, + rtol=1e-10, + err_msg="HingeAtPoint should pass through (xh, yh)", + ) + + +def test_hingeatpoint_func_evaluates_line2_region_correctly(): + """HingeAtPoint should follow line2 where m2*(x-x2) < m1*(x-x1). + + For x>5: -x+15 < 2x + """ + xh, yh, m1, m2 = 5.0, 10.0, 2.0, -1.0 + x1 = xh - yh / m1 # = 0 + x2 = xh - yh / m2 # = 15 + + x_test = np.array([6.0, 8.0, 10.0, 12.0, 15.0]) + expected = m2 * (x_test - x2) # [-1*(x-15)] = [9, 7, 5, 3, 0] + + x_dummy = np.linspace(0, 15, 50) + y_dummy = np.minimum(m1 * x_dummy, m2 * (x_dummy - x2)) + obj = HingeAtPoint(x_dummy, y_dummy, guess_xh=xh, guess_yh=yh) + result = obj.function(x_test, xh, yh, m1, m2) + + np.testing.assert_allclose( + result, + expected, + rtol=1e-10, + err_msg="HingeAtPoint should follow line2 in right region", + ) + + +# ----------------------------------------------------------------------------- +# HingeAtPoint E2. Parameter Recovery Tests (Clean Data) +# ----------------------------------------------------------------------------- + + +def test_hingeatpoint_fit_recovers_exact_parameters_from_clean_data( + clean_hingeatpoint_data, +): + """Fitting noise-free HingeAtPoint data should recover parameters within 2%.""" + x, y, w, true_params, derived = clean_hingeatpoint_data + + obj = HingeAtPoint(x, y, guess_xh=true_params["xh"], guess_yh=true_params["yh"]) + obj.make_fit() + + for param, true_val in true_params.items(): + fitted_val = obj.popt[param] + if abs(true_val) < 0.1: + assert ( + abs(fitted_val - true_val) < 0.05 + ), f"{param}: fitted={fitted_val:.4f}, true={true_val:.4f}" + else: + rel_error = abs(fitted_val - true_val) / abs(true_val) + assert rel_error < 0.02, ( + f"{param}: fitted={fitted_val:.4f}, true={true_val:.4f}, " + f"rel_error={rel_error:.2%} exceeds 2% tolerance" + ) + + +def test_hingeatpoint_fit_recovers_positive_slopes_parameters( + hingeatpoint_positive_slopes_data, +): + """Fitting HingeAtPoint with positive slopes should recover parameters within 2%.""" + x, y, w, true_params = hingeatpoint_positive_slopes_data + + obj = HingeAtPoint(x, y, guess_xh=true_params["xh"], guess_yh=true_params["yh"]) + obj.make_fit() + + for param, true_val in true_params.items(): + fitted_val = obj.popt[param] + if abs(true_val) < 0.1: + assert ( + abs(fitted_val - true_val) < 0.1 + ), f"{param}: fitted={fitted_val:.4f}, true={true_val:.4f}" + else: + rel_error = abs(fitted_val - true_val) / abs(true_val) + assert rel_error < 0.02, ( + f"{param}: fitted={fitted_val:.4f}, true={true_val:.4f}, " + f"rel_error={rel_error:.2%}" + ) + + +# ----------------------------------------------------------------------------- +# HingeAtPoint E3. Noisy Data Tests +# ----------------------------------------------------------------------------- + + +def test_hingeatpoint_fit_with_noise_recovers_parameters_within_2sigma( + noisy_hingeatpoint_data, +): + """Fitted HingeAtPoint parameters should be within 2sigma of true values.""" + x, y, w, true_params = noisy_hingeatpoint_data + + obj = HingeAtPoint(x, y, guess_xh=true_params["xh"], guess_yh=true_params["yh"]) + obj.make_fit() + + for param, true_val in true_params.items(): + fitted_val = obj.popt[param] + sigma = obj.psigma[param] + deviation = abs(fitted_val - true_val) + + assert deviation < 2 * sigma, ( + f"{param}: |fitted({fitted_val:.4f}) - true({true_val:.4f})| = " + f"{deviation:.4f} exceeds 2sigma = {2*sigma:.4f}" + ) + + +# ----------------------------------------------------------------------------- +# HingeAtPoint E4. Derived Property Tests +# ----------------------------------------------------------------------------- + + +def test_hingeatpoint_x_intercepts_property_has_correct_values(clean_hingeatpoint_data): + """Verify x_intercepts.x1 = xh - yh/m1 and x_intercepts.x2 = xh - yh/m2. + + True: x1 = 5 - 10/2 = 0, x2 = 5 - 10/(-1) = 15 + """ + x, y, w, true_params, derived = clean_hingeatpoint_data + + obj = HingeAtPoint(x, y, guess_xh=true_params["xh"], guess_yh=true_params["yh"]) + obj.make_fit() + + xh = obj.popt["xh"] + yh = obj.popt["yh"] + m1 = obj.popt["m1"] + m2 = obj.popt["m2"] + + x1_expected = xh - yh / m1 + x2_expected = xh - yh / m2 + + np.testing.assert_allclose( + obj.x_intercepts.x1, + x1_expected, + rtol=1e-6, + err_msg="x_intercepts.x1 should equal xh - yh/m1", + ) + np.testing.assert_allclose( + obj.x_intercepts.x2, + x2_expected, + rtol=1e-6, + err_msg="x_intercepts.x2 should equal xh - yh/m2", + ) + + +def test_hingeatpoint_x_intercepts_is_namedtuple(clean_hingeatpoint_data): + """Verify x_intercepts is a namedtuple with x1 and x2 attributes.""" + x, y, w, true_params, derived = clean_hingeatpoint_data + + obj = HingeAtPoint(x, y, guess_xh=true_params["xh"], guess_yh=true_params["yh"]) + obj.make_fit() + + assert hasattr(obj.x_intercepts, "x1"), "x_intercepts should have x1 attribute" + assert hasattr(obj.x_intercepts, "x2"), "x_intercepts should have x2 attribute" + + +def test_hingeatpoint_lines_pass_through_hinge(clean_hingeatpoint_data): + """Verify both lines pass through (xh, yh).""" + x, y, w, true_params, derived = clean_hingeatpoint_data + + obj = HingeAtPoint(x, y, guess_xh=true_params["xh"], guess_yh=true_params["yh"]) + obj.make_fit() + + xh = obj.popt["xh"] + yh = obj.popt["yh"] + m1 = obj.popt["m1"] + m2 = obj.popt["m2"] + x1 = obj.x_intercepts.x1 + x2 = obj.x_intercepts.x2 + + yh_from_line1 = m1 * (xh - x1) + yh_from_line2 = m2 * (xh - x2) + + np.testing.assert_allclose( + yh_from_line1, + yh, + rtol=1e-6, + err_msg="Line1 should pass through (xh, yh)", + ) + np.testing.assert_allclose( + yh_from_line2, + yh, + rtol=1e-6, + err_msg="Line2 should pass through (xh, yh)", + ) + + +# ----------------------------------------------------------------------------- +# HingeAtPoint E5. Edge Cases and Interface Tests +# ----------------------------------------------------------------------------- + + +def test_hingeatpoint_function_signature(): + """Test that HingeAtPoint function has correct parameter signature.""" + x = np.array([0.0, 5.0, 10.0]) + y = np.array([0.0, 10.0, 5.0]) + obj = HingeAtPoint(x, y, guess_xh=5.0, guess_yh=10.0) + + sig = inspect.signature(obj.function) + params = tuple(sig.parameters.keys()) + + assert params == ( + "x", + "xh", + "yh", + "m1", + "m2", + ), f"Function should have signature (x, xh, yh, m1, m2), got {params}" + + +def test_hingeatpoint_popt_has_correct_keys(clean_hingeatpoint_data): + """Test that HingeAtPoint popt contains expected parameter names.""" + x, y, w, true_params, derived = clean_hingeatpoint_data + + obj = HingeAtPoint(x, y, guess_xh=true_params["xh"], guess_yh=true_params["yh"]) + obj.make_fit() + + expected_keys = {"xh", "yh", "m1", "m2"} + actual_keys = set(obj.popt.keys()) + + assert ( + actual_keys == expected_keys + ), f"popt keys should be {expected_keys}, got {actual_keys}" + + +def test_hingeatpoint_callable_interface(clean_hingeatpoint_data): + """Test that fitted HingeAtPoint object is callable and returns correct shape.""" + x, y, w, true_params, derived = clean_hingeatpoint_data + + obj = HingeAtPoint(x, y, guess_xh=true_params["xh"], guess_yh=true_params["yh"]) + obj.make_fit() + + x_test = np.array([1.0, 5.0, 10.0]) + y_pred = obj(x_test) + + assert ( + y_pred.shape == x_test.shape + ), f"Predicted shape {y_pred.shape} should match input shape {x_test.shape}" + assert np.all(np.isfinite(y_pred)), "All predicted values should be finite" + + +def test_hingeatpoint_insufficient_data_raises_error(): + """Fitting HingeAtPoint with insufficient data should raise InsufficientDataError.""" + x = np.array([1.0, 2.0, 3.0]) + y = np.array([2.0, 4.0, 3.0]) + + obj = HingeAtPoint(x, y, guess_xh=2.0, guess_yh=4.0) + + with pytest.raises(InsufficientDataError): + obj.make_fit() + + +def test_hingeatpoint_psigma_values_are_positive(noisy_hingeatpoint_data): + """Test that all parameter uncertainties are positive.""" + x, y, w, true_params = noisy_hingeatpoint_data + + obj = HingeAtPoint(x, y, guess_xh=true_params["xh"], guess_yh=true_params["yh"]) + obj.make_fit() + + for param, sigma in obj.psigma.items(): + assert sigma > 0, f"psigma['{param}'] = {sigma} should be positive"