#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 import itertools as it # 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))) ax = fig.add_subplot(111) #end plotting variables def get_invvarweighted(x,y) : return sum(x/y/y)/sum(1./y/y); #initial variables case='coswmap33' #cos30 or #coswmap33 dire='/mn/stornext/d5/unnif/sindex_bp/coswmap23_coswmap33ab/combab_011-195/' 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) # plot some stuff for referance to the cosmoglobe results. Either WMAP of Planck. print(case) if case=='coswmap33': #plot the old wmap # this is the combined ttplot and maxlike, do not use this! #oldwmapbetafilnavn="/mn/stornext/d5/unnif/unnif-mirfakd1/spectral/w9_n64_1deg_kopi/ut_paper_finalfinalregion_alphamean.txt" # this is the old ttplot wmap result oldwmapbetafilnavn="/mn/stornext/d5/unnif/sindex_bp/wmap23_wmap33_bootn/ut_region_scatter_alphamean.txt" with open(oldwmapbetafilnavn) as w: rx=[] bx=[] ux=[] for line in w: data = line.split() rx.append(float(data[0])) #plot the point next to the first bx.append(float(data[1])) ux.append(float(data[2])) plt.errorbar(rx,bx, yerr=ux, marker='x', label=r"9-yr {\it WMAP}",linewidth=0.5, color='red',markersize=3,capsize=2,ls='none', alpha=0.4) #plt.axhline(y=-2.99, color='red',linestyle='--',linewidth=0.5) #-2.99 is the old inv var weigh mean plt.ylim(-4.,-2.2) #plt.yticks(-4., -2.2, step=0.2) plt.yticks([-4.,-3.8, -3.6, -3.4, -3.2, -3.0, -2.8, -2.6, -2.4, -2.2],[-4.0,-3.8, -3.6, -3.4, -3.2, -3.0, -2.8, -2.6, -2.4, -2.2]) elif case=='cos30': # this is the planck2018 result oldwmapbetafilnavn="/mn/stornext/d5/unnif/sindex_bp/wmap23_pl2018/ut_region_scatter_alphamean.txt" print(oldwmapbetafilnavn) with open(oldwmapbetafilnavn) as w: rx=[] bx=[] ux=[] for line in w: data = line.split() rx.append(float(data[0])-0.1) #plot the point next to the first bx.append(float(data[1])) ux.append(float(data[2])) plt.errorbar(rx,bx, yerr=ux, marker='x', label=r"9-yr {\it WMAP} + {\it Planck} 2018",linewidth=0.5, color='red',markersize=3,capsize=2,ls='none', alpha=0.4) w.close() #plot the horisontal line #planckbetafilnavn="/mn/stornext/d5/unnif/sindex_bp/wmap23_pl2018/ut_sample_betas_invvar.txt" #with open(planckbetafilnavn) as planckfil: # planckbeta=float(planckfil.read()) #plt.axhline(y=planckbeta, color='red',linestyle='--',linewidth=0.5) #-2.99 is the old inv var weigh mean #planckfil.close() # this is the planc npipe result oldwmapbetafilnavn="/mn/stornext/d5/unnif/sindex_bp/wmap23_plnpipe/ut_region_scatter_alphamean.txt" with open(oldwmapbetafilnavn) as w: rx=[] bx=[] ux=[] for line in w: data = line.split() rx.append(float(data[0])) #plot the point next to the first two bx.append(float(data[1])) ux.append(float(data[2])) plt.errorbar(rx,bx, yerr=ux, marker='x', label=r"9-yr {\it WMAP} + {\it Planck} DR4",linewidth=0.5,markersize=3,capsize=2,ls='none', color='blue', alpha=0.4) #plot the horisontal line #planckbetafilnavn="/mn/stornext/d5/unnif/sindex_bp/wmap23_plnpipe/ut_sample_betas_invvar.txt" #with open(planckbetafilnavn) as planckfil: # planckbeta=float(planckfil.read()) #plt.axhline(y=planckbeta, color='blue',linestyle='--',linewidth=0.5) #-2.99 is the old inv var weigh mean #planckfil.close() # PLOT THE LFI LIKELIHOOD VALUE with uncertainty #https://www.aanda.org/articles/aa/full_html/2020/09/aa36386-19/F22.html #alpha=0.058 +- 0.004 is the template scaling factor from 28.6 to 70.0 GHz, which can be translated into a spectral index #That's roughly the band you want: beta_mean = -3.26, beta_low = -3.34, beta_high = -3.18: #-3.25796758275099 #gnuplot> print log((0.058-0.004) / 1.1369 * 1.0214) / log(70.795/28.587) #-3.33676787269876 #gnuplot> print log((0.058+0.004) / 1.1369 * 1.0214) / log(70.795/28.587) #-3.18442467987845 # y_est=-3.25796758275099 y_min=-3.33676787269876 y_max=-3.18442467987845 #plt.axhline(y=y_est, linestyle='-', linewidth=0.8) plt.plot((0.8,12.2),(y_est,y_est), label=r"{\it Planck} 2018 likelihood analysis", linestyle='-', linewidth=0.8) plt.fill_between(np.linspace(0.8,12.2,11), y_min, y_max, alpha=0.3) #################################################### plt.ylim(-6.1,-1.7) #plt.ylim(-4.,-2.2) #plt.yticks([-4.,-3.8, -3.6, -3.4, -3.2, -3.0, -2.8, -2.6, -2.4, -2.2],[-4.0,-3.8, -3.6, -3.4, -3.2, -3.0, -2.8, -2.6, -2.4, -2.2]) meanfilnavn=dire+filpre+'_sample_betas_mean.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=[] alphabetas_mat=[] alphauncerts_mat=[] #the loop over samples (cosmoglobe) #for samp in range(11,25): #one range for samp in it.chain(range(11,66), range(111,196)): #two split ranges print('reg='+str(reg)+' samp='+str(samp) +' '+ str(samp).zfill(3)) #read the spectral index in region #datadir=dire+'/../s0'+str(samp)+'/'+pre+str(reg)+'/' datadir=dire+'/../s'+str(samp).zfill(3)+'/'+pre+str(reg)+'/' #have the samp with 3 digits 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=[] b=[] with open(alphafilnavn) as f: for line in f: data = line.split() b.append(float(data[1])) u.append(float(data[2])) uncert=min(u) print(uncert) #could write out the uncert to a file here uncerts.append(uncert) #write the sample alpha values to a array of arrays alphabetas_mat.append(b) alphauncerts_mat.append(u) #end of sample loop ############################################ # stuff for the (main) sample big alpha and scatter plot ############################################ #ut_sample_scatter_alpha_spectral_index.txt nalphabetas_mat=np.array(alphabetas_mat) nalphauncerts_mat=np.array(alphauncerts_mat) #print(nalphabetas_mat) #print(nalphabetas_mat.shape) #shape har formen (#samples,#alphavalues) # creating files samplealphafilnavn=dire+pre+str(reg)+'/'+filpre+'_sample_scatter_alpha_spectral_index.txt' samplebetavaluefilnavn=dire+pre+str(reg)+'/'+filpre+'_sample_value_beta_scatter.txt' sampleuncertvaluefilnavn=dire+pre+str(reg)+'/'+filpre+'_sample_value_uncert_scatter.txt' print(samplealphafilnavn) reg_alphafile = open(samplealphafilnavn, 'w') # looping over the 18 alpha values for a in range(0,18): print(a) reg_alpha_beta=get_invvarweighted(nalphabetas_mat[:,a],nalphauncerts_mat[:,a]) reg_alpha_uncert=min(nalphauncerts_mat[:,a]) +stat.stdev(nalphabetas_mat[:,a]) print(reg_alpha_beta) print(reg_alpha_uncert) print(min(nalphauncerts_mat[:,a]) ,stat.stdev(nalphabetas_mat[:,a])) #writeto file reg_alphafile.write(str(alphanum[a])+'\t'+str(reg_alpha_beta)+'\t'+str(reg_alpha_uncert)+'\n') reg_alphafile.close() ############################ # stuff for the region plot ############################ print(" region plot values") s_betas=np.array(betas) s_uncerts=np.array(uncerts) final_beta=sum(s_betas/s_uncerts/s_uncerts)/sum(1./s_uncerts/s_uncerts) 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 straight mean beta='+str(stat.mean(s_betas))+str(min(s_uncerts))) print('The inverse weighted beta='+str(final_beta)+'+-'+str(final_uncert)) #print('howmanyregions='+str(howmanyregions)) regmeanfile.write(str(reg)+'\t'+str(mean_beta)+'\t'+str(min(s_uncerts))+'\n') reginvvarfile.write(str(reg)+'\t'+str(final_beta)+'\t'+str(final_uncert)+'\n') #write beta to ut_sample_value_beta_scatter.txt sample_betafile=open(samplebetavaluefilnavn, 'w') sample_betafile.write(str(final_beta)) sample_betafile.close() #write uncert to ut_sample_value_uncert_scatter.txt sample_uncertfile=open(sampleuncertvaluefilnavn, 'w') sample_uncertfile.write(str(final_uncert)) sample_uncertfile.close() #add them to a region array, for plotting purposes reg_betas.append(final_beta) reg_uncerts.append(final_uncert) #plot them if reg==1: #if case=='coswmap33': # lab=r"Cosmoglobe 23-33 GHz {\it WMAP}" #else: # lab=r"Cosmoglobe 23-30 GHz" lab=r"Cosmoglobe" else: lab="" plt.errorbar(reg+0.2,final_beta, yerr=final_uncert, marker='x', label=lab,linewidth=0.5, color='black',markersize=3,capsize=2,ls='none') regmeanfile.close() reginvvarfile.close() print("LOOP OVER REGIONS FINISHED") print("All beta values:") print(reg_betas) print("All uncert values:") print(reg_uncerts) print("** FINISHED LOOPING OVER REGIONS **") #calculating the mean of cosmoglobe betas c_betas=np.array(reg_betas) c_uncerts=np.array(reg_uncerts) #cos_invvar_beta=sum(c_betas/c_uncerts/c_uncerts)/sum(1./c_uncerts/c_uncerts) cos_invvar_beta=get_invvarweighted(c_betas,c_uncerts) cos_mean_beta=stat.mean(c_betas) cos_min_uncert=min(c_uncerts) print('The straight mean beta='+str(cos_mean_beta)+' +-'+str(cos_min_uncert)) print('The inverse weighted beta='+str(cos_invvar_beta)+' +-'+str(cos_min_uncert)) #plt.axhline(y=cos_invvar_beta, color='black',linestyle='--',linewidth=0.5) # general plot parameters plt.xlim(0.2,24.8) plt.xticks([1,5, 10, 15, 20, 24]) #ax.set_xticks([2,3,4], minor=True) ax.minorticks_on() plt.ylabel(r"Spectral index, $\beta_\mathrm{s}$", fontsize=10);plt.xlabel(r"Region number", fontsize=10); #leg = plt.legend(frameon=True, loc='lower right', prop={'size':7}) leg = plt.legend(frameon=True, loc='upper left', prop={'size':7}) # remove box around legend leg.get_frame().set_edgecolor("white") leg.get_frame().set_alpha(.0) plt.show() plt.savefig("region_beta_cosmoglobe_vs_wmap.png", bbox_inches='tight', bbox_extra_artists=[],pad_inches=0.03) #ut_sample_scatter_alpha_spectral_index.txt