text:
import pandas as pd
import numpy as np
from glob import glob
from os.path import join
import os
import matplotlib.pyplot as plt
import seaborn as sns
import re
import sklearn.metrics as metrics

def plot_result(df_full, x="test_score", conf_ctrl=[], input_type='',
                no_confs=False, sorty_by=None, join=False, beautify_io=True):

    input_types_cnt = 0
    if isinstance(df_full, (list, tuple)):
        if isinstance(input_type, (list, tuple)) and len(input_type)==len(df_full):
            df_full = [dfi.assign(i_type=i_type) for dfi, i_type in zip(df_full, input_type)]
            input_types_cnt = len(input_type)
        else:
            df_full = [dfi.assign(i_type=input_type) for dfi in df_full]
        df = pd.concat(df_full)
        is_multi = True
    else:
        df = df_full.copy()
        df = df.assign(i_type=input_type)
        is_multi = False
    # drop confound-related predictions if requested
    if "i_is_conf" in df.columns and no_confs:
        df = df[~(df.i_is_conf) & ~(df.o_is_conf)]

    # calculate p_values if permutation tests were run
    if f"permuted_{x}" in df.columns:
        df = calc_p_values(df, x)

    # hack: for compatability between old MLpipeline and new DLpipeline
    if "technique" not in df and "conf_ctrl_tech" in df:
        df = df.rename(columns={"conf_ctrl_tech":"technique"})
    # if conf_ctrl is not configured then automatically determine it
    if not conf_ctrl:
        conf_ctrl = df["technique"].unique()
    else:
        df = df[df["technique"].isin(conf_ctrl)]

    ios = df["io"].unique()

    # if the x metric is not in the results, then try to compute this metric from sklearn.metrics
    if x not in df.columns:
        df = compute_metric(df, x)

    # make io labels latex formated
    if beautify_io:
        df.loc[:,"io"] = df.apply(remap_io, axis=1)
    df[x] = df[x].apply(lambda x: x*100)

    # setup the figure properties
    sns.set(style='whitegrid', context='paper')
    fig, axes = plt.subplots(1, len(conf_ctrl),
                             sharex=True, sharey=True,
                             dpi=120, figsize=(4*len(conf_ctrl), 1+0.4*(len(ios)+input_types_cnt)))
    if not isinstance(axes, np.ndarray): axes = [axes]
    fig.gca().set_xlim([0,100])

    y="io"
    if is_multi:
        # if sample_size are different for different io (i.e. they are from different data subsets)
        # then print this info in y_label
        df = combine_io_sample_size(df)
        y="io_n"
    y_order = get_y_order(df, y, sorty_by)
    for ((t, dfi), ax) in zip(df.groupby("technique"), axes):
        assert (len(set(y_order) - set(dfi[y].unique())))==0 or ('baseline-cb' in df["technique"].unique())
        all_models = dfi["model"].unique()
        # plotting details
        palette = sns.color_palette()
        ci, dodge, scale, errwidth, capsize = 95, 0.4, 0.4, 0.9, 0.08
        ML_MODELS = ["LR", "SVM-lin", "SVM-rbf", "GB"]
        hue_order= [m for m in ML_MODELS if m in all_models]
        # if none of the ML_MODELS are in this df then don't differentiate by model
        if len(hue_order)==0: hue_order=None

        ax = sns.pointplot(y=y, x=x, order=y_order,
                           hue="model", hue_order=hue_order,
                           join=join, data=dfi, ax=ax,
                           ci=ci, errwidth=errwidth, capsize=capsize,
                           dodge=dodge, scale=scale, palette=palette)# todo bugfix ci is sensible only if each trial score is statistically independant
        ax.legend_.remove()
        ax.set_title("{} (ave. n={:.0f})".format(t.upper(), dfi["n_samples"].mean()))
        ax.set_xlabel(x)
        ax.set_ylabel("")

        # Add significance stars if permutation scores are available
        if f"permuted_{x}" in dfi.columns:
            # collect p_values as a dict of form {mean_accuracy: p_value}
            p_dict = {g[x].mean(): g["p_value"].iloc[0]
                      for i, g in dfi.groupby(["io", "model"])}
            # filter out the error bars from all lines
            # and store the position of their right edges to print the p_values at
            err_bars_pos = [(l.get_xdata()[1], l.get_ydata()[0])
                            for l in ax.lines if l.get_ydata()[0]==l.get_ydata()[1]]
            # also collect the position of the mean 'accuracy' points
            points_pos = [tuple(p_pos)  for c in ax.collections for p_pos in c.get_offsets()]

            for i, pos in enumerate(err_bars_pos):

                p_pos = points_pos[i]
                # check that the y position of both the point and error bar are the same
                assert p_pos[1] == pos[1], f"y position of mean pt ({p_pos}) and error bar ({pos}) don't match"
                # choose the error bar with the same mean accuracy point
                p =  p_dict[p_pos[0]]
                ast = return_asterisks(p)

                ax.annotate(ast, np.array(pos) + (0.02, 0.06),
                            color=palette[i//len(df[y].unique())], fontsize=5)

        # draw the chance line in the legend
        chance = (len(dfi[y].unique()))*[50] #todo
        for z, ch in enumerate(chance):
            ax.axvline(x=ch, ymin=z/len(chance), ymax=(z+1)/(len(chance)),
                        label="chance", c='gray', ls='--', lw=1.5)

    # add legend: add models info and chance label
    handles, legends = ax.get_legend_handles_labels()
    leg1 = fig.legend(handles[len(chance):], legends[len(chance):],
                      title="Models",
                      bbox_to_anchor=[1.13, 0.89], loc="upper right", fancybox=True, frameon=True)
    # only choose the first 'chance' legend
    leg2 = fig.legend([handles[0]], [legends[0]], loc="center right")
    fig.add_artist(leg1)

    fig.tight_layout()

    return fig

# Calculate p-value
def calc_p_values(df, x="test_score", viz=False):

    df["p_value"] = np.nan
    grp_order = ["io", "technique", "model"]
    if 'i_type' in df and len(df['i_type'].unique())>1:
        grp_order.insert(0, 'i_type')
    groups = df.dropna(subset=[x, f"permuted_{x}"]).groupby(grp_order)
    n_models = len(df["model"].unique())

    if viz:
        sns.reset_orig()
        n_rows = len(groups)//n_models
#         fig, axes = plt.subplots(n_rows, n_models,
#                                  sharex=True, sharey=False,
#                                  figsize=(20, n_models*n_rows))
        ## paper plot jugaad
        fig, axes = plt.subplots(1, 3,
                                 sharex=True, sharey=True,
                                 figsize=(12, 4))
        axes = np.ravel(axes)
        plt.xlim([0,1])

    for i, (g, rows) in enumerate(groups):

        p_vals = []
        rows = rows.filter(like=x)
        if viz:
            permuted_scores = []
            true_scores = rows[x]

        for _, r in rows.iterrows():
            p_scores = np.array(eval(r[f'permuted_{x}']))
            # calc p_value = (C + 1)/(n_permutations + 1)
            # where C is permutations whose score >= true_score
            true_score = r[x]
            C = np.sum(p_scores >= r[x])
            pi = (C+1)/(len(p_scores)+1)
            p_vals.extend([pi])

            if viz: permuted_scores.extend([*p_scores])
#         if np.std(permuted_scores)>=1:
#             print("[WARN] the p-values for {} have high variance across each test-set (trial). \
# Simply averaging the p-values across trials in such a case is not recommended.".format(g))

        df.loc[(df[grp_order]==g).all(axis=1), ["p_value", ""]] = np.mean(p_vals)

        if viz:
            ax = axes[i]
#             ax.set_title("Model={}".format(g[-1]))
#             if i%n_models == 0:
#                 ax.set_ylabel("{} with {}".format(*g[-3:-1]))
            ax.hist(permuted_scores, bins='auto', alpha=0.8)
            for true_score in true_scores:
                ax.axvline(true_score, color='r')
            # draw chance lines
            if x in ["test_score"]:
                # chance is 14% for site prediction and 50% for y and sex predictions
                chance = 0.5 if g[0][-1]!='c' else (1/7)
                ax.axvline(chance, color='grey', linestyle="--", lw=1.5)
                ax.set_xlim(0.,1.)
            ## paper plot jugaad
            ax.set_xticklabels([str(item) for item in range(0,120, 120//len(ax.get_xticklabels()))])
            inp = "X_{{{}}}".format(["14yr", "19yr", "22yr"][i])
            out = "y_{{{binge}}}"
            ax.set_title(r"${} \longmapsto {}$".format(inp,out))
#             if i==0: ax.set_ylabel("distribution / counts")
            if i==1: ax.set_xlabel("Balanced accuracy (%)")
            if i==0:
                from matplotlib.lines import Line2D
                custom_lines = [Line2D([0], [0], color="tab:blue", markerfacecolor="tab:blue", marker='o', markersize=5, lw=0),
                                Line2D([0], [0], color="tab:red", lw=1, linestyle="--")]
                ax.legend(custom_lines, ['permuted score', 'model score'], loc='lower left')

    sns.set(style='whitegrid')
    return df

def return_asterisks(pval):
    if pd.isna(pval):
        ast = ""
    elif pval <= 0.001:
        ast = '***'
    elif pval <= 0.01:
        ast = '**'
    elif pval <= 0.05:
        ast = '*'
    else:
        ast = 'n.s.'
    return ast

def remap_io(r):
    if r.i_is_conf:
        i = "c_{{{}}}".format(r.i.replace("_", "\_"))
    elif r.i == 'X' and ("i_type" in r):
        i = "X_{{{}}}".format(r.i_type.replace("_", "\_"))
    elif r.i == 'X':
        i = "X"
    else: # then its a label
        i = "y_{{{}}}".format(r.i.replace("_", "\_"))

    if r.o_is_conf:
        o = "c_{{{}}}".format(r.o.replace("_", "\_"))
    else: # then its a label
        o = "y_{{{}}}".format(r.o.replace("_", "\_"))

    new_io = r"${} \longmapsto {}$".format(i, o)
    return new_io

def combine_io_sample_size(df):
    df["io_n"] = ''
    # only estimate sample_size for the first technique to be plotted
    dfi = df.query("technique == '{}'".format(df.technique.unique()[0]))
    for io, dfii in dfi.groupby("io"):
        # set the same sample_size in io across all techniques to avoid different y_labels in subplots
        df.loc[(df["io"]==io), "io_n"] = "{}\n(n={:.0f})".format(io, dfii["n_samples"].mean())
    return df

def get_y_order(df, y="io", sorty_by=None):
    if sorty_by is not None:
        dfi = df.query("technique == '{}'".format(df.technique.unique()[0]))[[y, sorty_by]]
        y_order = dfi.groupby(y).mean().sort_values(by=sorty_by, ascending=True).index.values
        return y_order
    else:
        return df[y].sort_values().unique()

## METRICS ##
def specificity_score(y_true, y_pred):
    tn, fp, fn, tp = metrics.confusion_matrix(y_true, y_pred).ravel()
    return tn/(tn+fp)

def compute_metric(df, metric_name):

    metrics_map = {"recall_score":    metrics.recall_score,
                   "precision_score": metrics.precision_score,
                   "sensitivity":     metrics.recall_score,
                   "specificity":     specificity_score,
                   "accuracy_score":  metrics.accuracy_score,
                   "f1_score":        metrics.f1_score,
                  }
    metric=None
    for k in metrics_map.keys():
        if metric_name.lower() in k:
            metric = metrics_map[k]
            break
    if metric is None:
        raise ValueError(
            "ERROR: Invalid 'x' metric requested. Allowed metric_names are {}".format(
                metrics_map.keys()))

    df['y_true'] = df["test_lbls"].apply(lambda lbls: np.array(eval(lbls), dtype=int))
    df['y_pred'] = df["test_probs"].apply(lambda probs: np.argmax(np.array(eval(probs)), axis=1))

    df[metric_name] = df.apply(lambda x: metric(y_true=x.y_true, y_pred=x.y_pred), axis=1)

    return df


text:
import pandas as pd
import sys
import time
import numpy as np
import scipy as sp
import scipy.optimize as opt
import scipy.signal as signal
import os, json
import pygama.dataset as ds
import pygama.analysis.histograms as pgh
import pygama.dsp.transforms as pgt
import pygama.utils as pgu
import pygama.analysis.peak_fitting as pga
import matplotlib.pyplot as plt
from matplotlib.lines import Line2D
import matplotlib.colors as mcolors
from matplotlib.colors import LogNorm
from tqdm import tqdm
plt.style.use('style.mplstyle')

# for this code to work you need the raw_to_dsp, tier2, and spectrum files for the run of interest.

def main():

    plot_wfs()
    #flip_through_wfs()
    #ADC_difference()
    #ADC_difference_cut()
    #superpulse()
    #samples_scatter()
    #samples_hist()
    #baseline_hist()

def plot_wfs():

    if(len(sys.argv) != 2):
        print('Usage: wfs.py [run number]')
        sys.exit()

    start = time.time()

    with open("runDB.json") as f:
        runDB = json.load(f)
    tier_dir = os.path.expandvars(runDB["tier_dir"])

    df = pd.read_hdf('{}/t1_run{}.h5'.format(tier_dir,sys.argv[1]), '/ORSIS3302DecoderForEnergy')

    runtime = ds.DataSet(run=int(sys.argv[1]), md='./runDB.json').get_runtime()
    counts = len(df)
    counts_per_second = (len(df))/runtime

    print('total runtime = {} seconds'.format(runtime))
    print('counts = {}'.format(counts))
    print('counting rate = {} counts/second'.format(counts_per_second))
    exit()

    df = df.reset_index(drop=True)
    del df['energy']
    del df['channel']
    del df['energy_first']
    del df['ievt']
    del df['packet_id']
    del df['timestamp']
    del df['ts_hi']
    del df['ts_lo']

    def bl_sub(wf):
        return df.loc[wf,:]-df.iloc[wf,0:850].mean()

    nsamp = 3000

    xvals = np.arange(0,nsamp)

    for i in range(0,11):
        plt.plot(xvals, bl_sub(i), lw=1)
    plt.xlabel('Sample Number', ha='right', x=1.0)
    plt.ylabel('ADC Value', ha='right', y=1.0)
    plt.tight_layout()
    print('python script time = {} seconds'.format(time.time() - start))
    plt.show()

def flip_through_wfs():

    if(len(sys.argv) != 4):
        print('Usage: wfs.py [run number] [lower energy limit (keV)] [upper energy limit (keV)]')
        sys.exit()

    with open("runDB.json") as f:
        runDB = json.load(f)
    tier_dir = os.path.expandvars(runDB["tier_dir"])
    meta_dir = os.path.expandvars(runDB["meta_dir"])

    df = pd.read_hdf('{}/t1_run{}.h5'.format(tier_dir,sys.argv[1]), '/ORSIS3302DecoderForEnergy')
    df['e_cal'] = pd.read_hdf("{}/Spectrum_{}.hdf5".format(meta_dir,sys.argv[1]), key="df")['e_cal']
    df['AoverE'] = (pd.read_hdf("{}/t2_run{}.h5".format(tier_dir,sys.argv[1]))['current_max'])/df['e_cal']

    df = df.reset_index(drop=True)
    del df['energy']
    del df['channel']
    del df['energy_first']
    del df['ievt']
    del df['packet_id']
    del df['timestamp']
    del df['ts_hi']
    del df['ts_lo']

    df = df.loc[(df.AoverE>=0.05)&(df.AoverE<=0.06)]
    df = df.loc[(df.e_cal>int(sys.argv[2]))&(df.e_cal<int(sys.argv[3]))]
    df = df.reset_index(drop=True)
    df_2 = pd.DataFrame(df['e_cal'])
    del df['e_cal']
    del df['AoverE']

    def bl_sub(wf):
        return df.loc[wf,:]-df.iloc[wf,0:500].mean()

    def savgol(wf):
        return signal.savgol_filter(df.loc[wf,:]-df.iloc[wf,0:500].mean(), 47, 2)

    nsamp = 3000

    xvals = np.arange(0,nsamp)

    i = -1
    while True:
        if i != -1:
            inp = input()
            if inp == "q": exit()
            if inp == "p": i -= 2
        i += 1
        print(i)

        plt.cla()
        plt.plot(xvals, bl_sub(i), color="black", lw=2, label="raw wf, run {}, E = {:.03f} keV".format(str(sys.argv[1]), df_2['e_cal'][i]))
        plt.plot(xvals, savgol(i), color="red", lw=1, label="Savitzky-Golay Filter")
        plt.xlabel('Sample Number', ha='right', x=1.0)
        plt.ylabel('ADC Value', ha='right', y=1.0)
        plt.legend(frameon=True, loc='upper left', fontsize='x-small')
        plt.tight_layout()
        plt.show(block=False)
        plt.pause(0.01)

def ADC_difference():

    if(len(sys.argv) != 2):
        print('Usage: wfs.py [run number]')
        sys.exit()

    start = time.time()

    with open("runDB.json") as f:
        runDB = json.load(f)
    tier_dir = os.path.expandvars(runDB["tier_dir"])
    meta_dir = os.path.expandvars(runDB["meta_dir"])

    df = pd.read_hdf('{}/t1_run{}.h5'.format(tier_dir,sys.argv[1]), '/ORSIS3302DecoderForEnergy')
    df_2 = pd.read_hdf("{}/Spectrum_{}.hdf5".format(meta_dir,sys.argv[1]), key="df")

    df_2 = df_2.reset_index(drop=True)
    df = df.reset_index(drop=True)
    del df['energy']
    del df['channel']
    del df['energy_first']
    del df['ievt']
    del df['packet_id']
    del df['timestamp']
    del df['ts_hi']
    del df['ts_lo']
    #df_2['diff'] = [0]*len(df_2)

    #def dADC(wf):
        #return df.iloc[wf,1499:3000].mean()-df.iloc[wf,0:500].mean()

    #for i in range(len(df)):
        #df_2['diff'][i] = dADC(i)

    df_2['diff'] = df.iloc[:,1499:3000].mean(axis=1) - df.iloc[:,0:500].mean(axis=1)
    df_2['ratio'] = df_2['diff']/df_2['e_cal']
    print('python script time = {:.0f} seconds'.format(time.time() - start))

    plt.hist2d(df_2['e_cal'], df_2['diff'], np.arange(-5,100,0.1), norm=LogNorm())
    plt.xlim(0,50)
    plt.ylim(-5,100)
    plt.xlabel('Energy (keV)', ha='right', x=1.0)
    plt.ylabel('dADC', ha='right', y=1.0)
    cbar = plt.colorbar()
    cbar.ax.set_ylabel('Counts')
    plt.tight_layout()
    plt.show()

def ADC_difference_cut():

    if(len(sys.argv) != 2):
        print('Usage: wfs.py [run number]')
        sys.exit()

    start = time.time()

    with open("runDB.json") as f:
        runDB = json.load(f)
    tier_dir = os.path.expandvars(runDB["tier_dir"])
    meta_dir = os.path.expandvars(runDB["meta_dir"])

    df = pd.read_hdf('{}/t1_run{}.h5'.format(tier_dir,sys.argv[1]), '/ORSIS3302DecoderForEnergy')
    df_2 = pd.read_hdf("{}/Spectrum_{}.hdf5".format(meta_dir,sys.argv[1]), key="df")

    df_2 = df_2.reset_index(drop=True)
    df = df.reset_index(drop=True)
    del df['energy']
    del df['channel']
    del df['energy_first']
    del df['ievt']
    del df['packet_id']
    del df['timestamp']
    del df['ts_hi']
    del df['ts_lo']
    #df_2['diff'] = [0]*len(df_2)

    #def dADC(wf):
        #return df.iloc[wf,1499:3000].mean()-df.iloc[wf,0:500].mean()

    #for i in range(len(df)):
        #df_2['diff'][i] = dADC(i)

    df_2['diff'] = df.iloc[:,1499:3000].mean(axis=1) - df.iloc[:,0:500].mean(axis=1)
    df_2['ratio'] = df_2['diff']/df_2['e_cal']

    df_2 = df_2.loc[(df_2['diff']>3)&(df_2.e_cal>5)]
    print('python script time = {:.0f} seconds'.format(time.time() - start))

    plt.hist2d(df_2['e_cal'], df_2['diff'], np.arange(-5,100,0.1), norm=LogNorm())
    plt.xlim(0,50)
    plt.ylim(-5,100)
    plt.xlabel('Energy (keV)', ha='right', x=1.0)
    plt.ylabel('dADC', ha='right', y=1.0)
    cbar = plt.colorbar()
    cbar.ax.set_ylabel('Counts')
    plt.tight_layout()
    plt.show()

def superpulse():

    if(len(sys.argv) != 5):
        print('Usage: wfs.py [run number 1] [run number 2] [lower energy cut (keV)] [upper energy cut (keV)]')
        sys.exit()

    start = time.time()

    def pulse(run_number, lower_energy_cut, upper_energy_cut):

        with open("runDB.json") as f:
            runDB = json.load(f)
        tier_dir = os.path.expandvars(runDB["tier_dir"])
        meta_dir = os.path.expandvars(runDB["meta_dir"])

        df = pd.read_hdf('{}/t1_run{}.h5'.format(tier_dir,int(run_number)), '/ORSIS3302DecoderForEnergy')
        df_2 = pd.read_hdf('{}/Spectrum_{}.hdf5'.format(meta_dir,int(run_number)))

        df['e_cal'] = df_2['e_cal']

        df = df.loc[(df.e_cal>=float(lower_energy_cut))&(df.e_cal<=float(upper_energy_cut))]

        df = df.reset_index(drop=True)
        del df['energy']
        del df['channel']
        del df['energy_first']
        del df['ievt']
        del df['packet_id']
        del df['timestamp']
        del df['ts_hi']
        del df['ts_lo']
        del df['e_cal']

        superpulse = pd.DataFrame(df[0:-1][0:len(df)].mean(), columns=['superpulse'])
        superpulse['superpulse'] = superpulse['superpulse'] - superpulse['superpulse'][0:500].mean()

        return superpulse['superpulse'].values

    a = pulse(sys.argv[2],sys.argv[3],sys.argv[4])
    #b = pulse(sys.argv[1],sys.argv[3],sys.argv[4])

    #b = b*sum(a[2500:2999])/sum(b[2500:2999])

    nsamp = 3000
    xvals = np.arange(0,nsamp)
    plt.plot(xvals, a, lw=1, color='black', label='superpulse -- run {}, {}<E<{} keV'.format(sys.argv[2],sys.argv[3],sys.argv[4]))
    #plt.plot(xvals, b, lw=1, color='purple', label='superpulse -- run {}, {}<E<{} keV'.format(sys.argv[1],sys.argv[3],sys.argv[4]))
    #plt.plot(xvals, signal.savgol_filter(b, 47, 2), lw=1, color='blue', label='savgol filtered superpulse -- run {}, {}<E<{} keV'.format(sys.argv[1],sys.argv[3],sys.argv[4]))
    plt.plot(xvals, signal.savgol_filter(a, 47, 2), lw=1, color='red', label='savgol filtered superpulse -- run {}, {}<E<{} keV'.format(sys.argv[2],sys.argv[3],sys.argv[4]))
    plt.xlabel('Sample Number', ha='right', x=1.0)
    plt.ylabel('ADC Value', ha='right', y=1.0)
    plt.legend(frameon=True, loc='best', fontsize='small')
    plt.tight_layout()
    plt.show()

def samples_scatter():

    if(len(sys.argv) != 2):
        print('Usage: wfs.py [run number]')
        sys.exit()

    with open("runDB.json") as f:
        runDB = json.load(f)
    tier_dir = os.path.expandvars(runDB["tier_dir"])

    df = pd.read_hdf('{}/t1_run{}.h5'.format(tier_dir,sys.argv[1]), '/ORSIS3302DecoderForEnergy')

    runtime = ds.DataSet(run=int(sys.argv[1]), md='./runDB.json').get_runtime()
    counts_per_second = (len(df))/runtime

    df = df.reset_index(drop=True)
    del df['energy']
    del df['channel']
    del df['energy_first']
    del df['ievt']
    del df['packet_id']
    del df['timestamp']
    del df['ts_hi']
    del df['ts_lo']

    xvals = np.arange(0,len(df))

    plt.scatter(xvals, df[df.columns[454]], s=18, color='red', label='sample 454, run {}'.format(sys.argv[1]))
    plt.scatter(xvals, df[df.columns[456]], s=12, color='yellow', label='sample 456, run {}'.format(sys.argv[1]))
    plt.scatter(xvals, df[df.columns[455]], s=0.5, color='black', label='sample 455, run {}'.format(sys.argv[1]))
    plt.xlabel('waveform number', ha='right', x=1.0)
    plt.ylabel('ADV Value', ha='right', y=1.0)
    plt.legend(frameon=True, loc='best', fontsize='x-small')
    plt.tight_layout()
    plt.show()

def samples_hist():

    if(len(sys.argv) != 2):
        print('Usage: wfs.py [run number]')
        sys.exit()

    with open("runDB.json") as f:
        runDB = json.load(f)
    tier_dir = os.path.expandvars(runDB["tier_dir"])

    df = pd.read_hdf('{}/t1_run{}.h5'.format(tier_dir,sys.argv[1]), '/ORSIS3302DecoderForEnergy')

    runtime = ds.DataSet(run=int(sys.argv[1]), md='./runDB.json').get_runtime()
    counts_per_second = (len(df))/runtime

    df = df.reset_index(drop=True)
    del df['energy']
    del df['channel']
    del df['energy_first']
    del df['ievt']
    del df['packet_id']
    del df['timestamp']
    del df['ts_hi']
    del df['ts_lo']

    plt.hist(df[df.columns[455]], np.arange(-8000,-1000,1), histtype='step', color='black', label='sample 455, run {}'.format(sys.argv[1]))
    plt.hist(df[df.columns[456]], np.arange(-8000,-1000,1), histtype='step', color='red', label='sample 456, run {}'.format(sys.argv[1]))
    plt.hist(df[df.columns[454]], np.arange(-8000,-1000,1), histtype='step', color='blue', label='sample 454, run {}'.format(sys.argv[1]))
    plt.ylim(0,plt.ylim()[1])
    plt.xlabel('ADC Value', ha='right', x=1.0)
    plt.ylabel('Counts', ha='right', y=1.0)
    plt.legend(frameon=True, loc='best', fontsize='small')
    plt.tight_layout()
    plt.show()

def baseline_hist():

    if(len(sys.argv) != 2):
        print('Usage: wfs.py [run number]')
        sys.exit()

    with open("runDB.json") as f:
       runDB = json.load(f)
    tier_dir = os.path.expandvars(runDB["tier_dir"])

    df = pd.read_hdf('{}/t1_run{}.h5'.format(tier_dir,sys.argv[1]), '/ORSIS3302DecoderForEnergy')
    df_2 = pd.read_hdf('{}/t2_run{}.h5'.format(tier_dir,sys.argv[1]))

    df_2 = df_2.reset_index(drop=True)
    df = df.reset_index(drop=True)
    del df['energy']
    del df['channel']
    del df['energy_first']
    del df['ievt']
    del df['packet_id']
    del df['timestamp']
    del df['ts_hi']
    del df['ts_lo']

    df[df.columns[0]] = df[df.columns[0]] - df_2['bl_p0']
    baseline_samples = df[df.columns[0]].values

    for i in tqdm(range(1,501)):
        df[df.columns[i]] = df[df.columns[i]] - df_2['bl_p0']
        a = df[df.columns[int(i)]].values
        baseline_samples = np.append(baseline_samples, a)

    plt.hist(signal.savgol_filter(baseline_samples,47,2), np.arange(-50,50,0.25),  histtype='step', color='black', label='filtered baseline samples of waveforms, run {}'.format(sys.argv[1]))
    #plt.ylim(0,plt.ylim()[1])
    plt.xlabel('ADC Value', ha='right', x=1.0)
    plt.ylabel('Counts', ha='right', y=1.0)
    plt.semilogy()
    plt.legend(frameon=True, loc='best', fontsize='small')
    plt.tight_layout()
    plt.show()

if __name__ == '__main__':
        main()

ext: py

max_stars_repo_name: sweigart/pygama

max_issues_repo_name: sweigart/pygama

max_forks_repo_name: sweigart/pygama

reason: import numpy,import scipy


text:
"""Methods for patching the ACAS Xu network with SyReNN."""
from timeit import default_timer as timer
import numpy as np
from pysyrenn import FullyConnectedLayer
from prdnn import DDNN, FTRepair
from experiments.acas_repair import ACASRepair

class ACASMFT(ACASRepair):
    """Experiment testing patching performance on an ACAS Xu model."""
    def do_repair(self, train_regions, train_syrenn, syrenn_time):
        n_unique = len(set({tuple(point)
                            for upolytope in train_syrenn
                            for pre_poly in upolytope
                            for point in pre_poly}))
        samples_per_plane = [n_unique // len(train_regions) for _ in range(len(train_regions))]

        _patcher = FTRepair.from_spec_function(
            self.network, train_regions, self.property,
            samples_per_plane=samples_per_plane)
        self.record_artifact(_patcher.inputs, "train_inputs", "pickle")
        self.record_artifact(_patcher.labels, "train_labels", "pickle")

        patchable = [i for i, layer in enumerate(self.network.layers)
                     if isinstance(layer, FullyConnectedLayer)]
        for layer in patchable:
            patcher = FTRepair(self.network, _patcher.inputs.copy(), _patcher.labels.copy())
            patcher.layer = layer
            patcher.norm_objective = True
            patcher.auto_stop = False
            patcher.make_holdout_set()
            patcher.epochs = 1000
            patcher.lr = 0.001
            patcher.momentum = 0.9
            patcher.batch_size = 16
            patched = patcher.compute()
            assert patched is not None
            self.record_artifact(patched, f"{layer}_patched", "network")
            timing = patcher.timing.copy()
            timing["syrenn"] = syrenn_time
            self.record_artifact(timing, f"{layer}_timing", "pickle")

    def analyze(self):
        unpatched = self.read_artifact("unpatched")
        self.network = unpatched

        train_syrenn = self.read_artifact("train_syrenn")
        test_syrenn = self.read_artifact("test_syrenn")
        train_inputs = self.read_artifact("train_inputs")
        train_labels = self.read_artifact("train_labels")

        patchable = [i for i, layer in enumerate(self.network.layers)
                     if isinstance(layer, FullyConnectedLayer)]
        for layer in patchable:
            print("~~~~~ Layer:", layer, "~~~~~")
            patched = self.read_artifact(f"{layer}_patched")
            timing = self.read_artifact(f"{layer}_timing")

            total_points = sum(len(pre_poly) for upolytope in train_syrenn
                               for pre_poly in upolytope)
            print("Size of repair set:", total_points)

            og_train_outputs = np.argmax(unpatched.compute(train_inputs), axis=1)
            print("Number of repair set points originally buggy:",
                  np.sum(og_train_outputs != train_labels))

            gen_set, drawdown_set = self.find_counterexamples(unpatched, test_syrenn)
            print("Size of generalization, drawdown sets:", len(gen_set), len(drawdown_set))

            print("Time (seconds):", timing["total"])

            train_outputs = np.argmax(patched.compute(train_inputs), axis=1)
            print("Number of repair set points buggy after repair:",
                  np.sum(train_outputs != train_labels))

            dd_desired = self.property(drawdown_set)
            dd_outputs = np.argmax(patched.compute(drawdown_set), axis=1)
            print("Drawdown-set counterexamples after repair:", np.sum(dd_desired != dd_outputs))

            gen_desired = self.property(gen_set)
            gen_outputs = np.argmax(patched.compute(gen_set), axis=1)
            print("Generalization-set counterexamples after repair:", np.sum(gen_desired != gen_outputs))

        return True

if __name__ == "__main__":
    ACASMFT("acas_mft").main()

reason: import numpy


text:
from .Node import Node
from .Detector import Detector
from . import components
import sympy as sp
import numpy as np
import timeit
from scipy.sparse.linalg import lsqr

class Model:
    """Defines the optical network to be modelled.

    The Model class holds the information that defines the structure of the
    optical network to be modelled with strapy, along with functions for
    building and evaluating the model. In general only the `wavelength`
    attribute should be accessed directly by the user - all other attributes
    should be set through the member functions.

    Attributes
    ----------
    wavelength : float
            The vacuum (n = 1) wavelength of illuminating light.
    components : dict
            Dictionary of optical components that have been added to the model.
    detectors : dict
            Dictionary of detectors that have been added to the model.
    nodes : dict
            Dictionary of nodes that have been added to the model.
    symbols : list
            List of sympy symbols used by the model.
    equations : list
            List of sympy equations used by the model.
    rhsVariables : list
            List of sympy symbols that are included in the right hand side
            vector of the matrix equation.
    matrixVariables : list
            List of sympy symbols that are included in the network matrix.
    updated : list
            List of optical components that have changed since the model was
            last evaluated.
    useLambdify : bool
            True if sympy's lambdify functionality should be used to set right
            hand side vector and matrix before evaluation. If false only
            components that are included in the `updated` list are updated in
            the matrix equation before solving. Lambdify is currently a
            fallback if there is more than one sympy symbol per matrix element.
    """

    def __init__(self):
        self.wavelength = 633e-9
        self.components = {}
        self.detectors = {}
        self.nodes = {}
        self.symbols = []
        self.equations = []
        self.rhsVariables = []
        self.matrixVariables = []
        self.updated = []
        self.useLambdify = False

    def add_component(self, component, name, nodes):
        """Adds component to model and updates the node list.

        Parameters
        ----------
        component : strapy.components.component
                Type of component to be added.
        name : str
                Unique name of component.
        nodes : str or tuple of str
                Node(s) to which the component is attached.
        """

        # allows for passing string instead of tuple of strings
        if not(isinstance(nodes, tuple)):
            nodes = (nodes,)

        self.components[name] = component(name, nodes, self)
        self.updated.append(name)

        if len(nodes) != self.components[name].node_number:
            raise Exception(
                'Specified nodes ({}) not equal to required node',
                ' number ({}).'.format(
                    len(nodes),
                    self.components[name].node_number))

        # keep nodes list up to date when adding a component
        for node in nodes:
            if not(node in self.nodes):
                self.nodes[node] = Node(node)
            if len(self.nodes[node].components) > 1:
                raise Exception(
                    'Nodes may not be linked to more that two',
                    ' components. Node {} is attached to {}.'.format(
                        node, len(self.nodes[node].components) + 1))
            else:
                self.nodes[node].components.append(self.components[name])

    def add_detector(self, name, node, properties=('amplitude',)):
        """Adds detector to model and updates list of nodes.

        For details on currently implemented properties see
        :py:class:`strapy.Detector()`.

        Parameters
        ----------
        name : str
                Unique name of the detector.
        nodes : str
                Node monitored by detector.
        properties : tuple of str
                Optical properties to be logged - see
                :py:class:`strapy.Detector()` for current options.
        """

        node = (node,)
        self.detectors[name] = Detector(name, node, properties)

        if not(node[0] in self.nodes):
            self.nodes[node[0]] = Node(node[0])

    def build(self, verbose=False):
        """Builds network matrix from defined components.

        The model must be built before `evaluate()` is called. Additionally, if
        components or detectors are added to the model between runs, the model
        must be rebuilt.

        Parameters
        ----------
        verbose : bool
                If true details of the constructed optical network will be
                printed during the build process.
        """

        sourceFlag = False

        # check for at least one source in the network
        for component in self.components.values():
            if isinstance(component, components.Source):
                sourceFlag = True
                break

        if not(sourceFlag):
            raise Exception('No sources present.')

        # Threadable.
        for node in self.nodes.values():
            if len(node.components) != 2:
                raise Exception(
                    'Nodes must be linked to exactly two',
                    ' components, node {} is attached to {}.'.format(
                        node.name, len(node.components)))

            self.symbols.extend(node.symbols)

            if verbose:
                print('Node {}'.format(node.name))
                print('\tComponent 1: {}'.format(node.components[0].name))
                print('\tComponent 2: {}'.format(node.components[1].name))

        # Should be threadable.
        for component in self.components.values():
            component.initEquation(self.nodes)
            component.symbolIdxs = [[] for _ in range(len(component.symbols))]

            for line in range(len(component.equation.lhs)):
                self.equations.append(sp.Eq(component.equation.lhs[line],
                                            component.equation.rhs[line]))
            if isinstance(component, components.Source):
                self.rhsVariables.extend(component.symbols)
            else:
                self.matrixVariables.extend(component.symbols)

        networkMatrix, rhsVector = sp.linear_eq_to_matrix(self.equations,
                                                          self.symbols)
        self.matrixShape = networkMatrix.shape

        self.matrix = np.zeros(self.matrixShape, dtype=np.complex)
        self.rhs = np.zeros(rhsVector.shape, dtype=np.complex)

        # These loops may be thread safe, other than the lambdify command, and
        # the loops could therefore be parfored. Unknown if access to a sympy
        # matrix is actually thread safe, needs testing.
        for i in range(networkMatrix.shape[0]):
            for j in range(networkMatrix.shape[1]):
                # isinstance() is significantly faster than .equals(0)
                if not(isinstance(networkMatrix[i, j], sp.core.numbers.Zero)):
                    if networkMatrix[i, j].is_constant():
                        self.matrix[i, j] = sp.N(networkMatrix[i, j])
                    else:
                        for component in self.components.values():
                            for k, symbol in enumerate(component.symbols):
                                if len(networkMatrix[i, j].free_symbols) > 1:
                                    print(
                                        "Multiple symbols per matrix element",
                                        ", falling back to lambdify matrix",
                                        " setting.")
                                    self.useLambdify = True
                                    self.setMatrix = sp.lambdify(
                                        self.matrixVariables,
                                        networkMatrix,
                                        modules=["numpy"])
                                if networkMatrix[i, j] == symbol:
                                    component.symbolIdxs[k].append((i, j, 1))
                                if networkMatrix[i, j] == -symbol:
                                    component.symbolIdxs[k].append((i, j, -1))

        # Again, may be threadable.
        for i in range(len(rhsVector)):
            # isinstance() is significantly faster than .equals(0)
            if not(isinstance(rhsVector[i], sp.core.numbers.Zero)):
                if rhsVector[i].is_constant():
                    self.rhs[i, j] = sp.N(rhsVector[i])
                else:
                    for component in self.components.values():
                        for k, symbol in enumerate(component.symbols):
                            if len(rhsVector[i].free_symbols) > 1:
                                print(
                                    "Multiple symbols per rhs vector element",
                                    ", falling back to lambdify rhs vector",
                                    " setting.")
                                self.useLambdify = True
                                self.setRhs = sp.lambdify(
                                    self.rhsVariables,
                                    rhsVector,
                                    modules=["numpy"])
                            if rhsVector[i] == symbol:
                                component.symbolIdxs[k].append((i, 0, 1))
                            if rhsVector[i] == -symbol:
                                component.symbolIdxs[k].append((i, 0, -1))

        # identify the location in the solution vector which the detector
        # should detect.
        for detector in self.detectors.values():
            detector.node_index = self.symbols.index(
                self.nodes[detector.node[0]].symbols[0])

        self.matrixPassVector = np.empty((len(self.matrixVariables),),
                                         dtype=np.complex)
        self.rhsPassVector = np.empty((len(self.rhsVariables),),
                                      dtype=np.complex)

        self.sparcity = len(self.matrixPassVector) \
            / (self.matrixShape[0] * self.matrixShape[1])

        # Should be threadable.
        for component in self.components.values():
            # print(component.name, component.symbols, component.symbolIdxs)
            if isinstance(component, components.Source):
                component.set_slice = slice(self.rhsVariables.index(
                    component.symbols[0]),
                    self.rhsVariables.index(component.symbols[0])
                    + len(component.symbols))
            elif not(isinstance(component, components.Dump)):
                component.set_slice = slice(
                    self.matrixVariables.index(component.symbols[0]),
                    self.matrixVariables.index(component.symbols[0])
                    + len(component.symbols))

        if verbose:
            print('\nNetwork built, {} matrix.\n'.format(networkMatrix.shape))

    def _update(self):
        """Updates numerical values in model matrix.

        Model must have been built first.
        Should not be called externally.
        """

        # Threadable
        for key in self.updated:
            if isinstance(self.components[key], components.Source):
                self.rhsPassVector[self.components[key].set_slice] = \
                    self.components[key].setVals()
            elif not(isinstance(self.components[key], components.Dump)):
                self.matrixPassVector[self.components[key].set_slice] = \
                    self.components[key].setVals()
        self.updated.clear()

    def evaluate(self, timing=False):
        """Solve the network matrix and log optical properties at detectors.

        `build()` must have been called before the model is evaluated.

        Parameters
        ----------
        timing : bool
                If true, the times taken to set values in the right hand side
                vector and network matrix (`set_time`), solve the network
                equation (`solve_time`) and pull out the detected values
                (`detector_time`) are returned."""

        if timing:
            set_time = timeit.default_timer()

            if not(self.useLambdify):
                for key in self.updated:
                    if isinstance(self.components[key], components.Source):
                        vals = self.components[key].setVals()
                        for i, index in enumerate(
                                self.components[key].symbolIdxs):
                            self.rhs[index[0]] = vals[i] * index[2]
                    elif not(isinstance(self.components[key],
                                        components.Dump)):
                        vals = self.components[key].setVals()
                        for i, index in enumerate(
                                self.components[key].symbolIdxs):
                            self.matrix[index[0], index[1]
                                        ] = vals[i] * index[2]
                self.updated.clear()
            else:
                self._update()
                self.matrix = self.setMatrix(*self.matrixPassVector)
                self.rhs = self.setRhs(*self.rhsPassVector)

            set_time = timeit.default_timer() - set_time

            solve_time = timeit.default_timer()
            self.solution_vector = np.linalg.solve(self.matrix, self.rhs)
            solve_time = timeit.default_timer() - solve_time

            detector_time = timeit.default_timer()
            for detector in self.detectors.values():
                detector.update(self.solution_vector)
            detector_time = timeit.default_timer() - detector_time

            return (set_time, solve_time, detector_time)
        else:
            if not(self.useLambdify):
                # Should be threadable.
                for key in self.updated:
                    if isinstance(self.components[key], components.Source):
                        vals = self.components[key].setVals()
                        for i, symbol in enumerate(
                                self.components[key].symbolIdxs):
                            for index in symbol:
                                self.rhs[index[0]] = vals[i] * index[2]
                    elif not(isinstance(self.components[key],
                                        components.Dump)):
                        vals = self.components[key].setVals()
                        for i, symbol in enumerate(
                                self.components[key].symbolIdxs):
                            for index in symbol:
                                self.matrix[index[0], index[1]] = \
                                    vals[i] * index[2]
                self.updated.clear()
            else:
                self._update()
                self.matrix = self.setMatrix(*self.matrixPassVector)
                self.rhs = self.setRhs(*self.rhsPassVector)

            self.solution_vector = np.linalg.solve(self.matrix, self.rhs)

            for detector in self.detectors.values():
                detector.update(self.solution_vector)

max_stars_repo_name: scampy-project/scampy

max_issues_repo_name: scampy-project/scampy

max_forks_repo_name: scampy-project/scampy

reason: import numpy,from scipy,import sympy


text:
import numpy as np

import util.img
import util.plt
import util.dec
def test_ds_size():
    image_size = (15, 15)
    kernel_size = (2, 2)
    stride = (1, 1)
    np.testing.assert_equal(util.img.ds_size(image_size, kernel_size, stride), (14, 14))

    kernel_size = (2, 2)
    stride = (3, 3)
    np.testing.assert_equal(util.img.ds_size(image_size, kernel_size, stride), (5, 5))

    kernel_size = (3, 3)
    stride = (2, 2)
    np.testing.assert_equal(util.img.ds_size(image_size, kernel_size, stride), (7, 7))

    kernel_size = (3, 2)
    stride = (3, 2)
    np.testing.assert_equal(util.img.ds_size(image_size, kernel_size, stride), (5, 7))

@util.dec.print_test
def test_blur():
    image_path = '~/dataset/ICDAR2015/Challenge4/ch4_training_images/img_150.jpg'
    image = util.img.imread(image_path)
    blurred_3 = util.img.average_blur(image, (3, 3))
    blurred_5 = util.img.average_blur(image, (5, 5))
    util.plt.show_images(images = [image, blurred_3, blurred_5], titles = ['origin', '3x3', '5x5'], bgr2rgb = True, share_axis = True)

@util.dec.print_test
def test_rect_perimeter():
    p = util.img.rect_perimeter((5, 5), (80, 80))
    np.testing.assert_equal(p, 300)

    p = util.img.rect_perimeter((0, 0), (80, 80))
    np.testing.assert_equal(p, 320)

@util.dec.print_test
def test_rect_perimeter():
    lu, rb = (5, 5), (80, 80)
    p = util.img.rect_area(lu, rb)
    mask =  util.img.black((100, 100))
    util.img.rectangle(mask, lu, rb, border_width = -1, color = 1)
    util.test.assert_equal(p, np.sum(mask))

@util.dec.print_test
def _test_resize():
    size = np.array((100, 100))
    img = util.img.black(size)
    util.img.rectangle(img, (10, 10), (90, 90), color = 1, border_width = -1)
    util.img.rectangle(img, (20, 20), (80, 80), color = 2, border_width = -1)
    util.img.rectangle(img, (30, 30), (70, 70), color = 3, border_width = -1)
    util.img.rectangle(img, (40, 40), (60, 60), color = 4, border_width = -1)

    img1 = util.img.resize(img, f = 0.5)
    util.test.assert_equal(util.img.get_shape(img1), 0.5 * size)

    img2 = util.img.resize(img, fx = 0.1, fy = 1)
    print util.img.get_wh(img2), (0.1 * size[0], size[1])
    util.test.assert_equal(util.img.get_wh(img2), (0.1 * size[0], size[1]))

    img3 = util.img.resize(img, size = (300, 400))
    util.test.assert_equal(util.img.get_wh(img3), (300, 400))
#    util.plt.show_images([img, img1, img2, img3])

@util.dec.print_test
def _test_translate():
    img = util.img.white((100, 100))
    dst1 = util.img.translate(img, delta_x = 10, delta_y = 20)
    dst2 = util.img.translate(img, delta_x = 0, delta_y = -20)
    util.plt.show_images([img, dst1, dst2], gray = True)

@util.dec.print_test
def _test_get_cnt_rect_bbox():
    img = util.img.black((100, 100))
    util.img.circle(img, center = (50, 50), r = 20, color = 255, border_width = -1)
    util.img.imshow('image', img)

    cnt = util.mask.find_contour(img)[0]
    x,y, w, h = util.img.get_contour_rect_box(cnt)
    util.img.rectangle(img, left_up = (x, y), right_bottom =  (x + w, y + h), color = 200, border_width = 1)
    util.img.imshow('image', img)

@util.dec.print_test
def _test_get_cnt_region_in_rect():
    img = util.img.black((100, 100))
    util.img.rectangle(img, left_up = (50, 50), right_bottom = (90, 80), color = 255, border_width = -1)
    util.img.imshow('Image', img)
    img = util.img.rotate_about_center(img, angle = 30)
    util.img.imshow('Image', img)
    cnt = util.mask.find_contour(img)[0]
    min_area_region = util.img.get_contour_region_in_min_area_rect(img, cnt)
    util.img.imshow('Region', min_area_region)

    min_area_region = util.img.get_contour_region_in_rect(img, cnt)
    util.img.imshow('Region', min_area_region)

@util.dec.print_test
def _test_get_contour_region_iou():
    img = util.img.black((100, 100))
    rect1 = util.img.get_rect_points((10, 30), (60, 80))
    rect2 = util.img.get_rect_points((35, 5), (85, 55))
    #util.img.rectangle(img, rect1[0], rect1[-2], 1, -1)
    cnt1 = util.img.points_to_contour(rect1)
    cnt2 = util.img.points_to_contour(rect2)
    util.test.assert_almost_equal(util.img.get_contour_region_iou(img, cnt1, cnt2), 1. / 7, 1)

@util.dec.print_test
def _test_rect_iou():
    img = util.img.black((100, 100))
    rect1 = util.img.get_rect_points((10, 30), (60, 80))
    rect2 = util.img.get_rect_points((35, 5), (85, 55))
    cnt1 = util.img.points_to_contour(rect1)
    cnt2 = util.img.points_to_contour(rect2)
    # util.test.assert_almost_equal(util.img.get_contour_region_iou(img, cnt1, cnt2), 1. / 7, 1)
    iou = util.img.get_contour_region_iou(img, cnt1, cnt2)
    rects1 = [[10, 30, 60, 80], [35, 5, 85, 55], [90, 90, 95, 95]]
    rects2 = [[10, 30, 60, 80], [35, 5, 85, 55]]
    iou_matrix = util.img.get_rect_iou(rects1, rects2)
    print iou_matrix

@util.dec.print_test
def _test_contour_convex():

    mask = util.img.black((100, 100))
    util.img.rectangle(img = mask, left_up = (30, 30), right_bottom = (99, 99), color = 1, border_width = -1)
    util.img.rectangle(img = mask, left_up = (50, 50), right_bottom = (60, 60), color = 0, border_width = -1)
    util.img.rectangle(img = mask, left_up = (40, 20), right_bottom = (50, 40), color = 0, border_width = -1)

    # util.plt.show_images(images=[mask], titles=['Mask', 'Hole Mask'])

    contours = util.img.find_contours(mask)
    vis = util.img.black((100, 100))
    util.img.draw_contours(vis, contours, idx = 0, color = 1, border_width = 1)

    hull = util.img.convex_hull(contours[0])
    util.img.draw_contours(mask, [hull], idx = -1, color = 1, border_width = 1)
    util.plt.show_images(images = [mask, vis])

if util.mod.is_main(__name__):
    util.init_logger()
    #test_rect_perimeter()
    #test_ds_size()
    #test_blur()
    #test_rect_perimeter()
    #_test_resize()
    #_test_translate()
    #_test_get_cnt_rect_bbox()
    #_test_get_cnt_region_in_rect()
    #_test_get_contour_region_iou()
    #_test_rect_iou()
    _test_contour_convex()


text:
# Copyright The PyTorch Lightning team.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
import inspect
from typing import Any, Callable, Dict, Optional, Sequence, Tuple, Type

import numpy as np

from flash.core.data.data_source import DefaultDataKeys, LabelsState
from flash.core.integrations.icevision.transforms import from_icevision_record
from flash.core.utilities.imports import _ICEVISION_AVAILABLE
from flash.image.data import ImagePathsDataSource

if _ICEVISION_AVAILABLE:
    from icevision.core.record import BaseRecord
    from icevision.core.record_components import ClassMapRecordComponent, FilepathRecordComponent, tasks
    from icevision.data.data_splitter import SingleSplitSplitter
    from icevision.parsers.parser import Parser

class IceVisionPathsDataSource(ImagePathsDataSource):
    def predict_load_data(self, data: Tuple[str, str], dataset: Optional[Any] = None) -> Sequence[Dict[str, Any]]:
        return super().predict_load_data(data, dataset)

    def load_sample(self, sample: Dict[str, Any]) -> Dict[str, Any]:
        record = sample[DefaultDataKeys.INPUT].load()
        return from_icevision_record(record)

    def predict_load_sample(self, sample: Dict[str, Any]) -> Dict[str, Any]:
        if isinstance(sample[DefaultDataKeys.INPUT], BaseRecord):
            return self.load_sample(sample)
        filepath = sample[DefaultDataKeys.INPUT]
        sample = super().load_sample(sample)
        image = np.array(sample[DefaultDataKeys.INPUT])

        record = BaseRecord([FilepathRecordComponent()])
        record.filepath = filepath
        record.set_img(image)
        record.add_component(ClassMapRecordComponent(task=tasks.detection))
        return from_icevision_record(record)

class IceVisionParserDataSource(IceVisionPathsDataSource):
    def __init__(self, parser: Optional[Type["Parser"]] = None):
        super().__init__()
        self.parser = parser

    def load_data(self, data: Tuple[str, str], dataset: Optional[Any] = None) -> Sequence[Dict[str, Any]]:
        if self.parser is not None:
            if inspect.isclass(self.parser) and issubclass(self.parser, Parser):
                root, ann_file = data
                parser = self.parser(ann_file, root)
            elif isinstance(self.parser, Callable):
                parser = self.parser(data)
            else:
                raise ValueError("The parser must be a callable or an IceVision Parser type.")
            dataset.num_classes = parser.class_map.num_classes
            self.set_state(LabelsState([parser.class_map.get_by_id(i) for i in range(dataset.num_classes)]))
            records = parser.parse(data_splitter=SingleSplitSplitter())
            return [{DefaultDataKeys.INPUT: record} for record in records[0]]
        raise ValueError("The parser argument must be provided.")

    def predict_load_data(self, data: Any, dataset: Optional[Any] = None) -> Sequence[Dict[str, Any]]:
        result = super().predict_load_data(data, dataset)
        if len(result) == 0:
            result = self.load_data(data, dataset)
        return result

max_stars_repo_name: Isaac-Flath/lightning-flash

max_issues_repo_name: Isaac-Flath/lightning-flash

max_forks_repo_name: Isaac-Flath/lightning-flash


text:
import numpy as np
from matplotlib import pyplot
from live2anime.io.load import load_images

# dataset path
path = 'data/png_files/'
# load dataset
images_live_action, images_anime = load_images(path)
print('Loaded: ', images_live_action.shape, images_anime.shape)

# save as compressed numpy array
filename = 'data/npz_files/pictures.npz'
np.savez_compressed(filename, images_live_action, images_anime)
print('Saved dataset: ', filename)

# load the prepared dataset
# load the dataset
data = np.load('data/npz_files/pictures.npz')
images_live_action, images_anime = data['arr_0'], data['arr_1']
print('Loaded: ', images_live_action.shape, images_anime.shape)
# plot source images
n_samples = 3
for i in range(n_samples):
	pyplot.subplot(2, n_samples, 1 + i)
	pyplot.axis('off')
	pyplot.imshow(images_live_action[i].astype('uint8'))
# plot target image
for i in range(n_samples):
	pyplot.subplot(2, n_samples, 1 + n_samples + i)
	pyplot.axis('off')
	pyplot.imshow(images_anime[i].astype('uint8'))
pyplot.show()


text:
import numpy as np
import matplotlib.pyplot as plt

def plot_concat(x,figsize=(20,10),cmap="terrain",interpolation="sinc",xlabel=None,ylabel=None,xlabels=[],ylabels=[],colorbar=True,compose=False,path=None):
    fig = plt.figure(figsize=figsize)
    ax = plt.gca()

    if compose:
        plot = plt.imshow(x,interpolation=interpolation,extent=[0,1,0,1],aspect="auto")
    else:
        plot = plt.imshow(x,interpolation=interpolation,extent=[0,1,0,1],cmap=cmap,aspect="auto")

    if len(xlabels):
        xtickpos = np.append(np.arange(0,1,1/(len(xlabels)-1)),[1])
        ax.set_xticks(xtickpos)
        ax.set_xticklabels(xlabels)

    if len(ylabels):
        ytickpos = np.append(np.arange(0,1,1/(len(ylabels)-1)),[1])
        ax.set_yticks(ytickpos)
        ax.set_yticklabels(ylabels)

    if not compose and colorbar:
        fig.colorbar(plot, ax=ax)

    plt.xticks(rotation=90)

    if xlabel:
        ax.set_xlabel(xlabel)

    if ylabel:
        ax.set_ylabel(ylabel)

    if path is not None:
        plt.savefig(path)

    return fig,ax,plot

def plot_aggr(x,figsize=(20,10),cmap="terrain",interpolation="sinc",xlabel=None,ylabel=None,xlabels=[],ylabels=[],aggr="mean",colorbar=True,compose=False,path=None):
    fig = plt.figure(figsize=figsize)
    ax = plt.gca()

    if compose:
        plot = plt.imshow(x,interpolation=interpolation,extent=[0,1,0,1],aspect="auto")
    else:
        plot = plt.imshow(x,interpolation=interpolation,extent=[0,1,0,1],cmap=cmap,aspect="auto")

    if len(xlabels):
        xtickpos = np.append(np.arange(0,1,1/(len(xlabels)-1)),[1])
        ax.set_xticks(xtickpos)
        ax.set_xticklabels(xlabels)

    if len(ylabels):
        ytickpos = np.append(np.arange(0,1,1/(len(ylabels)-1)),[1])
        ax.set_yticks(ytickpos)
        ax.set_yticklabels(ylabels)

    if not compose and colorbar:
        fig.colorbar(plot, ax=ax)

    if xlabel:
        ax.set_xlabel(xlabel)

    if ylabel:
        ax.set_ylabel(ylabel)

    plt.xticks(rotation=90)

    if path is not None:
        plt.savefig(path)

    return fig,ax,plot

max_stars_repo_name: CONABIO/yuntu

max_issues_repo_name: CONABIO/yuntu

max_forks_repo_name: CONABIO-audio/yuntu


text:
from collections import namedtuple, Counter
import os
import re

import matplotlib.pyplot as plt
import numpy as np
from PIL import Image
import tweepy
from wordcloud import WordCloud, STOPWORDS

Tweet = namedtuple('Tweet', 'id text created likes rts')

TWITTER_ACCOUNT = 'anthlis'

all_tweets_excl_rts_mentions = ' '.join([tw.text.lower() for tw in tweets
                                        if not tw.text.startswith('RT') and not tw.text.startswith('@')])

pb_mask = np.array(Image.open("twit_outline.png"))
stopwords = set(STOPWORDS)

stopwords.add('co')
stopwords.add('https')

wc = WordCloud(background_color="white", max_words=2000, mask=pb_mask,
               stopwords=stopwords)

wc.generate(all_tweets_excl_rts_mentions)

plt.figure(figsize=(15, 15))
plt.imshow(wc, interpolation="bilinear")
plt.margins(x=0, y=0)
plt.axis("off")


text:
from braket.circuits import circuit
import numpy

@circuit.subroutine(register=True)
def CUGate(theta, phi, lam, gam, dtype=None):
    theta, phi, lam, gam = float(theta), float(phi), float(lam), float(gam)
    cos = numpy.cos(theta/2)
    sin = numpy.sin(theta/2)
    exp1 = numpy.exp(1j*gam)
    exp2 = numpy.exp(1j*(gam+lam))
    exp3 = numpy.exp(1j*(gam+phi))
    exp4 = numpy.exp(1j*(gam+phi+lam))
    return numpy.array([[1, 0, 0, 0],
                        [0, 1, 0, 0],
                        [0, 0, exp1*cos, -exp2*sin],
                        [0, 0, exp3*sin, exp4*cos]], dtype=dtype)

max_stars_repo_name: Artificial-Brain/quantumcat

max_issues_repo_name: Artificial-Brain/quantumcat

max_forks_repo_name: Artificial-Brain/quantumcat


text:
#!/usr/bin/env python
# -*- coding: utf-8 -*-
from __future__ import (absolute_import, division,
                        print_function, unicode_literals)

__author__ = "Yuji Ikeda"

import numpy as np
try:
    from group.mathtools import lcm
except:
    from mathtools import lcm

class Group(object):
    def __init__(self, Cayley_table):
        """

        : 2-dim array-like
        """
        self.set_Cayley_table(Cayley_table)
        self._find_order_of_group()
        self._find_identity()

        self._create_conjugacy_classes()
        self._create_orders_of_elements()
        self._create_exponent()  # After creating orders_or_elements
        self.create_center()
        self._create_commutator_subgroup()
        self._create_crst()
        self._create_character_table()

    def set_Cayley_table(self, Cayley_table):
        self._Cayley_table = np.array(Cayley_table)

    def _find_order_of_group(self):
        self._order = self._Cayley_table.shape[0]

    def _find_identity(self):
        # TODO(ikeda): How about arbitrary symbols for the Cayley table?
        order = self._order
        Cayley_table = self._Cayley_table
        for i in range(order):
            is_identity = True
            for j in range(order):
                if Cayley_table[i, j] != j:
                    is_identity = False
                    break
            if is_identity:
                identity = i
                break
        self._identity = identity

    def _create_conjugacy_classes(self):
        Cayley_table = self._Cayley_table
        order = self._order

        conjugacy_classes = np.ones(order, dtype=int) * -1  # Initialized by -1
        count = 0
        for i in range(order):
            if conjugacy_classes[i] == -1:
                # "i" belongs to a new conjygacy class.
                conjugacy_classes[i] = count
                count += 1
            else:
                # The conjugacy class including "i" is already obtained.
                continue
            current_conjugacy_class = conjugacy_classes[i]
            for j in range(order):
                jinv = self.create_inverse_element(j)
                k = Cayley_table[Cayley_table[j, i], jinv]
                if conjugacy_classes[k] == -1:
                    conjugacy_classes[k] = current_conjugacy_class
                elif conjugacy_classes[k] != current_conjugacy_class:
                    raise ValueError(
                        "The conjugacy class of k is not equal to that of i.")
        self._conjugacy_classes = conjugacy_classes
        self._create_orders_of_conjugacy_classes()

    def _create_orders_of_conjugacy_classes(self):
        conjugacy_classes = self._conjugacy_classes
        self._orders_of_conjugacy_classes = np.bincount(conjugacy_classes)

    def _create_character_table(self):
        pass
        # crst = self._crst
        # for tmp_matrix in crst:
        #     # Get left eigenvectors
        #     eigvals, eigvecs = np.linalg.eig(tmp_matrix.T)
        #     eigvecs = eigvecs.T
        #     print("eigvals:")
        #     print(eigvals)
        #     print("eigvecs:")
        #     print(eigvecs)

    def create_center(self):
        """
        Create the center of the group G.

        The center is a normal abelian subgroup of G.
        """
        order = self._order
        Cayley_table = self._Cayley_table

        center = []
        for i in range(order):
            is_central = True
            for j in range(order):
                jinv = self.create_inverse_element(j)
                iconj = Cayley_table[Cayley_table[j, i], jinv]
                if i != iconj:
                    is_central = False
                    break
            if is_central:
                center.append(i)  # "i" is conjugate to all the group elements.
        self._center = np.array(center)

    def _create_crst(self):
        order = self._order
        Cayley_table = self._Cayley_table
        conjugacy_classes = self._conjugacy_classes

        num_classes = np.unique(conjugacy_classes).size
        crst = np.zeros((num_classes, num_classes, num_classes), dtype=int)
        for iclass in range(num_classes):
            # the representative
            i = np.where(conjugacy_classes == iclass)[0][0]
            for j in range(order):
                for k in range(order):
                    if Cayley_table[j][k] == i:
                        jclass = conjugacy_classes[j]
                        kclass = conjugacy_classes[k]
                        crst[jclass, kclass, iclass] += 1
        self._crst = crst

    def create_centralizer(self, subset):
        order = self._order
        Cayley_table = self._Cayley_table

        centralizer = []
        for i in range(order):
            is_in_centralizer = True
            iinv = self.create_inverse_element(i)
            for s in subset:
                sconj = Cayley_table[Cayley_table[i, s], iinv]
                if s != sconj:
                    is_in_centralizer = False
                    break
            if is_in_centralizer:
                centralizer.append(i)
        self._centralizer = np.array(centralizer)

    def create_normalizer(self, subset):
        order = self._order
        Cayley_table = self._Cayley_table

        normalizer = []
        for i in range(order):
            is_in_normalizer = True
            iinv = self.create_inverse_element(i)
            for s in subset:
                sconj = Cayley_table[Cayley_table[i, s], iinv]
                if sconj not in subset:
                    is_in_normalizer = False
                    break
            if is_in_normalizer:
                normalizer.append(i)
        self._normalizer = np.array(normalizer)

    def _create_commutator_subgroup(self):
        """

        https://en.wikipedia.org/wiki/Commutator_subgroup
        """
        order = self._order
        commutator_subgroup = []
        for i in range(order):
            for j in range(order):
                k = self.create_commutator(i, j)
                commutator_subgroup.append(k)
        self._commutator_subgroup = np.unique(commutator_subgroup)

    def create_commutator(self, i, j):
        Cayley_table = self._Cayley_table
        i_inv = self.create_inverse_element(i)
        j_inv = self.create_inverse_element(j)
        k = Cayley_table[i_inv, Cayley_table[j_inv, Cayley_table[i, j]]]
        return k

    def create_inverse_element(self, element):
        Cayley_table = self._Cayley_table
        identity = self._identity
        inverse = np.where(Cayley_table[element] == identity)[0][0]
        return inverse

    def _create_orders_of_elements(self):
        order = self._order
        orders_of_elements = np.zeros(order, dtype=int)
        for i in range(order):
            orders_of_elements[i] = self.create_order_of_element(i)
        self._orders_of_elements = orders_of_elements

    def create_order_of_element(self, i):
        Cayley_table = self._Cayley_table
        identity = self._identity
        j = i
        count = 0
        while True:
            count += 1
            if j == identity:
                order_of_element = count
                return order_of_element
            else:
                j = Cayley_table[j, i]

    def _create_exponent(self):
        orders_of_elements = self._orders_of_elements
        exponent = lcm(*orders_of_elements)
        self._exponent = exponent

    def get_identity(self):
        return self._identity

    def get_order_of_group(self):
        return self._order

    def get_orders_of_elements(self):
        return self._orders_of_elements

    def get_exponent(self):
        return self._exponent

    def get_conjugacy_classes(self):
        return self._conjugacy_classes

    def get_orders_of_conjugacy_classes(self):
        return self._orders_of_conjugacy_classes

    def get_center(self):
        return self._center

    def get_centralizer(self):
        return self._centralizer

    def get_normalizer(self):
        return self._normalizer

    def get_commutator_subgroup(self):
        return self._commutator_subgroup

    def get_crst(self):
        return self._crst

max_stars_repo_name: sheikhahnaf/upho

max_issues_repo_name: sheikhahnaf/upho

max_forks_repo_name: sheikhahnaf/upho


text:
import dxtbx
import numpy as np
from scipy.ndimage.filters import median_filter
import matplotlib.pyplot as plt
from matplotlib.gridspec import GridSpec
from scitbx.matrix import sqr, col
import scitbx.matrix as matrix

def radial_cbf_smooth(filename):
    frame = dxtbx.load(filename)
    data = frame.get_raw_data().as_numpy_array()
    detector = frame.get_detector()
    beam = frame.get_beam()

    for panel in detector:
        s0 = -1 * matrix.col(beam.get_direction())
        beam_center = col(panel.get_beam_centre_px(s0))
        pix_mm = panel.get_pixel_size()
        origin = matrix.col(panel.get_origin())

    detector_distance = s0.dot(origin)

    y, x = np.indices((data.shape))
    r = np.sqrt((x - beam_center[0])**2 + (y - beam_center[1])**2)
    r = r.astype(np.int)
    data_mask = np.array(data, dtype=bool)
    data_mask[data<1e-6]=False
    r=r[data_mask]
    data=data[data_mask]
    tbin = np.bincount(r.ravel(), data.ravel())
    nr = np.bincount(r.ravel())
    with np.errstate(invalid='ignore', divide='ignore'):
        radial_avg = tbin / nr
    radial_avg_smooth = median_filter(radial_avg, 50)
    print("detector distance = {}".format(detector_distance))
    # plt.plot(radial_avg_smooth)
    # plt.savefig("smoothed_rad_avg.png", dpi=300)

    return radial_avg

# from sys import argv

# script, filename = argv
# radial_cbf_smooth(filename)

reason: import numpy,from scipy


text:
import numpy as np
import random
from torch.utils.data.dataset import Dataset

from cutmix.utils import rand_bbox

from imagecorruptions import corrupt
from PIL import Image

import torchvision.datasets as datasets

class CutMix(Dataset):
    def __init__(self, dataset, num_class, num_mix=1, beta=1., prob=1.0, transform=None):
        self.dataset = dataset
        self.num_class = num_class
        self.num_mix = num_mix
        self.beta = beta
        self.prob = prob
        self.transform = transform

    def __getitem__(self, index):
        img, lb = self.dataset[index]

        # i have a transformed image which has fixed dimensions
        img = np.asarray(img)
        normal_img = Image.fromarray(np.copy(img))

        # generate mixed sample
        lam = np.random.beta(self.beta, self.beta)
        rand_index = random.choice(range(len(self)))
        img2, lb2 = self.dataset[rand_index]
        while lb2 == lb:
            rand_index = (rand_index + 1) % len(self)
            img2, lb2 = self.dataset[rand_index]
        img2 = np.asarray(img2)
        bbx1, bby1, bbx2, bby2 = rand_bbox(img.shape, lam)
        img[bbx1:bbx2, bby1:bby2, :] = img2[bbx1:bbx2, bby1:bby2, :]

        # corrupt image here
        corr_number = np.random.choice(15)
        severity_index = np.random.choice(5)
        corr_img = corrupt(img, severity=severity_index+1, corruption_number=corr_number)
        corr_img = Image.fromarray(corr_img)

        if self.transform:
            normal_img, corr_img = self.transform[0](normal_img), self.transform[1](corr_img)

        return normal_img, lb, corr_img, self.num_class

    def __len__(self):
        return len(self.dataset)

if __name__ == "__main__":
    traindir = "/home/jatin/Imagenet/ILSVRC/Data/CLS-LOC/train"
    train_set = datasets.ImageFolder(traindir)

    cnc_set = CutMix(train_set, 1000)

    img, lb, corr, lb2 = cnc_set[256]

    corr.save("cnc.jpg")
    img.save("normal.jpg")


text:
# Copyright 2021 Huawei Technologies Co., Ltd
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
# ============================================================================

""" export model for CRNN """
import os
import numpy as np
import mindspore as ms
from mindspore import Tensor, context, load_checkpoint, export
from src.crnn import crnn
from src.model_utils.moxing_adapter import moxing_wrapper
from src.model_utils.config import config
from src.model_utils.device_adapter import get_device_id

context.set_context(mode=context.GRAPH_MODE, device_target=config.device_target, save_graphs=False)

def modelarts_pre_process():
    config.file_name = os.path.join(config.output_path, config.file_name)

@moxing_wrapper(pre_process=modelarts_pre_process)
def model_export():
    if config.device_target == "Ascend":
        context.set_context(device_id=get_device_id())

    config.batch_size = 1
    net = crnn(config, full_precision=config.device_target != 'Ascend')

    load_checkpoint(config.ckpt_file, net=net)
    net.set_train(False)

    input_data = Tensor(np.zeros([1, 3, config.image_height, config.image_width]), ms.float32)

    export(net, input_data, file_name=config.file_name, file_format=config.file_format)

if __name__ == '__main__':
    model_export()

max_stars_repo_name: mindspore-ai/models

max_issues_repo_name: mindspore-ai/models

max_forks_repo_name: mindspore-ai/models


text:
# Accepted, but disappointingly only faster than 5% of Python submissions.
import numpy as np
import operator as op
from functools import reduce
from typing import List
from algorith.leetc.trees.tree_node import TreeNode

# Definition for a binary tree node.
# class TreeNode:
#     def __init__(self, val=0, left=None, right=None):
#         self.val = val
#         self.left = left
#         self.right = right

class TreeNode:
    def __init__(self, key):
        self.key = int(key)
        self.val = int(key)
        self.right = None
        self.left = None

class Solution:
    def generateTrees(self, n: int) -> List[TreeNode]:
        return print_all_trees(n)

def ncr(n, r):
    if r < 0:
        return 0
    r = min(r, n-r)
    numer = reduce(op.mul, range(n, n-r, -1), 1)
    denom = reduce(op.mul, range(1, r+1), 1)
    return numer // denom

def n_catalan2(h, t):
    if t < 0:
        return 0
    elif t == 0:
        return 1
    else:
        return abs(ncr(h+t-1, t) - ncr(h+t-1, t-1))

def ix_to_path(ix, h, t):
    """
    Given an index, returns the Dyck path.
    args:
        ix: The index to be converted.
        t: The number of tails in the path.
        h: The number of heads in the path.
    """
    v = ix
    path = np.zeros(t + h)
    for i in range(t+h):
        tail_paths = n_catalan2(h, t-1)
        if v >= tail_paths:
            path[i] = 1
            h -= 1
            v -= tail_paths
            if h == 0:
                break
        else:
            t -= 1
    return path

def get_tree(path, start, end):
    if end <= start:
        return None
    targ = path[start]
    x = find_tar(path, targ, start)
    root = TreeNode((x-start)/2)
    root.left = get_tree(path, start+1, x-1)
    root.right = get_tree(path, x, end)
    # Fix the labels of the keys.
    return root

def get_tree2(path):
    root = get_tree(path, 0, len(path)-1)
    Tr(root)
    return root

class Tr():
    def __init__(self, x):
        self.i = 1
        self.root = self.label_inorder(x)

    def label_inorder(self, x):
        if x is not None:
            self.label_inorder(x.left)
            x.key = self.i
            x.val = self.i
            self.i += 1
            self.label_inorder(x.right)

class Tr1():
    def __init__(self, root):
        self.root = root
        self.arr = []
        self.tree_to_array(root)

    def tree_to_array(self, tn):
        """
        Converts tree to leetcode array.
        """
        if tn is not None:
            self.arr.append(tn.key)
        else:
            self.arr.append(None)
        if tn is not None and\
           (tn.left is not None or tn.right is not None):
            self.tree_to_array(tn.left)
            self.tree_to_array(tn.right)

def print_all_trees(n=3):
    n_trees = n_catalan2(n, n+1)
    trees = []
    for ix in range(n_trees):
        tn = ix_to_tree(ix, n)
        trees.append(tn)
    return trees

def ix_to_tree(ix, n):
    """
    Converts index to lexicographic binary search tree.
    args:
        ix: The index.
        n: The number of nodes in the tree.
    """
    path = ix_to_path(ix, n, n+1)
    # The last entry is a redundant tail. So, remove it.
    tree_nd = print_tree(path[:len(path)-1])
    return tree_nd

def print_tree(path=[1, 1, 1, 0, 0, 0]):
    for i in range(len(path)):
        if path[i] == 0:
            path[i] = -1
    path = np.concatenate(([0], np.cumsum(path)))
    tn = get_tree2(path)
    return tn

def find_tar(path, targ, start):
    i = start+1
    while path[i] != targ:
        i += 1
    return i

if __name__=="__main__":
    sol = Solution()
    sol.generateTrees(3)


text:
# -*- coding: utf-8 -*-
"""
Created on Mon Apr 16 21:04:56 2018

@author: Poul
"""
#here the first table code
import random
from random import randint
import numpy
import pandas
import math
"""
I have added some comments here and there to help you understand the code. I hope it will be alright, sorry for the mess !
CRUCIAL HOW TO RUN :
1) First run the programm until the line with "checkpoint table 1" (around line 100)
3) complete carefully what the programm asks, this will influence a lot what the output is
2) you will get an excel in which you have the different mutants, in this excel put the experimental replicates values instead of the "X"
3) then run the rest of the program, you will get a second excel with the results !

"""
mutations_list = []
replicate_list = []
mutant_list = []
value_list=[]
your_study = input("Do you use selectivity or conversion values? Please answer with S (Selectivity) or C (Conversion): ")
mutation_number = int(input("Please indicate your mutation number: "))
replicate_number = int(input("Please indicate your replicate number (if some replicates are faulty, please fill the table with the average of the others otherwise the program might give unexpected results) : "))
your_data = input("Please enter the name of your replicate table (don't forget the file extension !): ")
your_data2 = input("Please enter the name of the file you want your results in (don't forget the file extension !): ")
#very important lines that determine a lot the output of the code. This gives flexibility to the code and intercations with the user.
for elt3 in range(1, replicate_number+1):
    replicate_list.append("Replicate n°%s" % (elt3))
for elt2 in range(1,mutation_number+1):
    mutations_list.append(input("Please indicate the mutation n°%s: " % (elt2)))
for elt4 in range(1,2**mutation_number+1):
    mutant_list.append("Mutant %s" % (elt4)) #these lines are very imortant to make a list of the mutations, a list of the replicates names and the mutants names. They will be used to make the table and combinations.

def create_combination(mutation_number):#this function creates the mutant combinations based on the number you indicated
    dic_list = []
    while len(mutant_list)>len(dic_list):
        for elt in mutant_list:
            elt = {}
            for elt2 in range(1, mutation_number+1): #here we attribute a number for + and - and roll the dice to obtain a random combination under the form of a dictionary !
                evolution_dice = random.randint(0,1)
                if evolution_dice == 0:
                    elt[mutations_list[elt2-1]] = "+"
                else:
                    elt[mutations_list[elt2-1]] = "-"
            count = 0
            for elt3 in dic_list:
                if elt.items() != elt3.items(): #this line will scan each combination of the list and compare it to the new combination
                        count+= 1
            if count == len(dic_list):
                dic_list.append(elt) #we add this combination to a new list. If this combination is already in the list, then we thrash it and do it again

    return dic_list

Mutant_number = len(create_combination(mutation_number)) #the number of mutants if equal to the number of combinations
box = ["X"]*Mutant_number*(mutation_number+replicate_number) #here we make a number of empty cases filled with "X"proportional to the number of mutants and mutations
final_table1 = numpy.reshape(box,((Mutant_number),(mutation_number+len(replicate_list))))  #here we make a matrix with the boxes above with the number of rows being equal to the number of mutants and the number of column being equal to the number of mutations + the number of replicates

for mutant in create_combination(mutation_number):
    sign_list = []
    for combinations in mutant:
        sign_list.append(mutant[combinations])
    value_list.append(sign_list) #this creates a list of the the combinations

def value_list_sorter(value_list): #this put the combinations of signs together based on the number of + they have
    sorted_values = []
    final = []
    ref_list = []
    count = 0
    value_dic = {}
    for value in value_list:
        count = value.count("+")
        ref_list.append(count)
        value_dic[str(value)]=count
    ref_list.sort()
    for num in ref_list:
        for item in value_dic.items():
            if num == item[1]:
                sorted_values.append(item[0])
                del value_dic[item[0]]
                break
    for elt in sorted_values:
        for elt2 in value_list:
            if elt == str(elt2):
                final.append(elt2)
    return final
final_value_list = value_list_sorter(value_list) #this is our sorted sign list !
def table_filler(final_table1, final_value_list, mutation_number):#this will fill the matrix with our ordered sign list
    i = 0
    while i < mutation_number:
        j = 0
        while j < len(mutant_list):
            final_table1[j][i] = final_value_list[j][i]
            j = j + 1
        i = i + 1
    return final_table1

excel_table1 = pandas.DataFrame(table_filler(final_table1, final_value_list, mutation_number), columns=mutations_list+replicate_list, index=mutant_list)
writer = pandas.ExcelWriter(your_data)
excel_table1.to_excel(writer, sheet_name="sheet_name", index = True)   #finally we write everything on a new excel, of which the name is given by the user
writer.close()
print("checkpoint after table1--------------------------------------------")

#from here the second table code

table1 = pandas.read_excel(your_data) #here we take the input given by the user in the excel tqble we previously generated
foundment = table1.iloc[:,:mutation_number]#foundment is the matrix of all the signs and mutants we obtained above
total_replicate_data = table1.iloc[:, mutation_number:]
replicate_matrix = total_replicate_data.values #this isolates the values that the user put in the input excel file. That way we obtain a matrix with the number of lines = number of mutants and number of columns = number of replicates
data_array = table1.values#This is all the data (signs/mutants and associated replicates)

def mean_and_sd_maker(data_array): #this function will look into the vqlues of each mutants and make an average and standard deviation out of it. In the final table those are called "experimental average" and "experimental standard deviation"
    data_dic = {}
    mean_and_sd = []
    for array in data_array:
        data = array[mutation_number:]
        data_float = numpy.array(data).astype(numpy.float64)
        mutant = str(array[:mutation_number])
        average = float(numpy.average(data_float))
        std = float(numpy.std(data_float))/math.sqrt(replicate_number)
        data_dic[mutant] = [average,std]
    for elt in data_dic.values():
        mean_and_sd.append(elt)
    return data_dic, mean_and_sd

#This function gives a tuple (dictionary of mutants associated with mean and std, array of mean and std)
mean_and_sd_dic = mean_and_sd_maker(data_array)[0] #here we just take the first element of the tuple, which is the dictionarry. I frankly don't even remember why I did a tuple and not just the dictionary but hey)
mean_and_sd_array = numpy.reshape(mean_and_sd_maker(data_array)[1], ((Mutant_number),2))
foundment_values = foundment.values #this line is to obtain readable values from the variable foundment, which will then be a matrix of signs.

def origin_finder(foundment_values): #this is the first function that will permit to find possible combinations between mutqnts. This one is useful to find double mutqnts. For exqmple [+ - + -] and [- + - +].
    additivity_list = []
    foundment_values[foundment_values == "+"] = 1 #here I change the + and - for 1 and 0. This is useful for calculations
    foundment_values[foundment_values == "-"] = 0
    i= 1
    while i < len(foundment_values)-1: #I go through the sign mqtrix
        j = i
        while j < len(foundment_values)-1:#and a second time, so I cqn isolqte two combinqtions qt q time qnd compare them
            res = foundment_values[i] + foundment_values[j+1] #so here we hqve this vqriqble "res". For example if the two combinations are [+ - +] and [- + -], res will be [1,1,1]. However, if the combinations are [+ + -] and [+ - +], res will be [2, 1, 1]
            for array in foundment_values: #we tqke this res and compare it to the mutants we have
                if numpy.array_equal(res, array) == True: #if res is equal to one of the mutant, we have found a combination !
                    additivity_list.append((list(res),(i+1,j+2))) #here we write the combination in a tuple with the combination and what mutants form it
            j = j + 1
        i = i + 1
    return additivity_list

origins = origin_finder(foundment_values)
def please_more_combinations(origins,foundment_values, mutations_list): #now is probably the trickiest function I had to do. The code above works for double mutants but not for triple, quadruple etc... The idea is that I use recurcivity to obtain new combinations
    final_comb_list = origins #we retake the combination list we obtain in the previous function
    cycle_number = len(mutations_list) #here we define the number of cycle the function will do. Everytime we have a new mutation the nuber of cycle increases by 1
    comb_comparator = []
    if cycle_number>1: #that is the recursivity condition. The function will stop after the number of cycles is down to one
        for comb in final_comb_list:#so the idea is to scan the comb list we obtained above. In that case we can make combinations of combinations to obtain more combinations !
            for array2 in foundment_values[1:]:
                res2 = numpy.array(comb[0] + array2)
                for array3 in foundment_values:
                    if numpy.array_equal(res2, array3) == True: #same principle as above
                        new_comb = list(comb[1])
                        new_comb.append(foundment_values.tolist().index(array2.tolist())+1)
                        count = 0
                        for elt in final_comb_list:
                            a_comb = list(elt[1])
                            a_comb.sort()
                            comb_comparator.append(a_comb)
                            new_comb.sort()
                        for elt2 in comb_comparator:#those lines make sure the newly formed combination has not been already made. This is probably a litle more complicated but I must admit I don't really recall everything
                            if elt2 == new_comb:
                                count += 1
                        if count == 0:
                            final_comb_list.append((list(res2), tuple(new_comb)))
        cycle_number = cycle_number - 1#we delete one cycle and go on
    else:
        please_more_combinations(final_comb_list)#and we repeat the function for the final cycle but with the final list as a variable
    return final_comb_list

all_combinations = please_more_combinations(origins,foundment_values, mutations_list)

#here will be made the combinations table
count_list = []
for elt in all_combinations:
    count_list.append((elt[0]).count(1))
count_list.sort() #this is just a variable coresponding to the number of combinations
ordered_combs = []
for elt in count_list:
    for elt2 in all_combinations:
        if list(elt2[0]).count(1) == elt:
            all_combinations.remove(elt2)
            ordered_combs.append(elt2)
#I think this was to remove any potential duplicate of combinations that somehow ended up in the list

comb_index=[]
for elt in range(1,len(ordered_combs)+1):
    comb_index.append("Combination n°%s" % (elt))
#this line is important for the final table, it gives a proper name to each combination

combs_only = []
for elt in ordered_combs:
    combs_only.append(elt[1])
#this gives a list of the mutant combinations only

signs_only = []
for elt in ordered_combs:
    signs_only.append(elt[0])
#same as above but for the signs only

reshaped_signs = numpy.reshape(signs_only, ((len(signs_only),(len(mutations_list)))))
reshaped_combs = numpy.reshape(combs_only, ((len(signs_only), 1)))
#reshqping everything to have a god format for the final table

final_comb_table = numpy.c_[reshaped_signs,reshaped_combs]
final_comb_table[final_comb_table == 1] = "+"
final_comb_table[final_comb_table == 0] = "-"
foundment_values[foundment_values == 1] = "+"
foundment_values[foundment_values == 0] = "-" #reconverting all 1 and 0 into + and -
foundment_values = numpy.c_[foundment_values, mean_and_sd_array] #we also add the averages and standard deviation (experimental) to the sign matrix

def theoretical_stats_selectivity(mean_and_sd_dic, final_comb_table, combs_only, mean_and_sd_array):
    grand_final = []
    all_of_it = []
    for elt in final_comb_table:
        for elt2 in mean_and_sd_dic.keys():
            if str(elt[:mutation_number]) == str(elt2):
                elt = numpy.append(elt, list(mean_and_sd_dic[elt2]))
        for elt3 in combs_only:
            if numpy.array_equal(elt[len(mutations_list)], elt3) == True:
                theor_mean = numpy.array([0])
                replicate_values = numpy.zeros((1,len(replicate_matrix[0])))
                for elt4 in elt3:
                    target = mean_and_sd_array[elt4-1][0]
                    theor_mean = numpy.add(theor_mean, target)
                    target2 =replicate_matrix[elt4-1]
                    replicate_values = numpy.add(replicate_values,target2)
                theor_sd = (numpy.std(replicate_values))/math.sqrt(replicate_number)
                elt =  numpy.append(elt, list(theor_mean))
                elt = numpy.append(elt, theor_sd)
                grand_final.append(elt)
    for elt5 in grand_final:
        at_last = (elt5[len(mutations_list)+1:][0])-(elt5[len(mutations_list)+1:][2])
        elt5 =  numpy.append(elt5, at_last)
        all_of_it.append(elt5)

    return numpy.array(all_of_it)
#the function above calculates the theoretical average and standard deviations based on the article that Carlos and his colleagues has written. This is for selectivity values

def theoretical_stats_conversion(mean_and_sd_dic, final_comb_table, combs_only, mean_and_sd_array):
    grand_final = []
    all_of_it = []
    keys = list(mean_and_sd_dic.keys())
    WT = keys[0]
    avgWT = mean_and_sd_dic[WT][0]
    for elt in final_comb_table:
        for elt2 in mean_and_sd_dic.keys():
            if str(elt[:mutation_number]) == str(elt2):
                elt = numpy.append(elt, list(mean_and_sd_dic[elt2]))
        for elt3 in combs_only:
            if numpy.array_equal(elt[len(mutations_list)], elt3) == True:
                theor_mean = numpy.array([0])
                replicate_values = numpy.zeros((1,len(replicate_matrix[0])))
                for elt4 in elt3:
                    new_target = []
                    target = mean_and_sd_array[elt4-1][0]-avgWT
                    theor_mean = numpy.add(theor_mean, target)
                    target2 =replicate_matrix[elt4-1]
                    for value in target2:
                        value = value - avgWT
                        new_target.append(value)
                    replicate_values = numpy.add(replicate_values,new_target)
                good_one = list(theor_mean)[0]
                good_one = avgWT+good_one
                theor_sd = (numpy.std(replicate_values))/math.sqrt(replicate_number)
                elt =  numpy.append(elt, good_one)
                elt = numpy.append(elt, theor_sd)
                grand_final.append(elt)
    for elt5 in grand_final:
        at_last = (elt5[len(mutations_list)+1:][0])-(elt5[len(mutations_list)+1:][2])
        elt5 =  numpy.append(elt5, at_last)
        all_of_it.append(elt5)

    return numpy.array(all_of_it)

#this time for conversion, which is a little different albeit very close.

if your_study == "S": #at the very beginning of the code we ask for "selectivity" or "conversion" from the USER so that the program actually adapt to what the user wants.
    all_of_it = theoretical_stats_selectivity(mean_and_sd_dic, final_comb_table, combs_only, mean_and_sd_array)
elif your_study == "C":
    all_of_it = theoretical_stats_conversion(mean_and_sd_dic, final_comb_table, combs_only, mean_and_sd_array)

def what_epistasis_sign_selectivity(all_of_it, replicate_matrix):
    sign = []
    epi_list = []
    what_epi = []
    i = 0
    for elt in all_of_it:
        noinspi = elt[len(mutations_list)+1:]
        Gexp = noinspi[0]
        Gexpstd = noinspi[1]
        Gcomb = noinspi[2]
        Gcombstd = noinspi[3]
        GexpES = Gexp - Gexpstd
        GcombES = Gcomb + Gcombstd
        GexpES2 = Gexp + Gexpstd
        GcombES2 = Gcomb - Gcombstd
        if GexpES < GcombES and Gexp > Gcomb:
            sign.append("Additive")
        elif GexpES2 > GcombES2 and Gexp < Gcomb :
            sign.append("Additive")
        elif Gexp < Gcomb :
            sign.append("- ")
        elif Gexp > Gcomb:
            sign.append("+ ")
    for elt2 in combs_only:
        combavg = []
        for lign in all_of_it:
            if lign[len(mutations_list)] == elt2:
                double_mutant_avg = lign[len(mutations_list)+1]
        for elt3 in elt2:
            mutant_avg = replicate_matrix[elt3-1]
            mutant_avg = float(numpy.average(mutant_avg))
            combavg.append(mutant_avg)
        count = 0
        for avg in combavg:
            if avg < 0:
                count = count - 1
            elif avg > 0 :
                count = count + 1
        if abs(count) == len(combavg):
            if count > 0 and double_mutant_avg > 0 or count < 0 and double_mutant_avg < 0:
                epi_list.append("Magnitude epistasis")
            elif count > 0 and double_mutant_avg < 0 or count < 0 and double_mutant_avg > 0:
                epi_list.append("Reciprocal sign epistasis")
        elif abs(count) != len(combavg):
            epi_list.append("Sign epistasis")

    while i < len(sign):
        if sign[i] != "Additive":
            what_epi.append(sign[i] + epi_list[i])
        else:
            what_epi.append(sign[i])
        i = i + 1
    return what_epi
#finally the great last function that also uses Carlos'equations to determine the nature of epistasis.
def what_epistasis_sign_conversion(all_of_it, replicate_matrix):
    sign = []
    epi_list = []
    what_epi = []
    i = 0
    keys = list(mean_and_sd_dic.keys())
    WT = keys[0]
    avgWT = mean_and_sd_dic[WT][0]
    for elt in all_of_it:
        noinspi = elt[len(mutations_list)+1:]
        Gexp = noinspi[0] - avgWT
        Gexpstd = noinspi[1]
        Gcomb = noinspi[2]- avgWT
        Gcombstd = noinspi[3]
        GexpES = Gexp - Gexpstd
        GcombES = Gcomb + Gcombstd
        GexpES2 = Gexp + Gexpstd
        GcombES2 = Gcomb - Gcombstd
        if GexpES < GcombES and Gexp > Gcomb:
            sign.append("Additive")
        elif GexpES2 > GcombES2 and Gexp < Gcomb :
            sign.append("Additive")
        elif Gexp < Gcomb :
            sign.append("- ")
        elif Gexp > Gcomb:
            sign.append("+ ")
    for elt2 in combs_only:
        combavg = []
        for lign in all_of_it:
            if lign[len(mutations_list)] == elt2:
                double_mutant_avg = lign[len(mutations_list)+1]
        for elt3 in elt2:
            mutant_avg = replicate_matrix[elt3-1] - avgWT
            mutant_avg = float(numpy.average(mutant_avg))
            combavg.append(mutant_avg)
        count = 0
        for avg in combavg:
            if avg < 0:
                count = count - 1
            elif avg > 0 :
                count = count + 1
        if abs(count) == len(combavg):
            if count > 0 and double_mutant_avg > 0 or count < 0 and double_mutant_avg < 0:
                epi_list.append("Magnitude epistasis")
            elif count > 0 and double_mutant_avg < 0 or count < 0 and double_mutant_avg > 0:
                epi_list.append("Reciprocal sign epistasis")
        elif abs(count) != len(combavg):
            epi_list.append("Sign epistasis")

    while i < len(sign):
        if sign[i] != "Additive":
            what_epi.append(sign[i] + epi_list[i])
        else:
            what_epi.append(sign[i])
        i = i + 1
    return what_epi
#similarily to before, conversion and sleectivity differ a little so I had to adapt the code
if your_study == "C":#same logic as before regarding adaptation of the code to what the user wants
    epistasis = what_epistasis_sign_conversion(all_of_it, replicate_matrix)
elif your_study == "S":
    epistasis = what_epistasis_sign_selectivity(all_of_it, replicate_matrix)
all_of_it = numpy.c_[all_of_it,epistasis]
#this all_of_it value is all the data we need, across the program we complete it as it goes

suppinfo = ["Combinations", "Experimental average", "Experimental standard deviation", "Theoretical average", "Theoretical standard deviation", "Exp.avg - Theor.avg", "Epistasis type"]

excel_table2= pandas.DataFrame(all_of_it, columns=mutations_list+suppinfo, index=comb_index)
excel_table3= pandas.DataFrame(foundment_values, columns=mutations_list+["Average", "Standard deviation"], index=mutant_list)
writer2 = pandas.ExcelWriter(your_data2)
excel_table2.to_excel(writer2, sheet_name="Theoretical results table", index = True)
excel_table3.to_excel(writer2, sheet_name="Experimental results table", index = True)
writer2.close()
#and here are the lines to write the final excel table ! THe final file has two sheet, one with all the values and combinations, and the other with the experimental values only and the single mutants.