diff --git a/analysis/fig_quiet.py b/analysis/fig_quiet.py new file mode 100644 index 0000000..56136fc --- /dev/null +++ b/analysis/fig_quiet.py @@ -0,0 +1,444 @@ +""" +fig_quiet.py + +M1 paper Figure 2 + +Contributors: salvadordura@gmail.com +""" + +from shared import * +from netpyne import analysis + +# ---------------------------------------------------------------- +def fig_quiet(): + ''' Figure showing spontaneous activity: + - raster plot + - spike histogram + - voltage traces + - boxplot stats model and experiment''' + + + # --------------------------------------------------------------------------------------------------------------- + # Config + + raster = 1 + histogram = 1 + traces = 1 + stats_boxplot = 1 # boxplot of rates + + fontsiz = 16 + + dataFolder = '../data/' + if raster or histogram or traces: # 2 sec N=1 + batchLabel = 'v56_batch18' + simLabel = 'v56_batch18_0_0_0' # same as 'v56_batch7_7_3_0_0' + + sim, data, out, root = loadSimData(dataFolder, batchLabel, simLabel) + + elif stats_boxplot: # 50 sec N=1 + + timeRange = [1000, 5000] + include = excpops+[SOMsubset,PVsubset] + xlim = [0, 40] + labelsModel = ['IT2', 'IT4', 'IT5A', 'IT5B', 'PT5B', 'IT6', 'CT6', 'SOM L2-6', 'PV L2-6'] + multipleTrials = True + loadAll = 1 + + # single trial + if not multipleTrials: + batchLabel = 'v56_batch7' + simLabel = 'v56_batch7_3_0_0' + sim, data, out, root = loadSimData(dataFolder, batchLabel, simLabel) + + # multiple trials + else: + batchLabel = 'v56_batch7' + simLabel = ['v56_batch7_3_'+str(iseed)+'_'+str(iconn) for iseed in range(5) for iconn in range(5)] + root = dataFolder+batchLabel+'/' + + if not loadAll: + + popGids, popNumCells = loadPopGids(dataFolder, labelsModel+['SOM2','SOM4','SOM5A','SOM5B','SOM6','PV2','PV4','PV5A','PV5B','PV6']) + popGids['SOM L2-6'] = popGids['SOM2']+popGids['SOM4']+popGids['SOM5A']+popGids['SOM5B']+popGids['SOM6'] + popGids['PV L2-6'] = popGids['PV2']+popGids['PV4']+popGids['PV5A']+popGids['PV5B']+popGids['PV6'] + + # load data from all sim files + statDataAll = [[]]*len(simLabel) + for isim,simLab in enumerate(simLabel): # get data from each sim + filename = root+simLab+'.json' + print(filename) + sim,data,out = utils.plotsFromFile(filename, raster=0, stats=0, rates=0, syncs=0, hist=0, psd=0, traces=0, grang=0, plotAll=0, popColors=popColors) + + statDataAll[isim] = [] + for isubset, subset in enumerate(labelsModel): + cellGids = popGids[subset] + dfSpks, _, _ = analysis.utils.getSpktSpkid(cellGids, timeRange, sim) # get quiet spikes + dfRates = dfSpks.groupby("spkid").count().div((timeRange[1]-timeRange[0])/1000.0) + dfRates = dfRates.reindex(cellGids, fill_value=0.0) + statDataAll[isim].append(list(dfRates.spkt)) + print(subset, len(dfRates.spkt), np.mean(dfRates.spkt), np.median(dfRates.spkt), scipy.stats.iqr(dfRates.spkt)) + + # save + statData = {'statDataAll': statDataAll} #, 'gidsData': gidsData, 'ynormsData': ynormsData} + with open(root+'%s_statDataAll_boxplot.json'%(simLabel[0][:-2]), 'w') as fileObj: + json.dump(statData, fileObj) + + # load All + else: + with open(root+'%s_statDataAll_boxplot.json'%(simLabel[0][:-2]), 'r') as fileObj: + statData = json.load(fileObj) + filename = root+simLabel[0]+'.json' + + + plt.style.use('seaborn-ticks') + + # --------------------------------------------------------------------------------------------------------------- + # raster + if raster: + timeRange = [2000, 4000] + include = allpops + orderBy = ['pop', 'y'] + + + from netpyne.analysis.spikes_legacy import plotRaster + fig1 = plotRaster(include=include, timeRange=timeRange, labels='overlay', + popRates=0, orderInverse=True, lw=0, markerSize=3.5, marker='.', popColors=popColors, + showFig=0, saveFig=0, figSize=(8.5*0.9*0.75, 7), orderBy=orderBy)# + ax = plt.gca() + + [i.set_linewidth(0.5) for i in ax.spines.values()] # make border thinner + plt.xticks([2000, 3000, 4000], ['2', '3', '4'], fontsize=fontsiz) + plt.yticks([0, 5000, 10000], [0, 5000, 10000], fontsize=fontsiz) + + plt.ylabel('Neuron ID', fontsize=fontsiz) #Neurons (ordered by NCD within each pop)') + plt.xlabel('Time (s)', fontsize=fontsiz) + + plt.title('') + filename='%s%s_raster_%d_%d_%s_small.png'%(root, simLabel, timeRange[0], timeRange[1], orderBy) + plt.savefig(filename, dpi=600) + + + # --------------------------------------------------------------------------------------------------------------- + # histogram + if histogram: + + timeRange = [2000, 4000] + include = allpops + orderBy = ['pop', 'y'] + + binSize=5 + measure = 'count' + graphType = 'bar' + + from netpyne.analysis.spikes_legacy import plotSpikeHist + fig_hist = plotSpikeHist(include=excpops, timeRange=timeRange, figSize=(8.5*0.9*0.75, 2.5), + popColors=popColors, legend=False, showFig=0, saveFig=0, linewidth=0.5, binSize=binSize, graphType='bar', + axis=True, measure=measure, scalebarLoc='upper center') + + ax=plt.gca() + ax.get_legend().remove() + ax.spines['right'].set_visible(False) + ax.spines['top'].set_visible(False) + ax.spines['left'].set_linewidth(0.5) + ax.spines['bottom'].set_linewidth(0.5) + + plt.tight_layout() + plt.subplots_adjust(bottom=0.15, top=1.0, right=0.9, left=0.15) + plt.yticks(fontsize=fontsiz) + plt.xticks([2000, 3000, 4000], ['2', '3', '4'], fontsize=fontsiz) + + filename='%s%s_spikehist_%d_%d_bin-%d_%s_%s.png'%(root, simLabel, timeRange[0], timeRange[1], binSize, measure, graphType) + plt.savefig(filename, dpi=600) + + + # --------------------------------------------------------------------------------------------------------------- + # traces + if traces: + # to plot all cells recorded + batchLabel = 'v56_batch23' + simLabel = 'v56_batch23_merged' + + sim, data, out, root = loadSimData(dataFolder, batchLabel, simLabel) + + sim.cfg.recordTraces = {'V_soma': {}} + + fontsiz = 20 + timeRange = [2000, 4000] + + + # plot all cells recorded + ''' + include = list(range(10073)) + timeRange = [2000, 4000] + recordStep = 0.1 + t = np.arange(timeRange[0], timeRange[1]+recordStep, recordStep) + for gid in include: + fullTrace = sim.allSimData['V_soma']['cell_'+str(gid)] + vtrace = np.array(fullTrace[int(timeRange[0]/recordStep):int(timeRange[1]/recordStep)]) + plt.figure(figsize=(12, 4)) + plt.plot(t[:len(vtrace)], vtrace, linewidth=1, color='black') + plt.ylim(-100,40) + plt.xlabel('time (ms)') + plt.ylabel('mV') + plt.savefig('%s/V_traces/v56_batch23_V_soma_time_%d_%d_cell_%d' % (root, timeRange[0], timeRange[1], gid)) + ''' + + + # plot subset of cells with the correct format + from netpyne.support.scalebar import add_scalebar + + def addScaleBar(timeRange=timeRange, loc=1): + ax = plt.gca() + sizex = (timeRange[1]-timeRange[0])/20.0 + add_scalebar(ax, hidex=False, hidey=True, matchx=False, matchy=True, sizex=sizex, sizey=None, unitsx='ms', unitsy='mV', scalex=1, scaley=1, loc=loc, pad=-1, borderpad=0.5, sep=4, prop=None, barcolor="black", barwidth=3) + plt.axis('off') + + # IT/PT cells for quiet + cells = [4504] + timeRange = [2000, 4000] + recordStep = 0.1 + + t = np.arange(timeRange[0], timeRange[1]+recordStep, recordStep) + + for gid in cells: + fullTrace = sim.allSimData['V_soma']['cell_'+str(gid)] + vtrace = np.array(fullTrace[int(timeRange[0]/recordStep):int(timeRange[1]/recordStep)]) + + plt.figure(figsize=(15*1.5*0.5, 2.5)) + plt.plot(t[:len(vtrace)], vtrace, linewidth=1, color=popColors['IT5B']) + baseline = min(vtrace) + plt.ylim(baseline, baseline+60) #TRUNCATE AT 60mV above baseline! + #addScaleBar() + plt.axis('off') + + ax = plt.gca() + plt.title('') + plt.tight_layout() + plt.savefig('%s/V_traces/%s_V_soma_cell_%d_%d_%d_x1.5_axis.png' % (root, simLabel, gid, timeRange[0], timeRange[1]), dpi=200) + + + # --------------------------------------------------------------------------------------------------------------- + # stats boxplot + if stats_boxplot: + + fontsiz = 20 + + # Model stats - single trial (N=1) + if not multipleTrials: + fig1,modelData = sim.analysis.plotSpikeStats(include=include, figSize=(8,4), timeRange=timeRange, xlim=xlim, + stats = ['rate'], legendLabels=labelsModel, fontSize = fontsiz, popColors=popColors, showFig=0, dpi=300, saveFig=False) + + # Model stats - multiple trials (N=25) + else: + + minValue = 0.01 + # combine data + statDataAll = statData['statDataAll'] + + statDataCombined = {} + for ipop in range(len(include)): + statDataCombined[ipop] = [] + for isim in range(len(simLabel)): + nonzeroData = statDataAll[isim][ipop] + nonzeroData = [x for x in nonzeroData if x>=minValue] + statDataCombined[ipop].extend(nonzeroData) + + modelData = {'statData': statDataCombined} + statData = modelData['statData'] + + # save combiend data + with open(root+'%s_statData_boxplot.json'%(batchLabel), 'w') as fileObj: + json.dump(statData, fileObj) + + + # Combine model stats with experimental data stats + expData = {} + from scipy.io import loadmat + + ## Este18 L5 RS and PV baseline + matData = loadmat(dataFolder+'/Este18/data_v2.mat') # [0][0][2][0] = RS; [0][0][2][0] = PV+ + + expData['Este18_L5_RS_baseline'] = list(matData['baseline'][0][0][0][0]) + \ + list(matData['sound_base'][0][0][0][0]) + \ + list(matData['vibration_base'][0][0][0][0]) + \ + list(matData['movement_base'][0][0][0][0]) + + expData['Este18_PV_baseline'] = list(matData['baseline'][0][0][2][0]) + \ + list(matData['sound_base'][0][0][2][0]) + \ + list(matData['vibration_base'][0][0][2][0]) + \ + list(matData['movement_base'][0][0][2][0]) + + ## Zagh15 L5 RS and PV baseline + with open(dataFolder+'/Zagh15/spike_rates.json', 'r') as f: # + jsdonData = json.load(f) + expData['Zagh15_L5_RS_baseline'] = jsdonData['pre-stimulus']['RSU'] + expData['Zagh15_PV_baseline'] = jsdonData['pre-stimulus']['FSU'] + + ## Schi15 + import schi15 + df = schi15.readExcelFiringRatesAndMetada() + df.dropna() + + expData['Schi15_L5B_RS_baseline'] = list(df.query("condition == 'main' and code!='23'").quietFR) # main, quiet (medium VL, medium ih) + expData['Schi15_IT2_RS_baseline'] = list(df.query("condition=='main' and code=='23'").quietFR) # main, quiet (medium VL, medium ih) + expData['Schi15_IT5B_RS_baseline'] = list(df.query("condition=='main' and cell_class=='IT'").quietFR) # main, quiet (medium VL, medium ih) + expData['Schi15_PT5B_RS_baseline'] = list(df.query("condition=='main' and cell_class=='PT'").quietFR) # main, quiet (medium VL, medium ih) N=3! remove + + + ## Dacr19 - waiting + + ## LiCh15 (ALM) L5 RS and PV baseline + matData = loadmat(dataFolder+'/LiCh15/ALM_compiled_all_data.mat') + expData['LiCh15_L5_RS_baseline'] = [] + expData['LiCh15_PV_baseline'] = [] + for i, (x1, x2) in enumerate(zip(list(matData['spk_count_yes_all'][:, 4]),list(matData['spk_count_no_all'][:, 4]))): + if not np.isnan(np.mean([x1,x2])): + if matData['cellType_all'][i, 0] == 1: # Pyr + expData['LiCh15_L5_RS_baseline'].append(np.mean([x1,x2])) + elif matData['cellType_all'][i,0] == 2: # FS + expData['LiCh15_PV_baseline'].append(np.mean([x1,x2])) + + + ## Econ18 (ALM) L5 PT baseline + matData = loadmat(dataFolder+'/Econ18/Econ18.mat') + expData['Econ18_PT5B_RS_baseline'] = [x[0] for x in list(matData['ratePTdown'])+list(matData['ratePTup'])] + + + # ensure exp data values are >= minValue + for k,v in expData.items(): + vmin = [x for x in v if x >= minValue] + expData[k] = vmin + + # combine model and experimental data + + ## statData[0] - IT2 + ## statData[1] - IT4 + ## statData[2] - IT5A + ## statData[3] - IT5B + ## statData[4] - PT5B + + + statData = [statData[0]] \ + + [expData['Schi15_IT2_RS_baseline']] \ + + [statData[1]]+[statData[2]]+[statData[3]] \ + + [expData['Schi15_IT5B_RS_baseline']] \ + + [statData[4]] \ + + [expData['Econ18_PT5B_RS_baseline']] \ + + [statData[3]+statData[4]] \ + + [expData['Schi15_L5B_RS_baseline']] \ + + [statData[2]+statData[3]+statData[4]] \ + + [expData['Este18_L5_RS_baseline'], expData['Zagh15_L5_RS_baseline'], expData['LiCh15_L5_RS_baseline']] \ + + [statData[5]] + [statData[6]] + + labels = [ 'L2/3 IT', + 'L2 IT \n(Schiemann 2015)', \ + 'L4 IT', 'L5A IT', 'L5B IT', \ + 'L5B IT \n(Schiemann 2015)', \ + 'L5B PT', \ + 'L5B PT \n(Economo 2018)', \ + 'L5B IT,PT', \ + 'L5B IT,PT \n(Schiemann 2015)', \ + 'L5 IT,PT', \ + 'L5 IT,PT \n(Estebanez 2018)', 'L5 IT,PT \n(Zagha 2015)', 'L5 IT,PT \n(Li 2015)',\ + 'L6 IT', 'L6 CT'] + + + labels_orig = [ 'IT2/3', + 'IT2 \n(Schiemann 2015)', \ + 'IT4', 'IT5A', 'IT5B', \ + 'IT5B \n(Schiemann 2015)', \ + 'PT5B', \ + 'PT5B \n(Economo 2018)', \ + 'IT5B+PT5B', \ + 'IT5B+PT5B \n(Schiemann 2015)', \ + 'IT5A+IT5B+PT5B', \ + 'IT5A+IT5B+PT5B \n(Estebanez 2018)', 'IT5A+IT5B+PT5B \n(Zagha 2015)', 'IT5A+IT5B+PT5B \n(Li 2015)',\ + 'IT6', 'CT6'] + + popColors['IT5A+IT5B+PT5B'] = '#BC01FF' + popColors['IT5B+PT5B'] = '#800080' + popColors['IT2/3'] = popColors['IT2'] + + # update pop color labels + for i,l in enumerate(labels_orig): + if '(' not in l: + popColors[labels[i]] = popColors[l] + + colors = [popColors[p] if p in popColors else '#999999' for p in labels] + + + # plot boxplot comparing model and experimental data + plt.figure(figsize=(10,14)) + meanpointprops = dict(marker = (5, 1, 0), markeredgecolor = 'black', markerfacecolor = 'white') + + flierprops = dict(marker='.', markerfacecolor='gray', markersize=2,linestyle='none', markeredgecolor='gray') + + bp=plt.boxplot(statData[::-1], labels=labels[::-1], notch=False, flierprops=flierprops, meanprops=meanpointprops, whis=1.5, widths=0.6, vert=False, showmeans=True, patch_artist=True) #labels[::-1] #positions=np.array(range(len(statData)))+0.4, + + # mean and std in errbar + capsize = 8 + lw = 2 + elw=3.0 + marker = 'o' + ms = 15 + mec = 'black' + + x = [np.mean(v) for v in statData[::-1]] + xerr = [np.std(v) for v in statData[::-1]] + y = np.array(range(len(x)))[::-1] + + plt.xlabel('Rate (Hz)', fontsize=fontsiz) + plt.ylabel('', fontsize = fontsiz) + plt.subplots_adjust(left=0.3,right=0.95, top=0.9, bottom=0.1) + + icolor=0 + borderColor = 'k' + for i in range(0, len(bp['boxes'])): + icolor = i + bp['boxes'][i].set_facecolor(colors[::-1][icolor]) + bp['boxes'][i].set_linewidth(2) + # we have two whiskers! + bp['whiskers'][i*2].set_color(borderColor) + bp['whiskers'][i*2 + 1].set_color(borderColor) + bp['whiskers'][i*2].set_linewidth(2) + bp['whiskers'][i*2 + 1].set_linewidth(2) + bp['medians'][i].set_color(borderColor) + bp['medians'][i].set_linewidth(3) + for c in bp['caps']: + c.set_color(borderColor) + c.set_linewidth(2) + + ax = plt.gca() + ax.spines['top'].set_visible(False) + ax.spines['right'].set_visible(False) + ax.spines['bottom'].set_visible(False) + ax.get_xaxis().tick_bottom() + ax.get_yaxis().tick_left() + ax.tick_params(axis='x', length=0) + #ax.tick_params(axis='y', direction='out') + ax.grid(axis='x', color="0.9", linestyle='-', linewidth=1) + ax.set_axisbelow(True) + if xlim: ax.set_xlim(xlim) + axisFontSize(ax, fontsiz-4) + + + filename = root+batchLabel+'_boxplot_%d_%d'%(timeRange[0], timeRange[1]) + #filename = root+batchLabel+'_errbar_%d_%d'%(timeRange[0], timeRange[1]) + plt.savefig(filename, dpi=600) + + + # stats and rank-test p-value + for i,stat in enumerate(statData): + print('%s: %.2f, %.2f, %.2f, %.2f' % (labels[i], np.mean(stat), np.std(stat), np.median(stat), scipy.stats.iqr(stat))) + + # L5 IT,PT exp1 vs exp2 + scipy.stats.mannwhitneyu(statData[11], statData[12]) + + # L5 IT,PT model vs exp2 + scipy.stats.mannwhitneyu(statData[10], statData[12]) + + + +# Main code +if __name__ == '__main__': + fig_quiet() diff --git a/analysis/shared.py b/analysis/shared.py new file mode 100644 index 0000000..a9877c7 --- /dev/null +++ b/analysis/shared.py @@ -0,0 +1,92 @@ +""" +paper.py + +Paper figures + +Contributors: salvadordura@gmail.com +""" + +import utils + +import json +import numpy as np +import scipy +from matplotlib import pyplot as plt +import pandas as pd +import seaborn as sb +import os +import pickle +from netpyne.support.scalebar import add_scalebar +from matplotlib import cm +from matplotlib.colors import ListedColormap + +# --------------------------------------------------------------------------------------------------------------- +# Population params +allpops = ['IT2','SOM2','PV2','IT4','IT5A','SOM5A','PV5A','IT5B','PT5B','SOM5B','PV5B', + 'IT6','CT6','SOM6','PV6'] +excpops = ['IT2','IT4','IT5A','IT5B','PT5B','IT6','CT6'] +inhpops = ['SOM2','PV2', 'SOM5A','PV5A', 'SOM5B','PV5B', 'SOM6', 'PV6'] +excpopsu = ['IT2','IT4','IT5A','PT5B'] +ITsubset = ('IT2', 'IT4', 'IT5A', 'IT6') +SOMsubset = ('SOM2', 'SOM5A','SOM5B', 'SOM6') +PVsubset = ('PV2', 'PV5A', 'PV5B', 'PV6') +try: + with open('../sim/cells/popColors.pkl', 'rb') as fileObj: + popColors = pickle.load(fileObj)['popColors'] + popColors['S2'] = [0.90,0.76,0.00] + popColors['TPO'] = [52/255.0, 138/255.0, 49/255.0] #[232/255.0, 37/255.0, 101/255.0] #'firebrick' #popColors['S2'] #[253/255.0, 102/255.0, 2/255.0] # + popColors['M2'] = [0.42,0.67,0.9] + popColors[ITsubset] = popColors['IT5A'] + popColors[SOMsubset] = popColors['SOM L2-6'] = popColors['SOM2-6'] = popColors['SOM5A'] + popColors[PVsubset] = popColors['PV L2-6'] = popColors['PV2-6'] = popColors['PV5A'] +except: + pass + +def loadSimData(dataFolder, batchLabel, simLabel): + ''' load sim file''' + root = dataFolder+batchLabel+'/' + sim,data,out = None, None, None + if isinstance(simLabel, str): + if os.path.exists(root+simLabel+'.json'): + filename = root + simLabel + '.json' + print(filename) + elif os.path.exists(root+simLabel+'.pkl'): + filename = root + simLabel + '.pkl' + print(filename) + + sim, data, out = utils.plotsFromFile(filename, raster=0, stats=0, rates=0, syncs=0, hist=0, psd=0, traces=0, grang=0, plotAll=0, + include = {'raster': allpops}, + popColors=popColors) + + return sim, data, out, root + +def axisFontSize(ax, fontsize): + for tick in ax.xaxis.get_major_ticks(): + tick.label.set_fontsize(fontsize) + for tick in ax.yaxis.get_major_ticks(): + tick.label.set_fontsize(fontsize) + +def loadPopGids(dataFolder, popLabels): + # load all data + # sim, data, out, root = loadSimData(dataFolder, batchLabel, simLabel) + from netpyne import sim + with open(dataFolder + 'v56_batch2_net/v56_batch2_net_0_0.pkl', 'rb') as f: + netLoad = pickle.load(f) + + allCells = netLoad['net']['cells'] + + L5Bmin=0.47 + L5Bmax=0.8 + L5Bmid = L5Bmin + (L5Bmax - L5Bmin) / 2 + + # set gids for each pop + popGids = {popLabel: [c['gid'] for c in allCells if c['tags']['pop'] == popLabel] for popLabel in popLabels if isinstance(popLabel, str)} + popGids['upperPT5B'] = tuple([c['gid'] for c in allCells \ + if L5Bmin <= c['tags']['ynorm'] <= L5Bmid and c['tags']['pop']=='PT5B']) + popGids['lowerPT5B'] = tuple([c['gid'] for c in allCells \ + if L5Bmid <= c['tags']['ynorm'] <= L5Bmax and c['tags']['pop']=='PT5B']) + popGids['L5B'] = popGids['IT5B'] + popGids['PT5B'] + + popNumCells = {popLabel: len(popGids[popLabel]) for popLabel in popGids} + + return popGids, popNumCells \ No newline at end of file diff --git a/analysis/utils.py b/analysis/utils.py new file mode 100644 index 0000000..af2f58a --- /dev/null +++ b/analysis/utils.py @@ -0,0 +1,563 @@ +""" +utils.py + +General functions to analyse simulation data + +Contributors: salvadordura@gmail.com +""" +import json +import pickle +import numpy as np +from pylab import * +from itertools import product +import pandas as pd +import seaborn as sb +from scipy import stats +from pprint import pprint +from collections import OrderedDict +try: + basestring +except NameError: + basestring = str + +from netpyne import specs + +def toPandas(params, data): + if 'simData' in data[data.keys()[0]]: + rows = [list(d['paramValues'])+[s for s in d['simData'].values()] for d in data.values()] + cols = [str(d['label']) for d in params]+[s for s in data[data.keys()[0]]['simData'].keys()] + else: + rows = [list(d['paramValues'])+[s for s in d.values()] for d in data.values()] + cols = [str(d['label']) for d in params]+[s for s in data[data.keys()[0]].keys()] + + df = pd.DataFrame(rows, columns=cols) + df['simLabel'] = data.keys() + + colRename=[] + for col in list(df.columns): + if col.startswith("[u'"): + colName = col.replace(", u'","_'").replace("[u","").replace("'","").replace("]","").replace(", ","_") + colRename.append(colName) + else: + colRename.append(col) + print(colRename) + df.columns = colRename + + return df + + +def setPlotFormat(numColors=8): + # plt.style.use('ggplot') + #plt.style.use(['dark_background', 'presentation']) + plt.style.use('seaborn-whitegrid') + + plt.rcParams['font.size'] = 12 + plt.rcParams['axes.titlesize'] = 14 + plt.rcParams['axes.labelsize'] = 12 + plt.rcParams['legend.fontsize'] = 'large' + + NUM_COLORS = numColors + colormap = plt.get_cmap('nipy_spectral') + colorlist = [colormap(1.*i/NUM_COLORS) for i in range(NUM_COLORS)] + + # add red and blue at the beginning + # colorlist.insert(0, [1, 0, 0]) + # colorlist.insert(0, [0, 0, 1]) + + plt.rc('axes', prop_cycle=(cycler('color', colorlist))) + + +def compare(source_file, target_file, source_key=None, target_key=None): + from deepdiff import DeepDiff + with open(source_file, 'r') as fileObj: + if source_file.endswith('.json'): + source = json.load(fileObj, object_pairs_hook=OrderedDict) + elif source_file.endswith('.pkl'): + source = pickle.load(fileObj) + if source_key: source = source[source_key] + + with open(target_file, 'r') as fileObj: + if target_file.endswith('.json'): + target = json.load(fileObj, object_pairs_hook=OrderedDict) + elif source_file.endswith('.pkl'): + target = pickle.load(fileObj) + if target_key: target = target[target_key] + + ddiff = DeepDiff(source, target) + pprint(ddiff) + return ddiff + + +def readBatchData(dataFolder, batchLabel, loadAll=False, saveAll=True, vars=None, maxCombs=None, listCombs=None): + # load from previously saved file with all data + if loadAll: + print('\nLoading single file with all data...') + filename = '%s/%s/%s_allData.json' % (dataFolder, batchLabel, batchLabel) + with open(filename, 'r') as fileObj: + dataLoad = json.load(fileObj, object_pairs_hook=OrderedDict) + params = dataLoad['params'] + data = dataLoad['data'] + return params, data + + if isinstance(listCombs, basestring): + filename = str(listCombs) + with open(filename, 'r') as fileObj: + dataLoad = json.load(fileObj) + listCombs = dataLoad['paramsMatch'] + + # read the batch file and cfg + batchFile = '%s/%s/%s_batch.json' % (dataFolder, batchLabel, batchLabel) + with open(batchFile, 'r') as fileObj: + b = json.load(fileObj)['batch'] + + # read params labels and ranges + params = b['params'] + + # reorder so grouped params come first + preorder = [p for p in params if 'group' in p and p['group']] + for p in params: + if p not in preorder: preorder.append(p) + params = preorder + + # read vars from all files - store in dict + if b['method'] == 'grid': + labelList, valuesList = zip(*[(p['label'], p['values']) for p in params]) + valueCombinations = product(*(valuesList)) + indexCombinations = product(*[range(len(x)) for x in valuesList]) + data = {} + print('Reading data...') + missing = 0 + for i,(iComb, pComb) in enumerate(zip(indexCombinations, valueCombinations)): + if (not maxCombs or i<= maxCombs) and (not listCombs or list(pComb) in listCombs): + print(i, iComb) + # read output file + iCombStr = ''.join([''.join('_'+str(i)) for i in iComb]) + simLabel = b['batchLabel']+iCombStr + outFile = b['saveFolder']+'/'+simLabel+'.json' + try: + with open(outFile, 'r') as fileObj: + output = json.load(fileObj, object_pairs_hook=OrderedDict) + # save output file in data dict + data[iCombStr] = {} + data[iCombStr]['paramValues'] = pComb # store param values + if not vars: vars = output.keys() + + for key in vars: + if isinstance(key, tuple): + container = output + for ikey in range(len(key)-1): + container = container[key[ikey]] + data[iCombStr][key[1]] = container[key[-1]] + + elif isinstance(key, basestring): + data[iCombStr][key] = output[key] + + except: + print('... file missing') + missing = missing + 1 + output = {} + else: + missing = missing + 1 + + print('%d files missing' % (missing)) + + # save + if saveAll: + print('Saving to single file with all data') + filename = '%s/%s/%s_allData.json' % (dataFolder, batchLabel, batchLabel) + dataSave = {'params': params, 'data': data} + with open(filename, 'w') as fileObj: + json.dump(dataSave, fileObj) + + return params, data + + + +def plotsFromFile(filename, raster=True, stats=True, rates=False, syncs=False, hist=True, psd=True, grang=True, traces=True, plotAll=False, + timeRange=None, textTop='', include={'raster':[]}, popColors=None, orderBy='gid'): + from netpyne import specs + from collections import OrderedDict + + # try: + if filename.endswith('.json'): + with open(filename, 'rb') as fileObj: + data = json.load(fileObj, object_pairs_hook=OrderedDict) + elif filename.endswith('.pkl'): + with open(filename, 'rb') as fileObj: + data = pickle.load(fileObj) + # except: + # print('Error opening file %s' % (filename)) + # return 0 + + sim,data, out = plotsFromData(data, raster=raster, stats=stats, rates=rates, syncs=syncs, hist=hist, psd=psd, traces=traces, grang=grang, + plotAll=plotAll, timeRange=timeRange, textTop=textTop, include=include, popColors=popColors, orderBy=orderBy) + return sim,data, out + + +def plotsFromData(data, textTop = '', raster=True, stats=False, rates=False, syncs=False, hist=True, psd=True, traces=True, grang=True, + plotAll=0, timeRange=None, include=None, popColors=None, orderBy='gid'): + import matplotlib.pyplot as plt + from netpyne import specs,sim + if 'simConfig' in data: + cfg = specs.SimConfig(data['simConfig']) + else: + cfg = specs.SimConfig() + cfg.createNEURONObj = False + + sim.initialize() # create network object and set cfg and net params + sim.loadAll('', data=data, instantiate=False) + sim.setSimCfg(cfg) + + import collections + temp = collections.OrderedDict() + if len(include['raster']) == 0: + include['raster'] = popColors.keys() + + # order pops by keys in popColors + for k in include['raster']: + if k in sim.net.params.popParams: + temp[k] = sim.net.params.popParams.pop(k) + + # add remaining pops at the end + for k in list(sim.net.params.popParams.keys()): + temp[k] = sim.net.params.popParams.pop(k) + + sim.net.params.popParams = temp + + try: + print('Cells created: ',len(sim.net.allCells)) + except: + sim.net.createPops() + sim.net.createCells() + sim.setupRecording() + sim.cfg.createNEURONObj = False + sim.gatherData() + + sim.allSimData = data['simData'] + + plt.style.use('seaborn-ticks') + + out = None + + if raster or plotAll: + fig1 = sim.analysis.plotRaster(include=include['raster'], timeRange=timeRange, labels='overlay', popRates=False, orderInverse=True, lw=0, markerSize=5, + marker='.', popColors=popColors, showFig=0, saveFig=0, figSize=(10,8), orderBy=orderBy) + ax = plt.gca() + plt.text(0.05, 1.07, textTop, transform=ax.transAxes, fontsize=12, color='k') + plt.ylabel('Neuron number') + plt.title('') + plt.savefig(cfg.filename+'_raster_%d_%d_%s.png'%(timeRange[0], timeRange[1], orderBy), dpi=600) + + if stats or plotAll: + filename = cfg.filename+'_%d_%d'%(timeRange[0], timeRange[1]) + fig1 = sim.analysis.plotSpikeStats(include=include['stats'], figSize=(4,8), timeRange=timeRange, stats = ['rate', 'isicv'], popColors=popColors, showFig=0, saveFig=filename) + + if rates or plotAll: + midpoint = (timeRange[1]+timeRange[0]) / 2.0 + timeRanges = [[timeRange[0], midpoint], [midpoint, timeRange[1]]] + colors = [popColors[inc] if inc in popColors else (0,0,0) for inc in include['rates']] + out = sim.analysis.plotRates(include=include['rates'], timeRanges=timeRanges, figSize=(4,2), timeRangeLabels=['Pre stim', 'Post stim'], colors=colors, showFig=0, saveFig=1) + + if syncs or plotAll: + midpoint = (timeRange[1]+timeRange[0]) / 2.0 + timeRanges = [[timeRange[0], midpoint], [midpoint, timeRange[1]]] + colors = [popColors[inc] if inc in popColors else (0,0,0) for inc in include['syncs']] + fig1 = sim.analysis.plotSyncs(include=include['syncs'], timeRanges=timeRanges, figSize=(5,4), timeRangeLabels=['Pre stim', 'Post stim'], colors=colors, showFig=0, saveFig=1) + + if hist or plotAll: + fig2 = sim.analysis.plotSpikeHist(include=include['hist'], yaxis='rate', binSize=5, graphType='bar', timeRange=timeRange, popColors=popColors, showFig=False, saveFig=0, figSize=(12,8)) + ax = plt.gca() + plt.text(0.05, 1.07, textTop, transform=ax.transAxes, fontsize=12, color='k') + plt.savefig(cfg.filename+'_spikeHist_replot_%d_%d.png'%(timeRange[0], timeRange[1])) + + if psd or plotAll: + popColors['allCells'] = 'k' + fig3 = sim.analysis.plotRatePSD(include=include['psd'], timeRange=timeRange, Fs=160, smooth=16 , showFig=0, saveFig=0, popColors=popColors, figSize=(11,5)) + # ylim=[-55, 0] + ax = plt.gca() + plt.text(0.05, 1.07, textTop, transform=ax.transAxes, fontsize=12, color='k') + plt.savefig(cfg.filename+'_spikePSD_%d_%d.png'%(timeRange[0], timeRange[1]), dpi=300) + + if traces or plotAll: + colors = popColors; colors=['r','b'] + fig4 = sim.analysis.plotTraces(include=include['traces'], timeRange=timeRange, overlay = True, oneFigPer = 'trace', rerun = False, ylim=[-90, 30], axis='on', figSize = (12,4), saveData = None, saveFig = 1, showFig = 0) + + return sim, data, out + + +def getSpksPop(data, timeRange=None): + from netpyne import specs,sim + cfg = specs.SimConfig(data['simConfig']) + cfg.createNEURONObj = False + + sim.initialize() # create network object and set cfg and net params + sim.loadAll('', data=data) + sim.setSimCfg(cfg) + sim.net.createPops() + sim.net.createCells() + sim.gatherData() + + sim.allSimData = data['simData'] + spkts = sim.allSimData['spkt'] + spkids = sim.allSimData['spkid'] + + popLabels = sim.net.allPops.keys() + gidPops = [cell['tags']['pop'] for cell in sim.net.allCells] + popNumCells = [float(gidPops.count(pop)) for pop in popLabels] + numCellsPop = {pop:float(gidPops.count(pop)) for pop in popLabels} + spktPop = {} + for pop, popNum in zip(popLabels, popNumCells): + if timeRange: + spktPop[pop] = [spkt for spkt,spkid in zip(spkts,spkids) if sim.net.allCells[int(spkid)]['tags']['pop']==pop if timeRange[0]<=spkt<=timeRange[1]] + else: + spktPop[pop] = [spkt for spkt,spkid in zip(spkts,spkids) if sim.net.allCells[int(spkid)]['tags']['pop']==pop] + + return spktPop, numCellsPop + + +def readBatchUpdatePopRates(dataFolder, batchLabel, timeRange=None, saveAll=True, loadAll=False): + # load from previously saved file with all data + if loadAll: + from netpyne import specs + print('\nLoading single file with all data...') + filename = '%s/%s/%s_allData_popRates.json' % (dataFolder, batchLabel, batchLabel) + with open(filename, 'r') as fileObj: + dataLoad = json.load(fileObj, object_pairs_hook=OrderedDict) + params = dataLoad['params'] + data = dataLoad['data'] + return params, data + + # load single sim to get netParams and simConfig + params, data = readBatchData(dataFolder, batchLabel, maxCombs=1) + + # recreate net to get cell+pop params + from netpyne import specs,sim + cfg = specs.SimConfig(data.values()[0]['simConfig']) + cfg.createNEURONObj = False + sim.initialize() # create network object and set cfg and net params + sim.loadAll('', data=data.values()[0]) + sim.setSimCfg(cfg) + sim.net.createPops() + sim.net.createCells() + sim.gatherData() + + # caclulate pop stats + if not timeRange: timeRange = [0, cfg.duration] + tsecs = (timeRange[1]-timeRange[0])/1e3 + popLabels = sim.net.allPops.keys() + gidPops = [cell['tags']['pop'] for cell in sim.net.allCells] + popNumCells = [float(gidPops.count(pop)) for pop in popLabels] + #numCellsPop = {pop:float(gidPops.count(pop)) for pop in popLabels} + popLabels = sim.net.allPops.keys() + + # load spkt and spkid + params, data = readBatchData(dataFolder, batchLabel, vars=[('simData','spkt'),('simData','spkid'),('simData','popRates')], saveAll=False) + for d in data.values(): + spkts = d['spkt'] + spkids = d['spkid'] + + spktPop = {} + for pop, popNum in zip(popLabels, popNumCells): + spktPop[pop] = [spkt for spkt,spkid in zip(spkts,spkids) if sim.net.allCells[int(spkid)]['tags']['pop']==pop if timeRange[0]<=spkt<=timeRange[1]] + d['popRates'][pop] = len(spktPop[pop]) / popNum / tsecs # calcualte new popRate + d['spkt'] = None # remove to save space + d['spkid'] = None + + # save + if saveAll: + print('Saving to single file with all data') + filename = '%s/%s/%s_allData_popRates.json' % (dataFolder, batchLabel, batchLabel) + dataSave = {'params': params, 'data': data} + with open(filename, 'w') as fileObj: + json.dump(dataSave, fileObj) + + return params, data + + + + + +# --------- +# Function to plot boxplot of model vs experiment +def stats_boxplot(dfStats, x, y, hue, color, overlayLine, quietExp, moveExp, quietModel, moveModel, figSize, fontsize, palette=None, order=None): + #my_pal = {'Model': color, 'Experiment': '#999999'} + if palette: + my_pal = palette + else: + my_pal = {'Model': 'royalblue', 'Experiment': 'orange'} + + fig=plt.figure(figsize=figSize) + ax = sb.boxplot(x=x, y=y, hue=hue, data=dfStats, palette=my_pal, showfliers=False, order=order)#, boxprops=dict(alpha=.75), linewidth=0.5, width=0.3) # + + #plt.title (title, fontsize=fontsiz) + ax.get_legend().remove() + plt.tight_layout() + plt.ylabel('Firing Rate (Hz)', fontsize=fontsize) + plt.xlabel('') + ax = plt.gca() + ax.spines['right'].set_visible(False) + ax.spines['top'].set_visible(False) + ax.spines['bottom'].set_visible(False) + ax.spines['left'].set_visible(False) + for tick in ax.xaxis.get_major_ticks(): + tick.label.set_fontsize(fontsize) + for tick in ax.yaxis.get_major_ticks(): + tick.label.set_fontsize(fontsize) + fig.subplots_adjust(left=0.13) + + n = len(my_pal) + width=0.8/n + + onlyStar = True + + if overlayLine: + for i,(title,color) in enumerate(my_pal.items()): + if 'Experiment' in title: + data = [np.mean(quietExp), np.mean(moveExp)] + elif isinstance(quietModel, dict): + data = [np.mean(quietModel[title]), np.mean(moveModel[title])] + else: + data = [np.mean(quietModel), np.mean(moveModel)] + + if onlyStar: + plt.scatter(-0.4+(i*width)+width/2, data[0], marker='*', s=80, color='white', edgecolors='black') + plt.scatter(1-0.4+(i*width)+width/2, data[1], marker='*', s=80, color='white', edgecolors='black') + else: + plt.plot([-0.4+(i*width)+width/2, 1-0.4+(i*width)+width/2], data, marker='*', markersize=10, linewidth=2, mec='black', mfc=color, color=color) + + # fixed axis scale + maxRate = 30 + plt.ylim(0, maxRate) + + + +# --------- +# Function to plot lineplot with errorbar of model vs experiment +# +# Note: individual model line plots are taken from statDataAllQuiet and statDataAllMove, which contain all data points +# ie. filtering of >0.01 Hz has not been applied to these raw variables, so may not correspond to quietExp and moveExp +def stats_lineplot(title, quietExp, moveExp, quietModel, moveModel, statDataAllQuiet, statDataAllMove, modelInd, figSize, fontsize, printOutput=False, linestyle='-', colors={}): + + fig=plt.figure(figsize=figSize) + + avgFunc = np.mean + varFunc = np.std + + capsize = 8 # set to 0 for legend + lw = 2 + elw=1.0 + marker = 'o' + ms = 15 # set to 10 for legend + mec = 'black' + labelExp = 'Experiment' + + + # plot model + if not isinstance(quietModel, dict) and not isinstance(moveModel, dict) : # convert to common dict format + quietModel = {'Model': quietModel} + moveModel = {'Model': moveModel} + colors = {'Model': 'blue'} + + for k in quietModel: + ax2 = plt.errorbar([0.2,0.8], [avgFunc(quietModel[k]), avgFunc(moveModel[k])], yerr=[varFunc(quietModel[k]), varFunc(moveModel[k])], + linewidth=lw, capsize=capsize, color=colors[k], marker=marker, markersize=ms, mec='black', elinewidth=elw,alpha=1.0) #mec='blue', mec='black')#, mfc='blue') + for cap in ax2[1]: + cap.set_markeredgewidth(elw) + + + # plot experiment (move to top for legend) + if linestyle != '-': # either show exp data in lower alpha; or not include exp data at all + pass + # ax1 = plt.errorbar([0.2,0.8], [avgFunc(quietExp), avgFunc(moveExp)], yerr=[varFunc(quietExp), varFunc(moveExp)], + # linewidth=lw, capsize=capsize, color='orange', alpha=0.5, marker=marker, markersize=ms, mec='black', linestyle=linestyle, elinewidth=elw) #mec='orange',)#, mfc='orange') + # ax1[-1][0].set_linestyle(linestyle) + + else: + ax1 = plt.errorbar([0.2,0.8], [avgFunc(quietExp), avgFunc(moveExp)], yerr=[varFunc(quietExp), varFunc(moveExp)], + linewidth=lw, capsize=capsize, color='orange', marker=marker, markersize=ms, mec='black', linestyle=linestyle, elinewidth=elw, alpha=1.0) #mec='orange',)#, mfc='orange') + + for cap in ax1[1]: + cap.set_markeredgewidth(elw) + + + ax = plt.gca() + ax.spines['right'].set_visible(False) + ax.spines['top'].set_visible(False) + + + trialsQuiet = [] + trialsMove = [] + + for statTrialQuiet, statTrialMove in zip(statDataAllQuiet, statDataAllMove): + + plt.plot([0.2,0.8], [avgFunc(statTrialQuiet[modelInd]), avgFunc(statTrialMove[modelInd])], + linewidth=0.25, color='royalblue', alpha=0.8)# capsize=capsize, elinewidth=0.25) #, capsize=capsize) + + trialsQuiet.append(avgFunc(statTrialQuiet[modelInd])) + trialsMove.append(avgFunc(statTrialMove[modelInd])) + + plt.xticks([0.2,0.8], ['Quiet', 'Movement'], fontsize=fontsize) + plt.yticks(fontsize=fontsize) + plt.yticks([0,5,10,15,20,25],[0,5,10,15,20,25]) + plt.xlim(0,1) + plt.ylim(0,25) + plt.ylabel('Firing Rate (Hz)', fontsize=fontsize) + plt.subplots_adjust(left=0.18, bottom=0.1) + + + if printOutput: + print(title) + + for k in quietModel: + if len(k) > 0: print(k) + print('MODEL: %s +- %s (across cells): quiet=%.1f+-%.1f , move=%.1f+-%.1f' + % (avgFunc.__name__, varFunc.__name__, avgFunc(quietModel[k]), varFunc(quietModel[k]), avgFunc(moveModel[k]), varFunc(moveModel[k]))) + + print('EXP: %s +- %s (across cells): quiet=%.1f+-%.1f , move=%.1f+-%.1f' + % (avgFunc.__name__, varFunc.__name__, avgFunc(quietExp), varFunc(quietExp), avgFunc(moveExp), varFunc(moveExp))) + + + +# ------------- +# statistical tests +def stats_tests(title, quietExp, moveExp, quietModel, moveModel, sameLength=True): + import scipy + + def pvalueStars(prob): + if prob <= 0.001: + probstar = '***' + elif prob <= 0.01: + probstar = '**' + elif prob <= 0.05: + probstar = '*' + else: + probstar = 'ns' + return probstar + + + if sameLength: + if len(quietModel) > len(moveModel): + quietModel = quietModel[:len(moveModel)] + else: + moveModel = moveModel[:len(quietModel)] + + # Mann-Whitney U + dquiet_mw, probquiet_mw = scipy.stats.mannwhitneyu(quietExp, quietModel) + dmove_mw, probmove_mw = scipy.stats.mannwhitneyu(moveExp, moveModel) + probquietstar_mw = pvalueStars(probquiet_mw) + probmovestar_mw = pvalueStars(probmove_mw) + + # confidence interval / significance level, alpha; typically 0.05 + alpha = 0.05 + c_alpha = {0.2: 1.07, 0.15: 1.14, 0.1: 1.22, 0.05: 1.36, 0.025: 1.48, 0.01: 1.63, 0.005: 1.73} + + dstarquiet = c_alpha[alpha] * np.sqrt((len(quietExp)+len(quietModel) / len(quietExp)*len(quietModel))) + dstarmove = c_alpha[alpha] * np.sqrt((len(moveExp)+len(moveModel) / len(moveExp)*len(moveModel))) + + # null hypothesis (same distribution) rejected if alpha >= prob + # i.e. same distribution NOT rejected if alpha < prob + + print('\nMann-Whitney U:') + print('%s QUIET, N_exp=%d, N_model=%d, D=%f, D*=%f, alpha=%f, prob=%g %s' % (title, len(quietExp), len(quietModel), dquiet_mw, dstarquiet, alpha, probquiet_mw, probquietstar_mw)) + print('%s MOVE, N_exp=%d, N_model=%d, D=%f, D*=%f, alpha=%f, prob=%g %s' % (title, len(moveExp), len(moveModel), dmove_mw, dstarmove, alpha, probmove_mw, probmovestar_mw)) + print('')