#! /usr/bin/env python

from numpy import *
from pylab import *
import csv

'''Code that makes plots from GEM sensitivity analysis'''

num_each=7

fnew=open('emulator_inputs1.txt','r')
files=open('emulator_outputs.txt','r')

f=csv.reader(fnew, delimiter='\t')
f.next()
totalrows=7*100*num_each

inputs=zeros((7,totalrows))
count=0
for row in f:
	for a in range(0,7):
		inputs[a,count]=row[a]
	count+=1

j=csv.reader(files, delimiter=',')
j.next()
outputs=zeros(totalrows)
count=0
for row in j:
	outputs[count]=float(row[1])
	count+=1

dpm_stuff=zeros((8,100*num_each))
mass_stuff=zeros((8,100*num_each))
sigma_stuff=zeros((8,100*num_each))
area_stuff=zeros((8,100*num_each))
vg_stuff=zeros((8,100*num_each))
time_stuff=zeros((8,100*num_each))
mixing_stuff=zeros((8,100*num_each))

tot_num=num_each*100

dpm_stuff[:7,:]=inputs[:,:tot_num]
dpm_stuff[-1,:]=outputs[:tot_num]
mass_stuff[:7,:]=inputs[:,tot_num:2*tot_num]
mass_stuff[-1,:]=outputs[tot_num:2*tot_num]
sigma_stuff[:7,:]=inputs[:,2*tot_num:3*tot_num]
sigma_stuff[-1,:]=outputs[2*tot_num:3*tot_num]
area_stuff[:7,:]=inputs[:,3*tot_num:4*tot_num]
area_stuff[-1,:]=outputs[3*tot_num:4*tot_num]
vg_stuff[:7,:]=inputs[:,4*tot_num:5*tot_num]
vg_stuff[-1,:]=outputs[4*tot_num:5*tot_num]
time_stuff[:7,:]=inputs[:,5*tot_num:6*tot_num]
time_stuff[-1,:]=outputs[5*tot_num:6*tot_num]
mixing_stuff[:7,:]=inputs[:,6*tot_num:]
mixing_stuff[-1,:]=outputs[6*tot_num:]

# Array shenanigans

#--------------
# Plots
#--------------
figure(figsize=(8,9))
j1=subplot(4,2,1)
semilogx()
xlabel('D$_{pm_0}$ [nm]')
ylabel('$\Delta \%_{\sigma}$')
for a in range(0,100):
	middle=dpm_stuff[-1,a*num_each+3]
	plot(10**dpm_stuff[0,a*num_each:(a+1)*num_each],(dpm_stuff[-1,a*num_each:(a+1)*num_each]-middle)/middle*100, 'g-')
axhline(y=0,c='k',ls='-',linewidth=2.0)
xlim(xmin=10**dpm_stuff[0].min(),xmax=10**dpm_stuff[0].max())
l=[20,50,100]
xticks(l,l)
ylim(ymax=50,ymin=-50)

j2=subplot(4,2,2)
semilogx()
for a in range(0,100):
	middle=mass_stuff[-1,a*num_each+3]
	plot(10**mass_stuff[1,a*num_each:(a+1)*num_each],(mass_stuff[-1,a*num_each:(a+1)*num_each]-middle)/middle*100,'g-')
xlabel('Mass Flux [$\mu$g m$^{-2}$ s$^{-1}$]')
ylabel('$\Delta \%_{\sigma}$')
axhline(y=0,c='k',ls='-',linewidth=2.0)
xlim(xmin=10**mass_stuff[1].min(),xmax=10**mass_stuff[1].max())

ylim(ymax=100)
ylim(ymax=50,ymin=-50)
j3=subplot(4,2,3)
for a in range(0,100):
	middle=sigma_stuff[-1,a*num_each+3]
	plot(sigma_stuff[2,a*num_each:(a+1)*num_each],(sigma_stuff[-1,a*num_each:(a+1)*num_each]-middle)/middle*100,'g-')
xlabel('$\sigma_0$')
ylabel('$\Delta \%_{\sigma}$')
axhline(y=0,c='k',ls='-',linewidth=2.0)
xlim(xmin=sigma_stuff[2].min(),xmax=sigma_stuff[2].max())
ylim(ymax=2*100)
ylim(ymax=50,ymin=-50)

