#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_wmap33_bootn/'
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') 

#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)

    #end of sample loop
    #s_betas=np.array(betas)
    #s_uncerts=np.array(uncerts)
    #final_beta=
    #invvar_uncert=math.sqrt(1./(sum(1./s_uncerts/s_uncerts)))
    #mean_beta=stat.mean(s_betas)
    #min_uncert=min(s_uncerts)
    #print('uncert(min)='+str(min_uncert))
    #std_beta=stat.stdev(s_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 beta='+str(beta)+str(uncert))
    #print('The inverse weighted beta='+str(final_beta)+'+-'+str(final_uncert))
    #print('howmanyregions='+str(howmanyregions))

    regmeanfile.write(str(reg)+'\t'+str(beta)+'\t'+str(uncert)+'\n')
    #reginvvarfile.write(str(reg)+'\t'+str(final_beta)+'\t'+str(final_uncert)+'\n')

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

    #plot them
    #if reg==1:
    #    lab=r"Cosmoglobe {\it WMAP}"
    #else:
    #    lab=""
    #plt.errorbar(reg,final_beta, yerr=final_uncert, marker='x', label=lab,linewidth=0.5, color='black',markersize=3,capsize=2,ls='none')

regmeanfile.close()