import numpy
import pickle
import matplotlib
import matplotlib.pyplot as plt
import pprint
from adjustText import adjust_text
import matplotlib._pylab_helpers
import matplotlib.patheffects as path_effects
import math
import scipy.stats as st
import statsmodels.stats.multitest as mt
import pandas as pd
import glob
import load_swc
import neuron
import seaborn as sns
from neuron import h, gui

filedir = 'HL5PN2.swc'
cell = load_swc.main(filedir)
morph = load_swc.load_swc(filedir,cell)

pp = pprint.PrettyPrinter(indent=4)

font = {'family' : 'normal',
		'weight' : 'normal',
		'size'   : 20}

matplotlib.rc('font', **font)

vars = ['g_pas', 'e_pas', 'gbar_NaTgsomatic', 'gbar_Napsomatic', 'gbar_K_Psomatic', 'gbar_K_Tsomatic', 'gbar_Kv3_1somatic', 'gbar_Imsomatic', 'gbar_SKsomatic', 'decay_CaDynamicssomatic', 'gbar_Ca_HVAsomatic', 'gbar_Ca_LVAsomatic', 'gbar_Ihall', 'gbar_NaTgaxonal', 'gbar_Napaxonal', 'gbar_K_Paxonal', 'gbar_K_Taxonal', 'gbar_Kv3_1axonal', 'gbar_Imaxonal', 'gbar_SKaxonal', 'decay_CaDynamicsaxonal', 'gbar_Ca_HVAaxonal', 'gbar_Ca_LVAaxonal']
labels = [r'$G_{pas}$',r'$E_{pas}$',r'$G_{NaT,s}$',r'$G_{NaP,s}$',r'$G_{K_P,s}$',r'$G_{K_T,s}$',r'$G_{Kv3.1,s}$',r'$G_{M,s}$',r'$G_{SK,s}$',r'$Ca_{Dyn._{\tau,s}}$',r'$G_{Ca_{HVA},s}$',r'$G_{Ca_{LVA},s}$',r'$G_{H}$',r'$G_{NaT,a}$',r'$G_{NaP,a}$',r'$G_{K_P,a}$',r'$G_{K_T,a}$',r'$G_{Kv3.1,a}$',r'$G_{M,a}$',r'$G_{SK,a}$',r'$Ca_{Dyn._{\tau,a}}$',r'$G_{Ca_{HVA},a}$',r'$G_{Ca_{LVA},a}$']
lowerlims = [0.00000001,-100,0,0,0,0,0,0,0,20,0,0,0,0,0,0,0,0,0,0,20,0,0]
upperlims = [0.0002,-72,1,1e-2,1.5,1.5,1.5,5e-4,1.5,1000,1e-4,1e-2,0.0001,1,1e-2,1.5,1.5,1.5,5e-4,1.5,1000,1e-4,1e-2]

dirpath1 = 'young/'
dirpath2 = 'old/'

finalpop_pathy = dirpath1+'finalpop.pkl'
finalpop_patho = dirpath2+'finalpop.pkl'

halloffame_pathy = dirpath1+'halloffame.pkl'
halloffame_patho = dirpath2+'halloffame.pkl'

hist_pathy = dirpath1+'hist.pkl'
hist_patho = dirpath2+'hist.pkl'

logs_pathy = dirpath1+'logs.pkl'
logs_patho = dirpath2+'logs.pkl'

objectives_pathy = dirpath1+'objectives.pkl'
objectives_patho = dirpath2+'objectives.pkl'

finalpopy = pickle.load(open(finalpop_pathy,"rb"))
finalpopo = pickle.load(open(finalpop_patho,"rb"))

halloffamey = pickle.load(open(halloffame_pathy,"rb"))
halloffameo = pickle.load(open(halloffame_patho,"rb"))

histy = pickle.load(open(hist_pathy,"rb"))
histo = pickle.load(open(hist_patho,"rb"))

logsy = pickle.load(open(logs_pathy,"rb"))
logso = pickle.load(open(logs_patho,"rb"))

objectivesy = pickle.load(open(objectives_pathy,"rb"))
objectiveso = pickle.load(open(objectives_patho,"rb"))

# Note: with new scinet code, the objectives are all out of order and not necessarily the same between runs, so need to reconstruct
objectivesy = [obname.name for obname in objectivesy]
objectiveso = [obname.name for obname in objectiveso]

maxBaseSD = 2 # 2 works for young but not all the old models features
maxRMPSD = 2 # same as base SD
maxSagSD = 1 # i.e. higher constraint for sag since that is the focus
voltage_base_Compensator = 2*maxRMPSD # Compensates for having divided the target SD by 2 during Optimization to get better fitness
sag_amplitude_Compensator = 10*maxSagSD # Compensates for having divided the target SD by 10 during Optimization to get better fitness
AP_height_Compensator = [2,6] # since fits are not as good in old
AHP_depth_abs_Compensator = 3 # since fits are not as good in old
AHP_depth_abs_slow_Compensator = 3 # since fits are not as good in old
AHP_slow_time_Compensator = 3 # since fits are not as good in old
AP_width_Compensator = 3 # since fits are not as good in old
Spikecount_Compensator = 2 # since fits are not as good in old
mean_frequency_Compensator = [0.5,0.8] # since fits are not as good in old

def weightvec(objectives,ag):
	weights = numpy.ones(len(objectives))
	for w in range(0,len(weights)):
		if objectives[w] == 'voltage_base':
			weights[w] = weights[w]*voltage_base_Compensator
		elif objectives[w] == 'sag_amplitude':
			weights[w] = weights[w]*sag_amplitude_Compensator
		elif objectives[w] == 'AP_height':
			weights[w] = weights[w]*AP_height_Compensator[0] if ag == 'y' else weights[w]*AP_height_Compensator[1]
		elif objectives[w] == 'AHP_depth_abs':
			weights[w] = weights[w]*AHP_depth_abs_Compensator
		elif objectives[w] == 'AHP_depth_abs_slow':
			weights[w] = weights[w]*AHP_depth_abs_slow_Compensator
		elif objectives[w] == 'AHP_slow_time':
			weights[w] = weights[w]*AHP_slow_time_Compensator
		elif objectives[w] == 'AP_width':
			weights[w] = weights[w]*AP_width_Compensator
		elif objectives[w] == 'Spikecount':
			weights[w] = weights[w]*Spikecount_Compensator
		elif objectives[w] == 'mean_frequency':
			weights[w] = weights[w]*mean_frequency_Compensator[0] if ag == 'y' else weights[w]*mean_frequency_Compensator[1]
		else:
			weights[w] = weights[w]*maxBaseSD
	
	return weights

weightsy = weightvec(objectivesy,'y')
weightso = weightvec(objectiveso,'o')

print(set(zip(objectivesy,weightsy)))
print(set(zip(objectiveso,weightso)))

bestIndexY = 0
bestIndexO = 0#2

# print params of chosen model
chosen = 0
for i2 in range(0,len(halloffameo[0])):
	print(vars[i2]+': '+ str(halloffameo[chosen][i2]))


totalpopy = []
totalpopo = []
qualityval_y = []
qualityval_o = []