j4=subplot(4,2,4)
semilogx()
for a in range(0,100):
	middle=area_stuff[-1,a*num_each+3]
	plot(10**area_stuff[3,a*num_each:(a+1)*num_each],(area_stuff[-1,a*num_each:(a+1)*num_each]-middle)/middle*100,'g-')
l=[5, 10, 50]
xlabel('Fire Area [km$^2$]')
ylabel('$\Delta \%_{\sigma}$')
axhline(y=0,c='k',ls='-',linewidth=2.0)
xlim(xmin=10**area_stuff[3].min(),xmax=10**area_stuff[3].max())
xticks(l,l)
ylim(ymax=85)
ylim(ymax=50,ymin=-50)

j5=subplot(4,2,5)
for a in range(0,100):
	middle=vg_stuff[-1,a*num_each+3]
	plot(vg_stuff[4,a*num_each:(a+1)*num_each],(vg_stuff[-1,a*num_each:(a+1)*num_each]-middle)/middle*100,'g-')

axhline(y=0,c='k',ls='-',linewidth=2.0)
xlabel('$v_g$ [m s$^{-1}$]')
ylabel('$\Delta \%_{\sigma}$')
ylim(ymax=1.2*100)
ylim(ymax=50,ymin=-50)
j6=subplot(4,2,6)

for a in range(0,100):
	middle=time_stuff[-1,a*num_each]
	plot(time_stuff[5,a*num_each:(a+1)*num_each],(time_stuff[-1,a*num_each:(a+1)*num_each]/middle-1)*100,'g-')
xlabel('Time [min]')
ylabel('$\Delta \%_{\sigma}$')
xlim(xmax=time_stuff[5].max())
axhline(y=0,c='k',ls='-',linewidth=2.0)
#subplots.tight_layout()
ylim(ymax=50,ymin=-50)

j7=subplot(4,2,7)

for a in range(0,100):
	semilogx()
	middle=mixing_stuff[-1,a*num_each+3]
	plot(10**mixing_stuff[6,a*num_each:(a+1)*num_each],(mixing_stuff[-1,a*num_each:(a+1)*num_each]-middle)/middle*100,'g-')
xlabel('d$_{mixing}$ [m]')
ylabel('$\Delta \%_{\sigma}$')
xlim(xmin=10**mixing_stuff[6].min(), xmax=2500.)
ylim(ymax=0.7*100)
ylim(ymax=50,ymin=-50)
axhline(y=0,c='k',ls='-',linewidth=2.0)
#l=arange(150,5400,num_each00)
l=[150, 500,2000]
#l=map(round,10**mixing_stuff[6,:5])
xticks(l,l)
ylim(ymax=50,ymin=-50)
tight_layout()
j1.annotate('a', xy=(0.1,0.8), xycoords='axes fraction',fontsize=18, xytext=(-5,5),textcoords='offset points',ha='right',va='bottom')
j2.annotate('b', xy=(0.1,0.8), xycoords='axes fraction',fontsize=18, xytext=(-5,5),textcoords='offset points',ha='right',va='bottom')
j3.annotate('c', xy=(0.1,0.8), xycoords='axes fraction',fontsize=18, xytext=(-5,5),textcoords='offset points',ha='right',va='bottom')
j4.annotate('d', xy=(0.1,0.8), xycoords='axes fraction',fontsize=18, xytext=(-5,5),textcoords='offset points',ha='right',va='bottom')
j5.annotate('e', xy=(0.1,0.8), xycoords='axes fraction',fontsize=18, xytext=(-5,5),textcoords='offset points',ha='right',va='bottom')
j6.annotate('f', xy=(0.1,0.8), xycoords='axes fraction',fontsize=18, xytext=(-5,5),textcoords='offset points',ha='right',va='bottom')
j7.annotate('g', xy=(0.1,0.8), xycoords='axes fraction',fontsize=18, xytext=(-5,5),textcoords='offset points',ha='right',va='bottom')
subplots_adjust(hspace=0.4)
savefig('green.pdf', format='pdf', background='transparent')
show()

