text:
"""
gc_maps.py contains functions useful for running larger parameter scans and
also for analysing their results. Results of these scans are saved to
'/map_data', as their format is more compressed than the complete simulation
data stored in '/raw_data'.
"""

from __future__ import division
from itertools import product

import numpy as np
from os import listdir
import pandas as pd
import pylab

import seaborn
import sys
import time
import matplotlib.pyplot as plt
from matplotlib import rc
from matplotlib.pyplot import cm
import matplotlib as mpl

import cf
import gc_memo

reload(cf)
reload(gc_memo)

from gc_memo import main

# general plot settings
cc = 30  # number of clours for cycling through
cmap = seaborn.color_palette('Spectral', cc)
colors = iter(cm.rainbow(np.linspace(0, 1, cc)))
seaborn.set_style('ticks')
seaborn.set_context('talk')

plt.rc('text', usetex=True)
rc('font', **{'family': 'sans-serif', 'sans-serif': ['Helvetica'],
              'style': 'normal'})
rc("text.latex", preamble=["\\usepackage{helvet}\\usepackage{sfmath}"])
pylab.ion()

def map_params(dose=[1], LFdecay=[10*12], nGCs=[1], nLFs=[25],
               naive_pool=[1000], nkey=[1, 2, 10], p_err=[0.003],
               tinf=[[0*12]],
               p_block=[0.5], repeats=1):

    """ Function for mapping out the effects of different parameter
    (combinations) on the standard TUCHMI protocol. All arguments are lists,
    either containing only the default value or a set of values, in which case
    all combinations of list arguments will be executed the given number
    of times (repeat).

    The paramter set used is stored together with the complete timecourse of
    of E_bind, SHM, entropy, memory number (means and std where applicable)
    and exported into a .h5 file.

    For running on a cluster, accepts an ID argument (e.g. job ID) for easier
    handling of errors etc.
    """

    # open file stamped with systemtime if no other ID was provided in the call
    try:
        sys.argv[1]
    except IndexError:
        filepath = 'map_data/data{}.h5'.format(int(time.time()*100))
    else:
        filepath = 'map_data/data{}.h5'.format(sys.argv[1])
    print(filepath)
    datafile = pd.HDFStore(filepath)
    # dict for collecting result series
    seriesdict = {}

    # set endtime and days for evaluation
    cf.endtime = 126*12  # run until challenge timepoint
    # evaluate pool every day
    evaldays = np.arange(126)

    # get parameter combinations
    paramsets = list(product(dose, LFdecay, nGCs, nLFs, naive_pool, nkey,
                             p_err, tinf, p_block))
    # for every parameter set, run the simulation repeat times and write
    # results to the file
    for p in paramsets:
        # set parameters
        cf.dose = [p[0] for i in range(len(p[7]))]
        cf.LFdecay = p[1]
        cf.nGCs = p[2]
        cf.nLFs = p[3]
        cf.naive_pool = p[4]*cf.nGCs
        cf.memory_pool = 100*cf.nGCs  # fixed! (or change manually)
        cf.nkey = p[5]
        cf.lAg = p[5]
        cf.lAb = 220 - p[5]
        cf.p_err_FWR = p[6]
        cf.p_err_CDR = p[6]
        cf.tinf = p[7]
        cf.p_block_FWR = p[8]

        for r in range(repeats):
            # get lists to store individual simulation results
            l_mems = []
            l_KDs = []
            s_KDs = []  # std
            l_SHMs = []
            s_SHMs = []  # std
            l_Entrs = []
            # run simulation and get filepath or dict
            evaldays, l_mems, l_KDs, s_KDs, l_SHMs, s_SHMs, l_Entrs = \
                main(store_export='minimal', evalperday=1)
            # write these lists to file together with parameters used.
            # Identifier system time.
            ID = 'ID_{}'.format(time.time())
            series = pd.Series([cf.dose, cf.LFdecay, cf.nGCs, cf.nLFs,
                                cf.naive_pool, cf.nkey, cf.p_block_FWR,
                                cf.p_err_CDR, cf.tinf, cf.memory_pool,
                                evaldays, np.array(l_mems),
                                np.array(l_KDs), np.array(s_KDs),
                                np.array(l_SHMs), np.array(s_SHMs),
                                np.array(l_Entrs)],
                               index=['dose', 'LFdecay', 'nGCs', 'nLFs',
                                      'naive_pool', 'nkey', 'p_block',
                                      'p_err', 'tinf', 'mem_pool',
                                      'evaldays', 'memcount', 'E_bind',
                                      'E_bind_std', 'SHM', 'SHM_std',
                                      'entropy'])
            seriesdict[ID] = series
    # make dataframe and store it
    df = pd.DataFrame(seriesdict)
    df = df.transpose()
    datafile['data'] = df
    # close datafile
    datafile.close()
    print('END')

def map_import(filepathlist):
    """ Imports data from all given filepaths and provides quick and dirty
    plots for testing purposes. Give filepaths as strings, please."""
    # import all frames and concatenate them
    frames = []
    for path in filepathlist:
        frames.append(pd.read_hdf(path, key='data'))
    df = pd.concat(frames)

    # plot affinities over time
    fig, ((ax0, ax1), (ax2, ax3)) = plt.subplots(2, 2, sharex=True,
                                                 figsize=(15, 10))
    lines = []
    labels = []
    ax0.set_title('Normalised binding energy')
    ax0.set_ylim([0.6, 1])
    for index, row in df.iterrows():
        l, = ax0.plot(row['evaldays'], row['E_bind'])
        ax0.fill_between(row['evaldays'], row['E_bind']-row['E_bind_std'],
                         row['E_bind']+row['E_bind_std'], color='lightgrey')
        lines.append(l)
        labels.append('dose={}, LFdecay={}, nGCs={}, nLFs={}, naives={}, '
                      'nkey={}, pblock={}'.format(row['dose'], row['LFdecay'],
                                                  row['nGCs'], row['nLFs'],
                                                  row['naive_pool'],
                                                  row['nkey'],
                                                  row['p_block']))

    ax1.set_title('SHM counts')
    ax1.set_ylim([0, 70])
    for index, row in df.iterrows():
        l, = ax1.plot(row['evaldays'], row['SHM'])
        ax1.fill_between(row['evaldays'], row['SHM']-row['SHM_std'],
                         row['SHM']+row['SHM_std'], color='lightgrey')

    ax2.set_title('Memory cell number')
    for index, row in df.iterrows():
        l, = ax2.plot(row['evaldays'], row['memcount'])

    ax3.set_title('Entropy')
    for index, row in df.iterrows():
        l, = ax3.plot(row['evaldays'], row['entropy'])

    plt.legend(lines, labels, loc=(-0.5, -0.5))
    plt.tight_layout()
    seaborn.despine()
    pylab.savefig('figures/overview.pdf', bbox_inches='tight')

def data_import(folder='map_data'):
    """ Imports all simulation results found in the given folder.

    Prints all values that where found for each parameter as well as a list
    of parameter combinations and the number of simulation results available
    for each of them.

    Args:
        folder (string): Name of the folder to parse for data

    Returns:
        pandas datafile with the concatenated data.
    """

    # get all filenames in the data folder
    filenames = listdir(folder)
    # import and merge into one dataframe
    frames = []
    for name in filenames:
        if name.startswith('data'):
            frames.append(pd.read_hdf(folder+'/'+name, key='data'))
    df = pd.concat(frames)

    # replace values which contain lists by proxies as this is not useful
    df.dose = df.dose.apply(np.nanmean)
    df.dose = df.dose.apply(np.round, decimals=3)
    df.tinf = df.tinf.apply(len)
    # replace overall naive number by number per GC
    df.naive_pool = df.naive_pool.divide(df.nGCs)

    # find unique parameter values and print them
    print('LFdecay: {}'.format(df.LFdecay.unique()))
    print('nGCs: {}'.format(df.nGCs.unique()))
    print('nLFs: {}'.format(df.nLFs.unique()))
    print('naive_pool: {}'.format(df.naive_pool.unique()))
    print('nkey: {}'.format(df.nkey.unique()))
    print('p_block: {}'.format(df.p_block.unique()))
    print('p_err: {}'.format(df.p_err.unique()))
    print('dose: {}'.format(df.dose.unique()))
    print('tinf: {}'.format(df.tinf.unique()))

    gf = df.groupby(['dose', 'LFdecay', 'nGCs', 'nLFs', 'naive_pool', 'nkey',
                    'p_block', 'p_err', 'tinf']).count().mem_pool
    print(gf)

    return df

def aff_time_plot_mean(df, paramvar='nkey', dose=[1], LFdecay=[10*12],
                       nGCs=[1], nLFs=[25], naive_pool=[1000],
                       nkey=[1, 2], p_err=[0.003], tinf=[1],
                       p_block=[0.5]):
    """ Plots the simulation results matching the requirements. Import the
    data using data_import and pass the resultung dataframe.

    If there are several results matching identical criteria, their mean is
    calculated and plotted together with one std shaded area.

    Args:
        dataframe (pd df): Contains simulation results as produced by
                            data_import().
        Other parameters: Specify which simulation results are to be plotted.
    """
    # get subdataframe of only results where all requirements are fulfilled.
    sdf = df[(np.round(df.dose, 5).isin(dose)) & (df.LFdecay.isin(LFdecay)) &
             (df.nGCs.isin(nGCs)) & (df.nLFs.isin(nLFs))
             & (df.naive_pool.isin(naive_pool)) & (df.nkey.isin(nkey)) &
             (df.p_err.isin(p_err)) & (df.tinf.isin(tinf)) &
             (df.p_block.isin(p_block))]

    print(len(sdf))
    # group this remainder by everything
    grouped = sdf.groupby(['dose', 'LFdecay', 'nGCs', 'nLFs', 'naive_pool',
                           'nkey', 'p_block', 'p_err', 'tinf'])

    fig, (ax0, ax1, ax2) = plt.subplots(3, 1, sharex=True,
                                        figsize=(7, 13))
    lines = []
    labels = []
    # fake line object for explanation label
    l, = ax0.plot([1, 2, 3], [0.6, 0.7, 0.8], alpha=0)
    lines.append(l)
    labels.append('$dose, LFdecay, n_{GCs}, n_{LFs}, naive pool, n_{key}, p_{block}, p_{err}, n_{boost}$')
    colors = iter(cm.rainbow(np.linspace(0, 1, cc)))
    for key, group in grouped:
        col = next(colors)
        ax0.set_title('Mean Affinity in Memory Pool')
        ax0.set_ylim([0.6, 1])
        ax0.set_yticks([0.6, 0.7, 0.8, 0.9, 1.0])
        ax0.set_yticklabels(['$10^{-5}$', '$10^{-6}$', '$10^{-7}$',
                             '$10^{-8}$', '$10^{-9}$'])
        ax0.set_ylabel('$K_D$ of mean $E_\mathrm{Ab|Ag}$')
        labels.append(key)
        ebindlist = []
        for index, row in group.iterrows():
            ebindlist.append(row['E_bind'])
            timelist = row['evaldays']/12.
        Ebind = np.mean(np.array(ebindlist), axis=0)
        Estd = np.nanstd(np.array(ebindlist), axis=0)
        l, = ax0.plot(timelist, Ebind, color=col)
        ax0.fill_between(timelist, Ebind-Estd, Ebind+Estd, color=col,
                         alpha=0.3)
        lines.append(l)

        ax1.set_title('Mean SHMs in Memory Pool')
        ax1.set_ylabel('mean SHM count')
        shmlist = []
        for index, row in group.iterrows():
            shmlist.append(row['SHM'])
            timelist = row['evaldays']/12.
        SHMmean = np.mean(np.array(shmlist), axis=0)
        SHMstd = np.std(np.array(shmlist), axis=0)
        l, = ax1.plot(timelist, SHMmean, color=col)
        ax1.fill_between(timelist, SHMmean-SHMstd, SHMmean+SHMstd, color=col,
                         alpha=0.3)

        ax2.set_title('Normalised Shannon Entropy of Memory Pool')
        ax2.set_xticks([0, 7, 25, 35, 50, 63, 75, 100, 126])
        ax2.set_xticklabels(['0', 'I', '25', 'II', '50', 'III', '75', '100',
                             'C'])
        ax2.set_xlim([0, 130])
        ax2.set_xlabel('time (days)')
        ax2.set_ylabel('normalised H (bits)')
        ax2.set_ylim([0, 1])
        shannon = []
        for index, row in group.iterrows():
            shannon.append(np.array(row['entropy'].tolist()) /
                           np.log2(np.array(row['memcount'].tolist())))
            timelist = row['evaldays']/12
        shanmean = np.mean(np.array(shannon), axis=0)
        shanstd = np.std(np.array(shannon), axis=0)
        l, ax2.plot(timelist, shanmean, color=col)
        ax2.fill_between(timelist, shanmean-shanstd, shanmean+shanstd,
                         color=col, alpha=0.3)

        plt.legend(lines, labels, loc=(-0.5, -1))
        plt.tight_layout()
        seaborn.despine()
        pylab.savefig('figures/aff_time_map_mean_{}.pdf'.format(paramvar),
                       bbox_inches='tight')
    plt.show()

def aff_param_plot(df, paramvar='nkey', day=126, dose=[1], LFdecay=[10*12],
                   nGCs=[100], nLFs=[25], naive_pool=[1000],
                   nkey=[1, 2], p_err=[0.003], tinf=[3],
                   p_block=[0.5]):
    """ For the given parameter to vary (paramvar), finds all simulation
    results matching the other given parameter values and calculates mean
    affinity, SHMs, number of memory cells and entropy at the given day and
    plots these over paramvar with one std error bars. """
    # get subdataframe of only results where all requirements are fulfilled.
    sdf = df[(np.round(df.dose, 5).isin(dose)) & (df.LFdecay.isin(LFdecay)) &
             (df.nGCs.isin(nGCs)) & (df.nLFs.isin(nLFs))
             & (df.naive_pool.isin(naive_pool)) & (df.nkey.isin(nkey)) &
             (df.p_err.isin(p_err)) & (df.tinf.isin(tinf)) &
             (df.p_block.isin(p_block))]

    # group by all parameters that are NOT to be varied here
    grouplist = ['dose', 'LFdecay', 'nGCs', 'nLFs', 'naive_pool',
                 'nkey', 'p_block', 'p_err', 'tinf']
    grouplist.remove(paramvar)
    grouped = sdf.groupby(grouplist)
    # open plots
    fig, ((ax0, ax1), (ax2, ax3)) = plt.subplots(2, 2, sharex=True,
                                                 figsize=(8, 5))
    ax0.set_title('Mean Affinity in Memory Pool at day {}'.format(day))
    ax0.set_ylim([0.6, 1])
    ax0.set_yticks([0.6, 0.7, 0.8, 0.9, 1.0])
    ax0.set_yticklabels(['$10^{-5}$', '$10^{-6}$', '$10^{-7}$', '$10^{-8}$',
                         '$10^{-9}$'])
    ax0.set_ylabel('$K_D$ of mean $E_\mathrm{Ab|Ag}$')

    ax1.set_title('Mean SHMs in Memory Pool at day {}'.format(day))
    ax1.set_ylabel('number')

    ax2.set_title('Size of Memory Pool at day {}'.format(day))
    ax2.set_ylabel('cell number')
    ax2.ticklabel_format(style='sci', axis='y', scilimits=(0, 0))
    ax2.set_xlabel('${}$'.format(paramvar))

    ax3.set_title('Normalised Entropy of Memory Pool at day {}'.format(day))
    ax3.set_xlabel('${}$'.format(paramvar))
    ax3.set_ylabel('normalised H (bits)')
    ax3.set_ylim([0, 1])
    ax3.set_xlabel(paramvar)
    # find unique values of paramvar in each group
    # find all values at each paramval, plot means
    paramvar_per_group = []
    for key, group in grouped:
        paramvar_per_group = group[paramvar].unique().tolist()
        # affinity
        affs = []
        shms = []
        mems = []
        nentrs = []

        # for each paramvar val, get info
        for pv in paramvar_per_group:
            # move through runs and get required day info
            laff = []
            lshm = []
            lmem = []
            lnentr = []
            for run in range(len(group[(group[paramvar] == pv)])):
                laff.append(group[(group[paramvar] == pv)].E_bind[run][day-1])
                lshm.append(group[(group[paramvar] == pv)].SHM[run][day-1])
                lmem.append(group[(group[paramvar] == pv)].memcount[run][day-1])
                lnentr.append(group[(group[paramvar] == pv)].entropy[run][day-1]/
                              np.log2(group[(group[paramvar] == pv)].memcount[run][day-1]))
            affs.append((np.nanmean(laff), np.nanstd(laff)))
            shms.append((np.nanmean(lshm), np.nanstd(lshm)))
            mems.append((np.nanmean(lmem), np.nanstd(lmem)))
            nentrs.append((np.nanmean(lnentr), np.nanstd(lnentr)))

        # plot stuff gathered here
        ax0.errorbar(paramvar_per_group, [i[0] for i in affs], yerr=[i[1] for i in affs], fmt='o', label=key)
        ax1.errorbar(paramvar_per_group, [i[0] for i in shms], yerr=[i[1] for i in shms], fmt='o', label=key)
        ax2.errorbar(paramvar_per_group, [i[0] for i in mems], yerr=[i[1] for i in mems], fmt='o', label=key)
        ax3.errorbar(paramvar_per_group, [i[0] for i in nentrs], yerr=[i[1] for i in nentrs], fmt='o', label=key)

    # # fake label
    # ll1 = ' '.join(grouplist)
    # ax0.errorbar([0], [0], yerr=0, alpha=0, label='test')
    plt.legend(loc=(-0.5, -1))
    plt.tight_layout()
    seaborn.despine()
    pylab.savefig('figures/aff_{}_map.pdf'.format(paramvar),
                  bbox_inches='tight')

def heatmap_2p(df, paramvar1='dose', paramvar2='p_err', day=126,
               dose=[0.02, 0.05, 0.1, 0.2, 0.5, 1], LFdecay=[10*12],
               nGCs=[50], nLFs=[25], naive_pool=[1000],
               nkey=[10], p_err=[0, 0.0003, 0.001, 0.003, 0.01, 0.03], tinf=[3],
               p_block=[0.55]):

    # get subdataframe of only results where all requirements are fulfilled.
    sdf = df[(np.round(df.dose,3).isin(dose)) & (df.LFdecay.isin(LFdecay)) &
             (df.nGCs.isin(nGCs)) & (df.nLFs.isin(nLFs))
             & (df.naive_pool.isin(naive_pool)) & (df.nkey.isin(nkey)) &
             (df.p_err.isin(p_err)) & (df.tinf.isin(tinf)) &
             (df.p_block.isin(p_block))]

    # get a parameter value pair string for each non mapped parameter
    paravals = zip([dose, np.array(LFdecay)/12, nGCs, nLFs, naive_pool, nkey, p_block,
                    p_err, tinf], [r'$dose=$ ', r'$LFdecay=$ ', r'$n_{GCs}=$ ', r'$n_{LFs}=$ ',
                    r'$naive\ pool=$ ', r'$n_{key}=$ ', r'$p_{block}=$ ', r'$p_{err}=$ ', r'$n_{inf}=$ '])
    textblock = r'Parameters \\'
    for pv in paravals:
        textblock += pv[1] + (r'${}$\\'.format(pv[0]))
    print(textblock)

    # get unique values of the two variable parameters and prepare plot labels
    param1 = sdf[paramvar1].unique().tolist()
    param1.sort()
    labels1 = map(str, param1)
    print(labels1)

    param2 = sdf[paramvar2].unique().tolist()
    param2.sort()
    labels2 = map(str, param2)
    print(labels2)

    # get empty arrays of size p1*p2 to collect means later
    AffM = np.zeros((len(param1), len(param2)))
    AffS = np.zeros((len(param1), len(param2)))
    SHMM = np.zeros((len(param1), len(param2)))
    MemM = np.zeros((len(param1), len(param2)))
    EntrM = np.zeros((len(param1), len(param2)))

    # go through all combinations and get energy means, put to arrays
    for i in range(len(param1)):
        for j in range(len(param2)):
            laff = []
            lshm = []
            lmem = []
            lnentr = []
            for run in range(len(sdf[(sdf[paramvar1] == param1[i]) &
                                     (sdf[paramvar2] == param2[j])])):
                laff.append(sdf[(sdf[paramvar1] == param1[i]) &
                        (sdf[paramvar2] == param2[j])].E_bind[run][day-1])
            # for run in range(len(sdf[(aae(sdf[paramvar1], param1[i])) &
            #                          (aae(sdf[paramvar2], param2[j]))])):
            #     laff.append(sdf[(np.allclose(sdf[paramvar1], param1[i])) &
            #                              (np.allclose(sdf[paramvar2], param2[j]))].E_bind[run][day-1])
                # lshm.append(sdf[(sdf[paramvar1] == param1[i]) &
                #                          (sdf[paramvar2] == param2[j])].SHM[run][day-1])
                # lmem.append(sdf[(sdf[paramvar1] == param1[i]) &
                #                          (sdf[paramvar2] == param2[j])].memcount[run][day-1])
                # lnentr.append(sdf[(sdf[paramvar1] == param1[i]) &
                #                          (sdf[paramvar2] == param2[j])].entropy[run][day-1]/
                #               np.log2(sdf[(sdf[paramvar1] == param1[i]) &
                #                                        (sdf[paramvar2] == param2[j])].memcount[run][day-1]))

            AffM[i][j] = np.nanmean(laff)
            AffS[i][j] = np.nanstd(laff)
            SHMM[i][j] = np.nanmean(lshm)
            MemM[i][j] = np.nanmean(lmem)
            EntrM[i][j] = np.nanmean(lnentr)

    print(AffM)
    print(AffS)
    # # open plots
    fig, ax0 = plt.subplots(1, 1, figsize=(7,7))
    ax0.set_title('Mean affinity of memory pool on day {}'.format(day))

    sth = ax0.pcolor(AffM, vmin=0.6, vmax=0.9, cmap=plt.cm.coolwarm)
    ax0.set_xticks(np.arange(AffM.shape[1])+0.5, minor=False)
    ax0.set_yticks(np.arange(AffM.shape[0])+0.5, minor=False)
    ax0.set_xticklabels(labels2, minor=False)
    labels = ax0.get_xticklabels()
    plt.setp(labels, rotation=90)
    ax0.set_yticklabels(labels1, minor=False)

    ax0.set_xlabel(paramvar2)
    ax0.set_ylabel(paramvar1)
    # colourbar for affinities
    cbar_ax = fig.add_axes([1.02, 0.1, 0.025, 0.8])
    normKd = mpl.colors.Normalize(vmin=0.6,vmax=0.9)
    cmap=plt.cm.coolwarm
    cb1 = mpl.colorbar.ColorbarBase(cbar_ax, cmap=cmap,
                                    norm=normKd,
                                    orientation='vertical',extend='min')
    cb1.set_label('$K_D$ of mean $E_\mathrm{Ab|Ag}$ \n in memory pool (mol/l)')
    cb1.set_ticks([0.6,0.7,0.8,0.9])
    cb1.set_ticklabels(['$10^{-6}$','$10^{-7}$','$10^{-8}$','$10^{-9}$'])

    seaborn.despine()
    plt.text(5,0.,textblock)
    pylab.savefig('paramaps/heatmap_{}_{}_day{}.pdf'.format(paramvar1,paramvar2,day), bbox_inches='tight')