qualitysag_y = []
qualitysag_o = []
qualityrmp_y = []
qualityrmp_o = []
sagidx_y = objectivesy.index("sag_amplitude")
rmpidx_y = objectivesy.index("voltage_base")
sagidx_o = objectiveso.index("sag_amplitude")
rmpidx_o = objectiveso.index("voltage_base")
reject_counts_y = [0 for _ in range(0,len(weightsy))]
reject_counts_o = [0 for _ in range(0,len(weightso))]
# iterate over all models across all generations
print(histy.__dict__['genealogy_history'][1].fitness)
for l in range(1,len(histy.__dict__['genealogy_history'])+1):
	f = histy.__dict__['genealogy_history'][l]
	count = 0
	for t in range(0,len(weightsy)):
		if f.fitness.values[t] < weightsy[t]:
			count += 1
		else:
			reject_counts_y[t] += 1
	
	if count == len(weightsy):
		qualityval_y.append(numpy.sum(f.fitness.values))
		qualitysag_y.append(f.fitness.values[sagidx_y]/10)
		qualityrmp_y.append(f.fitness.values[rmpidx_y]/2)
		totalpopy.append(f)
	else:
		continue

for l in range(1,len(histo.__dict__['genealogy_history'])+1):
	f = histo.__dict__['genealogy_history'][l]
	count = 0
	for t in range(0,len(weightso)):
		if f.fitness.values[t] < weightso[t]:
			count += 1
		else:
			reject_counts_o[t] += 1
	
	if count == len(weightso):
		qualityval_o.append(numpy.sum(f.fitness.values))
		qualitysag_o.append(f.fitness.values[sagidx_o]/10)
		qualityrmp_o.append(f.fitness.values[rmpidx_o]/2)
		totalpopo.append(f)
	else:
		continue

# Remove duplicate models
import itertools
totalpopy.sort()
totalpopy = list(totalpopy for totalpopy,_ in itertools.groupby(totalpopy))
totalpopo.sort()
totalpopo = list(totalpopo for totalpopo,_ in itertools.groupby(totalpopo))

print('Young % rejected:')
print(pp.pformat(list(zip(objectivesy,100*numpy.array(reject_counts_y)/len(histy.__dict__['genealogy_history'])))))
print('Old % rejected:')
print(pp.pformat(list(zip(objectiveso,100*numpy.array(reject_counts_o)/len(histo.__dict__['genealogy_history'])))))

print(len(qualityval_y))
print(len(qualityval_o))

summed_score_y = []
summed_score_o = []
for c, bi in enumerate(totalpopy):
	summed_score_y.append(bi.fitness.values)
for c, bi in enumerate(totalpopo):
	summed_score_o.append(bi.fitness.values)

sorted_totalpopy = [i for _,i in sorted(zip(summed_score_y,totalpopy))]
sorted_totalpopo = [i for _,i in sorted(zip(summed_score_o,totalpopo))]

num_models_to_generate = 100

for c, bi in enumerate(sorted_totalpopy):
	if c > num_models_to_generate-1: break
	txtf = open('highly_ranked_models/biophys_HL5PN1y' + str(c+1) + '.hoc', "w")
	txtf.write('proc biophys_HL5PN1(){\n')
	txtf.write('	forsec $o1.all {\n')
	txtf.write('		insert pas\n')
	txtf.write('		insert Ih\n')
	txtf.write('		Ra = 100\n')
	txtf.write('		cm = 0.9\n')
	txtf.write('		e_pas = ' + str(bi[1]) + '\n')
	txtf.write('		g_pas = ' + str(bi[0]) + '\n')
	txtf.write('		gbar_Ih = ' + str(bi[12]) + '\n')
	txtf.write('		shift1_Ih = 144.76545935424588\n')
	txtf.write('		shift2_Ih = 14.382865335237211\n')
	txtf.write('		shift3_Ih = -28.179477866349245\n')
	txtf.write('		shift4_Ih = 99.18311385307702\n')
	txtf.write('		shift5_Ih = 16.42000098505615\n')
	txtf.write('		shift6_Ih = 26.699880497099517\n')
	txtf.write('	}\n')
	txtf.write('	$o1.distribute_channels("apic","gbar_Ih",5,0.5,24,950,-285,$o1.soma.gbar_Ih)\n')
	txtf.write('	$o1.distribute_channels("dend","gbar_Ih",5,0.5,24,950,-285,$o1.soma.gbar_Ih)\n')
	txtf.write('	$o1.distribute_channels("axon","gbar_Ih",5,0.5,24,950,-285,$o1.soma.gbar_Ih)\n')
	txtf.write('	\n')
	txtf.write('	forsec $o1.somatic {\n')
	txtf.write('		insert NaTg\n')
	txtf.write('		insert Nap\n')
	txtf.write('		insert K_P\n')
	txtf.write('		insert K_T\n')
	txtf.write('		insert Kv3_1\n')
	txtf.write('		insert SK\n')
	txtf.write('		insert Im\n')
	txtf.write('		insert Ca_HVA\n')
	txtf.write('		insert Ca_LVA\n')
	txtf.write('		insert CaDynamics\n')
	txtf.write('		ek = -85\n')
	txtf.write('		ena = 50\n')
	txtf.write('		gbar_NaTg = ' + str(bi[2]) + '\n')
	txtf.write('		vshiftm_NaTg = 0\n')
	txtf.write('		vshifth_NaTg = 10\n')
	txtf.write('		slopem_NaTg = 9\n')
	txtf.write('		slopeh_NaTg = 6\n')
	txtf.write('		gbar_Nap = ' + str(bi[3]) + '\n')
	txtf.write('		gbar_K_P = ' + str(bi[4]) + '\n')
	txtf.write('		gbar_K_T = ' + str(bi[5]) + '\n')
	txtf.write('		gbar_Kv3_1 = ' + str(bi[6]) + '\n')
	txtf.write('		vshift_Kv3_1 = 0\n')
	txtf.write('		gbar_SK = ' + str(bi[8]) + '\n')
	txtf.write('		gbar_Im = ' + str(bi[7]) + '\n')
	txtf.write('		gbar_Ca_HVA = ' + str(bi[10]) + '\n')
	txtf.write('		gbar_Ca_LVA = ' + str(bi[11]) + '\n')
	txtf.write('		gamma_CaDynamics = 0.0005\n')
	txtf.write('		decay_CaDynamics = ' + str(bi[9]) + '\n')
	txtf.write('	}\n')
	txtf.write('	forsec $o1.axonal {\n')
	txtf.write('		insert NaTg\n')
	txtf.write('		insert Nap\n')
	txtf.write('		insert K_P\n')
	txtf.write('		insert K_T\n')
	txtf.write('		insert Kv3_1\n')
	txtf.write('		insert SK\n')
	txtf.write('		insert Im\n')
	txtf.write('		insert Ca_HVA\n')
	txtf.write('		insert Ca_LVA\n')
	txtf.write('		insert CaDynamics\n')
	txtf.write('		ek = -85\n')
	txtf.write('		ena = 50\n')
	txtf.write('		gbar_NaTg = ' + str(bi[13]) + '\n')
	txtf.write('		vshiftm_NaTg = 0\n')
	txtf.write('		vshifth_NaTg = 10\n')
	txtf.write('		slopem_NaTg = 9\n')
	txtf.write('		slopeh_NaTg = 6\n')
	txtf.write('		gbar_Nap = ' + str(bi[14]) + '\n')
	txtf.write('		gbar_K_P = ' + str(bi[15]) + '\n')
	txtf.write('		gbar_K_T = ' + str(bi[16]) + '\n')
	txtf.write('		gbar_Kv3_1 = ' + str(bi[17]) + '\n')
	txtf.write('		vshift_Kv3_1 = 0\n')
	txtf.write('		gbar_SK = ' + str(bi[19]) + '\n')
	txtf.write('		gbar_Im = ' + str(bi[18]) + '\n')
	txtf.write('		gbar_Ca_HVA = ' + str(bi[21]) + '\n')
	txtf.write('		gbar_Ca_LVA = ' + str(bi[22]) + '\n')
	txtf.write('		gamma_CaDynamics = 0.0005\n')
	txtf.write('		decay_CaDynamics = ' + str(bi[20]) + '\n')
	txtf.write('	}\n')
	txtf.write('}\n')
	txtf.write('\n')

for c, bi in enumerate(sorted_totalpopo):
	if c > num_models_to_generate-1: break
	txtf = open('highly_ranked_models/biophys_HL5PN1o' + str(c+1) + '.hoc', "w")
	txtf.write('proc biophys_HL5PN1(){\n')
	txtf.write('	forsec $o1.all {\n')
	txtf.write('		insert pas\n')
	txtf.write('		insert Ih\n')
	txtf.write('		Ra = 100\n')
	txtf.write('		cm = 0.9\n')
	txtf.write('		e_pas = ' + str(bi[1]) + '\n')
	txtf.write('		g_pas = ' + str(bi[0]) + '\n')
	txtf.write('		gbar_Ih = ' + str(bi[12]) + '\n')
	txtf.write('		shift1_Ih = 144.76545935424588\n')
	txtf.write('		shift2_Ih = 14.382865335237211\n')
	txtf.write('		shift3_Ih = -28.179477866349245\n')
	txtf.write('		shift4_Ih = 99.18311385307702\n')
	txtf.write('		shift5_Ih = 16.42000098505615\n')
	txtf.write('		shift6_Ih = 26.699880497099517\n')
	txtf.write('	}\n')
	txtf.write('	$o1.distribute_channels("apic","gbar_Ih",5,0.5,24,950,-285,$o1.soma.gbar_Ih)\n')
	txtf.write('	$o1.distribute_channels("dend","gbar_Ih",5,0.5,24,950,-285,$o1.soma.gbar_Ih)\n')
	txtf.write('	$o1.distribute_channels("axon","gbar_Ih",5,0.5,24,950,-285,$o1.soma.gbar_Ih)\n')
	txtf.write('	\n')
	txtf.write('	forsec $o1.somatic {\n')
	txtf.write('		insert NaTg\n')
	txtf.write('		insert Nap\n')
	txtf.write('		insert K_P\n')
	txtf.write('		insert K_T\n')
	txtf.write('		insert Kv3_1\n')
	txtf.write('		insert SK\n')
	txtf.write('		insert Im\n')
	txtf.write('		insert Ca_HVA\n')
	txtf.write('		insert Ca_LVA\n')
	txtf.write('		insert CaDynamics\n')
	txtf.write('		ek = -85\n')
	txtf.write('		ena = 50\n')
	txtf.write('		gbar_NaTg = ' + str(bi[2]) + '\n')
	txtf.write('		vshiftm_NaTg = 0\n')
	txtf.write('		vshifth_NaTg = 10\n')
	txtf.write('		slopem_NaTg = 9\n')
	txtf.write('		slopeh_NaTg = 6\n')
	txtf.write('		gbar_Nap = ' + str(bi[3]) + '\n')
	txtf.write('		gbar_K_P = ' + str(bi[4]) + '\n')
	txtf.write('		gbar_K_T = ' + str(bi[5]) + '\n')
	txtf.write('		gbar_Kv3_1 = ' + str(bi[6]) + '\n')
	txtf.write('		vshift_Kv3_1 = 0\n')
	txtf.write('		gbar_SK = ' + str(bi[8]) + '\n')
	txtf.write('		gbar_Im = ' + str(bi[7]) + '\n')
	txtf.write('		gbar_Ca_HVA = ' + str(bi[10]) + '\n')
	txtf.write('		gbar_Ca_LVA = ' + str(bi[11]) + '\n')
	txtf.write('		gamma_CaDynamics = 0.0005\n')
	txtf.write('		decay_CaDynamics = ' + str(bi[9]) + '\n')
	txtf.write('	}\n')
	txtf.write('	forsec $o1.axonal {\n')
	txtf.write('		insert NaTg\n')
	txtf.write('		insert Nap\n')
	txtf.write('		insert K_P\n')
	txtf.write('		insert K_T\n')
	txtf.write('		insert Kv3_1\n')
	txtf.write('		insert SK\n')
	txtf.write('		insert Im\n')
	txtf.write('		insert Ca_HVA\n')
	txtf.write('		insert Ca_LVA\n')
	txtf.write('		insert CaDynamics\n')
	txtf.write('		ek = -85\n')
	txtf.write('		ena = 50\n')
	txtf.write('		gbar_NaTg = ' + str(bi[13]) + '\n')
	txtf.write('		vshiftm_NaTg = 0\n')
	txtf.write('		vshifth_NaTg = 10\n')
	txtf.write('		slopem_NaTg = 9\n')
	txtf.write('		slopeh_NaTg = 6\n')
	txtf.write('		gbar_Nap = ' + str(bi[14]) + '\n')
	txtf.write('		gbar_K_P = ' + str(bi[15]) + '\n')
	txtf.write('		gbar_K_T = ' + str(bi[16]) + '\n')
	txtf.write('		gbar_Kv3_1 = ' + str(bi[17]) + '\n')
	txtf.write('		vshift_Kv3_1 = 0\n')
	txtf.write('		gbar_SK = ' + str(bi[19]) + '\n')
	txtf.write('		gbar_Im = ' + str(bi[18]) + '\n')
	txtf.write('		gbar_Ca_HVA = ' + str(bi[21]) + '\n')
	txtf.write('		gbar_Ca_LVA = ' + str(bi[22]) + '\n')
	txtf.write('		gamma_CaDynamics = 0.0005\n')
	txtf.write('		decay_CaDynamics = ' + str(bi[20]) + '\n')
	txtf.write('	}\n')
	txtf.write('}\n')
	txtf.write('\n')

qualityval_y = numpy.array(qualityval_y)
qualityval_o = numpy.array(qualityval_o)
fig, ax = plt.subplots(figsize=(8, 8))
ax.hist(qualityval_y,50,facecolor='k', alpha=0.5,label='Young')
ax.hist(qualityval_o,50,facecolor='darkred', alpha=0.5,label='Old')
bottom, top = ax.get_ylim()
ax.set_ylim((bottom,top))
ax.set_xlabel('Sum of Standard Deviations')
ax.set_ylabel('Model Count')
fig.tight_layout()
fig.savefig('PLOTfiles/Quality.pdf', bbox_inches='tight')
fig.savefig('PLOTfiles/Quality.png', bbox_inches='tight')
plt.close(fig)

print(halloffameo.__dict__['items'][bestIndexO].__dict__['fitness'].__dict__['wvalues'])
bestysag = abs(halloffamey.__dict__['items'][bestIndexY].__dict__['fitness'].__dict__['wvalues'][sagidx_y]/10)
bestyrmp = abs(halloffamey.__dict__['items'][bestIndexY].__dict__['fitness'].__dict__['wvalues'][rmpidx_y]/2)
bestosag = abs(halloffameo.__dict__['items'][bestIndexO].__dict__['fitness'].__dict__['wvalues'][sagidx_o]/10)
bestormp = abs(halloffameo.__dict__['items'][bestIndexO].__dict__['fitness'].__dict__['wvalues'][rmpidx_o]/2)

qualitysag_y = numpy.array(qualitysag_y)
qualitysag_o = numpy.array(qualitysag_o)
qualityrmp_y = numpy.array(qualityrmp_y)
qualityrmp_o = numpy.array(qualityrmp_o)
(y_r, y_p) = st.pearsonr(qualityrmp_y, qualitysag_y)
(o_r, o_p) = st.pearsonr(qualityrmp_o, qualitysag_o)
yfit = numpy.polyfit(qualityrmp_y, qualitysag_y, 1)
ofit = numpy.polyfit(qualityrmp_o, qualitysag_o, 1)
statsy = 'R = ' + str(numpy.around(y_r,2)) + '; p = ' + str(numpy.around(y_p,4))
statso = 'R = ' + str(numpy.around(o_r,2)) + '; p = ' + str(numpy.around(o_p,4))

font = {'family' : 'normal',
		'weight' : 'normal',
		'size'   : 20}

matplotlib.rc('font', **font)

fig, ax = plt.subplots(figsize=(7, 6))
ax.scatter(qualityrmp_y,qualitysag_y,color='k', alpha=0.7,label=statsy)
ax.scatter(qualityrmp_o,qualitysag_o,color='darkred', alpha=0.7,label=statso)
ax.scatter(bestyrmp,bestysag,s=500,linewidths=2,edgecolors='k',color='darkgray', alpha=1)
ax.scatter(bestormp,bestosag,s=500,linewidths=2,edgecolors='darkred',color='mistyrose', alpha=1)
for i2 in range(0,len(halloffamey.__dict__['items'])):
	ysag = abs(halloffamey.__dict__['items'][i2].__dict__['fitness'].__dict__['wvalues'][sagidx_y]/10)
	yrmp = abs(halloffamey.__dict__['items'][i2].__dict__['fitness'].__dict__['wvalues'][rmpidx_y]/2)
	osag = abs(halloffameo.__dict__['items'][i2].__dict__['fitness'].__dict__['wvalues'][sagidx_o]/10)
	ormp = abs(halloffameo.__dict__['items'][i2].__dict__['fitness'].__dict__['wvalues'][rmpidx_o]/2)
	text1 = ax.text(yrmp,ysag,str(i2+1),color='white',horizontalalignment='center',verticalalignment='center',fontsize=20)
	text2 = ax.text(ormp,osag,str(i2+1),color='mistyrose',horizontalalignment='center',verticalalignment='center',fontsize=20)
	text1.set_path_effects([path_effects.Stroke(linewidth=3, foreground='black'),
					   path_effects.Normal()])
	text2.set_path_effects([path_effects.Stroke(linewidth=3, foreground='darkred'),
					   path_effects.Normal()])
# plt.plot(qualityrmp_y,qualityrmp_y*yfit[0]+yfit[1],'k')
# plt.plot(qualityrmp_o,qualityrmp_o*ofit[0]+ofit[1],'darkred')
# leg = plt.legend(fontsize=10,loc='upper right', handlelength=0, handletextpad=0, fancybox=True)
bottom, top = ax.get_ylim()
ax.set_ylim((bottom,top))
ax.set_xlabel('Resting Potential Error (SD)')
ax.set_ylabel('Sag Amplitude Error (SD)')
fig.tight_layout()
fig.savefig('PLOTfiles/Quality_SagVsRMP.pdf', dpi=300, bbox_inches='tight',transparent=True)
fig.savefig('PLOTfiles/Quality_SagVsRMP.png', dpi=300, bbox_inches='tight',transparent=True)
plt.close(fig)

font = {'family' : 'normal',
		'weight' : 'normal',
		'size'   : 20}

matplotlib.rc('font', **font)

print(len(totalpopy))
print(len(totalpopo))

finalpopy = numpy.transpose(totalpopy)
finalpopo = numpy.transpose(totalpopo)

halloffamey = numpy.transpose(halloffamey)
halloffameo = numpy.transpose(halloffameo)

print('GHy = ' + str(halloffamey[12][bestIndexY]))
print('GHo = ' + str(halloffameo[12][bestIndexO]))

dists = numpy.linspace(0,1700,1000)
GHy_apics = halloffamey[12][bestIndexY]*(0.5+(24/(1+numpy.exp((dists-950)/(-285)))))
GHo_apics = halloffameo[12][bestIndexO]*(0.5+(24/(1+numpy.exp((dists-950)/(-285)))))
print('Y/O GH Ratio at soma = ' + str(GHo_apics[0]/GHy_apics[0]))
print('Y/O GH Ratio at distal apical = ' + str(GHo_apics[-1]/GHy_apics[-1]))
fig, ax = plt.subplots(1, figsize=(8, 8), facecolor='white')
ax.fill_between(dists,GHo_apics,color='r',alpha=0.7)
ax.fill_between(dists,GHy_apics,color='k',alpha=0.9)
ax.set_ylabel(r'Model G$_H$ (S/cm$^2$)')
ax.set_xlabel(r'Distance from Soma ($\mu$m)')
ax.set_xlim(0,1700)
ax.set_ylim(0,0.00065)
plt.savefig('PLOTfiles/GHvsDist.pdf', bbox_inches='tight', dpi=300)
plt.savefig('PLOTfiles/GHvsDist.png', bbox_inches='tight', dpi=300)
plt.gcf().clear()
plt.cla()
plt.clf()
plt.close()
# halloffamey = numpy.concatenate((halloffamey,[GHy_apic]))
# halloffameo = numpy.concatenate((halloffameo,[GHo_apic]))
# vars.append('gbar_Ihapic')
# labels.append(r'$G_{H,apic_{'+str(dist2test)+'\mum}}$')

gen_numbersy = logsy.select('gen')
gen_numberso = logso.select('gen')
min_fitnessy = numpy.array(logsy.select('min'))
min_fitnesso = numpy.array(logso.select('min'))
max_fitnessy = logsy.select('max')
max_fitnesso = logso.select('max')
mean_fitnessy = numpy.array(logsy.select('avg'))
mean_fitnesso = numpy.array(logso.select('avg'))
std_fitnessy = numpy.array(logsy.select('std'))
std_fitnesso = numpy.array(logso.select('std'))

fig, ax = plt.subplots(1, figsize=(8, 8), facecolor='white')

stdminus = mean_fitnessy - std_fitnessy
stdplus = mean_fitnessy + std_fitnessy

ax.plot(
	gen_numbersy,
	mean_fitnessy,
	color='black',
	linewidth=2,
	label='Population Average')

ax.fill_between(
	gen_numbersy,
	min_fitnessy,
	max_fitnessy,
	color='black',
	alpha=0.4,
	linewidth=2,
	label=r'Population Standard Deviation')

stdminus = mean_fitnesso - std_fitnesso
stdplus = mean_fitnesso + std_fitnesso

ax.plot(
	gen_numberso,
	mean_fitnesso,
	color='r',
	linewidth=2,
	label='Population Average')

ax.fill_between(
	gen_numberso,
	min_fitnesso,
	max_fitnesso,
	color='r',
	alpha=0.4,
	linewidth=2,
	label=r'Population Standard Deviation')

