from __future__ import annotations
import logging
from functools import cached_property
import lmfit
import numpy as np
from scipy.interpolate import pchip_interpolate
from ewoksxas.fit.data import ArrayF64, XYData
logger = logging.getLogger(__name__)
[docs]
class LCAModel:
"""A class for linear combination analysis of 1D X and Y data."""
def __init__(
self,
sample: XYData,
components: dict[str, XYData],
xmin: float = -np.inf,
xmax: float = np.inf,
) -> None:
"""Initialize an analysis for a single instance of sample data.
Args:
sample: x and y data for observations.
components: names and x and y data for known components.
xmin: The lower bound for the Region of Interest (ROI).
xmax: The upper bound for the Region of Interest (ROI).
Raises:
ValueError: If no components have been added to the model.
"""
if not components:
raise ValueError("No components added. Cannot perform fit.")
self.xmin: float = max(sample.x.min(), xmin)
self.xmax: float = min(sample.x.max(), xmax)
self.sample: XYData = sample
self.components: dict[str, XYData] = components
self.roi: XYData = self._set_roi()
self._ncomps: int = len(self.components)
def _set_roi(self) -> XYData:
"""Find the region of interest mask and return x and y values in range."""
if self.xmin == self.sample.x.min() and self.xmax == self.sample.x.max():
return self.sample
roi_mask = (self.sample.x >= self.xmin) & (self.sample.x <= self.xmax)
return XYData(x=self.sample.x[roi_mask], y=self.sample.y[roi_mask])
def _build_component_matrix(self) -> tuple[ArrayF64, list[str]]:
"""Interpolates component y-values onto the ROI x-grid to build a 2D matrix."""
matrix = np.empty((len(self.components), self.roi.x.size))
names: list[str] = []
for i, (name, data) in enumerate(self.components.items()):
matrix[i] = pchip_interpolate(xi=data.x, yi=data.y, x=self.roi.x)
names.append(name)
return matrix, names
def _objective(self, params: lmfit.Parameters, matrix: ArrayF64) -> ArrayF64:
"""Calculates the residual array for the minimization algorithm."""
weights = np.array([params[f"P{i}"].value for i in range(self._ncomps)])
model_y = np.dot(weights, matrix)
return model_y - self.roi.y
[docs]
def fit(self, **kwargs) -> FitResult:
"""Executes the least-squares minimization.
Args:
**kwargs: keyword arguments passed to cipy.optimize.least_squares.
"""
matrix, parameter_names = self._build_component_matrix()
params = lmfit.Parameters()
initial_guess = 1.0 / self._ncomps
for i in range(self._ncomps):
params.add(f"P{i}", value=initial_guess, min=0.0)
result = lmfit.minimize(
fcn=self._objective,
params=params,
args=(matrix,),
method=kwargs.pop("method", "least_squares"),
**kwargs,
)
return FitResult(
fit=result,
sample=self.roi,
component_matrix=matrix,
parameter_names=parameter_names,
)
[docs]
@classmethod
def from_dicts(
cls,
sample: dict[str, ArrayF64],
components: dict[str, dict[str, ArrayF64]],
parameters: dict[str, float] | None = None,
) -> LCAModel:
"""Constructs an LCAModel entirely from dictionaries.
Args:
sample: x and y data for observations.
components: component names mapped to x and y data dicts.
parameters: region of interest boundaries (xmin, xmax).
"""
parameters: dict[str, float] = parameters or {}
return cls(
sample=XYData(x=sample["x"], y=sample["y"]),
components={
name: XYData(x=data["x"], y=data["y"])
for name, data in components.items()
},
xmin=parameters.get("xmin", -np.inf),
xmax=parameters.get("xmax", np.inf),
)
[docs]
class FitResult:
"""Fitting statistics of scipy.optimize.least_squares outputs."""
def __init__(
self,
fit: lmfit.minimizer.MinimizerResult,
sample: XYData,
component_matrix: ArrayF64,
parameter_names: list[str],
):
"""Fitting statistics of scipy.optimize.least_squares outputs.
Args:
fit: output from the least_squares method from scipy.optimize.
sample: x and y data of the observed sample data.
component_matrix: ND matrix of component y-values.
parameter_names: names of the parameters (e.g. Parameter 1, Parameters 2).
"""
self.result: lmfit.minimizer.MinimizerResult = fit
self.names: list[str] = parameter_names
self._sample_sst = np.square(sample.y - sample.y.mean()).sum()
self._params: lmfit.Parameters = self.result.params # type: ignore[attr-defined]
self.model = XYData(x=sample.x, y=np.dot(self.solution, component_matrix))
@cached_property
def solution_norm(self) -> ArrayF64:
"""Normalized values of parameters at the minimized solution."""
return self.solution / self.solution.sum()
@cached_property
def sigma_norm(self) -> ArrayF64:
"""Normalized deviation of parameter values at the minimized solution."""
return self.sigma / self.solution.sum()
@cached_property
def solution(self) -> ArrayF64:
"""Values of parameters at the minimized least-squares solution."""
return np.array([self._params[f"P{i}"].value for i in range(len(self.names))])
@cached_property
def sigma(self) -> ArrayF64:
"""Deviation (standard error) of parameter values at the minimized solution."""
stderrs = [self._params[f"P{i}"].stderr for i in range(len(self.names))]
return np.array([np.nan if se is None else se for se in stderrs])
@cached_property
def rmse(self) -> float:
"""Root mean squared error of the minimized least-squares solution."""
return np.sqrt(float(self.result.redchi))
@cached_property
def r2(self) -> float:
"""Coefficient of determination of the minimized least-squares solution."""
return 1 - (self.result.chisqr / self._sample_sst)
[docs]
def summary(self) -> str:
"""Returns a string formatted fit report generated by lmfit."""
return lmfit.fit_report(self.result)
def __repr__(self):
np.set_printoptions(threshold=10)
ilen = len("Component Names:")
return (
f"{'-' * 120}\n"
f"Linear Combination Analysis Results\n"
f"{'-' * 120}\n"
f"{'Component Names:':>{ilen}} {self.names}\n"
f"{'Solution:':>{ilen}} {self.solution}\n"
f"{'Solution Sigma:':>{ilen}} {self.sigma}\n"
f"{'RMSE:':>{ilen}} {self.rmse}\n"
f"{'-' * 120}\n"
f"Fit Details\n"
f"{'-' * 120}\n"
f"{self.summary()}"
)