#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_pl2018/'
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=[]

#the loop over regions
for reg in range(1,25):
    betas=[]
    uncerts=[]
    #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)
        betas.append(beta)

        #read the uncert in region
        alphafilnavn=datadir+filpre+'_scatter_alpha_spectral_index.txt'
        u=[]
        with open(alphafilnavn) as f:
            for line in f:
                data = line.split()
                #x.append(float(data[0]))
                u.append(float(data[2]))
        uncert=min(u)
        print(uncert)
        #could write out the uncert to a file here
        uncerts.append(uncert)
    print('The beta='+str(beta)+str(uncert))
    regmeanfile.write(str(reg)+'\t'+str(beta)+'\t'+str(uncert)+'\n')

    ##add them to a region array, for plotting purposes                                                                             
    #reg_betas.append(beta)
    #reg_uncerts.append(uncert)

#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  #no adding in quadrature here
print('uncert(min)='+str(min_uncert))
#std_beta=stat.stdev(reg_betas)
#print('std='+str(std_beta))
#final_uncert=math.sqrt(min_uncert*min_uncert + std_beta*std_beta)
#print('final uncert (quarature)='+str(final_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()