Last active
November 19, 2017 17:47
-
-
Save astatide/234cbc6f23862346aa24666a923ff927 to your computer and use it in GitHub Desktop.
PlotFunctions using Pandas
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| # unbound to bound is 0,1 | |
| # bound to unbound is 1,0 | |
| # | |
| # # Due to the data structure in SIMS.p, each sim set is correct for both steady state simulations (that is, they include data from BOTH directions) | |
| # so we only bother correcting one of these sets, since in the plot I'm only ever pulling in from one simulation key (hell of a lot easier to plot that way) | |
| CH4.CH4: | |
| 05.EQ.COLOR.BOUND.LIGHT.2.WALKERS: | |
| (0,1): 8 | |
| (1,0): 4 | |
| 06.EQ.COLOR.BOUND.LIGHT.4.WALKERS: | |
| (0,1): 6 | |
| (1,0): 5 | |
| 07.ST.NOCOLOR.BOUND.LIGHT.2.WALKERS: | |
| (0,1): 18 # From this simulation; disregard and update. | |
| (0,1): 39 | |
| (1,0): 5 | |
| 08.ST.NOCOLOR.BOUND.LIGHT.4.WALKERS: | |
| (0,1): 13 # Again, from this simulation. It's not what we want. | |
| (0,1): 15 | |
| (1,0): 5 | |
| 09.ST.NOCOLOR.UNBOUND.LIGHT.2.WALKERS: | |
| (0,1): 39 | |
| (1,0): 'NONE' | |
| 10.ST.NOCOLOR.UNBOUND.LIGHT.4.WALKERS: | |
| (0,1): 15 | |
| (1,0): 38 | |
| NA.CL: | |
| 05.EQ.COLOR.BOUND.LIGHT.2.WALKERS: | |
| (0,1): 5 | |
| (1,0): 5 | |
| 06.EQ.COLOR.BOUND.LIGHT.4.WALKERS: | |
| (0,1): 5 | |
| (1,0): 4 | |
| 07.ST.NOCOLOR.BOUND.LIGHT.2.WALKERS: | |
| (0,1): 8 | |
| (0,1): 5 | |
| (1,0): 6 | |
| 08.ST.NOCOLOR.BOUND.LIGHT.4.WALKERS: | |
| (0,1): 7 | |
| (0,1): 6 | |
| (1,0): 6 | |
| 09.ST.NOCOLOR.UNBOUND.LIGHT.2.WALKERS: | |
| (0,1): 5 | |
| (1,0): 'NONE' | |
| 10.ST.NOCOLOR.UNBOUND.LIGHT.4.WALKERS: | |
| (0,1): 6 | |
| (1,0): 'NONE' | |
| K.CROWN.ETHER: | |
| 05.EQ.COLOR.BOUND.LIGHT.2.WALKERS: | |
| (0,1): 12 | |
| (1,0): 27 | |
| 06.EQ.COLOR.BOUND.LIGHT.4.WALKERS: | |
| (0,1): 7 | |
| (1,0): 29 | |
| 07.ST.NOCOLOR.BOUND.LIGHT.2.WALKERS: | |
| (0,1): 'NONE' | |
| (0,1): 20 | |
| (1,0): 'NONE' | |
| 08.ST.NOCOLOR.BOUND.LIGHT.4.WALKERS: | |
| (0,1): 'NONE' | |
| (0,1): 8 | |
| (1,0): 27 | |
| 09.ST.NOCOLOR.UNBOUND.LIGHT.2.WALKERS: | |
| (0,1): 20 | |
| (1,0): 'NONE' | |
| 10.ST.NOCOLOR.UNBOUND.LIGHT.4.WALKERS: | |
| (0,1): 8 | |
| (1,0): 'NONE' |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| #! /usr/bin/python | |
| # SOME SORT OF COPYRIGHT OR LICENCE, I'M SURE. | |
| # IMPORT STATEMENTS | |
| import numpy as np | |
| import pandas as pd | |
| import seaborn as sns | |
| import h5py | |
| import os | |
| from matplotlib import pyplot as plt | |
| import pickle | |
| # Begin functions. The general idea is a two to N layer pipeline (3, at a minimum, unless we're plotting directly from the file). | |
| # After each step, until we hit the final step (the plotting interface), we are always passing around | |
| # a PANDAS dataframe. That way, we can add, run calculations, etc. | |
| # The first step in the pipeline is loading the data from the HDF5 file. | |
| # Afterwards, we perform some calculation, if necessary. | |
| # Finally, we call the type of plotting function we want. | |
| # As an example, first we load the kinetics files. | |
| def __loadHDF5__(h5name, h5key, i=None, j=None, dt=None, div_dt=False): | |
| # First, load up the file. | |
| # Our data may be shaped in a very complex manner, which we'll want to accomodate for. | |
| # We then index by timepoint such that we know how long the dataset is, simulation wise. | |
| h5file = h5py.File(h5name, 'r') | |
| # We need to try for 0, 1, and 2 dimensional data. I'm leaving it up to the calling function to figure out what it wants, however. | |
| if i == None and j == None: | |
| df = pd.DataFrame(h5file[h5key][:], columns=h5file[h5key].dtype.names, index=range(0, h5file[h5key].shape[0])) | |
| elif j == None: | |
| df = pd.DataFrame(h5file[h5key][:,i], columns=h5file[h5key].dtype.names, index=range(0, h5file[h5key].shape[0])) | |
| else: | |
| df = pd.DataFrame(h5file[h5key][:,i,j], columns=h5file[h5key].dtype.names, index=range(0, h5file[h5key].shape[0])) | |
| # Now add in timedata, if necessary. | |
| # This is simply the end_point * dt, if available. We'll use this in future calculations. | |
| # If we don't specify it, we clearly don't care about it. | |
| if 'iter_stop' in df and dt is not None: | |
| df['timepoint'] = pd.Series((df['iter_stop']-1)*dt) | |
| # For kinetic datasets, we want to be sure we're using dt. Ergo, we may wish to convert to real time. | |
| if div_dt: | |
| if 'expected' in df: | |
| df['expected'] /= dt | |
| if 'ci_lbound' in df: | |
| df['ci_lbound'] /= dt | |
| if 'ci_ubound' in df: | |
| df['ci_ubound'] /= dt | |
| return df | |
| def __createDF__(n, columns=None): | |
| # Here, we're just creating an empty dataframe of size with the columns of expected, ci_lbound, ci_ubound | |
| # The idea is we can generate this and use it to place data into. | |
| if columns == None: | |
| columns = ['expected', 'ci_ubound', 'ci_lbound', 'timepoint'] | |
| df = pd.DataFrame(np.zeros((n,len(columns))), columns=columns, index=range(0, n)) | |
| return df | |
| def __placePointIntoNewSet__(df_in, df_out, oi, ii=-1): | |
| # The idea here is we just create a DF and put in the (typically final point, but not always) | |
| df_out[oi] = df_in[ii] | |
| return df_out | |
| def __WESTCalculateAggregateTime__(df_summary, df_kin, dt=1): | |
| # Here, we calculate the aggregate time from a WESTPA dataset. | |
| df_kin['aggregate_time'] = np.cumsum(df_summary['n_particles']*dt) | |
| return df_kin | |
| def __efficiency_function__(df_in, df_ref, ref_i=None): | |
| # Using the PANDAS dataset, calculate the efficiency in a vectorized fashion. | |
| # I should really rename all this crap. | |
| # ref_i is the timepoint that we're using to calculate the reference time. If None, well, then stop sucking. | |
| # We should have expected, ci_ubound, and ci_lbound, and timepoint. | |
| if ref_i == -1: | |
| # We use this as a shortcut to calculate it relative to the last timepoint. Not sure why pandas can't do this already. | |
| ref_i = df_ref.shape[0]-1 | |
| err_in = ((df_in['ci_ubound'] - df_in['ci_lbound'])/(2*(df_in['expected'])))**2 | |
| err_ref = (((df_ref['ci_ubound'] - df_ref['ci_lbound'])/(2*(df_ref['expected'])))**2)[ref_i] | |
| tim_in = df_in['aggregate_time'] | |
| tim_ref = df_ref['aggregate_time'][ref_i] | |
| return (tim_ref/tim_in)*(err_ref/err_in) | |
| def WESTCreateDataFrame(directory, dt, i, j, west, kin, h5key, div_dt = False): | |
| west = __loadHDF5__(os.path.join(directory, west), h5key='summary', i=None, j=None, div_dt=div_dt) | |
| df = __loadHDF5__(os.path.join(directory, kin), h5key, i=i, j=j, dt=dt, div_dt=div_dt) | |
| df = __WESTCalculateAggregateTime__(west, df, dt=dt) | |
| df.name = directory | |
| return df | |
| def WESTTimeEfficiency(s_df, r_df, name=None): | |
| # Take already existing dataframes, and calculate efficiencies from them. | |
| # Use WESTCreateDataFrame to load thes appropriately. | |
| if name == None: | |
| s_df['efficiency_{}'.format(r_df.name)] = __efficiency_function__(s_df, r_df, ref_i=-1) | |
| else: | |
| s_df['efficiency_{}'.format(name)] = __efficiency_function__(s_df, r_df, ref_i=-1) | |
| return s_df | |
| def WESTCalculateKD(df_d, df_a, df_out, kinetics = True, conc=1): | |
| # d_df = dissociation dataframe, a_df = association data_frame. | |
| # Given that working with very high dimensional dataframes is just... why, really? | |
| # We just split it into one dataframe per direction. This is more convenient with | |
| # steady state simulations anyway, so I don't have to do anything strange to handle them. | |
| length = min(df_d['expected'].shape[0], df_a['expected'].shape[0]) | |
| if kinetics == True: | |
| df_out['expected'] = (df_d['expected'][:length]/df_a['expected'][:length]) | |
| else: | |
| # This might work? | |
| #df_out['expected'] = (df_d['expected'][:length]**2)/(df_a['expected'][:length]) | |
| df_out['expected'] = ((df_d['expected'][:length]**2)/(df_a['expected'][:length]))*conc | |
| err_d = ((df_d['ci_ubound'][:] - df_d['ci_lbound'][:])/(2*df_d['expected'][:])) | |
| err_a = ((df_a['ci_ubound'][:] - df_a['ci_lbound'][:])/(2*df_a['expected'][:])) | |
| err = np.sqrt(err_d[:length]**2 + err_a[:length]**2) | |
| df_out['ci_ubound'] = df_out['expected'] + err | |
| df_out['ci_lbound'] = df_out['expected'] - err | |
| return df_out | |
| # Now some dictionary functions to load, create, and store data. | |
| def WESTDataFrames(directories, west, direct, reweight, dt, a=1, d=0, conc=1): | |
| # Directories is an nxn numpy array containing the directories of simulations we wish to analyze. | |
| n = directories.shape[0] | |
| DataDict = {} | |
| DataDict['reweight'] = {} | |
| DataDict['direct'] = {} | |
| # First, Kinetics and populations.. | |
| for i in range(0, n): | |
| DataDict['direct'][i] = WESTCreateDataFrame(directories[i,i], west=west, kin=direct, h5key='color_prob_evolution', i=i, j=None, dt=dt, div_dt=False) | |
| DataDict['reweight'][i] = WESTCreateDataFrame(directories[i,i], west=west, kin=reweight, h5key='color_prob_evolution', i=i, j=None, dt=dt, div_dt=False) | |
| for j in range(0, n): | |
| if i != j: | |
| # One thing we want to do is correct for any concentration effects, if applicable. That means dividing through by both dt and the concentration FOR THE ON RATE ONLY. | |
| if i == d and j == a: | |
| dc = dt*conc | |
| else: | |
| dc = dt | |
| DataDict['direct'][i,j] = WESTCreateDataFrame(directories[i,j], west=west, kin=direct, h5key='rate_evolution', i=i, j=j, dt=dc, div_dt=True) | |
| DataDict['reweight'][i,j] = WESTCreateDataFrame(directories[i,j], west=west, kin=reweight, h5key='rate_evolution', i=i, j=j, dt=dc, div_dt=True) | |
| # Now, we're going to calculate the KD using both kinetics and population. | |
| # State a is the association, state d is the dissociation, etc. | |
| # Use the direct for kinetics, and the reweighting for population. | |
| length = min(DataDict['direct'][a,d]['expected'].shape[0], DataDict['direct'][d,a].shape[0]) | |
| DataDict['KD_k'] = WESTCalculateKD(DataDict['direct'][a,d], DataDict['direct'][d,a], __createDF__(n=length), kinetics = True) | |
| DataDict['KD_p'] = WESTCalculateKD(DataDict['reweight'][d], DataDict['reweight'][a], __createDF__(n=length), kinetics = False, conc=conc) | |
| return DataDict | |
| def WESTPaper1SimulationsDict(): | |
| # Just a convenience function to do the work for me. Because why wouldn't I have it? | |
| LOC = '/home/varus/work/AdamLTC1Work/posttools/MOLECULAR.DISSOCIATION.SYSTEMS' | |
| SYSTEMS = ['CH4.CH4', 'NA.CL', 'K.CROWN.ETHER'] | |
| DT = { 'CH4.CH4': 0.5, 'NA.CL': 5.0, 'K.CROWN.ETHER': 0.5 } | |
| CONC = { 'CH4.CH4': 0.1192, 'NA.CL': 0.1196, 'K.CROWN.ETHER': 0.03709 } | |
| # We'll need to take the basis state time into account (not that it seems to have a large effect on anything). | |
| BF = { 'CH4.CH4': '/home/varus/work/4.SYSTEMS.BF/CH4/odld.normaltau', 'NA.CL': '/home/varus/work/4.SYSTEMS.BF/NACL/odld', 'K.CROWN.ETHER': '/home/varus/work/4.SYSTEMS.BF/KCROWN/odld' } | |
| BF_DT = { 'CH4.CH4': 2000.01, 'NA.CL': 1000.05, 'K.CROWN.ETHER': 200.001 } | |
| direct = 'ANALYSIS/FINESTATES/direct.h5' | |
| reweight = 'ANALYSIS/FINESTATES/reweight.h5' | |
| HW_E = ['01.EQ.NOCOLOR.BOUND.HEAVY', '02.EQ.COLOR.BOUND.HEAVY'] | |
| HW_S = ['03.ST.NOCOLOR.BOUND.HEAVY', '04.ST.NOCOLOR.UNBOUND.HEAVY'] | |
| LW_E = ['05.EQ.COLOR.BOUND.LIGHT.2.WALKERS', '06.EQ.COLOR.BOUND.LIGHT.4.WALKERS'] | |
| LW_S = ['07.ST.NOCOLOR.BOUND.LIGHT.2.WALKERS', '08.ST.NOCOLOR.BOUND.LIGHT.4.WALKERS', '09.ST.NOCOLOR.UNBOUND.LIGHT.2.WALKERS', '10.ST.NOCOLOR.UNBOUND.LIGHT.4.WALKERS'] | |
| CS_D = [ '34.ST.NOCOLOR.HEAVY', '79.ST.NOCOLOR.LIGHT.2.WALKERS', '810.ST.NOCOLOR.LIGHT.4.WALKERS'] | |
| MASTER = {} | |
| for sys in SYSTEMS: | |
| MASTER[sys] = {} | |
| # Brute Force, first. | |
| directories = np.array([[BF[sys],BF[sys]],[BF[sys],BF[sys]]]) | |
| MASTER[sys]['BF'] = WESTDataFrames(directories,west='west.h5', direct=direct,reweight=reweight,dt=BF_DT[sys], conc=CONC[sys]) | |
| # Now, do the heavyweight. | |
| for hw in HW_E: | |
| directories = np.array([[os.path.join(LOC,sys,hw),os.path.join(LOC,sys,hw)],[os.path.join(LOC,sys,hw),os.path.join(LOC,sys,hw)]]) | |
| MASTER[sys][hw] = WESTDataFrames(directories,west='west.h5', direct=direct,reweight=reweight,dt=DT[sys], conc=CONC[sys]) | |
| # Now, do the SS. | |
| directories = np.array([[os.path.join(LOC,sys,HW_S[0]),os.path.join(LOC,sys,HW_S[0])],[os.path.join(LOC,sys,HW_S[1]),os.path.join(LOC,sys,HW_S[1])]]) | |
| MASTER[sys][HW_S[0]] = WESTDataFrames(directories,west='west.h5', direct=direct,reweight=reweight,dt=DT[sys], conc=CONC[sys]) | |
| MASTER[sys][HW_S[1]] = WESTDataFrames(directories,west='west.h5', direct=direct,reweight=reweight,dt=DT[sys], conc=CONC[sys]) | |
| for lw in LW_E: | |
| directories = np.array([[os.path.join(LOC,sys,lw),os.path.join(LOC,sys,lw)],[os.path.join(LOC,sys,lw),os.path.join(LOC,sys,lw)]]) | |
| MASTER[sys][lw] = WESTDataFrames(directories,west='west.h5', direct=direct,reweight=reweight,dt=DT[sys], conc=CONC[sys]) | |
| # We'll do the cumulative, now. | |
| MASTER[sys][lw]['CUMULATIVE'] = {} | |
| for i in range(2,51): | |
| directories = np.array([[os.path.join(LOC,sys,lw,'CUMULATIVE',str(i).zfill(2)),os.path.join(LOC,sys,lw,'CUMULATIVE',str(i).zfill(2))],[os.path.join(LOC,sys,lw,'CUMULATIVE',str(i).zfill(2)),os.path.join(LOC,sys,lw,'CUMULATIVE',str(i).zfill(2))]]) | |
| MASTER[sys][lw]['CUMULATIVE'][i] = WESTDataFrames(directories,west='west.h5', direct=direct,reweight=reweight,dt=DT[sys], conc=CONC[sys]) | |
| for z in range(0,2): | |
| directories = np.array([[os.path.join(LOC,sys,LW_S[z]),os.path.join(LOC,sys,LW_S[z])],[os.path.join(LOC,sys,LW_S[z+2]),os.path.join(LOC,sys,LW_S[z+2])]]) | |
| MASTER[sys][LW_S[z]] = WESTDataFrames(directories,west='west.h5', direct=direct,reweight=reweight,dt=DT[sys], conc=CONC[sys]) | |
| MASTER[sys][LW_S[z+2]] = WESTDataFrames(directories,west='west.h5', direct=direct,reweight=reweight,dt=DT[sys], conc=CONC[sys]) | |
| # We'll do the cumulative, now. | |
| MASTER[sys][LW_S[z]]['CUMULATIVE'] = {} | |
| MASTER[sys][LW_S[z+2]]['CUMULATIVE'] = {} | |
| for i in range(2,51): | |
| directories = np.array([[os.path.join(LOC,sys,LW_S[z],'CUMULATIVE',str(i).zfill(2)),os.path.join(LOC,sys,LW_S[z],'CUMULATIVE',str(i).zfill(2))],[os.path.join(LOC,sys,LW_S[z+2],'CUMULATIVE',str(i).zfill(2)),os.path.join(LOC,sys,LW_S[z+2],'CUMULATIVE',str(i).zfill(2))]]) | |
| MASTER[sys][LW_S[z]]['CUMULATIVE'][i] = WESTDataFrames(directories,west='west.h5', direct=direct,reweight=reweight,dt=DT[sys], conc=CONC[sys]) | |
| MASTER[sys][LW_S[z+2]]['CUMULATIVE'][i] = WESTDataFrames(directories,west='west.h5', direct=direct,reweight=reweight,dt=DT[sys], conc=CONC[sys]) | |
| return MASTER | |
| #SIMS = WESTPaper1SimulationsDict() | |
| #pickle.dump(SIMS, open("SIMS.p", "wb")) | |
| #SIMS = pickle.load(open("SIMS.p", "rb")) | |
| #PAPER1Figure(SIMS) | |
| #kin = __loadHDF5__('/run/media/varus/AdamLTC1Work/MOLECULAR.DISSOCIATION.SYSTEMS/CH4.CH4/02.EQ.COLOR.BOUND.HEAVY/ANALYSIS/FINESTATES/direct.h5', 'rate_evolution', i=0, j=1, dt=0.5) | |
| #west = __loadHDF5__('/run/media/varus/AdamLTC1Work/MOLECULAR.DISSOCIATION.SYSTEMS/CH4.CH4/02.EQ.COLOR.BOUND.HEAVY/west.h5', 'summary') | |
| #kin = __WESTCalculateAggregateTime__(west, kin, dt=0.5) | |
| #print(__efficiency_function__(kin, kin, ref_i=-1)) | |
| #deff WESTNewTimeEfficiency(simulation, reference, s_dt, r_dt, i, j, s_w='west.h5', s_k='ANALYSIS/FINESTATES/direct.h5', r_w=None, r_k=None, h5key='rate_evolution'): | |
| #s_df = WESTCreateDataFrame(simulation, s_dt, i, j, s_w, s_k, h5key) | |
| # Uses the defaults and creates a dataframe. | |
| #def WESTCreateDataFrame(directory, dt, i, j, west, kin, h5key): | |
| #sim = WESTCreateDataFrame('/home/donimo/234cbc6f23862346aa24666a923ff927/02.EQ.COLOR.BOUND.HEAVY', west='west.h5', kin='ANALYSIS/FINESTATES/direct.h5', h5key='rate_evolution', i=0, j=1, dt=0.5) | |
| #directories = np.array([['/home/donimo/234cbc6f23862346aa24666a923ff927/02.EQ.COLOR.BOUND.HEAVY','/home/donimo/234cbc6f23862346aa24666a923ff927/02.EQ.COLOR.BOUND.HEAVY'],['/home/donimo/234cbc6f23862346aa24666a923ff927/02.EQ.COLOR.BOUND.HEAVY','/home/donimo/234cbc6f23862346aa24666a923ff927/02.EQ.COLOR.BOUND.HEAVY']]) | |
| #sim = WESTDataFrames(directories, west='west.h5', direct='ANALYSIS/FINESTATES/direct.h5', reweight='ANALYSIS/FINESTATES/reweight.h5', dt=0.5, conc=0.1192) | |
| #ref = WESTCreateDataFrame('/home/donimo/234cbc6f23862346aa24666a923ff927/02.EQ.COLOR.BOUND.HEAVY', west='west.h5', kin='ANALYSIS/FINESTATES/direct.h5', h5key='rate_evolution', i=0, j=1, dt=0.5) | |
| #sim = WESTTimeEfficiency(sim, ref) | |
| #print(sim['KD_k']) | |
| #print(sim['KD_p']) | |
| #print(sim['direct'][0,1]) | |
| #kin2 = __createDF__(50) | |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| #! /usr/bin/python | |
| # SOME SORT OF COPYRIGHT OR LICENCE, I'M SURE. | |
| # IMPORT STATEMENTS | |
| import numpy as np | |
| import pandas as pd | |
| import seaborn as sns | |
| import h5py | |
| import os | |
| from WESTData import * | |
| import matplotlib | |
| #matplotlib.use('Agg') | |
| #matplotlib.rc('text', usetex=True) | |
| #matplotlib.rcParams['text.latex.preamble']=[r"\usepackage{amsmath}"] | |
| from matplotlib import pyplot as plt | |
| plt.switch_backend('Agg') | |
| import pickle | |
| import yaml | |
| sns.set_style("white") | |
| sns.set_style('ticks') | |
| sns.set_context('paper') | |
| sns.axes_style({'font.family': ['monospace'], | |
| 'font.sans-serif': ['monospace'] | |
| }) | |
| sns.set(font='sans-serif', style='ticks') | |
| #sns.set_palette('husl') | |
| sns.set_palette('deep') | |
| class FigureManager(): | |
| def __init__(WEST): | |
| self.systems = ['CH4.CH4', 'NA.CL', 'K.CROWN.ETHER'] | |
| self.bf = { 'CH4.CH4': '/home/varus/work/4.SYSTEMS.BF/CH4/odld.normaltau', 'NA.CL': '/home/varus/work/4.SYSTEMS.BF/NACL/odld', 'K.CROWN.ETHER': '/home/varus/work/4.SYSTEMS.BF/KCROWN/odld' } | |
| self.WEST = WEST | |
| self.calculate_efficiency() | |
| def calculate_efficiency(self): | |
| # This does all the efficiencies. | |
| HW = ['01.EQ.NOCOLOR.BOUND.HEAVY', '02.EQ.COLOR.BOUND.HEAVY', '03.ST.NOCOLOR.BOUND.HEAVY', '04.ST.NOCOLOR.UNBOUND.HEAVY'] | |
| LW = ['05.EQ.COLOR.BOUND.LIGHT.2.WALKERS', '06.EQ.COLOR.BOUND.LIGHT.4.WALKERS', '07.ST.NOCOLOR.BOUND.LIGHT.2.WALKERS', '08.ST.NOCOLOR.BOUND.LIGHT.4.WALKERS'] | |
| BF = { 'CH4.CH4': '/home/varus/work/4.SYSTEMS.BF/CH4/odld.normaltau', 'NA.CL': '/home/varus/work/4.SYSTEMS.BF/NACL/odld', 'K.CROWN.ETHER': '/home/varus/work/4.SYSTEMS.BF/KCROWN/odld' } | |
| WEST = self.WEST | |
| for sys in self.systems: | |
| for lw in HW: | |
| for ikey, key in enumerate([(0,1),(1,0),'KD_k', 'KD_p']): | |
| #print(WEST[sys][lw].keys()) | |
| if key == (0,1) or key == (1,0): | |
| #print(dir(WEST[sys][lw]['direct'][key])) | |
| WESTTimeEfficiency(WEST[sys][lw]['direct'][key], WEST[sys]['BF']['direct'][key], name='BF') | |
| #print(WEST[sys][lw]['d=irect'][key]) | |
| else: | |
| WESTTimeEfficiency(WEST[sys][lw][key], WEST[sys]['BF'][key], name='BF') | |
| #print(WEST[sys][lw][key]) | |
| for lw in LW: | |
| for ikey, key in enumerate([(0,1),(1,0),'KD_k', 'KD_p']): | |
| #print(WEST[sys][lw].keys()) | |
| if key == (0,1) or key == (1,0): | |
| #print(dir(WEST[sys][lw]['direct'][key])) | |
| WESTTimeEfficiency(WEST[sys][lw]['direct'][key], WEST[sys]['BF']['direct'][key], name='BF') | |
| #print(WEST[sys][lw]['direct'][key]) | |
| else: | |
| WESTTimeEfficiency(WEST[sys][lw][key], WEST[sys]['BF'][key], name='BF') | |
| #print(WEST[sys][lw][key]) | |
| for sys in self.systems: | |
| for lw in LW: | |
| for ikey, key in enumerate([(0,1),(1,0),'KD_k', 'KD_p']): | |
| #print(WEST[sys][lw].keys()) | |
| for i in range(2,51): | |
| if key == (0,1) or key == (1,0): | |
| #print(dir(WEST[sys][lw]['direct'][key])) | |
| WESTTimeEfficiency(WEST[sys][lw]['CUMULATIVE'][i]['direct'][key], WEST[sys]['BF']['direct'][key], name='BF') | |
| else: | |
| WESTTimeEfficiency(WEST[sys][lw]['CUMULATIVE'][i][key], WEST[sys]['BF'][key], name='BF') | |
| self.WEST = WEST | |
| def PAPER1Figure1(WEST): | |
| # Okay, here we have the dictionary loaded (I think) and all ready to go. Do some analysis! | |
| # First, we want to run the efficiency for each relative to the brute force. | |
| SYSTEMS = ['CH4.CH4', 'NA.CL', 'K.CROWN.ETHER'] | |
| START = yaml.load(open("cumulative.yaml", 'rb')) | |
| # Either one is good, actually, due to the way we set this up. | |
| LW = ['05.EQ.COLOR.BOUND.LIGHT.2.WALKERS', '06.EQ.COLOR.BOUND.LIGHT.4.WALKERS', '07.ST.NOCOLOR.BOUND.LIGHT.2.WALKERS', '08.ST.NOCOLOR.BOUND.LIGHT.4.WALKERS'] | |
| BF = { 'CH4.CH4': '/home/varus/work/4.SYSTEMS.BF/CH4/odld.normaltau', 'NA.CL': '/home/varus/work/4.SYSTEMS.BF/NACL/odld', 'K.CROWN.ETHER': '/home/varus/work/4.SYSTEMS.BF/KCROWN/odld' } | |
| # Now, calculate the efficiencies. | |
| for sys in SYSTEMS: | |
| for lw in LW: | |
| for ikey, key in enumerate([(0,1),(1,0),'KD_k', 'KD_p']): | |
| #print(WEST[sys][lw].keys()) | |
| for i in range(2,51): | |
| if key == (0,1) or key == (1,0): | |
| #print(dir(WEST[sys][lw]['direct'][key])) | |
| WESTTimeEfficiency(WEST[sys][lw]['CUMULATIVE'][i]['direct'][key], WEST[sys]['BF']['direct'][key], name='BF') | |
| else: | |
| WESTTimeEfficiency(WEST[sys][lw]['CUMULATIVE'][i][key], WEST[sys]['BF'][key], name='BF') | |
| # Okay, done. Now we should be able to simply plot them, yes...? | |
| fig, ax = plt.subplots(3,4, sharey=False, figsize=(7,4)) | |
| handles = [] | |
| labels = ['EQ2', 'EQ4', 'SS2', 'SS4'] | |
| fontsize = 6 | |
| titlesize = 8 | |
| fig.suptitle('Simulation Efficiency vs. # of Simulations', fontsize=titlesize, fontweight='bold') | |
| #Ylabels = [r"$\boldsymbol{\mathrm{CH_4/CH_4}}$", 'NA.CL', 'K.CE'] | |
| #Ylabels = [r"$\mathrm{CH_4/CH_4}$", 'NA.CL', 'K.CE'] | |
| Ylabels = [r"$\mathbf{CH_4/CH_4}}$", '$\mathbf{Na^+/Cl^-}$', '$\mathbf{K^+/CE}$'] | |
| xtick_labels = ['LW2', 'LW4', 'HW50'] | |
| TITLES = [r"$\mathbf{K_{on}}$", r"$\mathbf{K_{off}}$", r"$\mathbf{K_{D} = k_{off}/k_{on}}$", r"$\mathbf{K_{D} = p_U^2/p_{B}}$"] | |
| for ipass in range(0, 2): | |
| for isys, sys in enumerate(SYSTEMS): | |
| for ilw, lw in enumerate(LW): | |
| for ikey, key in enumerate([(0,1),(1,0),'KD_k', 'KD_p']): | |
| ax[isys,ikey].tick_params(axis='both', labelsize=fontsize, length=3) | |
| # We're plotting the efficiency as a function of adding more simulations, so create the dataset. | |
| sims = np.zeros(51) | |
| # We have the key and the value. So... | |
| # The efficiency is done relative to the directory name of the system they're compared against. So, for BF, we need the folder name. | |
| for i in range(2,51): | |
| if key == (0,1) or key == (1,0): | |
| length = WEST[sys][lw]['CUMULATIVE'][i]['direct'][key]['efficiency_{}'.format('BF')].shape[0]-1 | |
| sims[i] = (WEST[sys][lw]['CUMULATIVE'][i]['direct'][key]['efficiency_{}'.format('BF')])[length] | |
| else: | |
| length = WEST[sys][lw]['CUMULATIVE'][i][key]['efficiency_{}'.format('BF')].shape[0]-1 | |
| sims[i] = (WEST[sys][lw]['CUMULATIVE'][i][key]['efficiency_{}'.format('BF')])[length] | |
| if isys == 0: | |
| ax[isys,ikey].set_title(TITLES[ikey], fontsize=titlesize, fontweight='bold') | |
| if ikey == 0: | |
| ax[isys,ikey].set_ylabel(Ylabels[isys], fontsize=titlesize, fontweight='bold') | |
| if key == (0,1) or key == (1,0): | |
| print(START[sys][lw]) | |
| POINT = START[sys][lw][str(key).replace(' ', '')] | |
| if POINT == 'NONE': | |
| POINT = 50 | |
| X = range(POINT, 51) | |
| if ipass == 0: | |
| ax[isys,ikey].plot(range(2,51), sims[2:], color='black', alpha=.3) | |
| else: | |
| handles.append(ax[isys,ikey].plot(X, sims[POINT:], label=labels[ilw])) | |
| #ax[isys,ikey].vlines(POINT, sims[POINT]-1, sims[POINT]+1, colors=handles[-1][0].get_color()) | |
| ax[isys,ikey].scatter(POINT, sims[POINT], s=20, color=handles[-1][0].get_color()) | |
| else: | |
| # If it's the KD, plot the maximum start point. | |
| POINT = max(START[sys][lw][str((0,1)).replace(' ', '')], START[sys][lw][str((1,0)).replace(' ', '')]) | |
| if POINT == 'NONE': | |
| POINT = 50 | |
| if ipass == 0: | |
| ax[isys,ikey].plot(range(2,51), sims[2:], color='black', alpha=.3) | |
| else: | |
| X = range(POINT, 51) | |
| handles.append(ax[isys,ikey].plot(X, sims[POINT:], label=labels[ilw])) | |
| #ax[isys,ikey].vlines(POINT, sims[POINT]-1, sims[POINT]+1, colors=handles[-1][0].get_color()) | |
| ax[isys,ikey].scatter(POINT, sims[POINT], s=20, color=handles[-1][0].get_color()) | |
| #if isys == 2: | |
| #ax[itype,ikey].set_xticks([0,1,2,4,5,6,8,9,10]) | |
| #ax[itype,ikey].set_xticklabels(xtick_labels*3, rotation=45) | |
| if isys != 2: | |
| ax[isys,ikey].set_xticks([]) | |
| #plt.figlegend(handles[:4], labels[:4], loc='upper right') | |
| #plt.legend(loc=(-1.75,-0.5), ncol=4, fontsize=fontsize) | |
| plt.figlegend(loc='lower center',ncol=4, fontsize=fontsize) | |
| plt.tight_layout() | |
| sns.despine() | |
| plt.savefig('LWvsBF.pdf') | |
| def PAPER1Figure2(WEST): | |
| # Okay, here we have the dictionary loaded (I think) and all ready to go. Do some analysis! | |
| # First, we want to run the efficiency for each relative to the brute force. | |
| SYSTEMS = ['CH4.CH4', 'NA.CL', 'K.CROWN.ETHER'] | |
| # Either one is good, actually, due to the way we set this up. | |
| HW = ['01.EQ.NOCOLOR.BOUND.HEAVY', '02.EQ.COLOR.BOUND.HEAVY', '03.ST.NOCOLOR.BOUND.HEAVY', '04.ST.NOCOLOR.UNBOUND.HEAVY'] | |
| LW = ['05.EQ.COLOR.BOUND.LIGHT.2.WALKERS', '06.EQ.COLOR.BOUND.LIGHT.4.WALKERS', '07.ST.NOCOLOR.BOUND.LIGHT.2.WALKERS', '08.ST.NOCOLOR.BOUND.LIGHT.4.WALKERS'] | |
| BF = { 'CH4.CH4': '/home/varus/work/4.SYSTEMS.BF/CH4/odld.normaltau', 'NA.CL': '/home/varus/work/4.SYSTEMS.BF/NACL/odld', 'K.CROWN.ETHER': '/home/varus/work/4.SYSTEMS.BF/KCROWN/odld' } | |
| # Now, calculate the efficiencies. | |
| # If we've run the first, this is fine to just do the HW. | |
| for sys in SYSTEMS: | |
| for lw in HW: | |
| for ikey, key in enumerate([(0,1),(1,0),'KD_k', 'KD_p']): | |
| #print(WEST[sys][lw].keys()) | |
| if key == (0,1) or key == (1,0): | |
| #print(dir(WEST[sys][lw]['direct'][key])) | |
| WESTTimeEfficiency(WEST[sys][lw]['direct'][key], WEST[sys]['BF']['direct'][key], name='BF') | |
| #print(WEST[sys][lw]['d=irect'][key]) | |
| else: | |
| WESTTimeEfficiency(WEST[sys][lw][key], WEST[sys]['BF'][key], name='BF') | |
| #print(WEST[sys][lw][key]) | |
| for lw in LW: | |
| for ikey, key in enumerate([(0,1),(1,0),'KD_k', 'KD_p']): | |
| #print(WEST[sys][lw].keys()) | |
| if key == (0,1) or key == (1,0): | |
| #print(dir(WEST[sys][lw]['direct'][key])) | |
| WESTTimeEfficiency(WEST[sys][lw]['direct'][key], WEST[sys]['BF']['direct'][key], name='BF') | |
| #print(WEST[sys][lw]['direct'][key]) | |
| else: | |
| WESTTimeEfficiency(WEST[sys][lw][key], WEST[sys]['BF'][key], name='BF') | |
| #print(WEST[sys][lw][key]) | |
| # Okay, done. Now we should be able to simply plot them, yes...? | |
| fig, ax = plt.subplots(2,4, sharey=False, figsize=(7,4)) | |
| handles = [] | |
| labels = [r"$\mathrm{CH_4/CH_4}}$", '$\mathrm{Na^+/Cl^-}$', '$\mathrm{K^+/CE}$'] | |
| xtick_labels = ['LW2', 'LW4', 'HW50'] | |
| fontsize = 6 | |
| titlesize = 8 | |
| fig.suptitle('Simulation Efficiency vs. Target Walker Count', fontsize=titlesize, fontweight='bold') | |
| TITLES = [r"$\mathbf{K_{on}}$", r"$\mathbf{K_{off}}$", r"$\mathbf{K_{D} = k_{off}/k_{on}}$", r"$\mathbf{K_{D} = p_U^2/p_{B}}$"] | |
| Ylabels = [r"$\mathbf{Equilibrium}$", r"Steady State"] | |
| EQ = ['05.EQ.COLOR.BOUND.LIGHT.2.WALKERS', '06.EQ.COLOR.BOUND.LIGHT.4.WALKERS', '02.EQ.COLOR.BOUND.HEAVY'] | |
| SS = ['07.ST.NOCOLOR.BOUND.LIGHT.2.WALKERS', '08.ST.NOCOLOR.BOUND.LIGHT.4.WALKERS', '03.ST.NOCOLOR.BOUND.HEAVY'] | |
| for itype in range(0, 2): | |
| for ikey, key in enumerate([(0,1),(1,0),'KD_k', 'KD_p']): | |
| # Just do this once. | |
| ax[itype,ikey].axhline(1, color='black', lw=0.5, alpha=0.25) | |
| for isys, sys in enumerate(SYSTEMS): | |
| X = [] | |
| VEC = [] | |
| ax[itype,ikey].set_yscale('log') | |
| # None of the fucking minor ticks. | |
| ax[itype,ikey].set_ylim([0.01,1000]) | |
| ax[itype,ikey].tick_params(axis='both', labelsize=fontsize, length=3) | |
| if itype == 0: | |
| SIMS = EQ | |
| else: | |
| SIMS = SS | |
| for ilw, lw in enumerate(SIMS): | |
| if key == (0,1) or key == (1,0): | |
| #handles.append(ax[isys,ikey].plot(sims, label=labels[ilw])) | |
| length = WEST[sys][lw]['direct'][key]['efficiency_BF'].shape[0]-1 | |
| VEC.append(WEST[sys][lw]['direct'][key]['efficiency_BF'][length]) | |
| else: | |
| length = WEST[sys][lw][key]['efficiency_BF'].shape[0]-1 | |
| VEC.append(WEST[sys][lw][key]['efficiency_BF'][length]) | |
| X.append((isys*(len(SIMS)+1))+ilw) | |
| #handles.append(ax[itype,ikey].bar(X, VEC, label=labels[isys])) | |
| # On some versions, it seems to need log=True, otherwise it just... borks out and won't do the bar properly | |
| # (seems to be that it doesn't understand that I'm plotting a log plot). | |
| ax[itype,ikey].bar(X, VEC, label=labels[isys], log=True) | |
| #ax[itype,ikey].scatter(X, VEC, label=labels[isys]) | |
| #sns.barplot(X, VEC, ax=ax[itype,ikey]) | |
| if itype == 1: | |
| ax[itype,ikey].set_xticks([0,1,2,4,5,6,8,9,10]) | |
| ax[itype,ikey].set_xticklabels(xtick_labels*3, rotation=45) | |
| else: | |
| ax[itype,ikey].set_xticks([]) | |
| ax[itype,ikey].set_title(TITLES[ikey], fontsize=titlesize, fontweight='bold') | |
| #ax[itype,ikey].set_xticklabels(xtick_labels*3, rotation=45) | |
| if ikey == 0: | |
| ax[itype,ikey].set_yticks([0.01, 0.1, 1, 25, 50, 100, 1000]) | |
| ax[itype,ikey].set_yticklabels([0.01, 0.1, 1, 25, 50, 100, 1000]) | |
| ax[itype,ikey].set_ylabel(Ylabels[itype], fontsize=titlesize, fontweight='bold') | |
| else: | |
| ax[itype,ikey].set_yticklabels([]) | |
| ax[itype,ikey].set_yticks([]) | |
| ax[itype,ikey].minorticks_off() | |
| plt.figlegend(loc='lower center',ncol=3, fontsize=fontsize) | |
| plt.tight_layout() | |
| sns.despine() | |
| plt.savefig('all_walkers.pdf') | |
| def PAPER1AggregateTime(WEST): | |
| # Okay, here we have the dictionary loaded (I think) and all ready to go. Do some analysis! | |
| # First, we want to run the efficiency for each relative to the brute force. | |
| SYSTEMS = ['CH4.CH4', 'NA.CL', 'K.CROWN.ETHER'] | |
| # Either one is good, actually, due to the way we set this up. | |
| HW = ['01.EQ.NOCOLOR.BOUND.HEAVY', '02.EQ.COLOR.BOUND.HEAVY', '03.ST.NOCOLOR.BOUND.HEAVY', '04.ST.NOCOLOR.UNBOUND.HEAVY'] | |
| LW = ['05.EQ.COLOR.BOUND.LIGHT.2.WALKERS', '06.EQ.COLOR.BOUND.LIGHT.4.WALKERS', '07.ST.NOCOLOR.BOUND.LIGHT.2.WALKERS', '08.ST.NOCOLOR.BOUND.LIGHT.4.WALKERS', '09.ST.NOCOLOR.UNBOUND.LIGHT.2.WALKERS', '10.ST.NOCOLOR.UNBOUND.LIGHT.4.WALKERS'] | |
| BF = { 'CH4.CH4': '/home/varus/work/4.SYSTEMS.BF/CH4/odld.normaltau', 'NA.CL': '/home/varus/work/4.SYSTEMS.BF/NACL/odld', 'K.CROWN.ETHER': '/home/varus/work/4.SYSTEMS.BF/KCROWN/odld' } | |
| # Now, calculate the efficiencies. | |
| # If we've run the first, this is fine to just do the HW. | |
| for isys, sys in enumerate(SYSTEMS): | |
| AGGREGATE = 0 | |
| for ilw, lw in enumerate(HW+LW): | |
| length = WEST[sys][lw]['direct'][(1,0)]['aggregate_time'].shape[0]-1 | |
| AGGREGATE += WEST[sys][lw]['direct'][(1,0)]['aggregate_time'][length] | |
| print('{}: {} us'.format(sys, AGGREGATE/1000/1000)) | |
| #TEMPLATE = yaml.load(open("template.yaml", 'rb')) | |
| #print(TEMPLATE) | |
| SIMS = pickle.load(open("SIMS.p", "rb")) | |
| PAPER1AggregateTime(SIMS) | |
| PAPER1Figure1(SIMS) | |
| PAPER1Figure2(SIMS) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment