update utils/ folder to make a python package called sddm

This commit adds an sddm python package to the utils/ folder. This allows me to
consolidate code used across all the various scripts. This package is now
installed by default to /home/tlatorre/local/lib/python2.7/site-packages so you
should add the following to your .bashrc file:

    export PYTHONPATH=$HOME/local/lib/python2.7/site-packages/:$PYTHONPATH

before using the scripts installed to ~/local/bin.

1 files changed, 44 insertions, 288 deletions

diff --git a/utils/plot-fit-results b/utils/plot-fit-results
index 18139ff..e154b21 100755
--- a/utils/plot-fit-results
+++ b/utils/plot-fit-results
@@ -15,265 +15,62 @@
 # this program. If not, see <https://www.gnu.org/licenses/>.
 
 from __future__ import print_function, division
-import numpy as np
-from scipy.stats import iqr, norm, beta
-from matplotlib.lines import Line2D
-import itertools
-
-IDP_E_MINUS  =    20
-IDP_MU_MINUS =    22
-
-SNOMAN_MASS = {
-    20: 0.511,
-    21: 0.511,
-    22: 105.658,
-    23: 105.658
-}
-
-def plot_hist(x, label=None):
-    # determine the bin width using the Freedman Diaconis rule
-    # see https://en.wikipedia.org/wiki/Freedman%E2%80%93Diaconis_rule
-    h = 2*iqr(x)/len(x)**(1/3)
-    n = max(int((np.max(x)-np.min(x))/h),10)
-    bins = np.linspace(np.min(x),np.max(x),n)
-    plt.hist(x, bins=bins, histtype='step', label=label)
-
-def plot_legend(n):
-    plt.figure(n)
-    ax = plt.gca()
-    handles, labels = ax.get_legend_handles_labels()
-    new_handles = [Line2D([],[],c=h.get_edgecolor()) for h in handles]
-    plt.legend(handles=new_handles,labels=labels)
-
-def get_stats(x):
-    """
-    Returns a tuple (mean, error mean, std, error std) for the values in x.
-
-    The formula for the standard error on the standard deviation comes from
-    https://stats.stackexchange.com/questions/156518.
-    """
-    mean = np.mean(x)
-    std = np.std(x)
-    n = len(x)
-    u4 = np.mean((x-mean)**4)
-    error = np.sqrt((u4-(n-3)*std**4/(n-1))/n)/(2*std)
-    return mean, std/np.sqrt(n), std, error
-
-def iqr_std_err(x):
-    """
-    Returns the approximate standard deviation assuming the central part of the
-    distribution is gaussian.
-    """
-    x = x.dropna()
-    n = len(x)
-    if n == 0:
-        return np.nan
-    # see https://stats.stackexchange.com/questions/110902/error-on-interquartile-range
-    std = iqr(x.values)/1.3489795
-    return 1.573*std/np.sqrt(n)
-
-def iqr_std(x):
-    """
-    Returns the approximate standard deviation assuming the central part of the
-    distribution is gaussian.
-    """
-    x = x.dropna()
-    n = len(x)
-    if n == 0:
-        return np.nan
-    return iqr(x.values)/1.3489795
-
-def quantile_error(x,q):
-    """
-    Returns the standard error for the qth quantile of `x`. The error is
-    computed using the Maritz-Jarrett method described here:
-    https://www.itl.nist.gov/div898/software/dataplot/refman2/auxillar/quantse.htm.
-    """
-    x = np.sort(x)
-    n = len(x)
-    m = int(q*n+0.5)
-    A = m - 1
-    B = n - m
-    i = np.arange(1,len(x)+1)
-    w = beta.cdf(i/n,A,B) - beta.cdf((i-1)/n,A,B)
-    return np.sqrt(np.sum(w*x**2)-np.sum(w*x)**2)
-
-def q90_err(x):
-    """
-    Returns the error on the 90th percentile for all the non NaN values in a
-    Series `x`.
-    """
-    x = x.dropna()
-    n = len(x)
-    if n == 0:
-        return np.nan
-    return quantile_error(x.values,0.9)
-
-def q90(x):
-    """
-    Returns the 90th percentile for all the non NaN values in a Series `x`.
-    """
-    x = x.dropna()
-    n = len(x)
-    if n == 0:
-        return np.nan
-    return np.percentile(x.values,90.0)
-
-def median(x):
-    """
-    Returns the median for all the non NaN values in a Series `x`.
-    """
-    x = x.dropna()
-    n = len(x)
-    if n == 0:
-        return np.nan
-    return np.median(x.values)
-
-def median_err(x):
-    """
-    Returns the approximate error on the median for all the non NaN values in a
-    Series `x`. The error on the median is approximated assuming the central
-    part of the distribution is gaussian.
-    """
-    x = x.dropna()
-    n = len(x)
-    if n == 0:
-        return np.nan
-    # First we estimate the standard deviation using the interquartile range.
-    # Here we are essentially assuming the central part of the distribution is
-    # gaussian.
-    std = iqr(x.values)/1.3489795
-    median = np.median(x.values)
-    # Now we estimate the error on the median for a gaussian
-    # See https://stats.stackexchange.com/questions/45124/central-limit-theorem-for-sample-medians.
-    return 1/(2*np.sqrt(n)*norm.pdf(median,median,std))
-
-def std_err(x):
-    x = x.dropna()
-    mean = np.mean(x)
-    std = np.std(x)
-    n = len(x)
-    if n == 0:
-        return np.nan
-    elif n == 1:
-        return 0.0
-    u4 = np.mean((x-mean)**4)
-    error = np.sqrt((u4-(n-3)*std**4/(n-1))/n)/(2*std)
-    return error
-
-# Taken from https://raw.githubusercontent.com/mwaskom/seaborn/c73055b2a9d9830c6fbbace07127c370389d04dd/seaborn/utils.py
-def despine(fig=None, ax=None, top=True, right=True, left=False,
-            bottom=False, offset=None, trim=False):
-    """Remove the top and right spines from plot(s).
-
-    fig : matplotlib figure, optional
-        Figure to despine all axes of, default uses current figure.
-    ax : matplotlib axes, optional
-        Specific axes object to despine.
-    top, right, left, bottom : boolean, optional
-        If True, remove that spine.
-    offset : int or dict, optional
-        Absolute distance, in points, spines should be moved away
-        from the axes (negative values move spines inward). A single value
-        applies to all spines; a dict can be used to set offset values per
-        side.
-    trim : bool, optional
-        If True, limit spines to the smallest and largest major tick
-        on each non-despined axis.
-
-    Returns
-    -------
-    None
-
-    """
-    # Get references to the axes we want
-    if fig is None and ax is None:
-        axes = plt.gcf().axes
-    elif fig is not None:
-        axes = fig.axes
-    elif ax is not None:
-        axes = [ax]
-
-    for ax_i in axes:
-        for side in ["top", "right", "left", "bottom"]:
-            # Toggle the spine objects
-            is_visible = not locals()[side]
-            ax_i.spines[side].set_visible(is_visible)
-            if offset is not None and is_visible:
-                try:
-                    val = offset.get(side, 0)
-                except AttributeError:
-                    val = offset
-                _set_spine_position(ax_i.spines[side], ('outward', val))
-
-        # Potentially move the ticks
-        if left and not right:
-            maj_on = any(
-                t.tick1line.get_visible()
-                for t in ax_i.yaxis.majorTicks
-            )
-            min_on = any(
-                t.tick1line.get_visible()
-                for t in ax_i.yaxis.minorTicks
-            )
-            ax_i.yaxis.set_ticks_position("right")
-            for t in ax_i.yaxis.majorTicks:
-                t.tick2line.set_visible(maj_on)
-            for t in ax_i.yaxis.minorTicks:
-                t.tick2line.set_visible(min_on)
-
-        if bottom and not top:
-            maj_on = any(
-                t.tick1line.get_visible()
-                for t in ax_i.xaxis.majorTicks
-            )
-            min_on = any(
-                t.tick1line.get_visible()
-                for t in ax_i.xaxis.minorTicks
-            )
-            ax_i.xaxis.set_ticks_position("top")
-            for t in ax_i.xaxis.majorTicks:
-                t.tick2line.set_visible(maj_on)
-            for t in ax_i.xaxis.minorTicks:
-                t.tick2line.set_visible(min_on)
-
-        if trim:
-            # clip off the parts of the spines that extend past major ticks
-            xticks = ax_i.get_xticks()
-            if xticks.size:
-                firsttick = np.compress(xticks >= min(ax_i.get_xlim()),
-                                        xticks)[0]
-                lasttick = np.compress(xticks <= max(ax_i.get_xlim()),
-                                       xticks)[-1]
-                ax_i.spines['bottom'].set_bounds(firsttick, lasttick)
-                ax_i.spines['top'].set_bounds(firsttick, lasttick)
-                newticks = xticks.compress(xticks <= lasttick)
-                newticks = newticks.compress(newticks >= firsttick)
-                ax_i.set_xticks(newticks)
-
-            yticks = ax_i.get_yticks()
-            if yticks.size:
-                firsttick = np.compress(yticks >= min(ax_i.get_ylim()),
-                                        yticks)[0]
-                lasttick = np.compress(yticks <= max(ax_i.get_ylim()),
-                                       yticks)[-1]
-                ax_i.spines['left'].set_bounds(firsttick, lasttick)
-                ax_i.spines['right'].set_bounds(firsttick, lasttick)
-                newticks = yticks.compress(yticks <= lasttick)
-                newticks = newticks.compress(newticks >= firsttick)
-                ax_i.set_yticks(newticks)
 
 if __name__ == '__main__':
     import argparse
     import numpy as np
     import h5py
     import pandas as pd