# ax.plot([min(gen_numbersy) - 1, max(gen_numbersy) + 1],[numpy.max(qualityval_y),numpy.max(qualityval_y)],'k',ls='dashed')
# ax.plot([min(gen_numberso) - 1, max(gen_numberso) + 1],[numpy.max(qualityval_o),numpy.max(qualityval_o)],'r',ls='dashed')
ax.set_xlim(min(gen_numbersy) - 1, max(gen_numbersy) + 1)
ax.set_xlabel('Generation #')
ax.set_ylabel('# Standard Deviations')
ax.set_ylim([10, 100000])
ax.set_yscale('log')
plt.savefig('PLOTfiles/Performance.pdf', bbox_inches='tight', dpi=300)
plt.savefig('PLOTfiles/Performance.png', bbox_inches='tight', dpi=300)
plt.gcf().clear()
plt.cla()
plt.clf()
plt.close()


texts1 = []
texts2 = []
xs = []
ys = []
ratios = []

fig, ax = plt.subplots(figsize=(15,14))
ax.set_xlabel('Old Model',fontsize=26)
ax.set_ylabel('Young Model',fontsize=26)
line1 = ax.plot(numpy.array([1e-7,1e+3]),numpy.array([1e-7,1e+3]),ls='dashed',color='k')
for i in range(0,len(halloffameo)):
	xs.append(halloffameo[i][bestIndexO])
	ys.append(halloffamey[i][bestIndexY])
	ratios.append(halloffamey[i][bestIndexY]/halloffameo[i][bestIndexO])
	if halloffameo[i][bestIndexO]>halloffamey[i][bestIndexY]:
		texts1.append(plt.text(halloffameo[i][bestIndexO], halloffamey[i][bestIndexY], labels[i]))
	else:
		texts2.append(plt.text(halloffameo[i][bestIndexO], halloffamey[i][bestIndexY], labels[i]))

scatter1 = ax.loglog(xs, ys, 'o', markersize=6, color='k')
ax.loglog(halloffameo[12][bestIndexO],halloffamey[12][bestIndexY],'o',fillstyle='none', markersize=14, color='g')
plots=[manager.canvas.figure for manager in matplotlib._pylab_helpers.Gcf.get_all_fig_managers()]
ax.grid(True)
ax.set_xlim(1e-7,1e+3)
ax.set_ylim(1e-7,1e+3)

expand = (1.2, 1.2)
forces = 0.4
adjust_text(texts1, precision=0.0000000001, ax=ax, autoalign='xy',
	expand_text=expand, expand_points=expand, expand_objects=expand, expand_align=expand,
	force_text=forces, force_points=forces, force_objects=forces,
	add_objects=line1, lim=99999,
	only_move={'points':'xy', 'texts':'xy', 'objects':'x'},
	arrowprops=dict(arrowstyle="wedge", color='dimgray', lw=1))
expand = (1.2, 1.2)
forces = 0.4
adjust_text(texts2, precision=0.0000000001, ax=ax, autoalign='xy',
	expand_text=expand, expand_points=expand, expand_objects=expand, expand_align=expand,
	force_text=forces, force_points=forces, force_objects=forces,
	add_objects=line1, lim=99999,
	only_move={'points':'xy', 'texts':'xy', 'objects':'y'},
	arrowprops=dict(arrowstyle="wedge", color='dimgray', lw=1))

fig.savefig('PLOTfiles/loglog_params.pdf')
fig.savefig('PLOTfiles/loglog_params.png')
plt.close(fig)

half = int(len(labels)/2)
fig = plt.figure(figsize=(15,7))
ax = fig.add_subplot(2,1,1)
ax.bar(labels[:half+1],ratios[:half+1],facecolor='k',edgecolor='k',linewidth=3,alpha=0.6)
ax.bar(labels[:half+1],ratios[:half+1],facecolor="None",edgecolor='k',linewidth=3,alpha=1)
ax.plot(numpy.array([-1,half+1]),numpy.array([1,1]),ls='dashed',color='r')
ax.set_ylabel('Y/O')
ax.set_yscale('log')
ax.set_xlim(-0.6,half+1-0.4)
ax2 = fig.add_subplot(2,1,2)
barlist2=ax2.bar(labels[half+1:],ratios[half+1:],facecolor='k',edgecolor='k',linewidth=3,alpha=0.6)
ax2.bar(labels[half+1:],ratios[half+1:],facecolor="None",edgecolor='k',linewidth=3,alpha=1)
barlist2[0].set_facecolor('g')
ax2.plot(numpy.array([-1,half]),numpy.array([1,1]),ls='dashed',color='r')
ax2.set_ylabel('Y/O')
ax2.set_yscale('log')
ax2.set_xlim(-0.6,half-0.4)
fig.savefig('PLOTfiles/ratios_params.pdf', bbox_inches='tight')
fig.savefig('PLOTfiles/ratios_params.png', bbox_inches='tight')
plt.close(fig)


# Create Total Conductance List
vars_Gtot = ['g_pas', 'gbar_NaTg', 'gbar_Nap', 'gbar_K_P', 'gbar_K_T', 'gbar_Kv3_1', 'gbar_Im', 'gbar_SK', 'gbar_Ca_HVA', 'gbar_Ca_LVA', 'gbar_Ihall']
labels_Gtot = [r'$G_{pas}$',r'$G_{NaT}$',r'$G_{NaP}$',r'$G_{K_P}$',r'$G_{K_T}$',r'$G_{Kv3.1}$',r'$G_{M}$',r'$G_{SK}$',r'$G_{Ca_{HVA}}$',r'$G_{Ca_{LVA}}$',r'$G_{H}$']

lowerlims2 = [0.00000001,0,0,0,0,0,0,0,0,0,0]
upperlims2 = [0.0002,1,1e-2,1.5,1.5,1.5,5e-4,1.5,1e-4,1e-2,0.0001]
lowerlims_Gtot = []
upperlims_Gtot = []
lowerlims_Gsoma = []
upperlims_Gsoma = []
# h('access soma')
h('distance()')
for i in range(0,len(lowerlims2)):
	h('Gval1 = 0')
	h('Gval2 = 0')
	h.Gval1 = lowerlims2[i]
	h.Gval2 = upperlims2[i]
	h('totG1 = 0')
	h('totG2 = 0')
	if (vars_Gtot[i] == 'g_pas'):
		h('forall for (x,0) totG1 += (area(x)*Gval1)*(1e+6/1e+8)')
		h('forall for (x,0) totG2 += (area(x)*Gval2)*(1e+6/1e+8)')
	elif (vars_Gtot[i] == 'gbar_Ihall'):
		h('forall for (x,0) totG1 += (area(x)*Gval1)*(1e+6/1e+8)*(0.5+(24/(1+exp((distance(x)-950)/(-285)))))')
		h('forall for (x,0) totG2 += (area(x)*Gval2)*(1e+6/1e+8)*(0.5+(24/(1+exp((distance(x)-950)/(-285)))))')
	else:
		h('forsec "soma" for (x,0) totG1 += (area(x)*Gval1)*(1e+6/1e+8)')
		h('forsec "axon" for (x,0) totG1 += (area(x)*Gval1)*(1e+6/1e+8)')
		h('forsec "soma" for (x,0) totG2 += (area(x)*Gval2)*(1e+6/1e+8)')
		h('forsec "axon" for (x,0) totG2 += (area(x)*Gval2)*(1e+6/1e+8)')
	lowerlims_Gtot.append(h.totG1)
	upperlims_Gtot.append(h.totG2)
	lowerlims_Gsoma.append(h.Gval1)
	upperlims_Gsoma.append(h.Gval2)

print(lowerlims_Gtot)
print(upperlims_Gtot)

