# -*- coding: utf-8 -*- 
# Nicolas Enfon 13/05/14
import numpy as np
from matplotlib.pyplot import *
import os

#Import data
data=np.loadtxt('part_file_day_h1_h2_db_ratio_ping_znb_zici_pnb_pici_moon',delimiter=',')
[part, file, day, h1, h2, db, ratio, ping, znb, zici, pnb, pici, moon] = data

#remove NaN
#sansNaN=[nb for nb in avecNaN if np.isnan(nb) == False]

#Make it boxplotable
def makeboxplotable(y,x):
    """Returns a array containing an array of y values, for each x value"""
    y_x = []
    for i in xrange( len( np.unique(x) ) ):
        y_x.append([ y[j] for j in xrange(len(y)) if x[j] == np.unique(x)[i]   ])
    return y_x

zici_day = makeboxplotable(zici,day)
znb_day = makeboxplotable(znb,day)
pici_day = makeboxplotable(pici,day)
pnb_day = makeboxplotable(pnb,day)
db_day = makeboxplotable(db,day)
ratio_day = makeboxplotable(ratio,day)

zici_moon = makeboxplotable(zici,moon)
znb_moon = makeboxplotable(znb,moon)
pici_moon = makeboxplotable(pici,moon)
pnb_moon = makeboxplotable(pnb,moon)
db_moon = makeboxplotable(db,moon)
ratio_moon = makeboxplotable(ratio,moon)

zici_h1 = makeboxplotable(zici,h1)
znb_h1 = makeboxplotable(znb,h1)
pici_h1 = makeboxplotable(pici,h1)
pnb_h1 = makeboxplotable(pnb,h1)
db_h1 = makeboxplotable(db,h1)
ratio_h1 = makeboxplotable(ratio,h1)

#Simple plots X vs Y
os.chdir('/net/nas-lsis-3/SABIOD/public_data/ANTARESV3_21/STAT_ANALYSIS/SIMPLE_PLOTS_X_VS_Y')
def var_name(var):
    """Careful, works only in some cases"""
    for name,value in globals().items() :
        if value is var :
            return name
all = [zici_day,znb_day,pici_day,pnb_day,db_day,ratio_day, zici_moon,znb_moon,pici_moon,pnb_moon,db_moon,ratio_moon, zici_h1,znb_h1,pici_h1,pnb_h1,db_h1,ratio_h1]
allname = ['zici_day','znb_day','pici_day','pnb_day','db_day','ratio_day', 'zici_moon','znb_moon','pici_moon','pnb_moon','db_moon','ratio_moon', 'zici_h1','znb_h1','pici_h1','pnb_h1','db_h1','ratio_h1']
for i in xrange(len(all)):
	figure(figsize=(20,5))
	boxplot(all[i],sym='')
	grid(True)
	title(allname[i].split('_')[0]+' VS '+allname[i].split('_')[1])
	savefig(  allname[i].split('_')[0]+'_vs_'+allname[i].split('_')[1]+'.png' ) 


#Simple histograms
os.chdir('/net/nas-lsis-3/SABIOD/public_data/ANTARESV3_21/STAT_ANALYSIS/HISTOGRAMS')
all = [db, ratio, znb, zici, pnb, pici]
allname = ['Energy_db', 'ratio_captured', 'ziphius_NB', 'ziphius_ICI', 'physeter_NB', 'physeter_ICI']
for i in xrange(len(all)):
	figure(figsize=(10,10))
	hist(all[i],bins=150)
	title('Histogram of '+allname[i])
	savefig('histo_'+allname[i]+'.png')


#2D histograms
os.chdir('/net/nas-lsis-3/SABIOD/public_data/ANTARESV3_21/STAT_ANALYSIS/HISTOGRAMS_2D')
all = [day, h1, db, ratio, znb, zici, pnb, pici, moon]
allname = ['day','h1','energy','ratio','ziphius_NB','ziphius_ICI','physeter_NB','physeter_ICI','moon']
for i in xrange(len(all)):
	for j in xrange(i+1, len(all)):
		figure(figsize=(10,10))
		h, x, y = np.histogram2d(all[j],all[i], bins=150)
		imshow(h, interpolation='nearest',extent=[y[0],y[-1],x[0],x[-1]],aspect='auto',cmap='bone',origin='lower')
		colorbar()
		xlabel(allname[i])
		ylabel(allname[j])
		title('Histogram 2D of '+allname[i]+' and '+allname[j])
		savefig('histo2d_'+allname[i]+'_'+allname[j])





"""
#R script
data=read.table('part_file_day_h1_h2_db_ratio_ping_znb_zici_pnb_pici_moon',sep=',')
data = t(data) #invert lines and rows
day = data[,3]
h1 = data[,4]
db = data[,6]
ratio = data[,7]
ping = data[,8]
znb = data[,9]
zici = data[,10]
pnb = data[,11]
pici = data[,12]
moon = data[,13]
#Caution: for a 'normal' ANOVA, we should have gaussian distribution and homoscedasticity:

"""









"""
#In case we need to reload the data
data = np.loadtxt('ToutReadmeAntaresV3_21houronly.txt',delimiter=',')
part = data[:,0]
file = data[:,1]
day_raw = data[:,2]
month = data[:,3]
h1 = data[:,5]
h2 = data[:,6]
wavsq = data[:,7]
ratio = data[:,8]
ping = data[:,9]  
znb = data[:,10]
zici = data[:,11]
pnb = data[:,14]
pici = data[:,15]
day = []
for i in xrange(len(month)):
    if month[i] == 8:
        day.append(0)
    else:
        day.append(day_raw[i] + (month[i] - 9) * 30 )
#conversions
db = np.log10(wavsq)
zici = 0.000512 * zici
pici = 0.000512 * pici
#Create the moon data: first 10 days are full moon, then nothing ...
moon = np.zeros(len(day))
i = day.index(9)
moon[:i] = 1
j = day.index(23)
k = day.index(38)
moon[j:k] = 1
l = day.index(52)
moon[l:] = 1

out = [part,file,day,h1,h2,db,ratio,ping,znb,zici,pnb,pici,moon]
np.savetxt('part_file_day_h1_h2_db_ratio_ping_znb_zici_pnb_pici_moon', out, delimiter=',')
"""