IPKM/nmreval
14
2
Fork 1
Code Issues 39 Pull Requests Actions Packages Projects Releases Wiki Activity

fit results are not complex for non-complex data

Browse Source
This commit is contained in:
Dominik Demuth 2024-02-01 19:09:17 +01:00
parent 4b6820af18
commit 7aecf962ba
3 changed files with 132 additions and 115 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

98
src/nmreval/fit/minimizer.py
View File

@ -361,7 +361,7 @@ class FitRoutine(object):
with np.errstate(all='ignore'): with np.errstate(all='ignore'):
res = optimize.least_squares(cost, p0, bounds=(lb, ub), max_nfev=500 * len(p0)) 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, self.make_results(data, res.x, var, data.para_keys, res.jac.shape,
err=err, corr=corr, partial_corr=partial_corr) err=err, corr=corr, partial_corr=partial_corr)
@ -375,7 +375,7 @@ class FitRoutine(object):
with np.errstate(all='ignore'): with np.errstate(all='ignore'):
res = optimize.least_squares(cost, p0, bounds=(lb, ub), max_nfev=500 * len(p0)) 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): for v, var_pars_k in zip(data, data_pars):
self.make_results(v, res.x, var, var_pars_k, res.jac.shape, self.make_results(v, res.x, var, var_pars_k, res.jac.shape,
err=err, corr=corr, partial_corr=partial_corr) 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)), 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) err=res.sd_beta, corr=corr, partial_corr=partial_corr)
def make_results(self, data, p, var_pars, used_pars, shape, def make_results(
err=None, corr=None, partial_corr=None): 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: if err is None:
err = [0] * len(p) err = [0] * len(p)
@ -498,52 +507,53 @@ class FitRoutine(object):
model = data.get_model() model = data.get_model()
self.result[idx] = FitResultCreator.make_with_model( self.result[idx] = FitResultCreator.make_with_model(
model, model=model,
data.x, x_orig=data.x,
data.y, y_orig=data.y,
actual_parameters, p=actual_parameters,
data.fun_kwargs, fun_kwargs=data.fun_kwargs,
data.we_string, we=data.we_string,
data.idx, idx=data.idx,
*shape, nobs=shape[0],
nvar=shape[1],
corr=actual_corr, corr=actual_corr,
pcorr=actual_pcorr, pcorr=actual_pcorr,
data_mode=data.complex_type,
) )
return self.result return self.result
@staticmethod def _calc_error(jac, chi, nobs, nvars):
def _calc_error(jac, chi, nobs, nvars): # copy of scipy.curve_fit to calculate covariance
# copy of scipy.curve_fit to calculate covariance # noinspection PyTupleAssignmentBalance
# noinspection PyTupleAssignmentBalance try:
try: _, s, vt = la.svd(jac, full_matrices=False)
_, s, vt = la.svd(jac, full_matrices=False) except ValueError as e:
except ValueError as e: # this may be issue #39: On entry to DGESSD parameter had an illegal value
# this may be issue #39: On entry to DGESSD parameter had an illegal value # catch this exception and ignore error calculation
# catch this exception and ignore error calculation logger.error(f'Error calculation failed with {e.args}')
logger.error(f'Error calculation failed with {e.args}') pcov = None
pcov = None else:
else: threshold = EPS * max(jac.shape) * s[0]
threshold = EPS * max(jac.shape) * s[0] s = s[s > threshold]
s = s[s > threshold] vt = vt[:s.size]
vt = vt[:s.size] pcov = np.dot(vt.T / s**2, vt) * chi / (nobs - nvars)
pcov = np.dot(vt.T / s**2, vt) * chi / (nobs - nvars)
if pcov is None: if pcov is None:
_err = np.zeros(nvars) _err = np.zeros(nvars)
corr = np.zeros((nvars, nvars)) corr = np.zeros((nvars, nvars))
else: else:
_err = np.sqrt(np.diag(pcov)) _err = np.sqrt(np.diag(pcov))
corr = pcov / (_err[:, None] * _err[None, :]) corr = pcov / (_err[:, None] * _err[None, :])
corr = corr.astype(np.float64) corr = corr.astype(np.float64)
try: try:
corr_inv = np.linalg.inv(corr) corr_inv = np.linalg.inv(corr)
corr_inv_diag = np.diag(np.sqrt(1 / np.diag(corr_inv))) 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 = -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. partial_corr[np.diag_indices_from(partial_corr)] = 1.
except np.linalg.LinAlgError: except np.linalg.LinAlgError:
partial_corr = corr partial_corr = corr
return _err, corr, partial_corr return _err, corr, partial_corr

2
src/nmreval/fit/model.py
View File

@ -9,7 +9,7 @@ from ._meta import MultiModel
from .parameter import Parameters, Parameter from .parameter import Parameters, Parameter
class Model(object): class Model:
def __init__(self, model, *args, **kwargs): def __init__(self, model, *args, **kwargs):
self.idx = kwargs.pop('idx', None) self.idx = kwargs.pop('idx', None)

147
src/nmreval/fit/result.py
View File

@ -11,6 +11,7 @@ from scipy.stats import f as fdist
from scipy.interpolate import interp1d from scipy.interpolate import interp1d
from ._meta import MultiModel from ._meta import MultiModel
from .model import Model
from .parameter import Parameter from .parameter import Parameter
from ..data.points import Points from ..data.points import Points
from ..data.signals import Signal from ..data.signals import Signal
@ -36,17 +37,30 @@ class FitResultCreator:
else: else:
resid = kwargs['y'] - y_orig 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, return FitResult(
kwargs['name'], stats, idx) 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 @staticmethod
def make_with_model( def make_with_model(
model: 'Model', model: 'Model',
x_orig: np.ndarray, x_orig: np.ndarray,
y_orig: np.ndarray, y_orig: np.ndarray,
p: 'Parameters', p: list,
fun_kwargs: dict, fun_kwargs: dict,
we: str, we: str,
idx: str | None, idx: str | None,
@ -54,7 +68,7 @@ class FitResultCreator:
nvar: int, nvar: int,
corr: np.ndarray, corr: np.ndarray,
pcorr: np.ndarray, pcorr: np.ndarray,
data_mode: int = 1, data_mode: int,
) -> FitResult: ) -> FitResult:
if np.all(x_orig > 0) and (np.max(x_orig) > 100 * np.min(x_orig)): if np.all(x_orig > 0) and (np.max(x_orig) > 100 * np.min(x_orig)):
islog = True islog = True
@ -84,20 +98,12 @@ class FitResultCreator:
actual_mode = fun_kwargs['complex_mode'] actual_mode = fun_kwargs['complex_mode']
fun_kwargs['complex_mode'] = 0 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: fun_kwargs['complex_mode'] = actual_mode
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 stats = calc_statistics(_y, resid, nobs, nvar)
stats = FitResultCreator.calc_statistics(_y, resid, nobs, nvar)
varied = [p.var for p in parameters.values()] varied = [p.var for p in parameters.values()]
if corr is None: if corr is None:
@ -138,38 +144,9 @@ class FitResultCreator:
pcorr=partial_correlation, pcorr=partial_correlation,
islog=islog, islog=islog,
func=model, 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): class FitResult(Points):
@ -192,7 +169,8 @@ class FitResult(Points):
pcorr: np.ndarray = None, pcorr: np.ndarray = None,
islog: bool = False, islog: bool = False,
func=None, func=None,
**kwargs data_complex: int = 1,
**kwargs,
): ):
self.parameter, name = self._prepare_names(params, name) self.parameter, name = self._prepare_names(params, name)
@ -214,6 +192,7 @@ class FitResult(Points):
self.y_data = y_data self.y_data = y_data
self._model_name = name self._model_name = name
self._func = func self._func = func
self._data_complex = data_complex
@staticmethod @staticmethod
def _prepare_names(parameter: dict, modelname: str): def _prepare_names(parameter: dict, modelname: str):
@ -422,20 +401,9 @@ class FitResult(Points):
if self.func is None: if self.func is None:
raise ValueError('no fit function available to calculate new y values') 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() new_fit = self.copy()
y_values = self.func.func(self.p_final, x_values, **self.fun_kwargs) y_values = self.func.func(self.p_final, x_values, **self.fun_kwargs)
if not actual_mode < 0: y_values = check_complex(y_values, self.fun_kwargs.get('complex_mode', -1), self._data_complex)
if actual_mode == 1:
y_values.imag = 0
elif actual_mode == 2:
y_values.real = 0
self.fun_kwargs['complex_mode'] = actual_mode
new_fit.set_data(x_values, y_values, y_err=0.0) 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') raise ValueError('no fit function available to calculate new y values')
part_functions = [] part_functions = []
actual_mode = -1 actual_mode = self.fun_kwargs.get('complex_mode', -1)
if 'complex_mode' in self.fun_kwargs:
actual_mode = self.fun_kwargs['complex_mode']
self.fun_kwargs['complex_mode'] = 0
for sub_name, sub_y in zip(self.func.sub_name(), self.func.sub(self.p_final, x_values, **self.fun_kwargs)): 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: sub_y = check_complex(sub_y, actual_mode, self._data_complex)
if actual_mode == 1:
sub_y.imag = 0
elif actual_mode == 2:
sub_y.real = 0
if np.iscomplexobj(sub_y):
part_functions.append(Signal(x_values, sub_y, name=sub_name)) part_functions.append(Signal(x_values, sub_y, name=sub_name))
else: else:
part_functions.append(Points(x_values, sub_y, name=sub_name)) 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 self.fun_kwargs['complex_mode'] = actual_mode
return part_functions 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