Gtotsy = []
Gtotso = []
Gsomay = []
Gsomao = []
count = 0
for i in range(0,len(finalpopy)-10):
	if ((vars[i] == 'e_pas') | (vars[i] == 'decay_CaDynamicssomatic') | (vars[i] == 'decay_CaDynamicsaxonal')):
		continue
	elif vars[i] == 'g_pas':
		Gsomay.append([])
		Gtotsy.append([])
		for Gval in finalpopy[i]:
			Gsomay[count].append(Gval)
			h('Gval = 0')
			h.Gval = Gval
			h('totG = 0')
			h('forall for (x,0) totG += (area(x)*Gval)*(1e+6/1e+8)')
			Gtotsy[count].append(h.totG)
		Gsomao.append([])
		Gtotso.append([])
		for Gval in finalpopo[i]:
			Gsomao[count].append(Gval)
			h('Gval = 0')
			h.Gval = Gval
			h('totG = 0')
			h('forall for (x,0) totG += (area(x)*Gval)*(1e+6/1e+8)')
			Gtotso[count].append(h.totG)
		count += 1
	elif vars[i] == 'gbar_Ihall':
		Gsomay.append([])
		Gtotsy.append([])
		for Gval in finalpopy[i]:
			Gsomay[count].append(Gval)
			h('Gval = 0')
			h.Gval = Gval
			h('totG = 0')
			h('forall for (x,0) totG += (area(x)*Gval)*(1e+6/1e+8)*(0.5+(24/(1+exp((distance(x)-950)/(-285)))))')
			Gtotsy[count].append(h.totG)
		Gsomao.append([])
		Gtotso.append([])
		for Gval in finalpopo[i]:
			Gsomao[count].append(Gval)
			h('Gval = 0')
			h.Gval = Gval
			h('totG = 0')
			h('forall for (x,0) totG += (area(x)*Gval)*(1e+6/1e+8)*(0.5+(24/(1+exp((distance(x)-950)/(-285)))))')
			Gtotso[count].append(h.totG)
		count += 1
	else:
		Gsomay.append([])
		Gtotsy.append([])
		for j in range(0,len(finalpopy[i])):
			Gsomay[count].append(finalpopy[i][j])
			Gval1 = finalpopy[i][j]
			Gval2 = finalpopy[i+11][j]
			h('Gval1 = 0')
			h('Gval2 = 0')
			h.Gval1 = Gval1
			h.Gval2 = Gval2
			h('totG = 0')
			h('forsec "soma" for (x,0) totG += (area(x)*Gval1)*(1e+6/1e+8)')
			h('forsec "axon" for (x,0) totG += (area(x)*Gval2)*(1e+6/1e+8)')
			Gtotsy[count].append(h.totG)
		Gsomao.append([])
		Gtotso.append([])
		for j in range(0,len(finalpopo[i])):
			Gsomao[count].append(finalpopo[i][j])
			Gval1 = finalpopo[i][j]
			Gval2 = finalpopo[i+11][j]
			h('Gval1 = 0')
			h('Gval2 = 0')
			h.Gval1 = Gval1
			h.Gval2 = Gval2
			h('totG = 0')
			h('forsec "soma" for (x,0) totG += (area(x)*Gval1)*(1e+6/1e+8)')
			h('forsec "axon" for (x,0) totG += (area(x)*Gval2)*(1e+6/1e+8)')
			Gtotso[count].append(h.totG)
		count += 1

Besty = []
Besto = []
Bestsomay = []
Bestsomao = []
for i in range(0,len(halloffamey)-10):
	if ((vars[i] == 'e_pas') | (vars[i] == 'decay_CaDynamicssomatic') | (vars[i] == 'decay_CaDynamicsaxonal')):
		continue
	elif vars[i] == 'g_pas':
		h('Gval = 0')
		h.Gval = halloffamey[i][bestIndexY]
		h('totG = 0')
		h('forall for (x,0) totG += (area(x)*Gval)*(1e+6/1e+8)')
		Bestsomay.append(h.Gval)
		Besty.append(h.totG)
		h('Gval = 0')
		h.Gval = halloffameo[i][bestIndexO]
		h('totG = 0')
		h('forall for (x,0) totG += (area(x)*Gval)*(1e+6/1e+8)')
		Bestsomao.append(h.Gval)
		Besto.append(h.totG)
	elif vars[i] == 'gbar_Ihall':
		h('Gval = 0')
		h.Gval = halloffamey[i][bestIndexY]
		h('totG = 0')
		h('forall for (x,0) totG += (area(x)*Gval)*(1e+6/1e+8)*(0.5+(24/(1+exp((distance(x)-950)/(-285)))))')
		Bestsomay.append(h.Gval)
		Besty.append(h.totG)
		h('Gval = 0')
		h.Gval = halloffameo[i][bestIndexO]
		h('totG = 0')
		h('forall for (x,0) totG += (area(x)*Gval)*(1e+6/1e+8)*(0.5+(24/(1+exp((distance(x)-950)/(-285)))))')
		Bestsomao.append(h.Gval)
		Besto.append(h.totG)
	else:
		Gval1 = halloffamey[i][bestIndexY]
		Gval2 = halloffamey[i+11][bestIndexY]
		h('Gval1 = 0')
		h('Gval2 = 0')
		h.Gval1 = Gval1
		h.Gval2 = Gval2
		h('totG = 0')
		h('forsec "soma" for (x,0) totG += (area(x)*Gval1)*(1e+6/1e+8)')
		h('forsec "axon" for (x,0) totG += (area(x)*Gval2)*(1e+6/1e+8)')
		Bestsomay.append(h.Gval1)
		Besty.append(h.totG)
		Gval1 = halloffameo[i][bestIndexO]
		Gval2 = halloffameo[i+11][bestIndexO]
		h('Gval1 = 0')
		h('Gval2 = 0')
		h.Gval1 = Gval1
		h.Gval2 = Gval2
		h('totG = 0')
		h('forsec "soma" for (x,0) totG += (area(x)*Gval1)*(1e+6/1e+8)')
		h('forsec "axon" for (x,0) totG += (area(x)*Gval2)*(1e+6/1e+8)')
		Bestsomao.append(h.Gval1)
		Besto.append(h.totG)

print(Besty)
print(Besto)

stats = []
xtickspoints = numpy.linspace(1,len(labels_Gtot),len(labels_Gtot))
fig = plt.figure(figsize=(10,4))
ax1 = fig.add_subplot(111)
for i in range(0,len(Gtotsy)):
	tstat, pval = st.ranksums(Gsomay[i],Gsomao[i])
	effectsize = (numpy.mean(Gsomay[i]) - numpy.mean(Gsomao[i]))/(numpy.std(numpy.concatenate((Gsomay[i],Gsomao[i]))))
	stats.append([vars_Gtot[i],tstat,pval,effectsize])
	
	fpy = numpy.array(Gsomay[i])
	fpo = numpy.array(Gsomao[i])
	normalizedy = (fpy-lowerlims_Gsoma[i])/(upperlims_Gsoma[i]-lowerlims_Gsoma[i])
	normalizedo = (fpo-lowerlims_Gsoma[i])/(upperlims_Gsoma[i]-lowerlims_Gsoma[i])
	
	min_normalizedy = numpy.percentile(normalizedy,5)
	min_normalizedo = numpy.percentile(normalizedo,5)
	max_normalizedy = numpy.percentile(normalizedy,95)
	max_normalizedo = numpy.percentile(normalizedo,95)
	
	median_normalizedy = numpy.median(normalizedy)
	median_normalizedo = numpy.median(normalizedo)
	
	best_normalizedy = (Bestsomay[i]-lowerlims_Gsoma[i])/(upperlims_Gsoma[i]-lowerlims_Gsoma[i]) # young
	best_normalizedo = (Bestsomao[i]-lowerlims_Gsoma[i])/(upperlims_Gsoma[i]-lowerlims_Gsoma[i]) # old
	
	asymmetric_errory = [[median_normalizedy-min_normalizedy], [max_normalizedy-median_normalizedy]]
	asymmetric_erroro = [[median_normalizedo-min_normalizedo], [max_normalizedo-median_normalizedo]]
	
	ax1.errorbar(xtickspoints[i]-0.1, median_normalizedy, yerr=asymmetric_errory, fmt='o', color='k')
	ax1.errorbar(xtickspoints[i]+0.1, median_normalizedo, yerr=asymmetric_erroro, fmt='o', color='r')

ax1.set_ylim(-0.01,1.01)
ax1.set_xticks(xtickspoints)
ax1.set_xticklabels(labels_Gtot, rotation=90)
fig.tight_layout()
fig.savefig('PLOTfiles/normalized_params_Gsoma.pdf')
fig.savefig('PLOTfiles/normalized_params_Gsoma.png')
plt.close(fig)

allpvals = [float(i) for i in numpy.transpose(stats)[2].tolist()]
pval_adj_BO = mt.multipletests(allpvals,method='bonferroni')
pval_adj_BH = mt.multipletests(allpvals,method='fdr_bh')
pval_adj_BO = pval_adj_BO[1]
pval_adj_BH = pval_adj_BH[1]

# Just ad list to dataframe now
df3 = pd.DataFrame(stats, columns=["Parameter", "t-stat", "p-value","Cohen's d"])
df3['p-values adjusted (Bonferroni)'] = pval_adj_BO
df3['p-values adjusted (BH)'] = pval_adj_BH
df3.to_csv('stats_somaG.csv')


stats = []
xtickspoints = numpy.linspace(1,len(labels_Gtot),len(labels_Gtot))
fig = plt.figure(figsize=(10,4))
ax1 = fig.add_subplot(111)
for i in range(0,len(Gtotsy)):
	tstat, pval = st.ranksums(Gtotsy[i],Gtotso[i])
	effectsize = (numpy.mean(Gtotsy[i]) - numpy.mean(Gtotso[i]))/(numpy.std(numpy.concatenate((Gtotsy[i],Gtotso[i]))))
	stats.append([vars_Gtot[i],tstat,pval,effectsize])
	
	fpy = numpy.array(Gtotsy[i])
	fpo = numpy.array(Gtotso[i])
	normalizedy = (fpy-lowerlims_Gtot[i])/(upperlims_Gtot[i]-lowerlims_Gtot[i])
	normalizedo = (fpo-lowerlims_Gtot[i])/(upperlims_Gtot[i]-lowerlims_Gtot[i])
	
	min_normalizedy = numpy.percentile(normalizedy,5)
	min_normalizedo = numpy.percentile(normalizedo,5)
	max_normalizedy = numpy.percentile(normalizedy,95)
	max_normalizedo = numpy.percentile(normalizedo,95)
	
	median_normalizedy = numpy.median(normalizedy)
	median_normalizedo = numpy.median(normalizedo)
	
	best_normalizedy = (Besty[i]-lowerlims_Gtot[i])/(upperlims_Gtot[i]-lowerlims_Gtot[i]) # young
	best_normalizedo = (Besto[i]-lowerlims_Gtot[i])/(upperlims_Gtot[i]-lowerlims_Gtot[i]) # old
	
	asymmetric_errory = [[median_normalizedy-min_normalizedy], [max_normalizedy-median_normalizedy]]
	asymmetric_erroro = [[median_normalizedo-min_normalizedo], [max_normalizedo-median_normalizedo]]
	
	ax1.errorbar(xtickspoints[i]-0.1, median_normalizedy, yerr=asymmetric_errory, fmt='o', color='k')
	ax1.errorbar(xtickspoints[i]+0.1, median_normalizedo, yerr=asymmetric_erroro, fmt='o', color='r')

ax1.set_ylim(-0.01,1.01)
ax1.set_xticks(xtickspoints)
ax1.set_xticklabels(labels_Gtot, rotation=90)
fig.tight_layout()
fig.savefig('PLOTfiles/normalized_params_Gtot.pdf')
fig.savefig('PLOTfiles/normalized_params_Gtot.png')
plt.close(fig)

allpvals = [float(i) for i in numpy.transpose(stats)[2].tolist()]
pval_adj_BO = mt.multipletests(allpvals,method='bonferroni')
pval_adj_BH = mt.multipletests(allpvals,method='fdr_bh')
pval_adj_BO = pval_adj_BO[1]
pval_adj_BH = pval_adj_BH[1]

# Just ad list to dataframe now
df2 = pd.DataFrame(stats, columns=["Parameter", "t-stat", "p-value","Cohen's d"])
df2['p-values adjusted (Bonferroni)'] = pval_adj_BO
df2['p-values adjusted (BH)'] = pval_adj_BH
df2.to_csv('stats_totalG.csv')

font = {'family' : 'normal',
		'weight' : 'normal',
		'size'   : 16}

matplotlib.rc('font', **font)

# Stats and normal tests
stats = []
for i in range(0,len(finalpopy)):
	
	tstat, pval = st.ranksums(finalpopy[i],finalpopo[i])
	effectsize = (numpy.mean(finalpopy[i]) - numpy.mean(finalpopo[i]))/(numpy.std(numpy.concatenate((finalpopy[i],finalpopo[i]))))
	stats.append([vars[i],tstat,pval,effectsize])
	
	plt.hist(finalpopy[i],50,facecolor='k', alpha=0.5,label='Young')
	plt.hist(finalpopo[i],50,facecolor='darkred', alpha=0.5,label='Old')
	bottom, top = plt.ylim()
	plt.plot([halloffamey[i][bestIndexY],halloffamey[i][bestIndexY]],[bottom,top],linestyle='dashed',color='k')
	plt.plot([halloffameo[i][bestIndexO],halloffameo[i][bestIndexO]],[bottom,top],linestyle='dashed',color='darkred')
	plt.ylim((bottom,top))
	# plt.xlim((lowerlims[i],upperlims[i]))
	if ((vars[i] == 'decay_CaDynamicssomatic') or (vars[i] == 'decay_CaDynamicsaxonal')):
		plt.xlabel(labels[i])
	elif (vars[i] == 'e_pas'):
		plt.xlabel(labels[i] + r' (mV)')
	else:
		plt.xlabel(labels[i] + r' (S/cm$^{2}$)')
		plt.ticklabel_format(style='sci', axis='x', scilimits=(0,0))
	
	# plt.ylabel('Count in Final Population')
	if (vars[i] == 'g_pas'):
		plt.legend()
	
	plt.tight_layout()
	plt.savefig('PLOTfiles/' + vars[i] + '.pdf', bbox_inches='tight')
	plt.savefig('PLOTfiles/' + vars[i] + '.png', bbox_inches='tight')
	plt.gcf().clear()
	plt.cla()
	plt.clf()
	plt.close()

allpvals = [float(i) for i in numpy.transpose(stats)[2].tolist()]
pval_adj_BO = mt.multipletests(allpvals,method='bonferroni')
pval_adj_BH = mt.multipletests(allpvals,method='fdr_bh')
pval_adj_BO = pval_adj_BO[1]
pval_adj_BH = pval_adj_BH[1]

