5: basic error tools in re #17

pymule/__init__.py

@@ -2,3 +2,5 @@ import numpy as np
 import matplotlib.pyplot as plt
 
 from vegas import importvegas
+from errortools import mergenumbers, plusnumbers, dividenumbers, timesnumbers,\
+                       mergeplots, addplots, combineNplots, scaleplot

pymule/errortools.py

+import numpy as np
+
+
+def mergenumbers(values, quiet=False):
+    if type(values) == list:
+        values = np.array(values)
+
+    weight = sum(1 / values[:,1]**2)
+    value = sum(values[:,0] / values[:,1]**2 / weight)
+    chi = 1./len(values) * sum(
+        ((values[:,0]-value) / values[:,1])**2
+    )
+    ans = np.array([value, np.sqrt(1/weight)])
+
+    if not quiet:
+        print chi
+        return ans
+    else:
+        return chi, ans
+
+
+def plusnumbers(*args):
+    if len(args) == 1:
+        lst = np.array(args[0])
+    else:
+        lst = np.array(args)
+    return np.array([
+        sum(lst[:,0]),
+        np.sqrt(sum(lst[:,1]**2))
+    ])
+
+
+def dividenumbers(a, b):
+    return np.array([
+        a[0]/b[0],
+        np.sqrt((b[1]**2*a[0]**2)/b[0]**4 + a[1]**2/b[0]**2)
+    ])
+
+
+def timesnumbers(a, b):
+    return np.array([
+        a[0]*b[0],
+        np.sqrt(b[0]**2*a[1]**2 + a[0]**2*b[1]**2)
+    ])
+
+
+def mergeplots(ps, returnchi=False):
+    if type(ps) == list:
+        ps = np.array(ps)
+
+    with np.errstate(divide='ignore', invalid='ignore'):
+        w = sum(1/ps[:,:,2]**2)
+        out = np.zeros(ps[0].shape)
+        out[:,0] = ps[0,:,0]
+
+        out[:,1] = sum(ps[:,:,1]/ps[:,:,2]**2 / w)
+        out[:,2] = np.sqrt(1/w)
+
+        # if one bin is empty, the zero causes trouble
+        out[~np.isfinite(out)] = 0
+
+        chi = 1./(len(ps)-1)*sum((ps[:,:,1]-out[:,1])**2/ps[:,:,2]**2)
+        chi = np.column_stack((out[:,0], chi))
+        chi[~np.isfinite(chi)] = 0
+
+    if returnchi:
+        return out, chi
+    else:
+        return out
+
+
+def addplots(a, b, sa=1., sb=1.):
+    # There must be a better way of doing this
+    maskA = [False]*len(a)
+    maskB = [False]*len(b)
+    for ia in range(len(a)):
+        for ib in range(len(b)):
+            if abs(1-a[ia, 0] / b[ib, 0]) < 1e-15:
+                maskA[ia] = True
+                maskB[ib] = True
+
+    x = a[maskA, 0]
+    y = sa*a[maskA, 1] + sb*b[maskB, 1]
+    e = np.sqrt(sa**2 * a[maskA, 2]**2 + sb**2 * b[maskB, 2]**2)
+
+    return np.column_stack((x, y, e))