Some tweaks to curve_fitting reporting

Author: gb119

Stoner/analysis/fitting/classes.py

@@ -172,7 +172,7 @@ class _Curve_Fit_Result:
     as a class to make handling easier.
     """
 
-    def __init__(self, popt=None, pcov=None, infodict=None, mesg=None, ier=None):
+    def __init__(self, infodict=None, mesg=None, ier=None, results=None):
         """Store the results of the curve fit full_output fit.
 
         Args:
@@ -187,30 +187,21 @@ def __init__(self, popt=None, pcov=None, infodict=None, mesg=None, ier=None):
             ier (int):
                 Numerical error message.
         """
-        self.popt = popt
-        self.pcov = pcov
-        if pcov is not None:
-            self.perr = np.sqrt(np.diag(pcov))
-        else:
-            self.perr = None
-        self.mesg = mesg
-        self.ier = ier
-        self.nfev = None
-        self.fvec = None
-        self.fjac = None
-        self.ipvt = None
-        self.qtf = None
-        self.func = None
+        self._mapping = {}
+        self.func = lambda *args: None
+        self.f_name = None
+        self.args = []
+        self.kwargs = {}
         self.p0 = None
         self._residual_vals = None
-        self.chisq = None
-        self.nfree = None
         self.f_name = None
-        self._infodict = infodict
-        if infodict:
-            for k in infodict:
-                setattr(self, k, infodict[k])
-        self.xdata = None
+        self._infodict = {}
+        self._results = {}
+        self.infodict = {k: None for k in ["mfev", "fvec", "fjac", "ipvt", "qtf"]} if infodict is None else infodict
+        self.results = (
+            {k: None for k in ["pop", "perr", "pcov", "mesg", "ier", "chisq", "nfree"]} if results is None else results
+        )
+        self.data = None
         self.labels = None
         self.units = None
 
@@ -250,14 +241,28 @@ def row(self):
     @property
     def infodict(self):
         """Wrapper infodict."""
-        return getattr(self, "_infodict", None)
+        return self._infodict
 
     @infodict.setter
     def infodict(self, value):
         """Wrapper for setting infodict and subkeys."""
+        self._infodict = value
         if isinstance(value, dict):
             for k in value:
-                setattr(self, k, value[k])
+                self._mapping[k] = "_infodict"
+
+    @property
+    def results(self):
+        """Wrapper for a results dictionary."""
+        return self._results
+
+    @results.setter
+    def results(self, value):
+        """Wrapper for setting results dictionary and nting keys."""
+        self._results = value
+        if isinstance(value, dict):
+            for k in value:
+                self._mapping[k] = "_results"
 
     @property
     def fit(self):
@@ -267,22 +272,18 @@ def fit(self):
     @property
     def fit_values(self):
         """Return the fit values if x data is set."""
-        if self.xdata is None or self.func is None or self.popt is None:
+        if self.data is None or self.func is None or self.popt is None:
             raise ValueError(
                 "Need to have some x-data, the fitting functions and optimal parameters  before calculating fit"
             )
-        return self.func(self.xdata, *self.popt)
+        return self.func(self.data.data[:, self.settings.columns.xcol], *self.popt)
 
     @property
-    def data(self):
-        """Return the data that was fitted."""
-        self._data = getattr(self, "_data", np.array([]))
-        return self._data
-
-    @data.setter
-    def data(self, data):
-        """Return the data that was fitted."""
-        self._data = data
+    def perr(self):
+        """Return standar error from covariance matrix."""
+        if "perr" in self.results and self.results["perr"] is not None:
+            return self.results["perr"]
+        return np.sqrt(np.diag(self.pcov))
 
     @property
     def report(self):
@@ -336,6 +337,27 @@ def params(self):
         """List the parameter class objects."""
         return get_func_params(self.func)
 
+    def __dir__(self):
+        """Extend the attribute directory."""
+        return super().__dir__() + list(self._mapping.keys())
+
+    def __getattr__(self, name):
+        """Defer to using things we're tracking in the mapping."""
+        try:
+            return super().__getattr__(name)
+        except AttributeError:
+            pass
+        if name in self._mapping:
+            return getattr(self, self._mapping[name]).get(name, None)
+        raise AttributeError(f"{name} is not an attribute or {self.__class__.__name__}.")
+
+    def __setattr__(self, name, value):
+        """Pass through if an atrbute is already set."""
+        if name != "_mapping" and name in self._mapping:
+            getattr(self, self._mapping[name])[name] = value
+            return
+        super().__setattr__(name, value)
+
     def fit_report(self):
         """Create a Fit report like lmfit does."""
         template = f"""[[ Model ]]

Stoner/analysis/fitting/functions.py

@@ -241,12 +241,12 @@ def _func(x, *beta):
     return result_obj
 
 
-def _normalise_fit_result(datafile, xcol, ycol, fit, result_obj):
+def _normalise_fit_result(datafile, settings, fit, result_obj):
     """Normalise the fit results based on the fit instance."""
     func = result_obj.func
     args = result_obj.args
-    result_obj.data = datafile.data[:, ycol]
-    result_obj.xdata = datafile.data[:, xcol]
+    result_obj.data = datafile
+    result_obj.settings = settings
 
     match fit:
         case _Curve_Fit_Result():
@@ -273,23 +273,19 @@ def _normalise_fit_result(datafile, xcol, ycol, fit, result_obj):
             perr = fit.perr
             nfev = fit.nfev
             nfree = len(datafile) - len(popt)
-            fit_data = func(result_obj.xdata, *popt)
-            chisq = np.sum((result_obj.data - fit_data) ** 2) / nfree
+            fit_data = func(datafile // settings.columns.xcol, *popt)
+            chisq = np.sum((datafile.data[:, settings.columns.ycol] - fit_data) ** 2) / nfree
         case _:
             raise RuntimeError("Unable to understand {type(fit)} as a fitting result")
-    result_obj.popt = popt
-    result_obj.perr = perr
-    result_obj.nfev = nfev
-    result_obj.chisq = chisq
-    result_obj.nfree = nfree
+    result_obj.results = {"popt": popt, "perr": perr, "nfev": nfev, "chisq": chisq, "nfree": nfree}
     return result_obj
 
 
 def _record_curve_fit_result(datafile, func, fit, settings):
     """Annotate the DataFile object with the curve_fit result."""
     result_obj = _Curve_Fit_Result()
     result_obj = _normalise_model_func(func, settings.prefix, result_obj)
-    result_obj = _normalise_fit_result(datafile, settings.columns.xcol, settings.columns.ycol, fit, result_obj)
+    result_obj = _normalise_fit_result(datafile, settings, fit, result_obj)
 
     result_obj.add_metadata(datafile)
 
@@ -721,14 +717,15 @@ def p0_func(ydata,x=xdata):
             sigma = None
         for var in ["xcol", "ycol", "zcol", "xerr", "yerr", "zerr", "scale_covar"]:
             kwargs.pop(var, None)
-        report = _Curve_Fit_Result(*_curve_fit(_func, xdat, ydat, **kwargs))
+        popt, pcov, infodict, msg, ier = _curve_fit(_func, xdat, ydat, **kwargs)
+        report = _Curve_Fit_Result()
+        report.settings = settings
         report.func = func
-        if p0 is None:
-            report.p0 = np.ones(len(report.popt))
-        else:
-            report.p0 = p0
+        report.results = {"popt": popt, "pcov": pcov, "mesg": msg, "ier": ier}
+        report.infodict = infodict
+        report.p0 = np.ones(len(report.popt)) if p0 is None else p0
         report.data = datafile
-        report.residual_vals = data.y - report.fvec
+        report.residual_vals = ydat - report.fvec
         report.chisq = (report.residual_vals**2).sum()
         report.nfree = len(datafile) - len(report.popt)
         report.chisq /= report.nfree
@@ -860,11 +857,15 @@ def differential_evolution(datafile, model, xcol=None, ycol=None, p0=None, sigma
     kwargs.pop("polish", None)
     kwargs["full_output"] = True
     kwargs["absolute_sigma"] = abs_sigma
-    polish = _Curve_Fit_Result(*_curve_fit(model.func, data.x, data.y[0], sigma=data.e[0], p0=fit.x, **kwargs))
+    popt, pcov, infodict, mesg, ier = _curve_fit(model.func, data.x, data.y[0], sigma=data.e[0], p0=fit.x, **kwargs)
+    polish = _Curve_Fit_Result()
+    polish.results = {"popt": popt, "pcov": pcov, "mesg": mesg, "ier": ier}
+    polish.infodict = infodict
+    polish.data = datafile
+    polish.settings = settings
 
     polish.func = model.func
     polish.p0 = p0
-    polish.data = datafile
     polish.residual_vals = data.y - polish.fvec
     polish.chisq = (polish.residual_vals**2).sum()
     polish.nfree = len(datafile) - len(polish.popt)

scripts/PCAR-chi^2.ini

@@ -21,6 +21,8 @@ rescale_v: False
 # Warning, this may go badly wrong with certain DC waveforms
 remove_offset: False
 # Be clever about annotating the result on the plots
+# Decomposes the data into symmetric and antisymmetric parts before fitting
+decompose: True
 fancy_result: True
 #
 # Can switch between a least-squares fitting algorithm based on the lmfit module, or othogonal distance regression "odr"
@@ -76,9 +78,9 @@ units: meV
 value: 0.5
 vary: False
 min: 0.2
-max: 0.7
+max: 0.54
 label: P
-step: 0.05
+step: 0.025
 
 [Z]
 value: 0.17