From 7aecf962bacbd7b6f9733f4459d57c19e87a2387 Mon Sep 17 00:00:00 2001 From: Dominik Demuth Date: Thu, 1 Feb 2024 19:09:17 +0100 Subject: [PATCH] fit results are not complex for non-complex data --- src/nmreval/fit/minimizer.py | 98 ++++++++++++----------- src/nmreval/fit/model.py | 2 +- src/nmreval/fit/result.py | 147 ++++++++++++++++++----------------- 3 files changed, 132 insertions(+), 115 deletions(-) diff --git a/src/nmreval/fit/minimizer.py b/src/nmreval/fit/minimizer.py index 27e7492..18c3f31 100644 --- a/src/nmreval/fit/minimizer.py +++ b/src/nmreval/fit/minimizer.py @@ -361,7 +361,7 @@ class FitRoutine(object): with np.errstate(all='ignore'): res = optimize.least_squares(cost, p0, bounds=(lb, ub), max_nfev=500 * len(p0)) - err, corr, partial_corr = self._calc_error(res.jac, np.sum(res.fun**2), *res.jac.shape) + err, corr, partial_corr = _calc_error(res.jac, np.sum(res.fun**2), *res.jac.shape) self.make_results(data, res.x, var, data.para_keys, res.jac.shape, err=err, corr=corr, partial_corr=partial_corr) @@ -375,7 +375,7 @@ class FitRoutine(object): with np.errstate(all='ignore'): res = optimize.least_squares(cost, p0, bounds=(lb, ub), max_nfev=500 * len(p0)) - err, corr, partial_corr = self._calc_error(res.jac, np.sum(res.fun**2), *res.jac.shape) + err, corr, partial_corr = _calc_error(res.jac, np.sum(res.fun**2), *res.jac.shape) for v, var_pars_k in zip(data, data_pars): self.make_results(v, res.x, var, var_pars_k, res.jac.shape, err=err, corr=corr, partial_corr=partial_corr) @@ -458,9 +458,18 @@ class FitRoutine(object): self.make_results(v, res.beta, var, var_pars_k, (sum(len(d) for d in data), len(p0)), err=res.sd_beta, corr=corr, partial_corr=partial_corr) - def make_results(self, data, p, var_pars, used_pars, shape, - err=None, corr=None, partial_corr=None): - + def make_results( + self, + data: Data, + p: list[float], + var_pars: list[str], + used_pars: list[str], + shape: tuple[int, int], + err: list[float] = None, + corr: np.ndarray = None, + partial_corr: np.ndarray = None, + ): + print(data.complex_type) if err is None: err = [0] * len(p) @@ -498,52 +507,53 @@ class FitRoutine(object): model = data.get_model() self.result[idx] = FitResultCreator.make_with_model( - model, - data.x, - data.y, - actual_parameters, - data.fun_kwargs, - data.we_string, - data.idx, - *shape, + model=model, + x_orig=data.x, + y_orig=data.y, + p=actual_parameters, + fun_kwargs=data.fun_kwargs, + we=data.we_string, + idx=data.idx, + nobs=shape[0], + nvar=shape[1], corr=actual_corr, pcorr=actual_pcorr, + data_mode=data.complex_type, ) return self.result - @staticmethod - def _calc_error(jac, chi, nobs, nvars): - # copy of scipy.curve_fit to calculate covariance - # noinspection PyTupleAssignmentBalance - try: - _, s, vt = la.svd(jac, full_matrices=False) - except ValueError as e: - # this may be issue #39: On entry to DGESSD parameter had an illegal value - # catch this exception and ignore error calculation - logger.error(f'Error calculation failed with {e.args}') - pcov = None - else: - threshold = EPS * max(jac.shape) * s[0] - s = s[s > threshold] - vt = vt[:s.size] - pcov = np.dot(vt.T / s**2, vt) * chi / (nobs - nvars) +def _calc_error(jac, chi, nobs, nvars): + # copy of scipy.curve_fit to calculate covariance + # noinspection PyTupleAssignmentBalance + try: + _, s, vt = la.svd(jac, full_matrices=False) + except ValueError as e: + # this may be issue #39: On entry to DGESSD parameter had an illegal value + # catch this exception and ignore error calculation + logger.error(f'Error calculation failed with {e.args}') + pcov = None + else: + threshold = EPS * max(jac.shape) * s[0] + s = s[s > threshold] + vt = vt[:s.size] + pcov = np.dot(vt.T / s**2, vt) * chi / (nobs - nvars) - if pcov is None: - _err = np.zeros(nvars) - corr = np.zeros((nvars, nvars)) - else: - _err = np.sqrt(np.diag(pcov)) - corr = pcov / (_err[:, None] * _err[None, :]) + if pcov is None: + _err = np.zeros(nvars) + corr = np.zeros((nvars, nvars)) + else: + _err = np.sqrt(np.diag(pcov)) + corr = pcov / (_err[:, None] * _err[None, :]) - corr = corr.astype(np.float64) - try: - corr_inv = np.linalg.inv(corr) - corr_inv_diag = np.diag(np.sqrt(1 / np.diag(corr_inv))) - partial_corr = -1. * np.dot(np.dot(corr_inv_diag, corr_inv), corr_inv_diag) # Partial correlation matrix - partial_corr[np.diag_indices_from(partial_corr)] = 1. - except np.linalg.LinAlgError: - partial_corr = corr + corr = corr.astype(np.float64) + try: + corr_inv = np.linalg.inv(corr) + corr_inv_diag = np.diag(np.sqrt(1 / np.diag(corr_inv))) + partial_corr = -1. * np.dot(np.dot(corr_inv_diag, corr_inv), corr_inv_diag) # Partial correlation matrix + partial_corr[np.diag_indices_from(partial_corr)] = 1. + except np.linalg.LinAlgError: + partial_corr = corr - return _err, corr, partial_corr + return _err, corr, partial_corr diff --git a/src/nmreval/fit/model.py b/src/nmreval/fit/model.py index c80a2a3..46aacee 100644 --- a/src/nmreval/fit/model.py +++ b/src/nmreval/fit/model.py @@ -9,7 +9,7 @@ from ._meta import MultiModel from .parameter import Parameters, Parameter -class Model(object): +class Model: def __init__(self, model, *args, **kwargs): self.idx = kwargs.pop('idx', None) diff --git a/src/nmreval/fit/result.py b/src/nmreval/fit/result.py index 463a8eb..e380421 100644 --- a/src/nmreval/fit/result.py +++ b/src/nmreval/fit/result.py @@ -11,6 +11,7 @@ from scipy.stats import f as fdist from scipy.interpolate import interp1d from ._meta import MultiModel +from .model import Model from .parameter import Parameter from ..data.points import Points from ..data.signals import Signal @@ -36,17 +37,30 @@ class FitResultCreator: else: resid = kwargs['y'] - y_orig - stats = FitResultCreator.calc_statistics(resid, _y) + stats = calc_statistics(resid, _y) - return FitResult(kwargs['x'], kwargs['y'], x_orig, y_orig, params, dict(kwargs['choice']), resid, 0, 0, - kwargs['name'], stats, idx) + return FitResult( + x=kwargs['x'], + y=kwargs['y'], + x_data=x_orig, + y_data=y_orig, + params=params, + fun_kwargs=dict(kwargs['choice']), + resid=resid, + nobs=0, + nvar=0, + we='', + name=kwargs['name'], + stats=stats, + idx=idx, + ) @staticmethod def make_with_model( model: 'Model', x_orig: np.ndarray, y_orig: np.ndarray, - p: 'Parameters', + p: list, fun_kwargs: dict, we: str, idx: str | None, @@ -54,7 +68,7 @@ class FitResultCreator: nvar: int, corr: np.ndarray, pcorr: np.ndarray, - data_mode: int = 1, + data_mode: int, ) -> FitResult: if np.all(x_orig > 0) and (np.max(x_orig) > 100 * np.min(x_orig)): islog = True @@ -84,20 +98,12 @@ class FitResultCreator: actual_mode = fun_kwargs['complex_mode'] fun_kwargs['complex_mode'] = 0 - _y = model.func(p_final, _x, **fun_kwargs) + _y = check_complex(model.func(p_final, _x, **fun_kwargs), actual_mode, data_mode) + print(_y) - if not actual_mode < 0: - if actual_mode == 1: - _Y = _y.real - # _y.imag = 0 - elif actual_mode == 2: - # if data_mode == - _y = _y.imag - # _y.real = 0 + fun_kwargs['complex_mode'] = actual_mode - fun_kwargs['complex_mode'] = actual_mode - - stats = FitResultCreator.calc_statistics(_y, resid, nobs, nvar) + stats = calc_statistics(_y, resid, nobs, nvar) varied = [p.var for p in parameters.values()] if corr is None: @@ -138,38 +144,9 @@ class FitResultCreator: pcorr=partial_correlation, islog=islog, func=model, + data_complex=data_mode, ) - @staticmethod - def calc_statistics(y, residual, nobs=None, nvar=None): - chi = (residual**2).sum() - try: - r = 1 - chi/((y-np.mean(y))**2).sum() - except RuntimeWarning: - r = -9999 - - if nobs is None: - nobs = 1 - - if nvar is None: - nvar = 0 - - dof = nobs - nvar - loglikehood = nobs * np.log(chi / nobs) - - stats = { - 'chi^2': chi, - 'R^2': r, - 'AIC': loglikehood + 2 * nvar, - 'BIC': loglikehood + np.log(nobs) * nvar, - 'adj. R^2': 1 - (nobs-1) / (dof+1e-13) * (1-r), - 'red. chi^2': chi / (dof + 1e-13), - } - - stats['AICc'] = stats['AIC'] + 2*(nvar+1)*nvar / (dof - 1 + 1e-13) - - return stats - class FitResult(Points): @@ -192,7 +169,8 @@ class FitResult(Points): pcorr: np.ndarray = None, islog: bool = False, func=None, - **kwargs + data_complex: int = 1, + **kwargs, ): self.parameter, name = self._prepare_names(params, name) @@ -214,6 +192,7 @@ class FitResult(Points): self.y_data = y_data self._model_name = name self._func = func + self._data_complex = data_complex @staticmethod def _prepare_names(parameter: dict, modelname: str): @@ -422,20 +401,9 @@ class FitResult(Points): if self.func is None: raise ValueError('no fit function available to calculate new y values') - actual_mode = -1 - if 'complex_mode' in self.fun_kwargs: - actual_mode = self.fun_kwargs['complex_mode'] - self.fun_kwargs['complex_mode'] = 0 - new_fit = self.copy() y_values = self.func.func(self.p_final, x_values, **self.fun_kwargs) - if not actual_mode < 0: - if actual_mode == 1: - y_values.imag = 0 - elif actual_mode == 2: - y_values.real = 0 - - self.fun_kwargs['complex_mode'] = actual_mode + y_values = check_complex(y_values, self.fun_kwargs.get('complex_mode', -1), self._data_complex) new_fit.set_data(x_values, y_values, y_err=0.0) @@ -446,20 +414,13 @@ class FitResult(Points): raise ValueError('no fit function available to calculate new y values') part_functions = [] - actual_mode = -1 - if 'complex_mode' in self.fun_kwargs: - actual_mode = self.fun_kwargs['complex_mode'] - self.fun_kwargs['complex_mode'] = 0 + actual_mode = self.fun_kwargs.get('complex_mode', -1) for sub_name, sub_y in zip(self.func.sub_name(), self.func.sub(self.p_final, x_values, **self.fun_kwargs)): - if not actual_mode < 0: - if actual_mode == 1: - sub_y.imag = 0 - elif actual_mode == 2: - sub_y.real = 0 + sub_y = check_complex(sub_y, actual_mode, self._data_complex) + if np.iscomplexobj(sub_y): part_functions.append(Signal(x_values, sub_y, name=sub_name)) - else: part_functions.append(Points(x_values, sub_y, name=sub_name)) @@ -467,3 +428,49 @@ class FitResult(Points): self.fun_kwargs['complex_mode'] = actual_mode return part_functions + + +def check_complex(y, model_complex, data_complex): + if not np.iscomplexobj(y): + return y + + if model_complex == 1: + y.imag = 0 + if data_complex == 1: + y = y.real + elif model_complex == 2: + y.real = 0 + if data_complex == 1: + y = y.imag + + return y + + +def calc_statistics(y, residual, nobs=None, nvar=None): + chi = (residual**2).sum() + try: + r = 1 - chi/((y-np.mean(y))**2).sum() + except RuntimeWarning: + r = -9999 + + if nobs is None: + nobs = 1 + + if nvar is None: + nvar = 0 + + dof = nobs - nvar + loglikehood = nobs * np.log(chi / nobs) + + stats = { + 'chi^2': chi, + 'R^2': r, + 'AIC': loglikehood + 2 * nvar, + 'BIC': loglikehood + np.log(nobs) * nvar, + 'adj. R^2': 1 - (nobs-1) / (dof+1e-13) * (1-r), + 'red. chi^2': chi / (dof + 1e-13), + } + + stats['AICc'] = stats['AIC'] + 2*(nvar+1)*nvar / (dof - 1 + 1e-13) + + return stats