207-noncomplex-fits (#244)

Co-authored-by: Dominik Demuth <dominik.demuth@physik.tu-darmstadt.de>
Reviewed-on: IPKM/nmreval#244

closes #207

src/gui_qt/main/management.py

@@ -503,11 +503,25 @@ class UpperManagement(QtCore.QObject):
                         we = we_option
 
                     if m_complex is None or m_complex == 1:
+                        # model is not complex: m_complex = None
+                        # model is complex, fit real part: m_complex = 1
                         _y = data_i.y.real
-                    elif m_complex == 2 and np.iscomplexobj(data_i.y):
+                        data_complex = 1
-                        _y = data_i.y.imag
+                    elif m_complex == 2:
+                        # model is complex, fit imag part: m_complex = 2
+                        if np.iscomplexobj(data_i.y):
+                            # data is complex, use imag part
+                            _y = data_i.y.imag
+                            data_complex = 2
+                        else:
+                            # data is real
+                            _y = data_i.y
+                            data_complex = 1
                     else:
+                        # model is complex, fit complex: m_complex = 0
+                        # use data as given (complex or not)
                         _y = data_i.y
+                        data_complex = 0
 
                     _x = data_i.x
 
@@ -524,9 +538,9 @@ class UpperManagement(QtCore.QObject):
 
                     try:
                         if isinstance(we, str):
-                            d = fit_d.Data(_x[inside], _y[inside], we=we, idx=set_id)
+                            d = fit_d.Data(_x[inside], _y[inside], we=we, idx=set_id, complex_type=data_complex)
                         else:
-                            d = fit_d.Data(_x[inside], _y[inside], we=we[inside], idx=set_id)
+                            d = fit_d.Data(_x[inside], _y[inside], we=we[inside], idx=set_id, complex_type=data_complex)
                     except Exception as e:
                         raise Exception(f'Setting data failed for {set_id}')
 

src/nmreval/fit/data.py

@@ -6,8 +6,8 @@ from .model import Model
 from .parameter import Parameters, Parameter
 
 
-class Data(object):
+class Data:
-    def __init__(self, x, y, we=None, idx=None):
+    def __init__(self, x, y, we=None, idx=None, complex_type: int = 0):
         self.x = np.asarray(x)
         self.y = np.asarray(y)
         if self.y.shape[0] != self.x.shape[0]:
@@ -20,6 +20,7 @@ class Data(object):
         self.parameter = Parameters()
         self.para_keys: list = []
         self.fun_kwargs = {}
+        self.complex_type = complex_type
 
     def __len__(self):
         return self.y.shape[0]

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,17 @@ 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,
+    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,
+    ):
         if err is None:
             err = [0] * len(p)
 
@@ -498,52 +506,54 @@ class FitRoutine(object):
         model = data.get_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,
             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)
 
-        if pcov is None:
+def _calc_error(jac, chi, nobs, nvars):
-            _err = np.zeros(nvars)
+    # copy of scipy.curve_fit to calculate covariance
-            corr = np.zeros((nvars, nvars))
+    # noinspection PyTupleAssignmentBalance
-        else:
+    try:
-            _err = np.sqrt(np.diag(pcov))
+        _, s, vt = la.svd(jac, full_matrices=False)
-            corr = pcov / (_err[:, None] * _err[None, :])
+    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)
 
-        corr = corr.astype(np.float64)
+    if pcov is None:
-        try:
+        _err = np.zeros(nvars)
-            corr_inv = np.linalg.inv(corr)
+        corr = np.zeros((nvars, nvars))
-            corr_inv_diag = np.diag(np.sqrt(1 / np.diag(corr_inv)))
+    else:
-            partial_corr = -1. * np.dot(np.dot(corr_inv_diag, corr_inv), corr_inv_diag)  # Partial correlation matrix
+        _err = np.sqrt(np.diag(pcov))
-            partial_corr[np.diag_indices_from(partial_corr)] = 1.
+        corr = pcov / (_err[:, None] * _err[None, :])
-        except np.linalg.LinAlgError:
-            partial_corr = corr
 
-        return _err, corr, partial_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
 

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)
 

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,
+        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
     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,6 +68,7 @@ class FitResultCreator:
             nvar: int,
             corr:  np.ndarray,
             pcorr: np.ndarray,
+            data_mode: int,
     ) -> FitResult:
         if np.all(x_orig > 0) and (np.max(x_orig) > 100 * np.min(x_orig)):
             islog = True
@@ -83,17 +98,11 @@ 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)
 
-        if not actual_mode < 0:
+        fun_kwargs['complex_mode'] = actual_mode
-            if actual_mode == 1:
-                _y.imag = 0
-            elif actual_mode == 2:
-                _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()]
 
         if corr is None:
@@ -134,38 +143,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):
 
@@ -188,7 +168,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)
@@ -210,6 +191,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):
@@ -418,20 +400,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:
+        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)
 
@@ -442,20 +413,13 @@ class FitResult(Points):
             raise ValueError('no fit function available to calculate new y values')
 
         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)):
-            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))
-
             else:
                 part_functions.append(Points(x_values, sub_y, name=sub_name))
 
@@ -463,3 +427,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