"""
template.py

    Updated 2026-08-30 for Python 3 and the current curve.fit numerical rules.

    Adapted from University of Toronto's odr_fit_to_data.py:

    https://www.physics.utoronto.ca/apl/python/odr_fit_to_data.py

    Copyright (c) 2011 University of Toronto
    Original Version   :   11 September 2011 by Michael Luzi
    Last Modification  :   11 January 2013 by David Bailey
    Contact: David Bailey <dbailey@physics.utoronto.ca>

    This is a small, editable example for fitting a user-defined model with
    scipy.odr. Ordinary least squares is used when x uncertainties are disabled;
    orthogonal distance regression is used when x uncertainties are enabled.

    Steps to use:
    (1) Change fitfunction(p, x), where p is the parameter vector and x is the
        independent-variable array.
    (2) Set data_file to a tab- or whitespace-delimited file containing
        x, dx, y, and dy columns.
    (3) Set p_guess and ifixb. An ifixb value of zero fixes that parameter.
    (4) Select which uncertainty columns to use with errorx and errory.

    For information on scipy.odr, see:
    https://docs.scipy.org/doc/scipy/reference/odr.html

    License: Released under the MIT License; the full terms are appended below.
"""

from __future__ import annotations

import inspect
from collections.abc import Callable, Sequence
from typing import Any

import numpy
import scipy.odr
import scipy.special
from matplotlib import pyplot


# Define the function you would like to fit. A Gaussian is provided as an example.
def fitfunction(p: Sequence[float], x: numpy.ndarray) -> numpy.ndarray:
    # Constant background, peak height, center, and standard-deviation width.
    return p[0] + p[1] * numpy.exp(-((x - p[2]) ** 2) / (2 * p[3] ** 2))


# Initial values for the parameters in fitfunction.
p_guess = (10.0, 200.0, 1173.9, 0.4)

# A zero fixes the corresponding parameter; a one allows it to vary.
ifixb = (1, 1, 1, 1)

# A tab- or whitespace-delimited file with x, dx, y, and dy columns.
data_file = "data.txt"

# Select whether the dx and dy columns participate in the fit.
errorx = True
errory = True


def _as_data_column(values: Sequence[float], name: str) -> numpy.ndarray:
    result = numpy.asarray(values, dtype=float)
    if result.ndim != 1 or not numpy.all(numpy.isfinite(result)):
        raise ValueError(f"{name} must be a one-dimensional array of finite numbers")
    return result


def _prepare_odr_data(
    x: numpy.ndarray,
    dx: numpy.ndarray,
    y: numpy.ndarray,
    dy: numpy.ndarray,
    *,
    use_x_errors: bool,
    use_y_errors: bool,
) -> tuple[scipy.odr.RealData, int]:
    """Apply the same four uncertainty-mode choices used by curve.fit."""

    if use_x_errors and use_y_errors:
        return scipy.odr.RealData(x=x, y=y, sx=dx, sy=dy), 0
    if use_y_errors:
        return scipy.odr.RealData(x=x, y=y, sy=dy), 2
    if use_x_errors:
        # ODR needs a finite y weight for the x-only case. This is the same
        # effectively-zero y uncertainty used by the curve.fit engine.
        residual_epsilon = numpy.sqrt(numpy.finfo(float).eps)
        return scipy.odr.RealData(x=x, y=y, sx=dx, sy=y * residual_epsilon), 0
    return scipy.odr.RealData(x=x, y=y), 2


def perform_fit(
    model_function: Callable[[Sequence[float], numpy.ndarray], numpy.ndarray],
    initial_parameters: Sequence[float],
    fixed_parameters: Sequence[int],
    x_values: Sequence[float],
    dx_values: Sequence[float],
    y_values: Sequence[float],
    dy_values: Sequence[float],
    *,
    use_x_errors: bool,
    use_y_errors: bool,
) -> dict[str, Any]:
    """Fit one dataset and return the values used by the example report."""

    x = _as_data_column(x_values, "x")
    dx = _as_data_column(dx_values, "dx")
    y = _as_data_column(y_values, "y")
    dy = _as_data_column(dy_values, "dy")
    if len(x) < 2 or not (len(x) == len(dx) == len(y) == len(dy)):
        raise ValueError("x, dx, y, and dy must contain the same two or more rows")
    if use_x_errors and numpy.any(dx == 0.0):
        raise ValueError("selected x uncertainties must be nonzero")
    if use_y_errors and numpy.any(dy == 0.0):
        raise ValueError("selected y uncertainties must be nonzero")

    beta0 = tuple(float(value) for value in initial_parameters)
    parameter_mask = tuple(int(value) for value in fixed_parameters)
    if len(beta0) == 0 or len(parameter_mask) != len(beta0):
        raise ValueError("p_guess and ifixb must contain the same nonzero number of values")
    if any(value not in {0, 1} for value in parameter_mask):
        raise ValueError("ifixb values must be zero (fixed) or one (free)")

    fixed_count = sum(value == 0 for value in parameter_mask)
    degrees_of_freedom = len(x) - (len(beta0) - fixed_count)
    if degrees_of_freedom < 0:
        raise ValueError("enter at least as many data points as free fit parameters")

    data, fit_type = _prepare_odr_data(
        x,
        dx,
        y,
        dy,
        use_x_errors=use_x_errors,
        use_y_errors=use_y_errors,
    )
    model = scipy.odr.Model(model_function)
    # Match curve.fit's residual stopping policy; this is not a parameter-error bound.
    odr = scipy.odr.ODR(
        data, model, beta0=beta0, maxit=10_000, sstol=1e-14, ifixb=parameter_mask
    )
    odr.set_job(fit_type=fit_type, deriv=1)
    output = odr.run()

    # scipy.odr exposes cov_beta before residual-variance scaling. curve.fit
    # reports one-standard-deviation parameter uncertainties from the scaled
    # covariance for every uncertainty mode.
    output.cov_beta = output.cov_beta * output.res_var
    covariance = numpy.asarray(output.cov_beta, dtype=float)
    uncertainties = numpy.sqrt(numpy.diagonal(covariance))

    delta = numpy.asarray(output.delta, dtype=float)
    eps = numpy.asarray(output.eps, dtype=float)
    if use_x_errors and not use_y_errors:
        # The artificial y weight is only an ODR implementation detail.
        eps = numpy.zeros_like(y)
        output.eps = eps

    if use_x_errors and use_y_errors:
        residual = -numpy.sign(eps) * numpy.sqrt(delta**2 / dx**2 + eps**2 / dy**2)
        residual_uncertainty = numpy.ones_like(y)
    elif use_y_errors:
        residual = -eps
        residual_uncertainty = dy
    elif use_x_errors:
        residual = -delta
        residual_uncertainty = dx
    else:
        residual = -eps
        residual_uncertainty = numpy.zeros_like(y)

    return {
        "output": output,
        "parameters": numpy.asarray(output.beta, dtype=float),
        "uncertainties": uncertainties,
        "covariance": covariance,
        "degrees_of_freedom": degrees_of_freedom,
        "residual": residual,
        "residual_uncertainty": residual_uncertainty,
        "x": x,
        "dx": dx,
        "y": y,
        "dy": dy,
    }


