In [1]:
%load_ext autoreload
#%autoreload 2
%matplotlib inline
execfile ("_ImportScript.py")
In [2]:
import time
import matplotlib.pyplot as plt
In [3]:
beatbox.You.create_original_Universe()
#beatbox.You.initiate_simulated_universe(truncated_nmax=20)
In [4]:
beatbox.You.initiate_simulated_universe(truncated_nmax=5)
In [ ]:
beatbox.You.all_simulated_universes[0].show_CMB_T_map()
ordered_inds_largenmax = beatbox.You.all_simulated_universes[0].get_ordered_fn_indices()
In [5]:
if not os.path.isfile('../data/covCyy_lmax%d_lmin%d.txt' % (beatbox.Multiverse.truncated_lmax, beatbox.Multiverse.truncated_lmin)):
beatbox.You.read_Planck_samples()
beatbox.You.calculate_covariance_matrix(filename='lmax%d_lmin%d' % (beatbox.Multiverse.truncated_lmax, beatbox.Multiverse.truncated_lmin))
else:
beatbox.You.load_covariance_matrix(filename='covCyy_lmax%d_lmin%d.txt' % (beatbox.Multiverse.truncated_lmax, beatbox.Multiverse.truncated_lmin))
# Calculate the inverse of the a_y covariance matrix
beatbox.You.calculate_sdv_Cyy_inverse()
datamap = beatbox.You.all_simulated_universes[0].ay2ayreal_for_inference(beatbox.You.all_simulated_universes[0].ay)+beatbox.You.generate_one_realization_of_noise()
MOCK = 1
#execfile ("_ReconstructionScript.py")
In [6]:
beatbox.You.solve_for_3D_potential(datamap.T , A=None, print_alpha=0)
In [ ]:
# np.savetxt( "/Users/LaurencePeanuts/Dropbox/KleineBar_KlaineTiger/MUSIC/2018/Acov_inv.txt", beatbox.You.inv_A)
# np.savetxt( "/Users/LaurencePeanuts/Dropbox/KleineBar_KlaineTiger/MUSIC/2018/Yash1.txt", beatbox.You.Yash1)
# np.savetxt( "/Users/LaurencePeanuts/Dropbox/KleineBar_KlaineTiger/MUSIC/2018/Yash2.txt", beatbox.You.Yash2)
In [ ]:
We = beatbox.Universe()
We.fn = beatbox.You.reconstrunct_fn
#We.fn = beatbox.You.all_simulated_universes[-1].fn * 0.5
We.transform_3D_potential_into_alm(usedefault=1, fn=1)
We.show_CMB_T_map(title = "Best Fit Model", from_perspective_of="observer", max=100)
We.rearrange_fn_from_vector_to_grid()
We.evaluate_potential_given_fourier_coefficients()
ordered_inds_smallnmax = We.get_ordered_fn_indices()
In [ ]:
beatbox.You.all_simulated_universes[0].show_CMB_T_map()
In [ ]:
In [ ]:
beatbox.You.inv_A = beatbox.You.calculate_pure_A_inverse()
In [ ]:
# beatbox.You.load_A_matrix('nmax14', inv = 'y')
plt.plot(We.fn[ordered_inds_smallnmax ], '*')
plt.errorbar(range(len(We.fn)), We.fn[ordered_inds_smallnmax ], yerr = np.sqrt(np.diag(beatbox.You.inv_A))[ordered_inds_smallnmax ])
plt.plot(beatbox.You.all_simulated_universes[0].fn[ordered_inds_largenmax], 'x')
plt.axis([0, 256, -10, 10])
In [ ]:
frac = np.abs(np.divide((We.fn - beatbox.You.all_simulated_universes[0].fn), np.sqrt(np.diag(beatbox.You.inv_A))))
print frac
In [ ]:
print np.double(np.sum(frac<=1))/np.double(len(frac))
In [ ]:
numreal = 100
for i in range(numreal):
beatbox.You.initiate_simulated_universe()
In [ ]:
numCp = np.zeros(len(beatbox.You.all_simulated_universes[0].fn))
for j in range(len(beatbox.You.all_simulated_universes[0].fn)):
fns = [beatbox.You.all_simulated_universes[i].fn[j] for i in range(numreal)]
numCp[j] = 1./(np.std(fns)**2)
In [ ]:
print np.mean(numCp)
In [ ]:
beatbox.You.solve_for_3D_potential(datamap.T , A=None, print_alpha=0)
In [ ]:
plt.plot( np.divide(numCp , np.diag(beatbox.You.inv_Cf) ))
In [ ]:
ind = np.where(beatbox.Universe.kfilter>0)
print ind[]
In [ ]:
# Check calculated k covariance from samples, against prior covariance matrix
numCp = np.zeros(len(beatbox.You.all_simulated_universes[0].fn))
for j in range(len(beatbox.You.all_simulated_universes[0].fn)):
fns = [beatbox.You.all_simulated_universes[i].fn[j] for i in range(numreal)]
numCp[j] = np.std(fns) #1./(np.std(fns)**2)
plt.semilogy(numCp[ordered_inds_smallnmax])
plt.semilogy(1/np.sqrt(np.diag(beatbox.You.inv_Cf)[ordered_inds_smallnmax]),'--')
In [7]:
#Calculate coverage prob for 1sig
numreal = 100
cov_frac = np.zeros((numreal,1))
for i in range(numreal):
print i
beatbox.You.initiate_simulated_universe(truncated_nmax=10)
if i == 0:
ordered_inds_largenmax = beatbox.You.all_simulated_universes[i].get_ordered_fn_indices()
datamap = beatbox.You.all_simulated_universes[i].ay2ayreal_for_inference(beatbox.You.all_simulated_universes[i].ay)+beatbox.You.generate_one_realization_of_noise()
beatbox.You.solve_for_3D_potential(datamap.T , A=1, print_alpha=0)
We=beatbox.Universe()
beatbox.You.all_data_universes = np.append(beatbox.You.all_data_universes,We)
beatbox.You.all_data_universes[i].fn = beatbox.You.reconstrunct_fn
num_fn = len(beatbox.You.reconstrunct_fn)
if i == 0:
ordered_inds_smallnmax = beatbox.You.all_data_universes[i].get_ordered_fn_indices()
abs_diff = np.abs(beatbox.You.all_data_universes[i].fn[ordered_inds_smallnmax].reshape(-1,1) - (beatbox.You.all_simulated_universes[i].fn[ordered_inds_largenmax[:num_fn]]).reshape(-1,1))
err_sig = np.sqrt(np.diag(beatbox.You.inv_A)[ordered_inds_smallnmax])
diff_sigma = np.divide(abs_diff , err_sig ).reshape(-1,1)
cov_frac[i,0] = np.double(np.sum(diff_sigma<=1)) / np.double(len(diff_sigma))
print diff_sigma
In [8]:
print np.mean(cov_frac)
In [ ]:
ind = np.where(beatbox.You.all_simulated_universes[0].kfilter>0)
k, theta, phi = beatbox.You.all_simulated_universes[0].k[ind], np.arctan2(beatbox.You.all_simulated_universes[0].ky[ind],beatbox.You.all_simulated_universes[0].kx[ind]), np.arccos(beatbox.You.all_simulated_universes[0].kz[ind]/beatbox.You.all_simulated_universes[0].k[ind])
kvec_long = np.zeros(2*len(ind[1]))
kvec_long[:len(ind[1])] = k
kvec_long[len(ind[1]):] = k
kvec = np.zeros(len(ind[1]))
kvec[:len(ind[1])/2] = kvec_long[:len(ind[1])/2]
kvec[len(ind[1])/2:] = kvec_long[len(ind[1]):3*len(ind[1])/2]
#ind_for_ordered_fn = np.argsort(kvec)
In [ ]:
print kvec[ordered_inds_largenmax[:50]].shape
In [ ]:
len(We.fn)
In [ ]: