#import astropy
#import h5py
import math
import numpy as np
#from scipy.constants import c 
#from astropy.io import fits
import subprocess
import statistics as stat
from operator import mul
from operator import truediv

# plotting variables
from setup_matplotlib import *
from matplotlib.ticker import MaxNLocator
width = 8.8
# Create the plot
fig = plt.figure(figsize=(1.25*cm2inch(width), 6./8.*cm2inch(width)))
#end plotting variables


#initial variables
dire='/mn/stornext/d5/unnif/sindex_bp/wmap23_plnpipe/'
pre='b_'     #b_
filpre='ut'  #ut
alphaarr=('a00000','a00005','a00010','a00015','a00020','a00025','a00030','a00035','a00040','a00045','a00050','a00055','a00060','a00065','a00070','a00075','a00080','a00085')
alphanum=[]
for alphas in range(0,90,5):
    alphanum.append(alphas)
print(alphanum)

meanfilnavn=dire+filpre+'_region_scatter_alphamean.txt'
regmeanfile = open(meanfilnavn, 'w') 
invvarfilnavn=dire+filpre+'_sample_betas_invvar.txt'
reginvvarfile = open(invvarfilnavn, 'w') 
uncertfilnavn=dire+filpre+'_sample_uncert.txt'
uncertfile = open(uncertfilnavn, 'w') 

#howmanyregions=0
reg_betas=[]
reg_uncerts=[]
betas=[]
uncerts=[]
systbeta=[]


#the loop over regions
for reg in range(1,25):
    print('reg '+str(reg))
    #for samp in range(1,1):
    #if 1>0:
    #read the spectral index in region
    datadir=dire+'/'+pre+str(reg)+'/'
    spfilnavn=datadir+filpre+'_value_beta_scatter.txt'
    #print(spfilnavn)
    with open(spfilnavn) as spfil:
        beta=float(spfil.read())
    print('beta= '+str(beta))
    betas.append(beta)

    #read the uncert in region
    alphafilnavn=datadir+filpre+'_scatter_alpha_spectral_index.txt'
    u=[]
    b=[]
    with open(alphafilnavn) as f:
        for line in f:
            data = line.split()
            b.append(float(data[1]))
            u.append(float(data[2]))
    stat_uncert=min(u)
    print('stat uncert= ' + str(stat_uncert))
    syst_uncert=(max(b)-min(b))/2.
    print('max beta= '+str(max(b)))
    print('min beta= '+str(min(b)))
    print('syst uncert= ' + str(syst_uncert))
    #could write out the uncert to a file here
    uncert=math.sqrt(stat_uncert*stat_uncert + syst_uncert*syst_uncert)
    uncerts.append(uncert)
    print('total uncert= ' + str(uncert))    


    print('The beta='+str(beta)+'+-  '+str(uncert))
    regmeanfile.write(str(reg)+'\t'+str(beta)+'\t'+str(uncert)+'\n')

    
#end of region loop
reg_betas=np.array(betas)
reg_uncerts=np.array(uncerts)
final_beta=sum(reg_betas/reg_uncerts/reg_uncerts)/sum(1./reg_uncerts/reg_uncerts)
#invvar_uncert=math.sqrt(1./(sum(1./s_uncerts/s_uncerts)))
#mean_beta=stat.mean(s_betas)
min_uncert=min(reg_uncerts)
final_uncert=min_uncert  
print('uncert(min)='+str(min_uncert))

print('The inverse weighted beta='+str(final_beta)+'+-'+str(final_uncert))
reginvvarfile.write(str(final_beta))
uncertfile.write(str(final_uncert))


regmeanfile.close()
reginvvarfile.close()
uncertfile.close()