+    import itertools
+    from sddm import IDP_E_MINUS, IDP_MU_MINUS, SNOMAN_MASS
+    from sddm.plot import plot_hist, plot_legend, get_stats, despine, iqr_std_err, iqr_std, quantile_error, q90_err, q90, median, median_err, std_err
 
     parser = argparse.ArgumentParser("plot fit results")
     parser.add_argument("filenames", nargs='+', help="input files")
     parser.add_argument("--save", action="store_true", default=False, help="save plots")
     args = parser.parse_args()
 
+    if args.save:
+        # default \textwidth for a fullpage article in Latex is 16.50764 cm.
+        # You can figure this out by compiling the following TeX document:
+        #
+        #    \documentclass{article}
+        #    \usepackage{fullpage}
+        #    \usepackage{layouts}
+        #    \begin{document}
+        #    textwidth in cm: \printinunitsof{cm}\prntlen{\textwidth}
+        #    \end{document}
+
+        width = 16.50764
+        width /= 2.54 # cm -> inches
+        # According to this page:
+        # http://www-personal.umich.edu/~jpboyd/eng403_chap2_tuftegospel.pdf,
+        # Tufte suggests an aspect ratio of 1.5 - 1.6.
+        height = width/1.5
+        FIGSIZE = (width,height)
+
+        import matplotlib.pyplot as plt
+
+        font = {'family':'serif', 'serif': ['computer modern roman']}
+        plt.rc('font',**font)
+
+        plt.rc('text', usetex=True)
+    else:
+        # on retina screens, the default plots are way too small
+        # by using Qt5 and setting QT_AUTO_SCREEN_SCALE_FACTOR=1
+        # Qt5 will scale everything using the dpi in ~/.Xresources
+        import matplotlib
+        matplotlib.use("Qt5Agg")
+
+        import matplotlib.pyplot as plt
+
+        # Default figure size. Currently set to my monitor width and height so that
+        # things are properly formatted
+        FIGSIZE = (13.78,7.48)
+
+        # Make the defalt font bigger
+        plt.rc('font', size=22)
+
     # Read in all the data.
     #
     # Note: We have to add the filename as a column to each dataset since this