def main() -> None:
    # Lines starting with '#' are ignored. Each remaining row is x, dx, y, dy.
    x, dx, y, dy = numpy.loadtxt(data_file, comments="#", unpack=True)
    result = perform_fit(
        fitfunction,
        p_guess,
        ifixb,
        x,
        dx,
        y,
        dy,
        use_x_errors=errorx,
        use_y_errors=errory,
    )

    output = result["output"]
    parameters = result["parameters"]
    uncertainties = result["uncertainties"]
    covariance = result["covariance"]
    degrees_of_freedom = result["degrees_of_freedom"]

    print("***********************************************************")
    print("                    CURVE FIT TO DATA")
    print("***********************************************************\n")
    print("ODR algorithm stop reason: " + output.stopreason[0])
    print(f"\nFit {len(x)} data points from file: {data_file}")
    print("To model:")
    print(inspect.getsource(fitfunction))
    print("Estimated parameters and one-standard-deviation uncertainties")
    for index, (parameter, uncertainty) in enumerate(zip(parameters, uncertainties)):
        print(
            f"   p[{index}] = {parameter:10.5g} +/- {uncertainty:10.5g}"
            f"          (Starting guess: {p_guess[index]:10.5g})"
        )

    print("\nCorrelation matrix:")
    with numpy.errstate(divide="ignore", invalid="ignore"):
        scale = numpy.sqrt(numpy.outer(numpy.diagonal(covariance), numpy.diagonal(covariance)))
        correlation = numpy.divide(covariance, scale)
    for row in correlation:
        print(" ".join(f"{value: 8.3g}" for value in row))

    if errorx or errory:
        if degrees_of_freedom > 0:
            cdf = 100.0 * float(
                scipy.special.chdtrc(
                    degrees_of_freedom,
                    degrees_of_freedom * output.res_var,
                )
            )
            print(f"\nReduced chi-squared = {output.res_var:10.5f}, chi-squared CDF = {cdf:10.5f}%")
        else:
            print("\nReduced chi-squared is undefined with zero degrees of freedom.")
    else:
        print(f"\nSum of squared residuals = {numpy.sum(output.eps**2):10.5f}")

    figure = pyplot.figure(facecolor="0.98")
    fit_axes = figure.add_subplot(211)
    fit_axes.tick_params(labelbottom=False)
    fit_axes.set_ylabel("y")
    fit_axes.set_title("Curve Fit to Data")

    x_model = numpy.linspace(min(x), max(x), 1_000)
    fit_axes.plot(x, y, "ro", label="Data")
    fit_axes.plot(x_model, fitfunction(parameters, x_model), label="Fit")
    fit_axes.errorbar(
        x,
        y,
        xerr=dx if errorx else None,
        yerr=dy if errory else None,
        fmt="none",
        ecolor="red",
    )
    fit_axes.plot(
        x_model,
        fitfunction(p_guess, x_model),
        "g--",
        label="Starting guess",
    )

    # Draw the displacement from each observation to ODR's fitted point.
    for observed_x, observed_y, fitted_x, fitted_y in zip(x, y, output.xplus, output.y):
        fit_axes.plot([fitted_x, observed_x], [fitted_y, observed_y], "k-")
    fit_axes.legend(loc="best")
    fit_axes.grid()

    residual_axes = figure.add_subplot(212)
    residual_axes.errorbar(
        x=x,
        y=result["residual"],
        yerr=result["residual_uncertainty"],
        fmt="r+",
        label="Residuals",
    )
    residual_axes.set_xlim(fit_axes.get_xlim())
    residual_axes.axhline(y=0, color="blue")
    residual_axes.set_xlabel("x")
    residual_axes.set_ylabel("Residuals")
    residual_axes.grid()
    pyplot.show()


if __name__ == "__main__":
    main()


"""
Full text of MIT License:

    Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in
all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
THE SOFTWARE.
"""