def return_mean_std(df, dose=[0.5], LFdecay=[120],
                    nGCs=[50], nLFs=[25], naive_pool=[1000],
                    nkey=[10], p_err=[0.003], tinf=[3],
                    p_block=[0.55]):
    """ Averages the simulation results matching the requirements, returns
    mean and std of energy, SHMs and normalised H. Take care to feed only
    a single set of paramters for meaningful results.

    Args:
        dataframe (pd df): Contains simulation results as produced by
                            data_import().
        Other parameters: Specify which simulation results are to be averaged.
    """
    # get subdataframe of only results where all requirements are fulfilled.
    sdf = df[(np.round(df.dose, 5).isin(dose)) & (df.LFdecay.isin(LFdecay)) &
             (df.nGCs.isin(nGCs)) & (df.nLFs.isin(nLFs))
             & (df.naive_pool.isin(naive_pool)) & (df.nkey.isin(nkey)) &
             (df.p_err.isin(p_err)) & (df.tinf.isin(tinf)) &
             (df.p_block.isin(p_block))]

    print(len(sdf))
    # group this remainder by everything
    grouped = sdf.groupby(['dose', 'LFdecay', 'nGCs', 'nLFs', 'naive_pool',
                           'nkey', 'p_block', 'p_err', 'tinf'])

    for key, group in grouped:
        ebindlist = []
        for index, row in group.iterrows():
            ebindlist.append(row['E_bind'])
            timelist = row['evaldays']/12.
        Ebind = np.mean(np.array(ebindlist), axis=0)
        Estd = np.nanstd(np.array(ebindlist), axis=0)

        shmlist = []
        for index, row in group.iterrows():
            shmlist.append(row['SHM'])
            timelist = row['evaldays']/12.
        SHMmean = np.mean(np.array(shmlist), axis=0)
        SHMstd = np.std(np.array(shmlist), axis=0)

        shannon = []
        for index, row in group.iterrows():
            shannon.append(np.array(row['entropy'].tolist()) /
                           np.log2(np.array(row['memcount'].tolist())))
        shanmean = np.mean(np.array(shannon), axis=0)
        shanstd = np.std(np.array(shannon), axis=0)

        return Ebind[-1], Estd[-1], shanmean[-1], shanstd[-1]
        # return Ebind, Estd, shanmean, shanstd, timelist


text:
import numpy
from distutils.core import setup
from distutils.extension import Extension
from Cython.Build import cythonize

extensions = [
    Extension("connections", ["connections.pyx"],
              include_dirs=[numpy.get_include()],
              extra_compile_args=["-O3"],
              extra_link_args=["-O3"]),
    Extension("skeletonize", ["skeletonize.pyx"],
              include_dirs=[numpy.get_include()],
              extra_compile_args=["-O3"],
              extra_link_args=["-O3"])
]

setup(
    ext_modules=cythonize(extensions, annotate=True)
)

ext: py

max_forks_repo_name: Chillee/CV-Maze-Solver

reason: import numpy


text:
import torch
import torchvision
import torchvision.transforms as transforms

import math
import time

# train three models, Federated learning
# each iteration over a subset of parameters: 1) average 2) pass back average to slaves 3) SGD step
# initialize with pre-trained models (better to use common initialization)
# loop order: loop 0: parameters/layers   {
#               loop 1 : {  averaging (part of the model)
#                loop 2: { epochs/databatches  { train; } } } }
# repeat this Nloop times

torch.manual_seed(69)
default_batch=512 # no. of batches (50000/3)/default_batch
batches_for_epoch=33#(50000/3)/default_batch
Nloop=12 # how many loops over the whole network
Nepoch=1 # how many epochs?
Nadmm=3 # how many FA iterations

# regularization
lambda1=0.0001 # L1 sweet spot 0.00031
lambda2=0.0001 # L2 sweet spot ?

load_model=False
init_model=True
save_model=True
check_results=True

transform=transforms.Compose(
   [transforms.ToTensor(),
     transforms.Normalize((0.5,0.5,0.5),(0.5,0.5,0.5))])

# split 50000 training data into three
subset1=range(0,16666)
subset2=range(16666,33333)
subset3=range(33333,50000)

trainset1=torchvision.datasets.CIFAR10(root='./torchdata', train=True,
    download=True, transform=transform)
trainloader1 = torch.utils.data.DataLoader(trainset1, batch_size=default_batch, shuffle=False, sampler=torch.utils.data.SubsetRandomSampler(subset1),num_workers=1)
trainloader2 = torch.utils.data.DataLoader(trainset1, batch_size=default_batch, shuffle=False, sampler=torch.utils.data.SubsetRandomSampler(subset2),num_workers=1)
trainloader3 = torch.utils.data.DataLoader(trainset1, batch_size=default_batch, shuffle=False, sampler=torch.utils.data.SubsetRandomSampler(subset3),num_workers=1)

testset=torchvision.datasets.CIFAR10(root='./torchdata', train=False,
    download=True, transform=transform)

testloader=torch.utils.data.DataLoader(testset, batch_size=default_batch,
    shuffle=False, num_workers=0)

classes=('plane', 'car', 'bird', 'cat',
  'deer', 'dog', 'frog', 'horse', 'ship', 'truck')

import numpy as np

# define a cnn
from simple_models import *

net1=Net()
net2=Net()
net3=Net()

# update from saved models
if load_model:
  checkpoint=torch.load('./s1.model')
  net1.load_state_dict(checkpoint['model_state_dict'])
  net1.train()
  checkpoint=torch.load('./s2.model')
  net2.load_state_dict(checkpoint['model_state_dict'])
  net2.train()
  checkpoint=torch.load('./s3.model')
  net3.load_state_dict(checkpoint['model_state_dict'])
  net3.train()

########################################################################### helper functions
def init_weights(m):
  if type(m)==nn.Linear or type(m)==nn.Conv2d:
    torch.nn.init.xavier_uniform_(m.weight)
    m.bias.data.fill_(0.01)

def unfreeze_one_layer(net,layer_id):
  ' set all layers to not-trainable except the layer given by layer_id (0,1,..)'
  for ci,param in enumerate(net.parameters(),0):
    if (ci == 2*layer_id) or (ci==2*layer_id+1):
       param.requires_grad=True
    else:
       param.requires_grad=False

def unfreeze_all_layers(net):
  ' reset all layers to trainable'
  for ci,param in enumerate(net.parameters(),0):
    param.requires_grad=True

def get_trainable_values(net):
  ' return trainable parameter values as a vector (only the first parameter set)'
  trainable=filter(lambda p: p.requires_grad, net.parameters())
  paramlist=list(trainable)
  N=0
  for params in paramlist:
    N+=params.numel()
  X=torch.empty(N,dtype=torch.float)
  X.fill_(0.0)
  offset=0
  for params in paramlist:
    numel=params.numel()
    with torch.no_grad():
      X[offset:offset+numel].copy_(params.data.view_as(X[offset:offset+numel].data))
    offset+=numel

  return X

def put_trainable_values(net,X):
  ' replace trainable parameter values by the given vector (only the first parameter set)'
  trainable=filter(lambda p: p.requires_grad, net.parameters())
  paramlist=list(trainable)
  offset=0
  for params in paramlist:
    numel=params.numel()
    with torch.no_grad():
     params.data.copy_(X[offset:offset+numel].data.view_as(params.data))
    offset+=numel

def number_of_layers(net):
  ' get total number of layers (note: each layers has weight and bias , so count as 2) '
  for ci,param in enumerate(net.parameters(),0):
   pass
  return int((ci+1)/2) # because weight+bias belong to one layer

def distance_of_layers(net1,net2,net3):
  'find Eculidean distance of each layer (from the mean) and return this as vector (normalized by size of layer)'
  L=number_of_layers(net1)
  W=np.zeros(L)
  for ci in range(L):
   unfreeze_one_layer(net1,ci)
   unfreeze_one_layer(net2,ci)
   unfreeze_one_layer(net3,ci)
   W1=get_trainable_values(net1)
   W2=get_trainable_values(net2)
   W3=get_trainable_values(net3)
   N=W1.numel()
   Wm=(W1+W2+W3)/3
   W[ci]=(Wm-W1).norm()/N
   W[ci]+=(Wm-W2).norm()/N
   W[ci]+=(Wm-W3).norm()/N
  return W

def sthreshold(z,sval):
  """soft threshold a tensor
    if element(z) > sval, element(z)=sval
    if element(z) < -sval, element(z)=-sval
  """
  with torch.no_grad():
    T=nn.Softshrink(sval) # if z_i < -sval, z_i -> z_i +sval , ...
    z=T(z)
  return z

def verification_error_check(net1,net2,net3):
  correct1=0
  correct2=0
  correct3=0
  total=0

  for data in testloader:
    images,labels=data
    outputs=net1(Variable(images))
    _,predicted=torch.max(outputs.data,1)
    correct1 += (predicted==labels).sum()
    outputs=net2(Variable(images))
    _,predicted=torch.max(outputs.data,1)
    correct2 += (predicted==labels).sum()
    outputs=net3(Variable(images))
    _,predicted=torch.max(outputs.data,1)
    correct3 += (predicted==labels).sum()
    total += labels.size(0)

  print('Accuracy of the network on the %d test images:%%%f %%%f %%%f'%
     (total,100*correct1/total,100*correct2/total,100*correct3/total))

##############################################################################################

if init_model:
  # note: use same seed for random number generation
  torch.manual_seed(0)
  net1.apply(init_weights)
  torch.manual_seed(0)
  net2.apply(init_weights)
  torch.manual_seed(0)
  net3.apply(init_weights)

criterion1=nn.CrossEntropyLoss()
criterion2=nn.CrossEntropyLoss()
criterion3=nn.CrossEntropyLoss()

L=number_of_layers(net1)
# get layer ids in given order 0..L-1 for selective training
np.random.seed(0)# get same list
Li=net1.train_order_layer_ids()
# make sure number of layers match
if L != len(Li):
  print("Warning, expected number of layers and given layer ids do not agree")
else:
  print(Li)