@@ -351,47 +148,6 @@ if __name__ == '__main__':
 
     markers = itertools.cycle(('o', 'v')) 
 
-    if args.save:
-        # default \textwidth for a fullpage article in Latex is 16.50764 cm.
-        # You can figure this out by compiling the following TeX document:
-        #
-        #    \documentclass{article}
-        #    \usepackage{fullpage}
-        #    \usepackage{layouts}
-        #    \begin{document}
-        #    textwidth in cm: \printinunitsof{cm}\prntlen{\textwidth}
-        #    \end{document}
-
-        width = 16.50764
-        width /= 2.54 # cm -> inches
-        # According to this page:
-        # http://www-personal.umich.edu/~jpboyd/eng403_chap2_tuftegospel.pdf,
-        # Tufte suggests an aspect ratio of 1.5 - 1.6.
-        height = width/1.5
-        FIGSIZE = (width,height)
-
-        import matplotlib.pyplot as plt
-
-        font = {'family':'serif', 'serif': ['computer modern roman']}
-        plt.rc('font',**font)
-
-        plt.rc('text', usetex=True)
-    else:
-        # on retina screens, the default plots are way too small
-        # by using Qt5 and setting QT_AUTO_SCREEN_SCALE_FACTOR=1
-        # Qt5 will scale everything using the dpi in ~/.Xresources
-        import matplotlib
-        matplotlib.use("Qt5Agg")
-
-        import matplotlib.pyplot as plt
-
-        # Default figure size. Currently set to my monitor width and height so that
-        # things are properly formatted
-        FIGSIZE = (13.78,7.48)
-
-        # Make the defalt font bigger
-        plt.rc('font', size=22)
-
     fig3, ax3 = plt.subplots(3,1,figsize=FIGSIZE,num=3,sharex=True)
     fig4, ax4 = plt.subplots(3,1,figsize=FIGSIZE,num=4,sharex=True)