# Just ad list to dataframe now
df = pd.DataFrame(stats, columns=["Parameter", "t-stat", "p-value","Cohen's d"])
df['p-values adjusted (Bonferroni)'] = pval_adj_BO
df['p-values adjusted (BH)'] = pval_adj_BH
df.to_csv('stats.csv')

font = {'family' : 'normal',
		'weight' : 'normal',
		'size'   : 20}

matplotlib.rc('font', **font)

xtickspoints = numpy.linspace(1,len(labels),len(labels))
fig = plt.figure(figsize=(10,4))
ax1 = fig.add_subplot(111)
for i in range(0,len(finalpopy)):
	fpy = finalpopy[i]
	fpo = finalpopo[i]
	normalizedy = (fpy-lowerlims[i])/(upperlims[i]-lowerlims[i])
	normalizedo = (fpo-lowerlims[i])/(upperlims[i]-lowerlims[i])
	
	min_normalizedy = numpy.percentile(normalizedy,5)
	min_normalizedo = numpy.percentile(normalizedo,5)
	max_normalizedy = numpy.percentile(normalizedy,95)
	max_normalizedo = numpy.percentile(normalizedo,95)
	
	median_normalizedy = numpy.median(normalizedy)
	median_normalizedo = numpy.median(normalizedo)
	
	best_normalizedy = (halloffamey[i][bestIndexY]-lowerlims[i])/(upperlims[i]-lowerlims[i]) # young
	best_normalizedo = (halloffameo[i][bestIndexO]-lowerlims[i])/(upperlims[i]-lowerlims[i]) # old
	
	asymmetric_errory = [[median_normalizedy-min_normalizedy], [max_normalizedy-median_normalizedy]]
	asymmetric_erroro = [[median_normalizedo-min_normalizedo], [max_normalizedo-median_normalizedo]]
	
	ax1.errorbar(xtickspoints[i]-0.1, median_normalizedy, yerr=asymmetric_errory, fmt='o', color='k')
	ax1.errorbar(xtickspoints[i]+0.1, median_normalizedo, yerr=asymmetric_erroro, fmt='o', color='r')

ax1.set_xticks(xtickspoints)
ax1.set_xticklabels(labels, rotation=90)
fig.tight_layout()
fig.savefig('PLOTfiles/normalized_params.pdf')
fig.savefig('PLOTfiles/normalized_params.png')
plt.close(fig)

font = {'family' : 'normal',
		'weight' : 'normal',
		'size'   : 14}

matplotlib.rc('font', **font)

xtickspoints = numpy.array([1])
fig = plt.figure(figsize=(4,4))
ax1 = fig.add_subplot(111)
i = 12

fpy = finalpopy[i]
fpo = finalpopo[i]
normalizedy = (fpy-lowerlims[i])/(upperlims[i]-lowerlims[i])
normalizedo = (fpo-lowerlims[i])/(upperlims[i]-lowerlims[i])

min_normalizedy = numpy.percentile(normalizedy,5)
min_normalizedo = numpy.percentile(normalizedo,5)
max_normalizedy = numpy.percentile(normalizedy,95)
max_normalizedo = numpy.percentile(normalizedo,95)

best_normalizedy = (halloffamey[i][bestIndexY]-lowerlims[i])/(upperlims[i]-lowerlims[i]) # young
best_normalizedo = (halloffameo[i][bestIndexO]-lowerlims[i])/(upperlims[i]-lowerlims[i]) # old

asymmetric_errory = [[best_normalizedy-min_normalizedy], [max_normalizedy-best_normalizedy]]
asymmetric_erroro = [[best_normalizedo-min_normalizedo], [max_normalizedo-best_normalizedo]]

ax1.errorbar(xtickspoints-0.1, best_normalizedy, yerr=asymmetric_errory, fmt='o', color='k')
ax1.errorbar(xtickspoints+0.1, best_normalizedo, yerr=asymmetric_erroro, fmt='o', color='r')

ax1.set_xticks([xtickspoints-0.1,xtickspoints-0.1])
ax1.set_xlim(0.8,1.2)
ax1.set_xticklabels('', rotation=90)
fig.tight_layout()
fig.savefig('PLOTfiles/normalizedIH_params.pdf', bbox_inches='tight', dpi=300, transparent=True)
fig.savefig('PLOTfiles/normalizedIH_params.png', bbox_inches='tight', dpi=300, transparent=True)
plt.close(fig)

xtickspoints = numpy.array([1])
fig = plt.figure(figsize=(6,4))
ax1 = fig.add_subplot(111)
i = 12

fpy = finalpopy[i]
fpo = finalpopo[i]
normalizedy = (fpy-lowerlims[i])/(upperlims[i]-lowerlims[i])
normalizedo = (fpo-lowerlims[i])/(upperlims[i]-lowerlims[i])

min_normalizedy = numpy.percentile(normalizedy,5)
min_normalizedo = numpy.percentile(normalizedo,5)
max_normalizedy = numpy.percentile(normalizedy,95)
max_normalizedo = numpy.percentile(normalizedo,95)

best_normalizedy = (halloffamey[i][bestIndexY]-lowerlims[i])/(upperlims[i]-lowerlims[i]) # young
best_normalizedo = (halloffameo[i][bestIndexO]-lowerlims[i])/(upperlims[i]-lowerlims[i]) # old

asymmetric_errory = [[best_normalizedy-min_normalizedy], [max_normalizedy-best_normalizedy]]
asymmetric_erroro = [[best_normalizedo-min_normalizedo], [max_normalizedo-best_normalizedo]]

yvals = normalizedy[(normalizedy>min_normalizedy) & (normalizedy<max_normalizedy)]
ovals = normalizedo[(normalizedo>min_normalizedo) & (normalizedo<max_normalizedo)]

parts1 = ax1.violinplot(yvals, positions=xtickspoints-0.25, showmeans=False, showmedians=False,
		showextrema=False)
parts2 = ax1.violinplot(ovals, positions=xtickspoints+0.25, showmeans=False, showmedians=False,
		showextrema=False)

for pc in parts1['bodies']:
	pc.set_facecolor('dimgray')
	pc.set_edgecolor('black')
	pc.set_alpha(1)
for pc in parts2['bodies']:
	pc.set_facecolor('r')
	pc.set_edgecolor('r')
	pc.set_alpha(0.7)

xsesy = numpy.random.rand(len(yvals))*0.4+0.55
xseso = numpy.random.rand(len(ovals))*0.4+1.05
ax1.scatter(xsesy,yvals,s=1,c='k',marker='o')
ax1.scatter(xseso,ovals,s=1,c='k',marker='o')
ax1.scatter(xtickspoints-0.25,best_normalizedy,s=60,c='gold',marker='o',linewidths=0.5,edgecolors='k')
ax1.scatter(xtickspoints+0.25,best_normalizedo,s=60,c='gold',marker='o',linewidths=0.5,edgecolors='k')

ax1.spines['right'].set_visible(False)
ax1.spines['top'].set_visible(False)
ax1.set_xticks([xtickspoints-0.25,xtickspoints+0.25])
ax1.set_xticklabels('', rotation=90)
ax1.set_xlim(0.4,1.6)
fig.tight_layout()
fig.savefig('PLOTfiles/normalizedIHviolin_params.pdf', bbox_inches='tight', dpi=300, transparent=True)
fig.savefig('PLOTfiles/normalizedIHviolin_params.png', bbox_inches='tight', dpi=300, transparent=True)
plt.close(fig)