diff --git a/CHANGELOG.md b/CHANGELOG.md index ea2d589..697ab3b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,6 +5,9 @@ The format is based on [Keep a Changelog](http://keepachangelog.com/en/1.0.0/) and this project adheres to [Semantic Versioning](http://semver.org/spec/v2.0.0.html). ## [Unreleased] +### Added +- `CurveFit` class +- `FitResult` object ### Changed - `_model_function` method refactored - `_model_function` method renamed to `model_function` diff --git a/drux/__init__.py b/drux/__init__.py index e6788d6..cef8075 100644 --- a/drux/__init__.py +++ b/drux/__init__.py @@ -7,5 +7,6 @@ from .first_order import FirstOrderModel, FirstOrderParameters from .weibull import WeibullModel, WeibullParameters from .hopfenberg import HopfenbergModel, HopfenbergParameters +from .curve_fit import CurveFit, FitResult __version__ = DRUX_VERSION diff --git a/drux/curve_fit.py b/drux/curve_fit.py new file mode 100644 index 0000000..6b62c8a --- /dev/null +++ b/drux/curve_fit.py @@ -0,0 +1,175 @@ +# -*- coding: utf-8 -*- +"""Drux curve fitting utilities.""" + +import inspect +from dataclasses import dataclass +from typing import Dict, List, Optional, Sequence, Tuple, Type + +import numpy as np +from scipy.optimize import curve_fit + +from .base_model import DrugReleaseModel +from .first_order import FirstOrderModel +from .higuchi import HiguchiModel +from .hopfenberg import HopfenbergModel +from .weibull import WeibullModel +from .zero_order import ZeroOrderModel +from .messages import ( + ERROR_UNKNOWN_MODEL, + ERROR_TIME_RELEASE_LENGTH_MISMATCH, + ERROR_INSUFFICIENT_DATA_POINTS, + ERROR_UNKNOWN_FIT_PARAMETER, + ERROR_NO_FREE_PARAMETERS, + ERROR_NO_FIT_RESULT, +) + +# Single point of registration mapping a model name to its class. +MODEL_CLASSES = { + "zero_order": ZeroOrderModel, + "first_order": FirstOrderModel, + "higuchi": HiguchiModel, + "weibull": WeibullModel, + "hopfenberg": HopfenbergModel, +} + + +def _model_parameter_names(model_class: Type[DrugReleaseModel]) -> List[str]: + """Return a model's parameter names, in declaration order, excluding time. + + :param model_class: model class exposing a `model_function` signature + """ + parameters = inspect.signature(model_class.model_function).parameters + return [name for name in parameters if name != "t"] + + +@dataclass +class FitResult: + """ + Result of fitting a drug release model to experimental data. + + Attributes: + model_name (str): name of the fitted model + parameters (Dict[str, float]): known and fitted parameter values, keyed by name + r_squared (float): coefficient of determination of the fit + model (DrugReleaseModel): model instance initialized with the fitted parameters + """ + + model_name: str + parameters: Dict[str, float] + r_squared: float + model: DrugReleaseModel + + +class CurveFit: + """Fit a registered drug release model to experimental time-series data.""" + + def __init__( + self, + model_name: str, + time: Sequence[float], + release_profile: Sequence[float], + known_parameters: Optional[Dict[str, float]] = None, + ) -> None: + """ + Initialize the curve fitter. + + :param model_name: name of a registered model (see `MODEL_CLASSES`) + :param time: time points of the experimental data (s) + :param release_profile: a sequence of measured drug release at each time point + :param known_parameters: parameter values to hold fixed instead of fitting + """ + if model_name not in MODEL_CLASSES: + raise ValueError(ERROR_UNKNOWN_MODEL.format(model_name=model_name, available_models=sorted(MODEL_CLASSES))) + + time = np.asarray(time, dtype=float) + release_profile = np.asarray(release_profile, dtype=float) + if time.shape != release_profile.shape: + raise ValueError(ERROR_TIME_RELEASE_LENGTH_MISMATCH) + if time.size < 2: + raise ValueError(ERROR_INSUFFICIENT_DATA_POINTS) + + self._model_name = model_name + self._model_class = MODEL_CLASSES[model_name] + self._time = time + self._release_profile = release_profile + self._known_parameters = known_parameters or {} + self._parameter_names = _model_parameter_names(self._model_class) + + unknown_keys = sorted(set(self._known_parameters) - set(self._parameter_names)) + if unknown_keys: + raise ValueError( + ERROR_UNKNOWN_FIT_PARAMETER.format(unknown_keys=unknown_keys, model_name=model_name, parameter_names=self._parameter_names) + ) + + self._free_parameters = [p for p in self._parameter_names if p not in self._known_parameters] + if not self._free_parameters: + raise ValueError(ERROR_NO_FREE_PARAMETERS) + + self._fit_result: Optional[FitResult] = None + + def _merge_parameters(self, free_values: Sequence[float]) -> Dict[str, float]: + """Merge fitted values for free parameters with known fixed parameters. + + :param free_values: fitted values in the model's declared free-parameter order + """ + free_values = iter(free_values) + return { + name: self._known_parameters[name] if name in self._known_parameters else next(free_values) + for name in self._parameter_names + } + + def _equation(self, t: np.ndarray, *free_values: float) -> np.ndarray: + """Evaluate the model equation for a set of free parameter values. + + :param t: time values where the model should be evaluated + :param free_values: free parameter values passed by the optimizer + """ + parameters = self._merge_parameters(free_values) + return np.vectorize(lambda ti: self._model_class.model_function(ti, **parameters))(t) + + def fit( + self, + initial_guess: Optional[Sequence[float]] = None, + bounds: Optional[Tuple[Sequence[float], Sequence[float]]] = None, + ) -> FitResult: + """ + Estimate the unknown parameters via non-linear least squares. + + :param initial_guess: initial guess for each free parameter in declared order (default: 1.0 for each) + :param bounds: `(lower, upper)` bounds for free parameters forwarded to `curve_fit` (default: (~0, inf)) + """ + n_free = len(self._free_parameters) + p0 = list(initial_guess) if initial_guess is not None else [1.0] * n_free + lower_upper = bounds or ([1e-10] * n_free, [np.inf] * n_free) + + fitted_values, _ = curve_fit( + self._equation, self._time, self._release_profile, p0=p0, bounds=lower_upper, maxfev=10000 + ) + parameters = self._merge_parameters(fitted_values) + + predicted = self._equation(self._time, *fitted_values) + residual_sum_of_squares = np.sum((self._release_profile - predicted) ** 2) + total_sum_of_squares = np.sum((self._release_profile - np.mean(self._release_profile)) ** 2) + r_squared = 1 - residual_sum_of_squares / total_sum_of_squares if total_sum_of_squares > 0 else 0.0 + + fit_result = FitResult( + model_name=self._model_name, + parameters=parameters, + r_squared=r_squared, + model=self._model_class(**parameters), + ) + self._fit_result = fit_result + return fit_result + + def get_result(self) -> FitResult: + """Return the most recent fit result.""" + if self._fit_result is None: + raise ValueError(ERROR_NO_FIT_RESULT) + return self._fit_result + + def __repr__(self) -> str: + """Return a string representation of the curve fitter.""" + if self._fit_result is None: + return f"drux.CurveFit({self._model_name}, not fitted)" + params = ", ".join(f"{k}={v:.4f}" for k, v in self._fit_result.parameters.items()) + return f"drux.CurveFit({self._model_name}: {params}, R^2={self._fit_result.r_squared:.4f})" diff --git a/drux/messages.py b/drux/messages.py index 9d7b189..965607e 100644 --- a/drux/messages.py +++ b/drux/messages.py @@ -44,3 +44,11 @@ ERROR_INVALID_EROSION_CONSTANT = "Erosion rate constant (k0) must be non-negative." ERROR_INVALID_INITIAL_RADIUS = "Initial radius or half-thickness (a0) must be positive." ERROR_INVALID_GEOMETRY_FACTOR = "Geometry factor (n) must be 1 (slab), 2 (cylinder), or 3 (sphere)." + +# Error messages for curve fitting +ERROR_UNKNOWN_MODEL = "Unknown model '{model_name}'. Available models: {available_models}." +ERROR_TIME_RELEASE_LENGTH_MISMATCH = "Time and release profile must have the same length." +ERROR_INSUFFICIENT_DATA_POINTS = "At least 2 data points are required for fitting." +ERROR_UNKNOWN_FIT_PARAMETER = "Unknown parameter(s) {unknown_keys} for model '{model_name}'. Valid parameters: {parameter_names}." +ERROR_NO_FREE_PARAMETERS = "All parameters are known; there is nothing left to fit." +ERROR_NO_FIT_RESULT = "No fit result available. Run fit() first." diff --git a/requirements-latest.txt b/requirements-latest.txt index 013c21d..30990cc 100644 --- a/requirements-latest.txt +++ b/requirements-latest.txt @@ -1,2 +1,3 @@ numpy==2.2.6 matplotlib==3.10.3 +scipy==1.16.3 diff --git a/setup.py b/setup.py index 2147dca..2388148 100644 --- a/setup.py +++ b/setup.py @@ -37,7 +37,8 @@ def read_description() -> str: }, install_requires=[ 'numpy>=1.21.0', - 'matplotlib>=3.2.2' + 'matplotlib>=3.2.2', + 'scipy>=1.7.3' ], python_requires='>=3.7', classifiers=[