python/sbp_fit.py: report results in JSON

1 files changed, 66 insertions, 44 deletions

diff --git a/python/sbp_fit.py b/python/sbp_fit.py
index 0b7d29a..0ce8d4d 100644
--- a/python/sbp_fit.py
+++ b/python/sbp_fit.py
@@ -3,9 +3,13 @@
 #
 # Aaron LI
 # Created: 2016-03-13
-# Updated: 2016-03-14
+# Updated: 2016-03-31
 #
 # Changelogs:
+# 2016-03-31:
+#   * Remove `ci_report()'
+#   * Add `make_results()' to orgnize all results as s Python dictionary
+#   * Report results as json string
 # 2016-03-28:
 #   * Add `main()', `make_model()'
 #   * Use `configobj' to handle configurations
@@ -25,6 +29,7 @@ or the double-beta model:
   s(r) = s01 * [1.0 + (r/rc1)^2] ^ (0.5-3*beta1) +
          s02 * [1.0 + (r/rc2)^2] ^ (0.5-3*beta2) + bkg
 
+
 Sample config file:
 -------------------------------------------------
 name    = <NAME>
@@ -61,8 +66,8 @@ bkg   = 1.0e-9, 0.0,  1.0e-7
 -------------------------------------------------
 """
 
-__version__ = "0.4.0"
-__date__    = "2016-03-28"
+__version__ = "0.5.0"
+__date__    = "2016-03-31"
 
 import numpy as np
 import lmfit
@@ -76,6 +81,8 @@ import os
 import sys
 import re
 import argparse
+import json
+from collections import OrderedDict
 
 
 plt.style.use("ggplot")
@@ -456,54 +463,69 @@ class SbpFit:
         if not hasattr(self, "f_residual") or self.f_residual is None:
             self.set_residual()
         self.fitter = lmfit.Minimizer(self.f_residual, self.model.params)
-        self.results = self.fitter.minimize(method=method)
-        self.fitted_model = lambda x: self.model.func(x, self.results.params)
+        self.fitted = self.fitter.minimize(method=method)
+        self.fitted_model = lambda x: self.model.func(x, self.fitted.params)
 
     def calc_ci(self, sigmas=[0.68, 0.90]):
         # `conf_interval' requires the fitted results have valid `stderr',
         # so we need to re-fit the model with method `leastsq'.
-        results = self.fitter.minimize(method="leastsq",
-                params=self.results.params)
-        self.ci, self.trace = lmfit.conf_interval(self.fitter, results,
+        fitted = self.fitter.minimize(method="leastsq",
+                params=self.fitted.params)
+        self.ci, self.trace = lmfit.conf_interval(self.fitter, fitted,
                 sigmas=sigmas, trace=True)
 
-    @staticmethod
-    def ci_report(ci, with_offset=True):
+    def make_results(self):
         """
