# -*- coding: utf-8 -*-
"""
Created on Tue Apr 18 11:32:12 2017

@author: camacho
"""
import Gedi as gedi
import RV_function as RVfunc

import numpy as np; #np.random.seed(13042017)
import matplotlib.pylab as pl; pl.close()
import astropy.table as Table

                ##### TESTS USING MCMC #####

##### Spots data preparation ##################################################
RV_spots_data=[] #to contain all data in the end

# data from ex_amp0.rdb
rdb_data=Table.Table.read('ex_amp0.rdb',format='ascii')
RV_spot=rdb_data['RV_tot'][1:101]
RV_spot=np.array(RV_spot)
RV_spot=RV_spot.astype('Float64')
RV_spotfinal=np.concatenate((RV_spot,RV_spot,RV_spot,RV_spot),axis=0)

#to organize the data into just 30 measurments
spots_yy=[]  
for i in np.arange(4,401,4):
    a=(RV_spotfinal[i-4]+RV_spotfinal[i-3]+RV_spotfinal[i-2]+RV_spotfinal[i-1])*1000/4.
    spots_yy.append(a)
spots_data=[] 
for j in np.arange(1,100,3.3):
       spots_data.append(spots_yy[int(round(j))])
RV_spots_data.append(spots_data) #final data

# data from ex_amp1.rdb
rdb_data=Table.Table.read('ex_amp1.rdb',format='ascii')
RV_spot=rdb_data['RV_tot'][1:101]
RV_spot=np.array(RV_spot)
RV_spot=RV_spot.astype('Float64')
RV_spotfinal=np.concatenate((RV_spot,RV_spot,RV_spot,RV_spot),axis=0)

#to organize the data into just 30 measurments
spots_yy=[]  
for i in np.arange(4,401,4):
    a=(RV_spotfinal[i-4]+RV_spotfinal[i-3]+RV_spotfinal[i-2]+RV_spotfinal[i-1])*1000/4.
    spots_yy.append(a)
spots_data=[] 
for j in np.arange(1,100,3.3):
       spots_data.append(spots_yy[int(round(j))])
RV_spots_data.append(spots_data) #final data

# data from ex_amp2.rdb
rdb_data=Table.Table.read('ex_amp2.rdb',format='ascii')
RV_spot=rdb_data['RV_tot'][1:101]
RV_spot=np.array(RV_spot)
RV_spot=RV_spot.astype('Float64')
RV_spotfinal=np.concatenate((RV_spot,RV_spot,RV_spot,RV_spot),axis=0)

#to organize the data into just 30 measurments
spots_yy=[]  
for i in np.arange(4,401,4):
    a=(RV_spotfinal[i-4]+RV_spotfinal[i-3]+RV_spotfinal[i-2]+RV_spotfinal[i-1])*1000/4.
    spots_yy.append(a)
spots_data=[] 
for j in np.arange(1,100,3.3):
       spots_data.append(spots_yy[int(round(j))])
RV_spots_data.append(spots_data) #final data

# data from ex_amp3.rdb
rdb_data=Table.Table.read('ex_amp3.rdb',format='ascii')
RV_spot=rdb_data['RV_tot'][1:101]
RV_spot=np.array(RV_spot)
RV_spot=RV_spot.astype('Float64')
RV_spotfinal=np.concatenate((RV_spot,RV_spot,RV_spot,RV_spot),axis=0)

#to organize the data into just 30 measurments
spots_yy=[]  
for i in np.arange(4,401,4):
    a=(RV_spotfinal[i-4]+RV_spotfinal[i-3]+RV_spotfinal[i-2]+RV_spotfinal[i-1])*1000/4.
    spots_yy.append(a)
spots_data=[] 
for j in np.arange(1,100,3.3):
       spots_data.append(spots_yy[int(round(j))])
RV_spots_data.append(spots_data) #final data

# data from ex_amp4.rdb
rdb_data=Table.Table.read('ex_amp4.rdb',format='ascii')
RV_spot=rdb_data['RV_tot'][1:101]
RV_spot=np.array(RV_spot)
RV_spot=RV_spot.astype('Float64')
RV_spotfinal=np.concatenate((RV_spot,RV_spot,RV_spot,RV_spot),axis=0)

#to organize the data into just 30 measurments
spots_yy=[]  
for i in np.arange(4,401,4):
    a=(RV_spotfinal[i-4]+RV_spotfinal[i-3]+RV_spotfinal[i-2]+RV_spotfinal[i-1])*1000/4.
    spots_yy.append(a)
spots_data=[] 
for j in np.arange(1,100,3.3):
       spots_data.append(spots_yy[int(round(j))])
RV_spots_data.append(spots_data) #final data

# data from ex_amp5.rdb
rdb_data=Table.Table.read('ex_amp5.rdb',format='ascii')
RV_spot=rdb_data['RV_tot'][1:101]
RV_spot=np.array(RV_spot)
RV_spot=RV_spot.astype('Float64')
RV_spotfinal=np.concatenate((RV_spot,RV_spot,RV_spot,RV_spot),axis=0)

#to organize the data into just 30 measurments
spots_yy=[]  
for i in np.arange(4,401,4):
    a=(RV_spotfinal[i-4]+RV_spotfinal[i-3]+RV_spotfinal[i-2]+RV_spotfinal[i-1])*1000/4.
    spots_yy.append(a)
spots_data=[] 
for j in np.arange(1,100,3.3):
       spots_data.append(spots_yy[int(round(j))])
RV_spots_data.append(spots_data) #final data

# data from ex_amp6.rdb
rdb_data=Table.Table.read('ex_amp6.rdb',format='ascii')
RV_spot=rdb_data['RV_tot'][1:101]
RV_spot=np.array(RV_spot)
RV_spot=RV_spot.astype('Float64')
RV_spotfinal=np.concatenate((RV_spot,RV_spot,RV_spot,RV_spot),axis=0)

#to organize the data into just 30 measurments
spots_yy=[]  
for i in np.arange(4,401,4):
    a=(RV_spotfinal[i-4]+RV_spotfinal[i-3]+RV_spotfinal[i-2]+RV_spotfinal[i-1])*1000/4.
    spots_yy.append(a)
spots_data=[] 
for j in np.arange(1,100,3.3):
       spots_data.append(spots_yy[int(round(j))])
RV_spots_data.append(spots_data) #final data

# data from ex_amp7.rdb
rdb_data=Table.Table.read('ex_amp7.rdb',format='ascii')
RV_spot=rdb_data['RV_tot'][1:101]
RV_spot=np.array(RV_spot)
RV_spot=RV_spot.astype('Float64')
RV_spotfinal=np.concatenate((RV_spot,RV_spot,RV_spot,RV_spot),axis=0)

#to organize the data into just 30 measurments
spots_yy=[]  
for i in np.arange(4,401,4):
    a=(RV_spotfinal[i-4]+RV_spotfinal[i-3]+RV_spotfinal[i-2]+RV_spotfinal[i-1])*1000/4.
    spots_yy.append(a)
spots_data=[] 
for j in np.arange(1,100,3.3):
       spots_data.append(spots_yy[int(round(j))])
RV_spots_data.append(spots_data) #final data

# data from ex_amp8.rdb
rdb_data=Table.Table.read('ex_amp8.rdb',format='ascii')
RV_spot=rdb_data['RV_tot'][1:101]
RV_spot=np.array(RV_spot)
RV_spot=RV_spot.astype('Float64')
RV_spotfinal=np.concatenate((RV_spot,RV_spot,RV_spot,RV_spot),axis=0)

#to organize the data into just 30 measurments
spots_yy=[]  
for i in np.arange(4,401,4):
    a=(RV_spotfinal[i-4]+RV_spotfinal[i-3]+RV_spotfinal[i-2]+RV_spotfinal[i-1])*1000/4.
    spots_yy.append(a)
spots_data=[] 
for j in np.arange(1,100,3.3):
       spots_data.append(spots_yy[int(round(j))])
RV_spots_data.append(spots_data) #final data

# data from ex_amp9.rdb
rdb_data=Table.Table.read('ex_amp9.rdb',format='ascii')
RV_spot=rdb_data['RV_tot'][1:101]
RV_spot=np.array(RV_spot)
RV_spot=RV_spot.astype('Float64')
RV_spotfinal=np.concatenate((RV_spot,RV_spot,RV_spot,RV_spot),axis=0)

#to organize the data into just 30 measurments
spots_yy=[]  
for i in np.arange(4,401,4):
    a=(RV_spotfinal[i-4]+RV_spotfinal[i-3]+RV_spotfinal[i-2]+RV_spotfinal[i-1])*1000/4.
    spots_yy.append(a)
spots_data=[] 
for j in np.arange(1,100,3.3):
       spots_data.append(spots_yy[int(round(j))])
RV_spots_data.append(spots_data) #final data

# data from ex_amp10.rdb
rdb_data=Table.Table.read('ex_amp10.rdb',format='ascii')
RV_spot=rdb_data['RV_tot'][1:101]
RV_spot=np.array(RV_spot)
RV_spot=RV_spot.astype('Float64')
RV_spotfinal=np.concatenate((RV_spot,RV_spot,RV_spot,RV_spot),axis=0)

#to organize the data into just 30 measurments
spots_yy=[]  
for i in np.arange(4,401,4):
    a=(RV_spotfinal[i-4]+RV_spotfinal[i-3]+RV_spotfinal[i-2]+RV_spotfinal[i-1])*1000/4.
    spots_yy.append(a)
spots_data=[] 
for j in np.arange(1,100,3.3):
       spots_data.append(spots_yy[int(round(j))])
RV_spots_data.append(spots_data) #final data


##### RV data #################################################################
# Period(P) ~ 20 e 50 days
# Observations(space) ~ every 4 days
# Error(yerr) ~ 0.20 a 0.50 m
test=RVfunc.RV_circular(P=25,K=17,T=0,gamma=0,time=100,space=30)

t=test[0] #time

y_rv=np.array(test[1])
yerr=np.array([np.random.uniform(0.2,0.5) for n in y_rv])
y=np.array([n+m for n,m in zip(y_rv,yerr)]) #rv measurements


###### Signal to noise ratios  ################################################
#for i in range(0,11):
#    y_var=np.var(y)
#    noise_var= np.var(RV_spots_data[i])
#    SNR_var= y_var/noise_var
#    print 'SNR_var_%i ='%i, SNR_var
#    
#    y_amp= np.max(y)
#    noise_amp= np.max(RV_spots_data[i])
#    SNR_amp= y_amp/noise_amp
#    print 'SNR_amp_%i ='%i, SNR_amp
#    print ''
    

##### RV signal + noise #######################################################
y_amp0=np.array([n+m for n,m in zip(y,RV_spots_data[0])])
y_amp1=np.array([n+m for n,m in zip(y,RV_spots_data[1])])
y_amp2=np.array([n+m for n,m in zip(y,RV_spots_data[2])])
y_amp3=np.array([n+m for n,m in zip(y,RV_spots_data[3])])
y_amp4=np.array([n+m for n,m in zip(y,RV_spots_data[4])])
y_amp5=np.array([n+m for n,m in zip(y,RV_spots_data[5])])
y_amp6=np.array([n+m for n,m in zip(y,RV_spots_data[6])])
y_amp7=np.array([n+m for n,m in zip(y,RV_spots_data[7])])
y_amp8=np.array([n+m for n,m in zip(y,RV_spots_data[8])])
y_amp9=np.array([n+m for n,m in zip(y,RV_spots_data[9])])
y_amp10=np.array([n+m for n,m in zip(y,RV_spots_data[10])])


##### MCMC preparation ########################################################
runs=100000
#burns=round(runs/2)
burns=round(runs*0.75)

theta_means=[];theta_errors=[]
P_means=[];P_errors=[]

graphs=0 #0 graphs off, 1 graphs on

kernel= gedi.kernel.ExpSineSquared(1,1,1) + gedi.kernel.WhiteNoise(1)
parameters= [[15.0,20.0], [5.0,10.0], [20.0,30.0], [0.1,2.0]]


##### Optimization preparation ################################################
theta_means=[];theta_errors=[]
P_means=[];P_errors=[]


##### Amplitude ratio 0.5 #####################################################
print '### SNR 0.5 ###'

thetaESS_final= [];lESS_final=[];PESS_final=[];thetaWN_final=[]
for i in range(1,11):
    #initial conditions
    trial=gedi.kernel_mcmc.MCMC(kernel,t,y_amp0,yerr,parameters,runs,burns)
    print trial[0], trial[1]
    
    #graphics
    if graphs==1:
        xt=np.linspace(burns,runs,runs-burns)
        pl.figure()
        pl.title('log marginal likelihood')
        pl.plot(xt,trial[2],'k-')
        
        f, axarr = pl.subplots(2, 2)
        axarr[0, 0].plot(xt, trial[3][0])
        axarr[0, 0].set_title('amplitude')
        axarr[0, 1].plot(xt, trial[3][1])
        axarr[0, 1].set_title('lenght scale')
        axarr[1, 0].plot(xt, trial[3][2])
        axarr[1, 0].set_title('period')
        axarr[1, 1].plot(xt, trial[3][3])
        axarr[1, 1].set_title('white noise')
        pl.setp([a.get_xticklabels() for a in axarr[0, :]], visible=False)
    if graphs==0:
        pass
    
    #save results    
    thetaESS_final.append(trial[0].pars[0])
    lESS_final.append(trial[0].pars[1])
    PESS_final.append(trial[0].pars[2])
    thetaWN_final.append(trial[0].pars[3])
    
amp05_thetaM=np.mean(thetaESS_final); theta_means.append(amp05_thetaM)
amp05_thetaE=np.std(thetaESS_final); theta_errors.append(amp05_thetaE)
amp05_PM=np.mean(PESS_final); P_means.append(amp05_PM)
amp05_PE=np.std(PESS_final); P_errors.append(amp05_PE) 


##### Amplitude ratio 1 #######################################################
print '### SNR 1 ###'

thetaESS_final= [];lESS_final=[];PESS_final=[];thetaWN_final=[]
for i in range(1,11):
    #initial conditions
    trial=gedi.kernel_mcmc.MCMC(kernel,t,y_amp1,yerr,parameters,runs,burns)
    print trial[0], trial[1]
    
    #graphics
    if graphs==1:
        xt=np.linspace(burns,runs,runs-burns)
        pl.figure()
        pl.title('log marginal likelihood')
        pl.plot(xt,trial[2],'k-')
        
        f, axarr = pl.subplots(2, 2)
        axarr[0, 0].plot(xt, trial[3][0])
        axarr[0, 0].set_title('amplitude')
        axarr[0, 1].plot(xt, trial[3][1])
        axarr[0, 1].set_title('lenght scale')
        axarr[1, 0].plot(xt, trial[3][2])
        axarr[1, 0].set_title('period')
        axarr[1, 1].plot(xt, trial[3][3])
        axarr[1, 1].set_title('white noise')
        pl.setp([a.get_xticklabels() for a in axarr[0, :]], visible=False)
    if graphs==0:
        pass
    
    #save results    
    thetaESS_final.append(trial[0].pars[0])
    lESS_final.append(trial[0].pars[1])
    PESS_final.append(trial[0].pars[2])
    thetaWN_final.append(trial[0].pars[3])
    
amp1_thetaM=np.mean(thetaESS_final); theta_means.append(amp1_thetaM)
amp1_thetaE=np.std(thetaESS_final); theta_errors.append(amp1_thetaE)
amp1_PM=np.mean(PESS_final); P_means.append(amp1_PM)
amp1_PE=np.std(PESS_final); P_errors.append(amp1_PE)


##### Amplitude ratio 2 #######################################################
print '### SNR 2 ###'

thetaESS_final= [];lESS_final=[];PESS_final=[];thetaWN_final=[]
for i in range(1,11):
    #initial conditions
    trial=gedi.kernel_mcmc.MCMC(kernel,t,y_amp2,yerr,parameters,runs,burns)
    print trial[0], trial[1]
    
    #graphics
    if graphs==1:
        xt=np.linspace(burns,runs,runs-burns)
        pl.figure()
        pl.title('log marginal likelihood')
        pl.plot(xt,trial[2],'k-')
        
        f, axarr = pl.subplots(2, 2)
        axarr[0, 0].plot(xt, trial[3][0])
        axarr[0, 0].set_title('amplitude')
        axarr[0, 1].plot(xt, trial[3][1])
        axarr[0, 1].set_title('lenght scale')
        axarr[1, 0].plot(xt, trial[3][2])
        axarr[1, 0].set_title('period')
        axarr[1, 1].plot(xt, trial[3][3])
        axarr[1, 1].set_title('white noise')
        pl.setp([a.get_xticklabels() for a in axarr[0, :]], visible=False)
    if graphs==0:
        pass
    
    #save results    
    thetaESS_final.append(trial[0].pars[0])
    lESS_final.append(trial[0].pars[1])
    PESS_final.append(trial[0].pars[2])
    thetaWN_final.append(trial[0].pars[3])
    
amp2_thetaM=np.mean(thetaESS_final); theta_means.append(amp2_thetaM)
amp2_thetaE=np.std(thetaESS_final); theta_errors.append(amp2_thetaE)
amp2_PM=np.mean(PESS_final); P_means.append(amp2_PM)
amp2_PE=np.std(PESS_final); P_errors.append(amp2_PE)     


##### Amplitude ratio 3 #######################################################
print '### SNR 3 ###'

thetaESS_final= [];lESS_final=[];PESS_final=[];thetaWN_final=[]
for i in range(1,11):
    #initial conditions
    trial=gedi.kernel_mcmc.MCMC(kernel,t,y_amp3,yerr,parameters,runs,burns)
    print trial[0], trial[1]
    
    #graphics
    if graphs==1:
        xt=np.linspace(burns,runs,runs-burns)
        pl.figure()
        pl.title('log marginal likelihood')
        pl.plot(xt,trial[2],'k-')
        
        f, axarr = pl.subplots(2, 2)
        axarr[0, 0].plot(xt, trial[3][0])
        axarr[0, 0].set_title('amplitude')
        axarr[0, 1].plot(xt, trial[3][1])
        axarr[0, 1].set_title('lenght scale')
        axarr[1, 0].plot(xt, trial[3][2])
        axarr[1, 0].set_title('period')
        axarr[1, 1].plot(xt, trial[3][3])
        axarr[1, 1].set_title('white noise')
        pl.setp([a.get_xticklabels() for a in axarr[0, :]], visible=False)
    if graphs==0:
        pass
    
    #save results    
    thetaESS_final.append(trial[0].pars[0])
    lESS_final.append(trial[0].pars[1])
    PESS_final.append(trial[0].pars[2])
    thetaWN_final.append(trial[0].pars[3])
    
amp3_thetaM=np.mean(thetaESS_final); theta_means.append(amp3_thetaM)
amp3_thetaE=np.std(thetaESS_final); theta_errors.append(amp3_thetaE)
amp3_PM=np.mean(PESS_final); P_means.append(amp3_PM)
amp3_PE=np.std(PESS_final); P_errors.append(amp3_PE)     


##### Amplitude ratio 4 #######################################################
print '### SNR 4 ###'

thetaESS_final= [];lESS_final=[];PESS_final=[];thetaWN_final=[]
for i in range(1,11):
    #initial conditions
    trial=gedi.kernel_mcmc.MCMC(kernel,t,y_amp4,yerr,parameters,runs,burns)
    print trial[0], trial[1]
    
    #graphics
    if graphs==1:
        xt=np.linspace(burns,runs,runs-burns)
        pl.figure()
        pl.title('log marginal likelihood')
        pl.plot(xt,trial[2],'k-')
        
        f, axarr = pl.subplots(2, 2)
        axarr[0, 0].plot(xt, trial[3][0])
        axarr[0, 0].set_title('amplitude')
        axarr[0, 1].plot(xt, trial[3][1])
        axarr[0, 1].set_title('lenght scale')
        axarr[1, 0].plot(xt, trial[3][2])
        axarr[1, 0].set_title('period')
        axarr[1, 1].plot(xt, trial[3][3])
        axarr[1, 1].set_title('white noise')
        pl.setp([a.get_xticklabels() for a in axarr[0, :]], visible=False)
    if graphs==0:
        pass
    
    #save results    
    thetaESS_final.append(trial[0].pars[0])
    lESS_final.append(trial[0].pars[1])
    PESS_final.append(trial[0].pars[2])
    thetaWN_final.append(trial[0].pars[3])
    
amp4_thetaM=np.mean(thetaESS_final); theta_means.append(amp4_thetaM)
amp4_thetaE=np.std(thetaESS_final); theta_errors.append(amp4_thetaE)
amp4_PM=np.mean(PESS_final); P_means.append(amp4_PM)
amp4_PE=np.std(PESS_final); P_errors.append(amp4_PE)     


##### Amplitude ratio 5 #######################################################
print '### SNR 5 ###'

thetaESS_final= [];lESS_final=[];PESS_final=[];thetaWN_final=[]
for i in range(1,11):
    #initial conditions
    trial=gedi.kernel_mcmc.MCMC(kernel,t,y_amp5,yerr,parameters,runs,burns)
    print trial[0], trial[1]
    
    #graphics
    if graphs==1:
        xt=np.linspace(burns,runs,runs-burns)
        pl.figure()
        pl.title('log marginal likelihood')
        pl.plot(xt,trial[2],'k-')
        
        f, axarr = pl.subplots(2, 2)
        axarr[0, 0].plot(xt, trial[3][0])
        axarr[0, 0].set_title('amplitude')
        axarr[0, 1].plot(xt, trial[3][1])
        axarr[0, 1].set_title('lenght scale')
        axarr[1, 0].plot(xt, trial[3][2])
        axarr[1, 0].set_title('period')
        axarr[1, 1].plot(xt, trial[3][3])
        axarr[1, 1].set_title('white noise')
        pl.setp([a.get_xticklabels() for a in axarr[0, :]], visible=False)
    if graphs==0:
        pass
    
    #save results    
    thetaESS_final.append(trial[0].pars[0])
    lESS_final.append(trial[0].pars[1])
    PESS_final.append(trial[0].pars[2])
    thetaWN_final.append(trial[0].pars[3])
    
amp5_thetaM=np.mean(thetaESS_final); theta_means.append(amp5_thetaM)
amp5_thetaE=np.std(thetaESS_final); theta_errors.append(amp5_thetaE)
amp5_PM=np.mean(PESS_final); P_means.append(amp5_PM)
amp5_PE=np.std(PESS_final); P_errors.append(amp5_PE)


##### Amplitude ratio 6 #######################################################
print '### SNR 6 ###'

thetaESS_final= [];lESS_final=[];PESS_final=[];thetaWN_final=[]
for i in range(1,11):
    #initial conditions
    trial=gedi.kernel_mcmc.MCMC(kernel,t,y_amp6,yerr,parameters,runs,burns)
    print trial[0], trial[1]
    
    #graphics
    if graphs==1:
        xt=np.linspace(burns,runs,runs-burns)
        pl.figure()
        pl.title('log marginal likelihood')
        pl.plot(xt,trial[2],'k-')
        
        f, axarr = pl.subplots(2, 2)
        axarr[0, 0].plot(xt, trial[3][0])
        axarr[0, 0].set_title('amplitude')
        axarr[0, 1].plot(xt, trial[3][1])
        axarr[0, 1].set_title('lenght scale')
        axarr[1, 0].plot(xt, trial[3][2])
        axarr[1, 0].set_title('period')
        axarr[1, 1].plot(xt, trial[3][3])
        axarr[1, 1].set_title('white noise')
        pl.setp([a.get_xticklabels() for a in axarr[0, :]], visible=False)
    if graphs==0:
        pass
    
    #save results    
    thetaESS_final.append(trial[0].pars[0])
    lESS_final.append(trial[0].pars[1])
    PESS_final.append(trial[0].pars[2])
    thetaWN_final.append(trial[0].pars[3])
    
amp6_thetaM=np.mean(thetaESS_final); theta_means.append(amp6_thetaM)
amp6_thetaE=np.std(thetaESS_final); theta_errors.append(amp6_thetaE)
amp6_PM=np.mean(PESS_final); P_means.append(amp6_PM)
amp6_PE=np.std(PESS_final); P_errors.append(amp6_PE)     


##### Amplitude ratio 7 #######################################################
print '### SNR 7 ###'

thetaESS_final= [];lESS_final=[];PESS_final=[];thetaWN_final=[]
for i in range(1,11):
    #initial conditions
    trial=gedi.kernel_mcmc.MCMC(kernel,t,y_amp7,yerr,parameters,runs,burns)
    print trial[0], trial[1]
    
    #graphics
    if graphs==1:
        xt=np.linspace(burns,runs,runs-burns)
        pl.figure()
        pl.title('log marginal likelihood')
        pl.plot(xt,trial[2],'k-')
        
        f, axarr = pl.subplots(2, 2)
        axarr[0, 0].plot(xt, trial[3][0])
        axarr[0, 0].set_title('amplitude')
        axarr[0, 1].plot(xt, trial[3][1])
        axarr[0, 1].set_title('lenght scale')
        axarr[1, 0].plot(xt, trial[3][2])
        axarr[1, 0].set_title('period')
        axarr[1, 1].plot(xt, trial[3][3])
        axarr[1, 1].set_title('white noise')
        pl.setp([a.get_xticklabels() for a in axarr[0, :]], visible=False)
    if graphs==0:
        pass
    
    #save results    
    thetaESS_final.append(trial[0].pars[0])
    lESS_final.append(trial[0].pars[1])
    PESS_final.append(trial[0].pars[2])
    thetaWN_final.append(trial[0].pars[3])
    
amp7_thetaM=np.mean(thetaESS_final); theta_means.append(amp7_thetaM)
amp7_thetaE=np.std(thetaESS_final); theta_errors.append(amp7_thetaE)
amp7_PM=np.mean(PESS_final); P_means.append(amp7_PM)
amp7_PE=np.std(PESS_final); P_errors.append(amp7_PE)     


##### Amplitude ratio 8 #######################################################
print '### SNR 8 ###'

thetaESS_final= [];lESS_final=[];PESS_final=[];thetaWN_final=[]
for i in range(1,11):
    #initial conditions
    trial=gedi.kernel_mcmc.MCMC(kernel,t,y_amp8,yerr,parameters,runs,burns)
    print trial[0], trial[1]
    
    #graphics
    if graphs==1:
        xt=np.linspace(burns,runs,runs-burns)
        pl.figure()
        pl.title('log marginal likelihood')
        pl.plot(xt,trial[2],'k-')
        
        f, axarr = pl.subplots(2, 2)
        axarr[0, 0].plot(xt, trial[3][0])
        axarr[0, 0].set_title('amplitude')
        axarr[0, 1].plot(xt, trial[3][1])
        axarr[0, 1].set_title('lenght scale')
        axarr[1, 0].plot(xt, trial[3][2])
        axarr[1, 0].set_title('period')
        axarr[1, 1].plot(xt, trial[3][3])
        axarr[1, 1].set_title('white noise')
        pl.setp([a.get_xticklabels() for a in axarr[0, :]], visible=False)
    if graphs==0:
        pass
    
    #save results    
    thetaESS_final.append(trial[0].pars[0])
    lESS_final.append(trial[0].pars[1])
    PESS_final.append(trial[0].pars[2])
    thetaWN_final.append(trial[0].pars[3])
    
amp8_thetaM=np.mean(thetaESS_final); theta_means.append(amp8_thetaM)
amp8_thetaE=np.std(thetaESS_final); theta_errors.append(amp8_thetaE)
amp8_PM=np.mean(PESS_final); P_means.append(amp8_PM)
amp8_PE=np.std(PESS_final); P_errors.append(amp8_PE)     


##### Amplitude ratio 9 #######################################################
print '### SNR 9 ###'

thetaESS_final= [];lESS_final=[];PESS_final=[];thetaWN_final=[]
for i in range(1,11):
    #initial conditions
    trial=gedi.kernel_mcmc.MCMC(kernel,t,y_amp9,yerr,parameters,runs,burns)
    print trial[0], trial[1]
    
    #graphics
    if graphs==1:
        xt=np.linspace(burns,runs,runs-burns)
        pl.figure()
        pl.title('log marginal likelihood')
        pl.plot(xt,trial[2],'k-')
        
        f, axarr = pl.subplots(2, 2)
        axarr[0, 0].plot(xt, trial[3][0])
        axarr[0, 0].set_title('amplitude')
        axarr[0, 1].plot(xt, trial[3][1])
        axarr[0, 1].set_title('lenght scale')
        axarr[1, 0].plot(xt, trial[3][2])
        axarr[1, 0].set_title('period')
        axarr[1, 1].plot(xt, trial[3][3])
        axarr[1, 1].set_title('white noise')
        pl.setp([a.get_xticklabels() for a in axarr[0, :]], visible=False)
    if graphs==0:
        pass
    
    #save results    
    thetaESS_final.append(trial[0].pars[0])
    lESS_final.append(trial[0].pars[1])
    PESS_final.append(trial[0].pars[2])
    thetaWN_final.append(trial[0].pars[3])
    