from lbfgsnew import LBFGSNew # custom optimizer
import torch.optim as optim
############### loop 00 (over the full net)
for nloop in range(Nloop):
  ############ loop 0 (over layers of the network)
  for ci in Li:
   unfreeze_one_layer(net1,ci)
   unfreeze_one_layer(net2,ci)
   unfreeze_one_layer(net3,ci)
   trainable=filter(lambda p: p.requires_grad, net1.parameters())
   params_vec1=torch.cat([x.view(-1) for x in list(trainable)])

   # number of parameters trained
   N=params_vec1.numel()
   z=torch.empty(N,dtype=torch.float,requires_grad=False)
   z.fill_(0.0)

   #opt1=optim.Adam(filter(lambda p: p.requires_grad, net1.parameters()),lr=0.001)
   #opt2=optim.Adam(filter(lambda p: p.requires_grad, net2.parameters()),lr=0.001)
   #opt3=optim.Adam(filter(lambda p: p.requires_grad, net3.parameters()),lr=0.001)
   opt1 =LBFGSNew(filter(lambda p: p.requires_grad, net1.parameters()), history_size=10, max_iter=4, line_search_fn=True,batch_mode=True)
   opt2 =LBFGSNew(filter(lambda p: p.requires_grad, net2.parameters()), history_size=10, max_iter=4, line_search_fn=True,batch_mode=True)
   opt3 =LBFGSNew(filter(lambda p: p.requires_grad, net3.parameters()), history_size=10, max_iter=4, line_search_fn=True,batch_mode=True)

   ############# loop 1 (Federated avaraging for subset of model)
   for nadmm in range(Nadmm):
     ##### loop 2 (data)
     for epoch in range(Nepoch):
        running_loss1=0.0
        running_loss2=0.0
        running_loss3=0.0

        for i,(data1,data2,data3) in enumerate(zip(trainloader1,trainloader2,trainloader3),0):
           # get the inputs
           inputs1,labels1=data1
           inputs2,labels2=data2
           inputs3,labels3=data3
           # wrap them in variable
           inputs1,labels1=Variable(inputs1),Variable(labels1)
           inputs2,labels2=Variable(inputs2),Variable(labels2)
           inputs3,labels3=Variable(inputs3),Variable(labels3)

           trainable=filter(lambda p: p.requires_grad, net1.parameters())
           params_vec1=torch.cat([x.view(-1) for x in list(trainable)])
           trainable=filter(lambda p: p.requires_grad, net2.parameters())
           params_vec2=torch.cat([x.view(-1) for x in list(trainable)])
           trainable=filter(lambda p: p.requires_grad, net3.parameters())
           params_vec3=torch.cat([x.view(-1) for x in list(trainable)])

           # fc1 and fc3 layers have L1 and L2 regularization
           def closure1():
                 if torch.is_grad_enabled():
                    opt1.zero_grad()
                 outputs=net1(inputs1)
                 loss=criterion1(outputs,labels1)
                 if ci in net1.linear_layer_ids():
                    loss+=lambda1*torch.norm(params_vec1,1)+lambda2*(torch.norm(params_vec1,2)**2)
                 if loss.requires_grad:
                    loss.backward()
                 return loss
           def closure2():
                 if torch.is_grad_enabled():
                    opt2.zero_grad()
                 outputs=net2(inputs2)
                 loss=criterion2(outputs,labels2)
                 if ci in net2.linear_layer_ids():
                    loss+=lambda1*torch.norm(params_vec2,1)+lambda2*(torch.norm(params_vec2,2)**2)
                 if loss.requires_grad:
                    loss.backward()
                 return loss
           def closure3():
                 if torch.is_grad_enabled():
                    opt3.zero_grad()
                 outputs=net3(inputs3)
                 loss=criterion3(outputs,labels3)
                 if ci in net3.linear_layer_ids():
                    loss+=lambda1*torch.norm(params_vec3,1)+lambda2*(torch.norm(params_vec3,2)**2)
                 if loss.requires_grad:
                    loss.backward()
                 return loss

           # ADMM step 1
           opt1.step(closure1)
           opt2.step(closure2)
           opt3.step(closure3)

           # only for diagnostics
           outputs1=net1(inputs1)
           loss1=criterion1(outputs1,labels1).data.item()
           running_loss1 +=loss1
           outputs2=net2(inputs2)
           loss2=criterion2(outputs2,labels2).data.item()
           running_loss2 +=loss2
           outputs3=net3(inputs3)
           loss3=criterion3(outputs3,labels3).data.item()
           running_loss3 +=loss3

           print('layer=%d %d(%d) minibatch=%d epoch=%d losses %e,%e,%e'%(ci,nloop,N,i,epoch,loss1,loss2,loss3))
     # Federated averaging
     x1=get_trainable_values(net1)
     x2=get_trainable_values(net2)
     x3=get_trainable_values(net3)
     znew=(x1+x2+x3)/3
     dual_residual=torch.norm(z-znew).item()/N # per parameter
     print('dual (loop=%d,layer=%d,avg=%d)=%e'%(nloop,ci,nadmm,dual_residual))
     z=znew
     put_trainable_values(net1,z)
     put_trainable_values(net2,z)
     put_trainable_values(net3,z)

     if check_results:
       verification_error_check(net1,net2,net3)

print('Finished Training')

if save_model:
 torch.save({
     'model_state_dict':net1.state_dict(),
     'epoch':epoch,
     'optimizer_state_dict':opt1.state_dict(),
     'running_loss':running_loss1,
     },'./s1.model')
 torch.save({
     'model_state_dict':net2.state_dict(),
     'epoch':epoch,
     'optimizer_state_dict':opt2.state_dict(),
     'running_loss':running_loss2,
     },'./s2.model')
 torch.save({
     'model_state_dict':net3.state_dict(),
     'epoch':epoch,
     'optimizer_state_dict':opt3.state_dict(),
     'running_loss':running_loss3,
     },'./s3.model')

max_stars_repo_name: niklausliu/federated-pytorch-test

max_issues_repo_name: niklausliu/federated-pytorch-test

max_forks_repo_name: niklausliu/federated-pytorch-test


text:
# Copyright 2020-present NAVER Corp. Under BSD 3-clause license

import os
import os.path as path
from random import randint
import numpy as np
from PIL import Image, ImageDraw

from .features import image_keypoints_from_file

def image_keypoints_to_image(
        image: Image,
        keypoints: np.ndarray,
        radius: int = 2,
        filled: bool = True) -> Image:
    """
    Displays keypoints on top of the image.

    :param image: an image
    :param keypoints: the keypoints
    :param radius: radius of the drawn circles
    :param filled: True: draw discs, False: draw circles
    :return: a new Image
    """
    draw = ImageDraw.Draw(image)
    for coords in keypoints[:, 0:2]:
        coords_tl = (coords - radius).astype(int)
        coords_br = (coords + radius).astype(int)
        color = tuple([randint(0, 255) for _ in range(3)])
        draw.ellipse((coords_tl[0], coords_tl[1], coords_br[0], coords_br[1]), outline=color,
                     fill=color if filled else None)
    return image

def image_keypoints_to_image_file(
        output_filepath: str,
        image_filepath: str,
        keypoints_filepath: str,
        keypoint_dtype,
        keypoint_dsize,
        radius: int = 2) -> None:
    """
    Displays keypoints on top of the image and save it to image file.

    :param output_filepath: input path to output image of keypoints.
    :param image_filepath: input path to input image.
    :param keypoints_filepath: input path to keypoints file.
    :param keypoint_dtype: input data type of keypoints data (cf. binary).
    :param keypoint_dsize: input data size of keypoints data (cf. binary)
    :param radius: radius of the keypoint in image (in pixel).
    :return:
    """
    os.makedirs(path.dirname(output_filepath), exist_ok=True)
    keypoints = image_keypoints_from_file(keypoints_filepath, keypoint_dtype, keypoint_dsize)
    image = Image.open(image_filepath)
    image_of_keypoints = image_keypoints_to_image(image, keypoints, radius)
    image_of_keypoints.save(output_filepath)

max_stars_repo_name: v-mehta/kapture

max_issues_repo_name: v-mehta/kapture

max_forks_repo_name: v-mehta/kapture


text:
import numpy as np