-        Format and return the report for confidence intervals.
-        This function is intended to replace the `lmfit.ci_report()`
-        which do not use scitenfic notation for small values.
+        Make the `self.results' dictionary which contains all the fitting
+        results as well as the confidence intervals.
         """
-        pnames = list(ci.keys())
-        plen = max(map(len, pnames))
-        report = []
-        # header
-        sigmas = [ v[0] for v in ci.get(pnames[0]) if v[0] > 0.0 ]
-        fmt = "{:<%d}" % (plen+1) + \
-                "  {:11.2%}  {:11.2%}     _BEST_    {:11.2%}  {:11.2%}"
-        report.append(fmt.format(" ", *sigmas))
-        # parameters ci
-        for par in pnames:
-            values = [ v[1] for v in ci.get(par) ]
-            if with_offset:
-                N = len(values)
-                offsets = [ values[N//2] for i in range(N) ]
-                offsets[N//2] = 0.0
-                values = list(map(lambda a,b: a-b, values, offsets))
-            fmt = "{:<%d}" % (plen+1) + \
-                    "  {:+11.5g}  {:+11.5g}  {:11.5g}  {:+11.5g}  {:+11.5g}"
-            report.append(fmt.format(par+":", *values))
-        return "\n".join(report)
+        fitted = self.fitted
+        self.results = OrderedDict()
+        ## fitting results
+        self.results.update(
+                nfev   = fitted.nfev,
+                ndata  = fitted.ndata,
+                nvarys = fitted.nvarys,  # number of varible paramters
+                nfree  = fitted.nfree,  # degree of freem
+                chisqr = fitted.chisqr,
+                redchi = fitted.redchi,
+                aic    = fitted.aic,
+                bic    = fitted.bic)
+        params = fitted.params
+        pnames = list(params.keys())
+        pvalues = OrderedDict()
+        for pn in pnames:
+            par = params.get(pn)
+            pvalues[pn] = [par.value, par.min, par.max, par.vary]
+        self.results["params"] = pvalues
+        ## confidence intervals
+        if hasattr(self, "ci") and self.ci is not None:
+            ci = self.ci
+            ci_values = OrderedDict()
+            ci_sigmas = [ "ci%02d" % (v[0]*100) for v in ci.get(pnames[0]) ]
+            ci_names = sorted(list(set(ci_sigmas)))
+            ci_idx = { k: [] for k in ci_names }
+            for cn, idx in zip(ci_sigmas, range(len(ci_sigmas))):
+                ci_idx[cn].append(idx)
+            # parameters ci
+            for pn in pnames:
+                ci_pv = OrderedDict()
+                pv = [ v[1] for v in ci.get(pn) ]
+                # best
+                pv_best = pv[ ci_idx["ci00"][0] ]
+                ci_pv["best"] = pv_best
+                # ci of each sigma
+                pv2 = [ v-pv_best for v in pv ]
+                for cn in ci_names[1:]:
+                    ci_pv[cn] = [ pv2[idx] for idx in ci_idx[cn] ]
+                ci_values[pn] = ci_pv
+            self.results["ci"] = ci_values
 
     def report(self, outfile=sys.stdout):
-        if hasattr(self, "results") and self.results is not None:
-            p = lmfit.fit_report(self.results)
-            print(p, file=outfile)
-        if hasattr(self, "ci") and self.ci is not None:
-            print("-----------------------------------------------",
-                    file=outfile)
-            p = self.ci_report(self.ci)
-            print(p, file=outfile)
+        if not hasattr(self, "results") or self.results is None:
+            self.make_results()
+        jd = json.dumps(self.results, indent=2)
+        print(jd, file=outfile)
 
     def plot(self, ax=None, fig=None):
         if ax is None:
@@ -525,7 +547,7 @@ class SbpFit:
         ypred = self.fitted_model(xpred)
         ymin = min(min(self.ydata), min(ypred))
         ymax = max(max(self.ydata), max(ypred))
-        self.model.plot(params=self.results.params, xdata=xpred, ax=ax)
+        self.model.plot(params=self.fitted.params, xdata=xpred, ax=ax)
         ax.set_xscale("log")
         ax.set_yscale("log")
         ax.set_xlim(1.0, xmax)
@@ -537,8 +559,8 @@ class SbpFit:
         ax.set_xlabel("Radius (pixel)")
         ax.set_ylabel(r"Surface Brightness (photons/cm$^2$/pixel$^2$/s)")
         ax.text(x=xmax, y=ymax,
-                s="redchi: %.2f / %.2f = %.2f" % (self.results.chisqr,
-                    self.results.nfree, self.results.chisqr/self.results.nfree),
+                s="redchi: %.2f / %.2f = %.2f" % (self.fitted.chisqr,
+                    self.fitted.nfree, self.fitted.chisqr/self.fitted.nfree),
                 horizontalalignment="right", verticalalignment="top")
         return (fig, ax)