amp9_thetaM=np.mean(thetaESS_final); theta_means.append(amp9_thetaM)
amp9_thetaE=np.std(thetaESS_final); theta_errors.append(amp9_thetaE)
amp9_PM=np.mean(PESS_final); P_means.append(amp9_PM)
amp9_PE=np.std(PESS_final); P_errors.append(amp9_PE)     


##### Amplitude ratio 10 ######################################################
print '### SNR 10 ###'

thetaESS_final= [];lESS_final=[];PESS_final=[];thetaWN_final=[]
for i in range(1,11):
    #initial conditions
    trial=gedi.kernel_mcmc.MCMC(kernel,t,y_amp10,yerr,parameters,runs,burns)
    print trial[0], trial[1]
    
    #graphics
    if graphs==1:
        xt=np.linspace(burns,runs,runs-burns)
        pl.figure()
        pl.title('log marginal likelihood')
        pl.plot(xt,trial[2],'k-')
        
        f, axarr = pl.subplots(2, 2)
        axarr[0, 0].plot(xt, trial[3][0])
        axarr[0, 0].set_title('amplitude')
        axarr[0, 1].plot(xt, trial[3][1])
        axarr[0, 1].set_title('lenght scale')
        axarr[1, 0].plot(xt, trial[3][2])
        axarr[1, 0].set_title('period')
        axarr[1, 1].plot(xt, trial[3][3])
        axarr[1, 1].set_title('white noise')
        pl.setp([a.get_xticklabels() for a in axarr[0, :]], visible=False)
    if graphs==0:
        pass
    
    #save results    
    thetaESS_final.append(trial[0].pars[0])
    lESS_final.append(trial[0].pars[1])
    PESS_final.append(trial[0].pars[2])
    thetaWN_final.append(trial[0].pars[3])
    
amp10_thetaM=np.mean(thetaESS_final); theta_means.append(amp10_thetaM)
amp10_thetaE=np.std(thetaESS_final); theta_errors.append(amp10_thetaE)
amp10_PM=np.mean(PESS_final); P_means.append(amp10_PM)
amp10_PE=np.std(PESS_final); P_errors.append(amp10_PE)


##### Mass determination ######################################################
mass_means=[(n/28.4)*((m/365.2425)**(1./3)) for n,m in zip(theta_means,P_means)]

error_first=[n*m for n,m in zip(P_means,theta_errors)]
error_second=[n*m for n,m in zip(theta_means,P_errors)]
mass_errors=[n+m for n,m in zip(error_first,error_second)]


##### computation of the final graphic ######################################## 
X=[0.5,1,2,3,4,5,6,7,8,9,10]
Xexample=[0,0.5,1,2,3,4,5,6,7,8,9,10,11,12]
thetaexample=[17,17,17,17,17,17,17,17,17,17,17,17,17,17]
Pexample=[25,25,25,25,25,25,25,25,25,25,25,25,25,25]
Massexample=[0.245,0.245,0.245,0.245,0.245,0.245,0.245,0.245,0.245 \
                ,0.245,0.245,0.245,0.245,0.245]

pl.figure()
pl.plot(Xexample,thetaexample,'k:')
pl.errorbar(X,theta_means,yerr=theta_errors, fmt=".k", capsize=0)
pl.title('Radial velocity semi-amplitude')
pl.xlabel("S/N")
pl.ylabel("m/s")

pl.figure()
pl.plot(Xexample,Pexample,'k:')
pl.errorbar(X,P_means,yerr=P_errors, fmt=".k", capsize=0)
pl.title('Periods')
pl.xlabel("S/N")
pl.ylabel("days")

pl.figure()
pl.plot(Xexample,Massexample,'k:')
pl.errorbar(X,mass_means,yerr=mass_errors, fmt=".k", capsize=0)
pl.title('Masses')
pl.xlabel("S/N")
pl.ylabel("Jupiter")