31.4 xi in #17: produce plot

pymule/xicut.py

 from loader import multiintersect, addsets, mergeset
 from scipy import optimize
 import scipy.stats
+import matplotlib.pyplot as plt
 import numpy as np
 from errortools import *
 
@@ -58,6 +59,20 @@ def mergefkswithplot(sets):
     psets = [addkeyedsets(i) for i in sets]
     ans = mergeset(addkeyedsets(psets).values())
 
+    fig, (ax1, ax2) = plt.subplots(1, 2)
+    ax1.set_xscale('log')
+    ax2.set_xscale('log')
+    ax1.set_xlabel('$\\xi_c$')
+    ax2.set_xlabel('$\\xi_c$')
+    ax1.set_ylabel('$\\sigma^{(%d)}_i$' % (n-1))
+    ax2.set_ylabel(
+        r'$\tfrac{\sigma^{(%d)}}{\langle\sigma^{(%d)}\rangle}-1$' % (n-1, n-1)
+    )
+    ax2.ticklabel_format(style='sci',axis='y')
+    ax2.yaxis.major.formatter.set_powerlimits([1,1])
+
+    xval = np.exp(np.linspace(-7,0))
+
     totcoeff = np.zeros(n)
     totcovar = np.zeros((n,n))
     pvalues = []
@@ -72,4 +87,45 @@ def mergefkswithplot(sets):
         totcoeff += coeff
         totcovar += covar
 
-    return ans
+        ax1.plot(
+            xval, get_val(xval, coeff),
+            c='C%d' % ind
+        )
+        ax1.errorbar(
+            values[:,0], values[:,1], values[:,2],
+            linestyle='none', marker='.',
+            c='C%d' % ind
+        )
+
+    norm = ans['value'][0]
+
+    tot = combineNplots(addplots,pvalues)
+    tot0 = [1,1/norm,1/norm] * tot - [0,1,0]
+
+    band = get_errorbands(xval, totcoeff, totcovar)
+    ax2.plot(band[:,0], band[:,1]/norm-1)
+    ax2.fill_between(
+        band[:,0],
+        band[:,2]/norm-1, band[:,3]/norm-1,
+        alpha=0.1, color='C0'
+    )
+    ax2.errorbar(
+        tot0[:,0], tot0[:,1], tot0[:,2],
+        linestyle='none', marker='.',
+        c='C0'
+    )
+
+    ax1.plot(band[:,0], band[:,1])
+    ax1.errorbar(
+        tot[:,0], tot[:,1], tot[:,2],
+        linestyle='none', marker='.',
+        c='C0'
+    )
+
+    fig.legend(
+        ['$\\sigma^{(%d)}_{n}$' % (n-1)] +
+        ['$\\sigma^{(%d)}_{n+%d}$' % (i, n-1) for i in range(1, n)] +
+        ['$\\sigma^{(%d)}$' % (n-1)]
+    )
+
+    return fig, ans