def calib_err(confidence, correct, p='2', beta=100):
    # beta is target bin size
    idxs = np.argsort(confidence)
    confidence = confidence[idxs]
    correct = correct[idxs]
    bins = [[i * beta, (i + 1) * beta] for i in range(len(confidence) // beta)]
    bins[-1] = [bins[-1][0], len(confidence)]

    cerr = 0
    total_examples = len(confidence)
    for i in range(len(bins) - 1):
        bin_confidence = confidence[bins[i][0]:bins[i][1]]
        bin_correct = correct[bins[i][0]:bins[i][1]]
        num_examples_in_bin = len(bin_confidence)

        if num_examples_in_bin > 0:
            difference = np.abs(np.nanmean(bin_confidence) - np.nanmean(bin_correct))

            if p == '2':
                cerr += num_examples_in_bin / total_examples * np.square(difference)
            elif p == '1':
                cerr += num_examples_in_bin / total_examples * difference
            elif p == 'infty' or p == 'infinity' or p == 'max':
                cerr = np.maximum(cerr, difference)
            else:
                assert False, "p must be '1', '2', or 'infty'"

    if p == '2':
        cerr = np.sqrt(cerr)

    return cerr

def soft_f1(confidence, correct):
    wrong = 1 - correct

    # # the incorrectly classified samples are our interest
    # # so they make the positive class
    # tp_soft = np.sum((1 - confidence) * wrong)
    # fp_soft = np.sum((1 - confidence) * correct)
    # fn_soft = np.sum(confidence * wrong)

    # return 2 * tp_soft / (2 * tp_soft + fn_soft + fp_soft)
    return 2 * ((1 - confidence) * wrong).sum()/(1 - confidence + wrong).sum()

def tune_temp(logits, labels, binary_search=True, lower=0.2, upper=5.0, eps=0.0001):
    logits = np.array(logits)

    if binary_search:
        import torch
        import torch.nn.functional as F

        logits = torch.FloatTensor(logits)
        labels = torch.LongTensor(labels)
        t_guess = torch.FloatTensor([0.5*(lower + upper)]).requires_grad_()

        while upper - lower > eps:
            if torch.autograd.grad(F.cross_entropy(logits / t_guess, labels), t_guess)[0] > 0:
                upper = 0.5 * (lower + upper)
            else:
                lower = 0.5 * (lower + upper)
            t_guess = t_guess * 0 + 0.5 * (lower + upper)

        t = min([lower, 0.5 * (lower + upper), upper], key=lambda x: float(F.cross_entropy(logits / x, labels)))
    else:
        import cvxpy as cx

        set_size = np.array(logits).shape[0]

        t = cx.Variable()

        expr = sum((cx.Minimize(cx.log_sum_exp(logits[i, :] * t) - logits[i, labels[i]] * t)
                    for i in range(set_size)))
        p = cx.Problem(expr, [lower <= t, t <= upper])

        p.solve()   # p.solve(solver=cx.SCS)
        t = 1 / t.value

    return t

def get_measures(confidence, correct):
    rms = calib_err(confidence, correct, p='2')
    mad = calib_err(confidence, correct, p='1')
    sf1 = soft_f1(confidence, correct)

    return rms, mad, sf1

def print_measures(rms, mad, sf1, method_name='Baseline'):
    print('\t\t\t\t\t\t\t' + method_name)
    print('RMS Calib Error (%): \t\t{:.2f}'.format(100 * rms))
    print('MAD Calib Error (%): \t\t{:.2f}'.format(100 * mad))
    print('Soft F1 Score (%):   \t\t{:.2f}'.format(100 * sf1))

def print_measures_with_std(rmss, mads, sf1s, method_name='Baseline'):
    print('\t\t\t\t\t\t\t' + method_name)
    print('RMS Calib Error (%): \t\t{:.2f}\t+/- {:.2f}'.format(100 * np.mean(rmss), 100 * np.std(rmss)))
    print('MAD Calib Error (%): \t\t{:.2f}\t+/- {:.2f}'.format(100 * np.mean(mads), 100 * np.std(mads)))
    print('Soft F1 Score (%):   \t\t{:.2f}\t+/- {:.2f}'.format(100 * np.mean(sf1s), 100 * np.std(sf1s)))

def show_calibration_results(confidence, correct, method_name='Baseline'):

    print('\t\t\t\t' + method_name)
    print('RMS Calib Error (%): \t\t{:.2f}'.format(
        100 * calib_err(confidence, correct, p='2')))

    print('MAD Calib Error (%): \t\t{:.2f}'.format(
        100 * calib_err(confidence, correct, p='1')))

    # print('Max Calib Error (%): \t\t{:.2f}'.format(
    #     100 * calib_err(confidence, correct, p='infty')))

    print('Soft F1-Score (%): \t\t{:.2f}'.format(
        100 * soft_f1(confidence, correct)))

    # add error detection measures?

reason: import numpy,import cvxpy


text:
from astropy import units as u
from json_to_dict import constants
import numpy as np
pi = np.pi

######################################################################################

def get_K(Ms, Mp, P, i, e=0*u.one, G=constants["G"]):
    top = (2*pi*G)**(1/3) * Mp * np.sin(i)
    bottom = P**(1/3) * (Ms + Mp)**(2/3) * np.sqrt(1-e**2)
    K = top/bottom
    K = K.decompose()
    return K

def get_Mp(Ms, K, P, e=0*u.one, G=constants["G"], i=0*u.deg, units=True):
    if units:
        Ms = (Ms.to(u.kg)).value
        K = (K.to(u.m/u.s)).value
        P = (P.to(u.s)).value
        e = e.value
        G = (G.to(u.Unit("m**3/(kg*s*s)"))).value
        i = (i.to(u.rad)).value

    A = ( K * P**(1/3)/(2*pi*G)**(1/3))**1.5 * (np.sqrt(1-e**2))**1.5
    Mp = A*Ms/(np.sin(i)-A)

    if units:
        Mp = Mp*u.kg

    return Mp

def get_Mp_2(Ms, K_true, P, minMp=0, maxMp=3, iters=100, i=90*u.deg):

    Mps = np.linspace(minMp,maxMp,iters) * constants["MJupiter"]
    Kdiffs = []
    for Mp in Mps:
        K = get_K(Ms, Mp, P, i)
        Kdiffs.append((K-K_true).value)
    Kdiffs = np.abs(np.array(Kdiffs))
    mindex = np.argmin((Kdiffs))

    return Kdiffs[mindex], Mps[mindex]/ constants["MJupiter"]

def get_a(mu, P, units=True, unit_return=u.AU):
    inner = mu*(P/(2*pi))**2
    a = inner**(1/3)
    if units:
        a = a.decompose()
        a = a.to(unit_return)

    return a

Ms = constants["MSun"]
K1 = 32*u.m/u.s
P1 = 150 * constants["TEarth"]
K2 = 88*u.m/u.s
P2 = 1200 * constants["TEarth"]
i=90*u.deg

#######################################################################################################################
# Part D
print("=====================================Part D============================================")
calc1=True
if calc1:
    Kdiff1, Mp1 = get_Mp_2(Ms, K1, P1, minMp=0, maxMp=1, iters=10000, i=i)
    print("(Minimum) mass of planet 1 (in Jupiter masses), Mp1: %s \nWith K error: %s m/s \n" %(Mp1, Kdiff1))

    Kdiff2, Mp2 = get_Mp_2(Ms, K2, P2, minMp=4, maxMp=5, iters=10000, i=i)
    print("(Minimum) mass of planet 2 (in Jupiter masses), Mp2: %s \nWith K error: %s m/s \n" %(Mp2, Kdiff2))

########################################################################################################################
# Part G
print("=====================================Part G============================================")
iG = 30*u.deg
calc2=True
if calc2:
    Kdiff1G, Mp1G = get_Mp_2(Ms, K1, P1, minMp=1, maxMp=2, iters=10000, i=iG)
    print("Mass of planet 1 (in Jupiter masses), for i=30 deg  Mp1G: %s \nWith K error: %s m/s \n" %(Mp1G, Kdiff1G))

    Kdiff2G, Mp2G = get_Mp_2(Ms, K2, P2, minMp=0, maxMp=100, iters=10000, i=iG)
    print("Mass of planet 2 (in Jupiter masses), for i=30 deg  Mp2G: %s \nWith K error: %s m/s \n" %(Mp2G, Kdiff2G))

########################################################################################################################
# Part H
print("=====================================Part H============================================")
mu = constants["G"]*constants["MSun"]
a1 = get_a(mu, P1)
print("Semi-major axis of planet 1, a1 (ie D1): %s \n" %a1)

a2 = get_a(mu, P2)
print("Semi-major axis of planet 2, a2 (ie D2): %s \n" %a2)

########################################################################################################################
# Part I
print("=====================================Part I============================================")

Mi = 0.8*Ms
mu_i = Mi*constants["G"]

calc3=True
if calc3:
    Kdiff1_I, Mp1_I = get_Mp_2(Mi, K1, P1, minMp=0, maxMp=1, iters=10000, i=i)
    print("(Minimum) mass of planet 1 (in Jupiter masses), for 0.8 solar mass star, Mp1_I: %s \nWith K error: %s m/s \n" %(Mp1_I, Kdiff1_I))

    Kdiff2_I, Mp2_I = get_Mp_2(Mi, K2, P2, minMp=3, maxMp=5, iters=10000, i=i)
    print("(Minimum) mass of planet 2 (in Jupiter masses), for 0.8 solar mass star, Mp2_I: %s \nWith K error: %s m/s \n" %(Mp2_I, Kdiff2_I))

a1_I = get_a(mu_i, P1)
print("Semi-major axis of planet 1, for a 0.8 solar mass star, a1_I (ie D1_I): %s \n" %a1_I)

a2_I = get_a(mu_i, P2)
print("Semi-major axis of planet 2, for a 0.8 solar mass star, a2_I (ie D2_I): %s \n" %a2_I)

max_stars_repo_name: MattTurnock/PlanetarySciencesMatt

max_issues_repo_name: MattTurnock/PlanetarySciencesMatt

max_forks_repo_name: MattTurnock/PlanetarySciencesMatt

reason: import numpy,from astropy


text:
from matplotlib import rc
rc('font',**{'family':'serif','serif':['Palatino']})
rc('text', usetex=True)
#import colormaps
from matplotlib import cm
from math import sqrt
import matplotlib.pyplot as plt
import matplotlib.path as mplPath
import numpy as np
from calculateforce import loadforces
Fel_, Fdrag_ = loadforces()
X_How_2d = np.array([[1.,4.5],[2.5,4.5],[2.5,1.1],[10.,1.1],[10.,-1.1],[2.5,-1.1],[2.5,-4.5],[1.,-4.5]])
def Fel(x,y,z):
    if x==0. and y==0.:
        return [0.,0.,Fel_(np.array([0,z]))[1]]
    else:
        rad=sqrt(x**2+y**2)
        x0=x/rad
        y0=y/rad
        return [x0*Fel_(rad,z)[0],y0*Fel_(rad,z)[0],Fel_(rad,z)[1]]
def Fdrag(x,y,z):
    if x==0. and y==0.:
        return [0.,0.,Fdrag_(np.array([0,z]))[1]]
    else:
        rad=sqrt(x**2+y**2)
        x0=x/rad
        y0=y/rad
        return [x0*Fdrag_(rad,z)[0],y0*Fdrag_(rad,z)[0],Fdrag_(rad,z)[1]]

#plt.ion()
#fig1=plt.figure(figsize=(18,12))
#fig=fig1.add_subplot()
#bar=fig1.add_subplot()
bar, fig = plt.subplots(figsize=(12,8))
ax=plt.axes()
def argument(x,y,z):
    return np.array([float(x),float(y),float(z)])
#def F(vec):
#    return [0.,0.,-2e-13]
def radius(x,y):
    return sqrt(x**2+y**2)
def sgn(x):
    if x<0:
        return -1
    elif x>0:
        return 1
    else:
        return 0

leftend=15.
#x_mem=np.linspace(X_How_2d[18][0],leftend,100)
#y_mem=np.zeros(x_mem.shape[0])+X_How_2d[18][1]
#x_mem_2=-x_mem
size=X_How_2d.shape[0]
X=np.zeros(size+1)
Y=np.zeros(size+1)
for index in range(size):
	X[index]=X_How_2d[index][0]
	Y[index]=X_How_2d[index][1]
X[size]=X[0]
Y[size]=Y[0]
X_2=-X

# whole domain: fac=0.1,p2=[

axes=plt.gca()
axes.set_ylim([-5,10])
axes.set_xlim([-10,10])

bbPath=mplPath.Path(X_How_2d)
Ny = 40
Nx = 50
plt.plot(X,Y,linewidth=3,color='black')
plt.plot(X_2,Y,linewidth=3,color='black')
#plt.plot(x_mem,y_mem,color='black',linewidth=1)
#plt.plot(x_mem_2,y_mem,color='black',linewidth=1)

Y, X = np.mgrid[-5:10:Ny*1j, -10:10:Nx*1j] #-5:30 and -30:30
U = np.zeros((Ny,Nx))
V = np.zeros((Ny,Nx))
for y in range(Ny):
    for x in range(Nx):
        if bbPath.contains_point((X[y][x],Y[y][x])) or bbPath.contains_point((-X[y][x],Y[y][x])):
            U[y][x] = 0
            V[y][x] = 0
        else:
            F=Fdrag(X[y][x],0.,Y[y][x])
            if sqrt(F[0]**2+F[1]**2+F[2]**2)<1e-16:
                U[y][x] = 0.
                V[y][x] = 0.
            else:
                U[y][x] = F[0]
                V[y][x] = F[2]
#for y in range(Ny):
#    for x in range(Nx):
#        if bbPath.contains_point((X[y][x],Y[y][x])) or bbPath.contains_point((-X[y][x],Y[y][x])):
#            U[y][x] = 0
#            V[y][x] = 0
#        else:
#            if Y[y][x]<-5.5 and (X[y][x]<-1.2 or X[y][x]>1.2):
#                U[y][x] = 0
#                V[y][x] = 0
#            else:
#                F=Fdrag(argument(X[y][x],0,Y[y][x]))
#                U[y][x] = F[0]
#                V[y][x] = F[2]

strm = plt.streamplot(X,Y,U,V,arrowsize=3, linewidth=2., density=2.5, cmap=cm.viridis, color=np.log10(np.sqrt(U*U+V*V)))
bar.colorbar(strm.lines)
ax.set_xlabel('X coordinate [nm]', fontsize=20)
ax.set_ylabel('Z coordinate [nm]', fontsize=20)
plt.title(r'$F_{\mathrm{drag}}(x)$', fontsize=25)
#plt.show()
plt.savefig('drag.eps')

max_stars_repo_name: jhwnkim/nanopores

max_issues_repo_name: jhwnkim/nanopores

max_forks_repo_name: jhwnkim/nanopores


text:
# Copyright (c) 2020, codenamewei
#
# 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 json
import time
import numpy as np
from konduit.client import Client
from konduit.load import client_from_file

current_milli_time = lambda: int(round(time.time() * 1000))

client = client_from_file("config.yaml")

total_count = 5
iteration = range(total_count)

before_milli_time = current_milli_time()

arr_input = [[0.272702493273322,0.0201936700818061]] #label: 0
#arr_input = [[0.867855051798607,0.597829895646543]] #label: 1

print("Get REST API Prediction")
for i in iteration:

    prediction = client.predict({client.input_names[0]: np.array(arr_input)})

    index = int(np.argmax(prediction))
    probabilities = np.max(prediction)

    print("Output: {}".format(prediction))
    print("Class Index: {}".format(index))
    print("Probabilities {:0.1f}".format(probabilities))

after_milli_time = current_milli_time()

total_time = after_milli_time - before_milli_time
milli_average = total_time / float(total_count)

print("Total time taken: {} ms".format(total_time))
print("Average time per request: {0:.2f} ms/request".format(milli_average))
print("Average time per request: {0:.2f} s/request".format(milli_average * 0.001))

max_stars_repo_name: codenamewei/konduit-serving-use-cases

max_issues_repo_name: codenamewei/konduit-serving-use-cases

max_forks_repo_name: codenamewei/konduit-serving-use-cases


text:
import pandas as pd
import numpy as np

from model.simulator import RunSingleTactic
from model.strategy_arrays import (
    EqualResFreqBreakdownArray, EqualSelectionArray,
)

# TOC
# RunAlongContourDF

# ThisStratDetailedDF
# ThisStratSummaryDF

class RunAlongContourDFs:
    def __init__(self, rand_pars, grid_output,
                 n_cont_points, strat_name) -> None:

        print(f"Running contour method: {strat_name}")

        if strat_name == "RFB":
            strat_obj = EqualResFreqBreakdownArray(grid_output)
        elif strat_name == "EqSel":
            strat_obj = EqualSelectionArray(grid_output)
        else:
            raise Exception(f"Invalid strat_name: {strat_name}")

        max_grid_EL = np.amax(grid_output.FY)

        self.df = ThisStratDetailedDF(rand_pars,
                                      n_cont_points,
                                      strat_obj,
                                      max_grid_EL).df

        self.summary = ThisStratSummaryDF(self.df,
                                          strat_obj.name,
                                          strat_obj.level,
                                          max_grid_EL).df

# End of RunAlongContourDF

class ThisStratDetailedDF:
    def __init__(self, rand_pars, n_cont_points, strat_obj, max_grid_EL) -> None:

        self.rand_pars = rand_pars
        self.n_cont_points = n_cont_points

        self.strat_name = strat_obj.name
        self.level = strat_obj.level

        self.strat_obj = strat_obj
        self.max_grid_EL = max_grid_EL

        self.df = self._get_df()

    def _get_df(self):

        cntr = self.strat_obj.find_contours(self.rand_pars,
                                            self.n_cont_points)
        out = self._find_df(cntr)

        return out

    def _find_df(self, contour):

        df_this_run = pd.DataFrame()

        if not contour:
            return df_this_run

        for dose1, dose2, cont_quant in zip(
                contour['x'],
                contour['y'],
                contour['cont_vals'],
        ):

            EL = self._get_EL_this_dose_combo(dose1, dose2)

            data = {"cont_quant": cont_quant,
                    "dose1": dose1,
                    "dose2": dose2,
                    "DS": dose1 + dose2,
                    "EL": EL}

            df_this_run = df_this_run.append(data, ignore_index=True)

        df_this_run['worked'] = max(df_this_run['EL']) >= self.max_grid_EL
        df_this_run['max_grid_EL'] = self.max_grid_EL

        return df_this_run

    def _get_EL_this_dose_combo(self, dose1, dose2):

        rp = self.rand_pars

        this_dose_conf = rp.get_single_conf(dose1, dose2)

        sing_run = RunSingleTactic(rp.fung_parms).run(this_dose_conf)

        EL = sing_run.failure_year

        return EL

# End of ThisStratDetailedDF

class ThisStratSummaryDF:

    def __init__(self, df, strat_name, level, max_grid_EL) -> None:
        self.df = df
        self.strat_name = strat_name
        self.level = level
        self.max_grid_EL = max_grid_EL

        self.LOW_THRESH = 1/3
        self.HIGH_THRESH = 2/3

        self.df = self.get_summary_df()

    def get_summary_df(self):
        df = self.df
        strat_name = self.strat_name

        if not df.shape[0]:
            return pd.DataFrame()

        lowDS_val = self._get_low_dose_max_EL()
        medDS_val = self._get_med_dose_max_EL()
        highDS_val = self._get_high_dose_max_EL()

        min_opt_dist_from_cntr = self._get_min_opt_dist_from_contour()

        data = {
            f"c_{strat_name[0]}_minContEL": min(df['EL']),
            f"c_{strat_name[0]}_maxContEL": max(df['EL']),

            f"c_{strat_name[0]}_min_opt_dist_from_contour": min_opt_dist_from_cntr,

            f"c_{strat_name[0]}_minDS": min(df['DS']),
            f"c_{strat_name[0]}_maxDS": max(df['DS']),

            f"c_{strat_name[0]}_worked_geq": self.max_grid_EL <= max(df['EL']),
            f"c_{strat_name[0]}_worked_equal": self.max_grid_EL == max(df['EL']),

            f"c_{strat_name[0]}_lowDoseMaxEL": lowDS_val,
            f"c_{strat_name[0]}_medDoseMaxEL": medDS_val,
            f"c_{strat_name[0]}_highDoseMaxEL": highDS_val,
        }

        return pd.DataFrame([data])

    def _get_low_dose_max_EL(self):
        df = self.df

        DS_thres_low = min(df['DS']) + self.LOW_THRESH * \
            (max(df['DS']) - min(df['DS']))

        filt = df[df['DS'] < DS_thres_low]

        return self._get_max_if_df_non_empty(filt)

    def _get_med_dose_max_EL(self):
        df = self.df

        DS_thres_low = min(df['DS']) + self.LOW_THRESH * \
            (max(df['DS']) - min(df['DS']))
        DS_thres_high = min(df['DS']) + self.HIGH_THRESH * \
            (max(df['DS']) - min(df['DS']))

        filt = df[((df['DS'] >= DS_thres_low) & (df['DS'] <= DS_thres_high))]

        return self._get_max_if_df_non_empty(filt)

    def _get_high_dose_max_EL(self):
        df = self.df

        DS_thres_high = min(df['DS']) + self.HIGH_THRESH * \
            (max(df['DS']) - min(df['DS']))

        filt = df[df['DS'] > DS_thres_high]

        return self._get_max_if_df_non_empty(filt)

    def _get_max_if_df_non_empty(self, df):
        if df.shape[0]:
            return max(df['EL'])
        else:
            return "NA"

    def _get_min_opt_dist_from_contour(self):
        df = self.df

        opt_df = df[df['EL'] == self.max_grid_EL]

        vec = abs(opt_df["cont_quant"] - self.level)

        if len(vec):
            out = min(vec)
            return out
        else:
            return "NA"

# End of ThisStratSummaryDF


text:
import sys

import numpy as np
from scipy.spatial.distance import cdist

from numpy.typing import ArrayLike

# try:
#     from enum import StrEnum
# except ImportError:
#     from backports.strenum import StrEnum

# class Strategy(StrEnum):
#     AFFINE = "affine"
#     # SIMILARITY = "similarity"
#     # RIGID = "rigid"

ORDER = "F"

def reshape_points(control_points: np.ndarray) -> np.ndarray:
    """Reshape NxD array into 1xDxNx1.

    Where D is the dimensionality, and N the number of points.

    Parameters
    ----------
    locations

    Returns
    -------
    reshaped array
    """
    n_locations, n_dimensions = control_points.shape
    return control_points.ravel().reshape(1, n_dimensions, n_locations, 1, order=ORDER)

class Transformer:
    def __init__(self, control_points: ArrayLike, deformed_control_points: ArrayLike):
        """Class for transforming points using Moving Least Squares.

        Given control point arrays must both be of same shape NxD,
        where N is the number of points,
        and D the dimensionality.

        Parameters
        ----------
        control_points : ArrayLike
        deformed_control_points : ArrayLike

        Raises
        ------
        ValueError
            Invalid control point array(s)
        """
        from_arr = np.asarray(control_points)
        to_arr = np.asarray(deformed_control_points)

        if from_arr.shape != to_arr.shape:
            raise ValueError("Control points must have the same shape")
        if from_arr.ndim != 2:
            raise ValueError("Control points must be 2D array")

        self.n_landmarks, self.ndim = from_arr.shape

        self.control_points: np.ndarray = reshape_points(from_arr)
        self.deformed_control_points: np.ndarray = reshape_points(to_arr)

    def _transform_affine(
        self, locs: np.ndarray, orig_cp: np.ndarray, deformed_cp: np.ndarray
    ) -> np.ndarray:
        n_locs = locs.shape[-1]

        # Pairwise distances between original control points and locations to transform
        # reshaped to 1,1,N_cp,N_l
        # jittered to avoid 0s
        distances = (
            cdist(
                orig_cp.ravel(ORDER).reshape(self.n_landmarks, self.ndim),
                locs.ravel(ORDER).reshape(n_locs, self.ndim),
                "sqeuclidean",
            ).reshape(1, 1, self.n_landmarks, n_locs)
            + sys.float_info.epsilon
        )
        weights = 1 / distances

        weights_inverse_norm = 1 / np.sum(weights, axis=2)

        # weighted centroids
        orig_star = np.einsum(
            "ijl,ijkl,ijkl->ijl", weights_inverse_norm, weights, orig_cp
        ).reshape(1, self.ndim, 1, n_locs)
        deformed_star = np.einsum(
            "ijl,ijkl,ijkl->ijl", weights_inverse_norm, weights, deformed_cp
        ).reshape(1, self.ndim, 1, n_locs)

        # distance to weighted centroids
        orig_hat = orig_cp - orig_star
        deformed_hat = deformed_cp - deformed_star

        Y = np.einsum("ijkl,mikl,mjkl->ijl", weights, orig_hat, orig_hat).reshape(
            self.ndim, self.ndim, 1, n_locs
        )

        Y_inv = np.zeros((self.ndim, self.ndim, 1, n_locs))
        for i in range(n_locs):
            Y_inv[:, :, 0, i] = np.linalg.inv(Y[:, :, 0, i])

        Z = np.einsum("ijkl,mikl,mjkl->ijl", weights, orig_hat, deformed_hat).reshape(
            self.ndim, self.ndim, 1, n_locs
        )

        vprime = (
            np.einsum("iakl,abkl,bjkl->ijkl", locs - orig_star, Y_inv, Z)
            + deformed_star
        )
        vprime = vprime.ravel(ORDER).reshape(n_locs, self.ndim)

        return vprime

    def transform(
        self,
        locations: ArrayLike,
        reverse=False,
        # strategy=Strategy.AFFINE,
    ) -> np.ndarray:
        """Transform some locations using the given control points.

        Uses the affine form of the MLS algorithm.

        Parameters
        ----------
        locations : ArrayLike
            NxD array of N locations in D dimensions to transform.
        reverse : bool, optional
            Transform from deformed space to original space, by default False

        Returns
        -------
        Deformed points
        """
        locs = np.asarray(locations)
        if locs.ndim != 2 or locs.shape[-1] != self.ndim:
            raise ValueError(
                "Locations must be 2D array of same width as control points, "
                f"got {locs.shape}"
            )

        orig_cp = self.control_points
        deformed_cp = self.deformed_control_points

        if reverse:
            orig_cp, deformed_cp = deformed_cp, orig_cp

        # Reshape into arrays with consistent indices
        locs = locs.ravel().reshape(1, self.ndim, 1, len(locs), order=ORDER)

        # if strategy == Strategy.AFFINE:
        #     return self._transform_affine(locs, orig_cp, deformed_cp)
        # else:
        #     raise ValueError(f"Unimplemented/ unknown strategy {strategy}")
        return self._transform_affine(locs, orig_cp, deformed_cp)

max_stars_repo_name: bdpedigo/molesq

max_issues_repo_name: bdpedigo/molesq

max_forks_repo_name: bdpedigo/molesq

reason: import numpy,from numpy,from scipy


text:
import argparse
import numpy as np
import tensorflow as tf
import time
import pickle

import maddpg.common.tf_util as U
from maddpg.trainer.maddpg import MADDPGAgentTrainer
import tensorflow.contrib.layers as layers

with U.single_threaded_session():

    episode_rewards = [0.0]  # sum of rewards for all agents
    agent_rewards = [[0.0] for _ in range(env.n)]  # individual agent reward
    final_ep_rewards = []  # sum of rewards for training curve
    final_ep_ag_rewards = []  # agent rewards for training curve
    agent_info = [[[]]]  # placeholder for benchmarking info
    saver = tf.train.Saver()
    obs_n = env.reset()
    episode_step = 0
    train_step = 0
    t_start = time.time()

    U.initialize()
    U.load_state("/tmp/policy/")
    saver = tf.train.Saver()

max_stars_repo_name: hamidosooli/MADDPG

max_issues_repo_name: hamidosooli/MADDPG

max_forks_repo_name: hamidosooli/MADDPG


text:
#!/usr/bin/env python3

# LICENSE
#
# _This file is Copyright 2018 by the Image Processing and Analysis Group (BioImage Suite Team). Dept. of Radiology & Biomedical Imaging, Yale School of Medicine._
#
# BioImage Suite Web is licensed under the Apache License, Version 2.0 (the "License");
#
# - you may not use this software except in compliance with the License.
# - You may obtain a copy of the License at [http://www.apache.org/licenses/LICENSE-2.0](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.__
#
# ENDLICENSE

import bis_path
import sys
import math
import numpy as np
import argparse
import bis_basemodule
import bis_objects
import modules_desc;
import biswrapper as libbis;

class approximateField(bis_basemodule.baseModule):

    def __init__(self):
        super().__init__();
        self.name='approximateField';

    def createDescription(self):
        return modules_desc.descriptions['approximateField'];

    def directInvokeAlgorithm(self,vals):
        print('oooo invoking: approximateField with vals', vals);
        input = self.inputs['input'];

        print('This far',libbis);
        try:
            self.outputs['output'] = libbis.approximateDisplacementFieldWASM2(input, {
                "spacing": vals['spacing'],
                "steps":  vals['steps'],
                "stepsize": vals['stepsize'],
                "lambda": vals['lambda'],
                "iterations": vals['iterations'],
                "tolerance": vals['tolerance'],
                'optimization': vals['optimization'],
                "windowsize": vals['windowsize'],
                "levels": vals['levels'],
                "resolution": vals['resolution'],
                "inverse": self.parseBoolean(vals['inverse']),
            }, self.parseBoolean(vals['debug']));
        except:
            e = sys.exc_info()[0]
            print('---- Failed to invoke algorithm',e);
            return False

        return True

if __name__ == '__main__':
    import bis_commandline; sys.exit(bis_commandline.loadParse(approximateField(),sys.argv,False));

#!/usr/bin/env python3

# LICENSE
#
# _This file is Copyright 2018 by the Image Processing and Analysis Group (BioImage Suite Team). Dept. of Radiology & Biomedical Imaging, Yale School of Medicine._
#
# BioImage Suite Web is licensed under the Apache License, Version 2.0 (the "License");
#
# - you may not use this software except in compliance with the License.
# - You may obtain a copy of the License at [http://www.apache.org/licenses/LICENSE-2.0](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.__
#
# ENDLICENSE

import bis_path
import sys
import math
import numpy as np
import argparse
import bis_basemodule
import bis_baseutils
import bis_objects
import modules_desc
import biswrapper as libbis;

class approximateField(bis_basemodule.baseModule):

    def __init__(self):
        super().__init__();
        self.name='approximateField';

    def createDescription(self):
        return modules_desc.descriptions['approximateField'];

    def directInvokeAlgorithm(self,vals):
        print('oooo invoking: approximateField', (vals));
        inputImg = self.inputs['input'];

        try:
            self.outputs['output'] = libbis.approximateDisplacementFieldWASM2(inputImg, {
                "spacing": (vals['spacing']),
                "steps": (vals['steps']),
                "stepsize": (vals['stepsize']),
                "lambda": (vals['lambda']),
                "iterations": (vals['iterations']),
                "tolerance": (vals['tolerance']),
                'optimization': bis_baseutils.getOptimizationCode(vals['optimization']),
                "windowsize": (vals['windowsize']),
                "levels": (vals['levels']),
                "resolution": (vals['resolution']),
                "inverse": self.parseBoolean(vals['inverse']),
            }, self.parseBoolean(vals['debug']));
        except:
            print('---- Failed to invoke algorithm nonlinear component of module');
            return False

        return True

if __name__ == '__main__':
    import bis_commandline; sys.exit(bis_commandline.loadParse(approximateField(),sys.argv,False));

max_stars_repo_name: KolCrooks/bisweb

max_issues_repo_name: KolCrooks/bisweb

max_forks_repo_name: KolCrooks/bisweb


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you 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 mxnet as mx
from mxnet.gluon import nn
from mxnet import amp
import numpy as np
import pytest

@pytest.fixture
def np_shape_array():
    flags = mx.npx.is_np_shape(), mx.npx.is_np_array(), mx.npx.is_np_default_dtype()
    mx.npx.set_np()
    yield
    mx.npx.set_np(*flags)

@pytest.fixture(scope='module')
def amp_init():
    amp.init()

def test_npi_concatenate_multicast(np_shape_array, amp_init):
    class Foo(nn.HybridBlock):
        def __init__(self, **kwargs):
            super().__init__(**kwargs)
            self.dense0 = nn.Dense(16, in_units=8)

        def forward(self, x):
            y = self.dense0(x)
            return mx.np.concatenate([y, x], axis=-1)

    foo = Foo()
    foo.initialize(ctx=mx.gpu())

    data = mx.np.ones((32, 8), ctx=mx.gpu())
    out = foo(data)
    assert out.dtype == np.float32

max_stars_repo_name: PawelGlomski-Intel/incubator-mxnet

max_issues_repo_name: PawelGlomski-Intel/incubator-mxnet

max_forks_repo_name: PawelGlomski-Intel/incubator-mxnet


text:
import numpy as np
from smbus2 import SMBus, i2c_msg
# from smbus import SMBus
import wiringpi

class SerialLink:
    """
    Class to communicate with castle serial link, takes Pyserial Serial class
    in constructor
    """

    def __init__(self, serial, device_id=0):
        """
        :param serial: Pyserial Serial class
        :param device_id: number between 0-63 associate with serial link address

        """

        assert isinstance(device_id, int)
        assert 0 <= device_id < 64

        # assign class
        # default baudrate is 9600 and device id is 0 on the serial link
        self.serial_link = serial

        self.device_id = device_id

        # Addresses for registers, taken from serial link manual
        self.register_dictionary = {
            "voltage": 0, "ripple": 1, "current": 2, "throttle": 3, "power": 4,
            "speed": 5, "temp": 6, "BEC voltage": 7, "BEC current": 8,
            "Raw NTC": 9, "RAW linear": 10, "link live": 25, "fail safe": 26,
            "e stop": 27, "packet in": 28, "packet out": 29, "check bad": 30,
            "packet bad": 31, "write throttle": 128, "write fail safe": 129,
            "write e stop": 130, "write packet in": 131,
            "write packet out": 132, "write check bad": 133,
            "write packet bad": 134
        }

        # Scales for converting from bytes to integers,
        # taken from serial link manual. -1 means do not scale
        self.scale_dictionary = {
            "voltage": 20.0, "ripple": 4.0, "current": 50.0, "throttle": 1.0,
            "power": 0.2502, "speed": 20416.66, "temp": 30.0,
            "BEC voltage": 4.0, "BEC current": 4.0, "Raw NTC": 63.8125,
            "RAW linear": 30.0, "link live": -1, "fail safe": -1, "e stop": -1,
            "packet in": -1, "packet out": -1, "check bad": -1, "packet bad": -1
        }

    def calc_checksum(self, rsp):
        """
        Verify the incoming checksum from the serial link
        :param rsp: a byte array, has shape (3,). The first two values are read,
         and the -1 value is the checksum
        :return: True if checksum is valid, false if not
        """

        intlist = [int(a) for a in rsp]

        # Calculate the checksum
        if np.sum(intlist) % 256 == 0:
            return True
        else:
            return False

    def append_checksum(self, array):
        """
        array is list of size 4. Appends the checksum of the first 4 bytes to
        the end of the list
        :param array: a list of size 4, contains the command we are writing the
        the serial link
        :return: array with  the checksum
        """

        array.append((0 - np.sum(array)) % 256)
        return array

    def convert_response(self, var, rsp):
        """
        :param var: The variable we requested, e.g. "voltage"
        :param rsp: the 3-byte response of the serial link. The last byte is a
        checksum
        :return: The converted response
        """

        # extract the values from the byte array
        values = rsp[:2]

        # get the the scale to multiply the values by
        scale = self.scale_dictionary[var]

        # if scale = -1, there is no scale associated with it in the manual
        if scale == -1:
            return values[0] * 256 + values[1]
        else:
            return (values[0] * 256 + values[1]) * scale / 2042.0

    def write_var(self, var, s):
        """
        Sends a command to the serial link
        :param var: the register to write to, e.g. "speed"
        :param s: integer between 0 and 65535 to write to the register
        :return: the response of the serial link, which is just confirmation of
        the command we sent
        with a checksum appended. This does not contain new information
        """
        assert 0 <= s < 256**2
        assert isinstance(s, int)

        # the list of bytes (as integers) to write to the serial link
        # byte 1 is the device address, always starts with 1xxxxxxx hence a 128
        # byte 2 is the register address
        # bytes 3-4 are the values to be written to the register
        # byte 5 is a checksum
        write_array = [128 + self.device_id,
                     self.register_dictionary[var],
                     int(s // 256), s % 256]

        write_array = self.append_checksum(write_array)

        self.serial_link.write(write_array)

        # 3 byte response from serial link confirming command was sent
        rsp = self.serial_link.read(3)
        return rsp

    def read_var(self, var):
        """
        Reads a value from the serial link register
        :param var: the register to read from, e.g. "speed"
        :return: the 3-byte response of the serial link
        """

        # byte 1 is the device address, always starts with 1xxxxxxx hence a 128
        # byte 2 is the register address
        # bytes 3-4 are ignored by the serial link since we are reading
        # byte 5 is a checksum
        read_array = [128 + self.device_id,
                      self.register_dictionary[var],
                      0, 0]
        read_array = self.append_checksum(read_array)

        self.serial_link.write(read_array)
        rsp = self.serial_link.read(3)

        # check that the response is indeed 3 bytes and the checksum is correct
        if len(rsp) == 3 and self.calc_checksum(rsp):
            return self.convert_response(var, rsp)
        else:
            arg = "Invalid response from serial link. " \
                  "The response was: {}".format(rsp)
            raise Exception(arg)

class SerialLinkI2C(SerialLink):
    def simple_test_wiringpi(self):
        bus = SMBus(1)
        # 7 bit address (will be left shifted to add the read write bit)
        address = 0x0b
        reg_read_voltage = self.register_dictionary['voltage']
        reg_write_throttle = self.register_dictionary['write throttle']

        #          76543210
        # add = int('00001011', 2)  # address 11
        # add = int('00010110', 2)  # address 11 + write bit
        # add = int('00010111', 2)  # address 11 + read bit
        # add = int('00001100', 2)  # address 12

        # var = 'voltage'
        var = 'write throttle'
        reg = self.register_dictionary[var]
        wlen = 5
        rlen = 3

        # b = bus.read_byte_data(i2c_addr=80, register=0)
        # print(b)

        # byte 1 is the 7bit device address
        # byte 2 is the register address
        # bytes 3-4 are ignored by the serial link since we are reading
        # byte 5 is a checksum

        # 10000 = 00100111 00010000
        b1 = int('00100111', 2)
        b0 = int('00010000', 2)

        # Write an array of registers
        data = [b1, b0]
        data = self.append_checksum(data)
        data = list(data)
        print('data %s' % data)
        bus.write_i2c_block_data(address, reg_write_throttle, data)
        # block = bus.read_i2c_block_data(address, reg_write_throttle, rlen)
        block = bus.read_byte(address)
        print('block %s' % block)
        block = bus.read_byte(address)
        print('block %s' % block)
        block = bus.read_byte(address)
        print('block %s' % block)

        bus.close()

    def simple_test_smbus(self):
        bus = SMBus(1)
        # 7 bit address (will be left shifted to add the read write bit)
        address = 0x0b
        reg_read_voltage = self.register_dictionary['voltage']
        reg_write_throttle = self.register_dictionary['write throttle']

        #          76543210
        # add = int('00001011', 2)  # address 11
        # add = int('00010110', 2)  # address 11 + write bit
        # add = int('00010111', 2)  # address 11 + read bit
        # add = int('00001100', 2)  # address 12

        # var = 'voltage'
        var = 'write throttle'
        reg = self.register_dictionary[var]
        wlen = 5
        rlen = 3

        # b = bus.read_byte_data(i2c_addr=80, register=0)
        # print(b)

        # byte 1 is the 7bit device address
        # byte 2 is the register address
        # bytes 3-4 are ignored by the serial link since we are reading
        # byte 5 is a checksum

        # 10000 = 00100111 00010000
        b1 = int('00100111', 2)
        b0 = int('00010000', 2)

        # Write an array of registers
        data = [b1, b0]
        data = self.append_checksum(data)
        data = list(data)
        print('data %s' % data)
        bus.write_i2c_block_data(address, reg_write_throttle, data)
        # block = bus.read_i2c_block_data(address, reg_write_throttle, rlen)
        block = bus.read_byte(address)
        print('block %s' % block)
        block = bus.read_byte(address)
        print('block %s' % block)
        block = bus.read_byte(address)
        print('block %s' % block)

        bus.close()

    def simple_test_smbus2(self):
        # Open i2c bus 1 and read one byte from address 80, offset 0
        bus = SMBus(1)
        # bus.open(1)
        add = int('0b', 16)
        # var = 'voltage'
        var = 'write throttle'
        reg = self.register_dictionary[var]
        wlen = 5
        rlen = 3

        # b = bus.read_byte_data(i2c_addr=80, register=0)
        # print(b)

        # byte 1 is the 7bit device address
        # byte 2 is the register address
        # bytes 3-4 are ignored by the serial link since we are reading
        # byte 5 is a checksum
        b1 = int('00100111', 2)    # d10000 = 00100111 00010000 = d39 d16
        b0 = int('00010000', 2)
        read_array = [add,
                      self.register_dictionary[var],
                      b1, b0]
        read_array = self.append_checksum(read_array)

        # # Write a block of 8 bytes to address 80 from offset 0
        data = [b0, b1]
        bus.write_i2c_block_data(i2c_addr=add, register=reg, data=data)
        print('data %s' % data)

        # Read a block of 3 bytes from address 80, offset 0
        # block = bus.read_i2c_block_data(i2c_addr=add, register=reg, length=3)
        # Returned value is a list of 16 bytes
        # print('block %s' % block)

        msg = i2c_msg.read(address=add, length=rlen)
        bus.i2c_rdwr(msg)
        # block = msg.data.contents.block[1:rlen + 1]
        for value in msg:
            print(value)
        # msg.data.contents.block[1:length + 1]

        bus.close()


text:
import sys,os
import random
import numpy as np
import matplotlib.pyplot as plt
sys.path.append('../')
from mrina.mri_utils import extractFluidMask
from mrina.recover import recover_vel
import argparse

home = os.getenv('HOME')

# Parameters
EPSILON = 0.5

def circle(shape, center, dist, useFluidMask=False, fluidMask=None):
  # probably a more efficient way to determine this
  pts = []
  for x in range(max(0,center[0]-dist), min(center[0]+dist, shape[0])):
    for y in range(max(0,center[1]-dist), min(center[1]+dist, shape[1])):
      if abs((x-center[0])**2 + (y-center[1])**2 - dist**2) < EPSILON**2:
        if(useFluidMask):
          if(fluidMask[x,y]):
            pts.append([x, y])
        else:
          pts.append([x, y])
  return pts

def pick_random_point(imsz, useFluidMask=False, fluidMask=None):
  if(useFluidMask):
    found = False
    while(not found):
      pt1 = (random.randint(0, imsz[0]-1), random.randint(0, imsz[1]-1))
      found = fluidMask[pt1[0],pt1[1]]
  else:
    pt1 = (random.randint(0, imsz[0]-1), random.randint(0, imsz[1]-1))
  return pt1

def select_points(dist, imsz, useFluidMask=False, fluidMask=None):
  # Select center at random
  pt1 = pick_random_point(imsz, useFluidMask, fluidMask)
  choices = circle(imsz, pt1, dist, useFluidMask, fluidMask)
  # Not sure why this loop is here...
  while not choices: # make sure there were points within range for distance
    pt1 = pick_random_point(imsz, useFluidMask, fluidMask)
    choices = circle(imsz, pt1, dist, useFluidMask, fluidMask)
  # Select pt2 at random from the list
  pt2 = tuple(random.choice(choices))
  if pt1 <= pt2: # pt1[0] < pt2[0] or (pt1[0] == pt2[0] and pt1[1] <= pt2[1]):
    return pt1,pt2
  return pt2, pt1

def get_points(size, num_pts, imsz, ptsdir, pointFileName, useFluidMask=False, fluidMask=None):
  print('Getting points of size ' + str(size) + ' and ' + str(num_pts) + ' points...')
  points = np.zeros((size, num_pts, 2, 2), dtype=int)
  for k in range(1,size+1):
    if(k % 10 == 0):
      print(k)
    for j in range(num_pts):
      pt1, pt2 = select_points(k,imsz,useFluidMask,fluidMask)
      while pt1[0] in points[k-1,:,0,0] and pt1[1] in points[k-1,:,0,1] and pt2[0] in points[k-1,:,1,0] and pt2[0] in points[k-1,:,1,1]:
        # Select new points if we've already selected these
        pt1, pt2 = select_points(k,imsz,useFluidMask,fluidMask)
      points[k-1, j, 0, 0] = pt1[0]
      points[k-1, j, 0, 1] = pt1[1]
      points[k-1, j, 1, 0] = pt2[0]
      points[k-1, j, 1, 1] = pt2[1]
  # Save Point File
  np.save(pointFileName,points)

def get_samples(recfile, singlechannel, vencdir):
  # Load reconstructed images
  samples = np.load(recfile)
  # Get Vencoding
  vencfile = vencdir + 'venc_n1' + '.npy'
  if os.path.exists(vencfile):
    venc = np.load(vencfile)
  else:
    venc = None
  # Recover the velocities
  if(singlechannel):
    samples = np.absolute(samples)
  else:
    samples = recover_vel(samples, venc)
  # Compute average reconstructed image
  avgimg = np.mean(samples, axis=0)
  # Compute reconstruction noise and return
  for k in range(samples.shape[0]):
    samples[k] = samples[k] - avgimg
  return samples

def get_saved_points(samples, size, num_pts, maindir, ptsdir, useFluidMask, fluidMaskFile):
  imsz = samples.shape[3:]

  # Get File Name
  fluidMask = None
  if(useFluidMask):
    if(fluidMaskFile == ''):
      print('Using fluid mask...')
    else:
      print('Using fluid mask from file: '+maindir+fluidMaskFile+'...')

    # Get File with Original Image
    orig_file = maindir + 'imgs_n1.npy'

    # Check if files exist and open the images
    if(os.path.exists(orig_file)):
      orig = np.load(orig_file).astype(np.complex)
    else:
      print('ERROR: Cannot find file with original image.')
      exit(-1)

    if(fluidMaskFile == ''):
      # Extract fluid mask for simple geometries
      fluidMask = extractFluidMask(orig)
    else:
      # Read mask from file
      fluidMask = np.load(maindir + fluidMaskFile)

    if(False):
      # Plot the fluid mask for debug
      plt.imshow(fluidMask,cmap='gray')
      plt.show()
      exit()

    # Set File Name for Point in the Fluid Mask
    pointFileName = ptsdir+'points_s' + str(size) + '_n' + str(num_pts) + '_fluidmask.npy'
  else:

    # Set the file name for points without fluid mask
    pointFileName = ptsdir+'points_s' + str(size) + '_n' + str(num_pts) + '.npy'

  # Read or generate the points
  if not os.path.isfile(pointFileName):
    get_points(size, num_pts, imsz, ptsdir, pointFileName, useFluidMask, fluidMask)
  points = np.load(pointFileName)
  return points

def get_coeff(size, num_pts, samples, points):
  print('in coeff', samples.shape)
  coeff = np.zeros((samples.shape[1], size, num_pts))
  for k in range(1, size+1):
    for j in range(num_pts):
      pt1 = (points[k-1, j, 0,0], points[k-1,j, 0,1])
      pt2 = (points[k-1, j, 1,0], points[k-1,j, 1,1])
      for v in range(samples.shape[1]):
        var1 = np.abs(samples[:, v, 0, pt1[0], pt1[1]])
        var2 = np.abs(samples[:, v, 0, pt2[0], pt2[1]])
        # If both sequences contain zeros, this generates a nan.
        # We replace nan with ones. This is only a problem when evaluating the
        # correlation between pixels in the background.
        coeff[v, k-1, j] = np.corrcoef(np.asarray([var1, var2]))[0,1] # correlation
  if(np.isnan(coeff).any()):
    print('WARNING: Found some NAN, probably due to sequences of zeros in the background. Replacing with correlation of 1.0')
    coeff[np.isnan(coeff)] = 1.0
  return coeff

def get_vals(recfile, savefile, maxcorrpixeldist, numpts, singlechannel, maindir, recdir, ptsdir, useFluidMask, fluidMaskFile, save_numpy=True):
  samples = get_samples(recfile, singlechannel, args.vencdir)
  print('Use Fluid Mask: ',useFluidMask)
  points  = get_saved_points(samples, maxcorrpixeldist, numpts, maindir, ptsdir, useFluidMask, fluidMaskFile)
  coeff   = get_coeff(maxcorrpixeldist, numpts, samples, points)
  if save_numpy:
    np.save(savefile, coeff)
  return coeff

def get_all(args):
  # size is max. distance to retrieve correlations for
  # num_pts is the number of points to average correlation across
  for noise_percent in [0.0, 0.01, 0.05, 0.1, 0.3]:
    for p in [0.25, 0.5, 0.75, 0.80, 0.85, 0.90, 0.95]:
      for wavetype in ['HAAR','DB8']:
        for algtype in ['CS','CSDEB','OMP']:
          for samptype in ['bernoulli', 'vardengauss']:
            recfile = args.recdir + 'rec_noise' + str(int(noise_percent*100)) + \
                                    '_p' + str(int(p*100)) + samptype + \
                                    '_n' + str(args.numsamples) + \
                                    '_w' + str(wavetype) + \
                                    '_a' + str(algtype) + '.npy'

            savefile = args.recdir + 'corrcoeff' + str(args.numpts) + \
                                     '_noise' + str(int(noise_percent*100)) + \
                                     '_p' + str(int(p*100)) + samptype + \
                                     '_n'+str(args.numsamples)+ \
                                     '_w' + str(wavetype) + \
                                     '_a' + str(algtype)

            # Change name for fluid Mask
            if(args.usefluidmask):
              savefile += '_fluid.npy'
            else:
              savefile += '.npy'

            if(os.path.exists(recfile)):
              get_vals(recfile, savefile, args.maxcorrpixeldist, args.numpts, args.singlechannel, args.maindir, args.recdir, args.ptsdir, args.usefluidmask, args.fluidmaskfile)
              print('Saved! Noise: ', noise_percent,', undersampling prob: ', p, ', undersampling mask: ', samptype, ', wavetype: ', wavetype, ', algorithm: ', algtype)

# MAIN
if __name__ == '__main__':

  # Init parser
  parser = argparse.ArgumentParser(description='Generate result images.')

  # numsamples
  parser.add_argument('-n', '--numsamples',
                      action=None,
                      # nargs='+',
                      const=None,
                      default=100,
                      type=int,
                      choices=None,
                      required=False,
                      help='number of repetitions',
                      metavar='',
                      dest='numsamples')
  # numpts
  parser.add_argument('-m', '--numpts',
                      action=None,
                      # nargs='+',
                      const=None,
                      default=50,
                      type=int,
                      choices=None,
                      required=False,
                      help='number of point pairs for computing correlations',
                      metavar='',
                      dest='numpts')
  # size
  parser.add_argument('-c', '--maxcorrpixeldist',
                      action=None,
                      # nargs='+',
                      const=None,
                      default=100,
                      type=int,
                      choices=None,
                      required=False,
                      help='Maximum integer pixel distance to use for correlation computation',
                      metavar='',
                      dest='maxcorrpixeldist')
  # recdir
  parser.add_argument('-d', '--maindir',
                      action=None,
                      # nargs='+',
                      const=None,
                      default='./',
                      type=str,
                      choices=None,
                      required=False,
                      help='folder containing the original image',
                      metavar='',
                      dest='maindir')
  # recdir
  parser.add_argument('-r', '--recdir',
                      action=None,
                      # nargs='+',
                      const=None,
                      default='./',
                      type=str,
                      choices=None,
                      required=False,
                      help='folder containing the non linear reconstructions',
                      metavar='',
                      dest='recdir')
  # ptsdir
  parser.add_argument('-t', '--ptsdir',
                      action=None,
                      # nargs='+',
                      const=None,
                      default='./',
                      type=str,
                      choices=None,
                      required=False,
                      help='folder containing the saved point coordinates',
                      metavar='',
                      dest='ptsdir')
  # vencdir
  parser.add_argument('-v', '--vencdir',
                      action=None,
                      # nargs='+',
                      const=None,
                      default='./',
                      type=str,
                      choices=None,
                      required=False,
                      help='folder containing the velocity encoding',
                      metavar='',
                      dest='vencdir')
  # singlechannel
  parser.add_argument('--singlechannel',
                      action='store_true',
                      default=False,
                      required=False,
                      help='treat the image as single-channel, without velocity components',
                      dest='singlechannel')
  # usefluidmask
  parser.add_argument('--usefluidmask',
                      action='store_true',
                      default=False,
                      required=False,
                      help='compute correlation only within fluid region',
                      dest='usefluidmask')

  # fluidmaskfile
  parser.add_argument('--fluidmaskfile',
                      action=None,
                      # nargs='+',
                      const=None,
                      default='',
                      type=str,
                      choices=None,
                      required=False,
                      help='name of the npy file containing the binary mask',
                      metavar='',
                      dest='fluidmaskfile')

  # Print Level
  parser.add_argument('-p', '--printlevel',
                      action=None,
                      # nargs='+',
                      const=None,
                      default=0,
                      type=float,
                      choices=None,
                      required=False,
                      help='print level, 0 - no print, >0 increasingly more information ',
                      metavar='',
                      dest='printlevel')

  # Parse Commandline Arguments
  args = parser.parse_args()

  # Save all the images
  get_all(args)

  # Completed!
  if(args.printlevel > 0):
    print('Completed!!!')

max_stars_repo_name: desResLab/mrina

max_issues_repo_name: desResLab/mrina

max_forks_repo_name: desResLab/mrina