In [1]:
import cmbfittest as fit
In [2]:
from scipy.ndimage.filters import gaussian_filter1d as gaussian_filter
from numpy.random import multivariate_normal as normal
from numpy.random import normal as N
from numpy import dot, sqrt, zeros, array, arange, save, mean, std
from numpy.linalg import inv, eigh
from scipy.linalg import sqrtm
from matplotlib.colors import LogNorm
In [3]:
### Instancing models
primordial = fit.model.model('primordial', ellmax = 5000, force_calculation = False)
lcdm = fit.model.model('lcdm', ellmax = 5000, force_calculation = False)
primordial.get_dcldphi()
lcdm.get_dcldphi()
primordial.get_covariance()
lcdm.get_covariance()
primordial.save_dcldphi()
lcdm.save_dcldphi()
In [4]:
#lcdm.dcldphi = gaussian_filter(lcdm.dcldphi, sigma = 20, axis = 0)
primordial.dcldphi= gaussian_filter(primordial.dcldphi, sigma = 20, axis = 0)
In [5]:
figure(figsize = (11,8.5), dpi = 400)
plot(primordial.dcldphi[:3000,:3000], linewidth = 2);
xlabel('multipole $l$', fontsize = 24)
ylabel(r'$\delta C_l \left[\mu {\rm k}^2\right]$', fontsize = 24)
title('Primordial Perturbation Responses', fontsize = 24)
savefig('figures/responses.png')
In [6]:
### SPT
spt = fit.data.data('spt', ellmax = 5000)
spt.get_covariance()
spt.get_bandpowers()
spt.get_window_functions()
errorbar(arange(47), spt.observation, yerr = sqrt(diag(spt.covariance)))
figure()
errorbar(arange(35,47), spt.observation[35:47], yerr = sqrt(diag(spt.covariance))[35:47])
Out[6]:
In [7]:
### WMAP
wmap = fit.data.data('wmap', ellmax = 5000)
wmap.get_covariance()
wmap.get_bandpowers()
wmap.get_window_functions()
errorbar(wmap.lslice, wmap.observation, yerr = sqrt(diag(wmap.covariance)))
Out[7]:
In [8]:
fisher_lcdm = fit.model.model('lcdm', ellmax = 5000)
parameter_prior_cov = inv(
dot(
dot(lcdm.dcldphi.T, wmap.windows.T),
dot(
inv(wmap.covariance),
dot(wmap.windows,lcdm.dcldphi)
)
)
+ dot(
dot(lcdm.dcldphi.T, spt.windows.T),
dot(
inv(spt.covariance),
dot(spt.windows,lcdm.dcldphi)
)
)
)
fisher_lcdm.dcldphi = dot(lcdm.dcldphi, parameter_prior_cov)
fisher_lcdm.get_covariance()
imshow(fisher_lcdm.covariance)
colorbar()
Out[8]:
In [9]:
#imshow(
# abs(wmap.covariance),
# norm = LogNorm(vmin = 10, vmax = 1e4)
# )
#colorbar()
#figure()
#imshow(abs((spt.covariance)),
# norm = LogNorm(vmin = 10, vmax=1e4)
# )
#colorbar()
#figure()
#plot(sqrt(diag(wmap.covariance)))
#figure()
#plot(sqrt(diag(spt.covariance)))
In [10]:
lda=fit.scheme.compression_scheme(
'LDA',
dcldphi = primordial.dcldphi,
noise_cov = [spt.covariance,wmap.covariance],
windows = [spt.windows, wmap.windows],
supermodel_cov = primordial.covariance,
kappa = 1e6,
SM_cov = fisher_lcdm.covariance,
nummodes = 20,
ellmax = 5000
)
In [11]:
lda.get_compression_modes()
In [12]:
from matplotlib.pyplot import cm
spt_ells = arange(650, 3000, 50)
nummodes = 4
figure(figsize = (11,8.5), dpi = 400)
colormap = cm.RdYlBu
#colors = [colormap(i) for i in np.linspace(0, 0.9, nummodes)]
colors = ['b', 'g', 'r', 'k']
for i in arange(nummodes):
plot(spt_ells, lda.modes[:47,i]/50, color = colors[i], linestyle='dashed', linewidth = 3);
plot(wmap.lslice, lda.modes[47:,i], color = colors[i], label = 'Mode %i'%(i+1), linewidth = 3);
legend(loc = 2)
xlabel('Multipole $l$', fontsize = 24)
ylabel('Filter Function $h_i(l)$', fontsize = 24)
savefig('figures/modes.png')
figure(figsize = (11,8.5), dpi = 400)
colors = ['b', 'g', 'r', 'k']
for i in arange(nummodes):
semilogx(spt_ells, lda.modes[:47,i]/50, color = colors[i], linestyle='dashed', linewidth = 3);
semilogx(wmap.lslice, lda.modes[47:,i], color = colors[i], label = 'Mode %i'%(i+1), linewidth = 3);
legend(loc = 2)
xlabel('Multipole $l$', fontsize = 24)
ylabel('Filter Function $h_i(l)$', fontsize = 24)
savefig('figures/modes_logscale.png')
In [13]:
print lcdm.dcldphi.T.shape
print spt.windows.shape
In [14]:
lda.modes.shape
for i in arange(6):
windowed_lcdm = []
windowed_lcdm.append(dot(lcdm.dcldphi.T[i,:], spt.windows.T))
windowed_lcdm.append(dot(lcdm.dcldphi.T[i,:], wmap.windows.T))
windowed_lcdm = [item for sublist in windowed_lcdm for item in sublist]
windowed_lcdm = array(windowed_lcdm)
print dot(real(lda.modes.T), windowed_lcdm)/sqrt(dot(windowed_lcdm, windowed_lcdm)* dot(real(lda.modes).T, real(lda.modes)));
#print dot(real(windowed_lcdm), windowed_lcdm)/sqrt(dot(windowed_lcdm, windowed_lcdm)* dot(windowed_lcdm, windowed_lcdm));
In [16]:
n = 1
chisqlist = zeros(n)
maxchisqlist = zeros(n)
globalz = zeros(n)
from numpy.random import normal as N
def randN(scale = 0):
if scale ==0:
return 0
else:
return N(scale = scale)
def add_noise(spectrum, sigma, orthogonal_modes):
return spectrum + dot(
orthogonal_modes,
array([randN(scale = sigma[j]) for j in arange(spectrum.size)])
)
for ii in arange(n):
cloud = fit.sim.gauss_approx(
observation = [spt.observation, wmap.observation],
covariance = [spt.covariance, wmap.covariance],
windows = [spt.windows, wmap.windows],
deltaCl = lcdm.dcldphi,
ellmax = 5000
)
cloud.get_bestfit()
#cloud.best_fit = spt_bestfit[:5001]
#print cloud.best_fit.shape
#plot(cloud.best_fit)
#variance, ortho_modes = eigh(spt.covariance[:3000, :3000])
#sigma = sqrt(variance)
#noisydraws = array(map(lambda x:add_noise(x, sigma, orthogonal_modes), draws.tolist()))
In [17]:
numsamples = 1000
draws = array([cloud.sample() for i in arange(numsamples)])
In [18]:
from numpy.random import normal
### Make fake WMAP experiments, using our biased posterior estimate
fake_wmap = dot(draws, wmap.windows.T)
##Diagonalize the wmap covariance matrix to add noise to each bin
variances, vecs = eigh(wmap.covariance)
for i, var in enumerate(variances):
for j in arange(numsamples):
fake_wmap[j,:] += normal(loc=0, scale = sqrt(var)) * vecs[:,i]
In [19]:
### Make fake SPT experiments, using our biased posterior estimate
fake_spt = dot(draws, spt.windows.T)
##Diagonalize the covariance matrix to add noise to each bin
variances, vecs = eigh(spt.covariance)
for i, var in enumerate(variances):
for j in arange(numsamples):
fake_spt[j,:] += normal(loc=0, scale = sqrt(var)) * vecs[:,i]
#print vecs[:,i].shape
In [20]:
plot((wmap.observation-dot(cloud.best_fit,wmap.windows.T))[:600], alpha = 1, color = 'black');
plot(fake_wmap[::100,:600].T, alpha = 0.05, color = 'green');
plot(draws[::10,:600].T, alpha = 0.01, color = 'blue');
xlim(0,600);
In [21]:
rc('text', usetex=True)
ax2 = gca()
plot(spt.observation - dot(cloud.best_fit,spt.windows.T), alpha = 1, color = 'black', linewidth = 2, label = 'True bandpowers' );
p1 = plot(fake_spt[1,:], alpha = 1, color = 'green', label = 'simmed bandpowers');
p2 = plot(dot(draws[1,:], spt.windows.T).T, alpha = 1, color = 'blue', label = 'simmed skies');
handles, labels = ax2.get_legend_handles_labels()
figure(figsize = (11,8.5), dpi = 400)
ax = gca()
plot(spt_ells, (spt.observation - dot(cloud.best_fit,spt.windows.T)), alpha = 1, color = 'black', linewidth = 4);
#plot(wmap.lslice[:650], (wmap.observation-dot(cloud.best_fit,wmap.windows.T))[:650], alpha = 1, color = 'black', linewidth = 1);
for i in arange(numsamples/10):
plot(spt_ells, fake_spt[10*i,:], alpha = 0.02, color = 'green');
plot(spt_ells, dot(draws[10*i,:], spt.windows.T), alpha = 0.02, color = 'blue');
#plot(wmap.lslice[:650], fake_wmap[10*i,:650].T, alpha = 0.01, color = 'green');
#plot(wmap.lslice[:650], dot(draws[10*i,:].T, wmap.windows.T)[:650], alpha = 0.01, color = 'blue');
legend(handles, labels);
xlim(650,3000)
ylim(-200, 300)
xlabel('Multipole $l$', fontsize = 24)
ylabel(r'$C_l - \langle C_l \rangle$[\mu {\rm K}]^2$', fontsize = 24)
savefig('figures/simmed_experiments.png')
In [22]:
labels
Out[22]:
In [23]:
plot(draws[::1000,:].T);
In [24]:
binned_truth = dot(cloud.best_fit, spt.windows.T)
#plot(dot(cloud.best_fit[:5001] - spt_bestfit[:5001], spt.windows.T), color = 'red', linewidth = 2)
plot(binned_truth - spt.observation, color = 'black', linewidth = 2)
plot(fake_spt[::50,:].T, alpha = .2);
In [25]:
### Concatenate the spt and wmap experiments and fake data:
observation = []
observation.append((spt.observation - dot(cloud.best_fit, spt.windows.T)).tolist())
observation.append((wmap.observation- dot(cloud.best_fit, wmap.windows.T)).tolist())
observation = [item for sublist in observation for item in sublist]
In [26]:
fake_observations = zeros([numsamples, len(observation)])
no_noise_observations = zeros([numsamples, len(observation)])
for i in arange(numsamples):
obs = []
obs.append(fake_spt[i,:])
obs.append(fake_wmap[i,:])
obs = [item for sublist in obs for item in sublist]
fake_observations[i,:] = obs
obs = []
obs.append(dot(draws[i,:], spt.windows.T))
obs.append(dot(draws[i,:], wmap.windows.T))
obs = [item for sublist in obs for item in sublist]
no_noise_observations[i,:] = obs
In [27]:
measured_dof = dot(observation, real(lda.modes))
cloud_dof = dot(fake_observations, real(lda.modes))
lcdm_dof = dot(no_noise_observations, real(lda.modes))
In [28]:
#reload(fit)
#globalchi2, global_zscore = fit.test.global_statistic(cloud_dof, measured_dof)
#maxchi2, max_location, max_zscore = fit.test.max_statistic(cloud_dof, measured_dof)
In [29]:
def distance_modulus(cloud_amps, observed_amps):
covariance = cov(cloud_amps.T)
means = mean(cloud_amps, axis = 0)
try:
distance_modulus = dot((observed_amps - means).T , dot(inv(covariance), (observed_amps - means)))
except:
distance_modulus = (observed_amps - means)**2/covariance
return distance_modulus
In [30]:
def global_statistic(cloud_amps, observed_amps):
numsamples, numparams = cloud_amps.shape
cloud_distances = array([distance_modulus(cloud_amps, cloud_amps[i,:]) for i in arange(numsamples)])
observed_distance = distance_modulus(cloud_amps, observed_amps)
expectation = mean(cloud_distances)
sigma = std(cloud_distances)
#chi2 = (observed_distance-expectation)**2/sigma**2
#expected_chi2 = subtract(cloud_distances, expectation)**2/sigma**2
#zscore = (chi2-mean(expected_chi2))/std(expected_chi2)
zscore = (observed_distance - expectation)/sigma
return observed_distance, zscore, expectation
In [31]:
def max_statistic(cloud_amps, observed_amps):
zscore = -1.0e-30
numsamples, numparams = cloud_amps.shape
for i in arange(1,numparams):
cloud_distances = array([distance_modulus(cloud_amps[:,:i], cloud_amps[j, :i]) for j in arange(numsamples)])
observed_distance = distance_modulus(cloud_amps[:,:i], observed_amps[:i])
#newchi2 = (observed_distance - expectation)**2/sigma**2
expectation = mean(cloud_distances)
sigma = std(cloud_distances)
#expected_chi2 = subtract(cloud_distances , expectation)**2/sigma**2
newscore = (observed_distance - expectation)/sigma
if abs(newscore) > abs(zscore):
chi2 = observed_distance
max_location = i
zscore = newscore
return observed_distance, max_location, zscore
In [32]:
x = global_statistic(cloud_dof, measured_dof)
y = max_statistic(cloud_dof, measured_dof)
print x, y
In [33]:
from numpy import std
numdof = 10
figure(figsize = (11,8.5), dpi = 400)
subplots_adjust(wspace = 0, hspace = 0)
for i in arange(numdof):
for j in arange(i+1):
if j == i:
ax = subplot(numdof, numdof, i + j*numdof + 1)
stdev = std(cloud_dof[:,i])
expected = mean(cloud_dof[:,i])
bins = linspace(-5 * stdev+expected, 5*stdev+expected, 50)
ax.hist(cloud_dof[:,i], bins, color = 'k')
ax.plot([measured_dof[i], measured_dof[i]], [0, sqrt(numsamples)], color = 'red', linewidth = 2)
ax.autoscale_view('tight')
setp( ax.get_yticklabels(), visible=False)
setp( ax.get_xticklabels(), visible=False)
# for tick in ax.xaxis.get_major_ticks():
# tick.label.set_fontsize(12)
# tick.label.set_rotation('vertical')
else:
ax=subplot(numdof,numdof, i+j*numdof+1)
ax.scatter(cloud_dof[:,i], cloud_dof[:,j], s = 1, alpha = .2, linewidths = 0, color = 'blue')
ax.scatter(lcdm_dof[:,i], lcdm_dof[:,j], s = 1, alpha = .2, linewidths = 0, color = 'k')
ax.scatter(measured_dof[i], measured_dof[j], s = 10, color = 'red')
setp( ax.get_xticklabels(), visible=False)
setp( ax.get_yticklabels(), visible=False)
ax.autoscale_view('tight')
savefig('figures/triangle_plot.png')
In [31]:
plot(lda.modes[:47,[6,16]])
figure()
semilogx(lda.modes[47:,[6,16]])
Out[31]:
In [32]:
spt_bestfit = np.load('spt_bestfit.npy')
In [33]:
plot(cloud.best_fit[:3000]- spt_bestfit[:3000])
Out[33]:
In [34]:
n = 1
chisqlist = zeros(n)
maxchisqlist = zeros(n)
globalz = zeros(n)
from numpy.random import normal as N
def randN(scale = 0):
if scale ==0:
return 0
else:
return N(scale = scale)
def add_noise(spectrum, sigma, orthogonal_modes):
return spectrum + dot(
orthogonal_modes,
array([randN(scale = sigma[j]) for j in arange(spectrum.size)])
)
for ii in arange(n):
cloud = fit.sim.gauss_approx(
observation = [spt.observation, wmap.observation],
covariance = [spt.covariance, wmap.covariance],
windows = [spt.windows, wmap.windows],
deltaCl = lcdm.dcldphi,
ellmax = 5000
)
#cloud.get_bestfit()
cloud.best_fit = spt_bestfit[:5001]
#print cloud.best_fit.shape
#plot(cloud.best_fit)
#variance, ortho_modes = eigh(spt.covariance[:3000, :3000])
#sigma = sqrt(variance)
#noisydraws = array(map(lambda x:add_noise(x, sigma, orthogonal_modes), draws.tolist()))
In [35]:
draws = array([cloud.sample() for i in arange(10000)])
In [36]:
from numpy.random import normal
### Make fake WMAP experiments, using our biased posterior estimate
fake_wmap = dot(draws, wmap.windows.T)
##Diagonalize the wmap covariance matrix to add noise to each bin
variances, vecs = eigh(wmap.covariance)
for i, var in enumerate(variances):
for j in arange(10000):
fake_wmap[j,:] += normal(loc=0, scale = sqrt(var)) * vecs[:,i]
In [37]:
### Make fake SPT experiments, using our biased posterior estimate
fake_spt = dot(draws, spt.windows.T)
##Diagonalize the covariance matrix to add noise to each bin
variances, vecs = eigh(spt.covariance)
for i, var in enumerate(variances):
for j in arange(10000):
fake_spt[j,:] += normal(loc=0, scale = sqrt(var)) * vecs[:,i]
In [38]:
plot((wmap.observation-dot(cloud.best_fit,wmap.windows.T))[:600], alpha = 1, color = 'black');
plot(draws[::10,:600].T, alpha = 0.01, color = 'blue');
plot(fake_wmap[::100,:600].T, alpha = 0.05, color = 'green');
xlim(0,600);
In [39]:
plot(spt.observation-dot(cloud.best_fit, spt.windows.T), alpha = 1, color = 'black', linewidth = 2);
plot(dot(draws[::10,:], spt.windows.T).T, alpha = 0.01, color = 'blue');
plot(fake_spt[::100,:].T, alpha = 0.1, color = 'green');
In [40]:
plot(draws[::1000,:3000].T);
In [41]:
### Concatenate the spt and wmap experiments and fake data:
observation = []
observation.append(spt.observation.tolist() - dot(cloud.best_fit, spt.windows.T.tolist()))
observation.append(wmap.observation.tolist() - dot(cloud.best_fit, wmap.windows.T.tolist()))
observation = [item for sublist in observation for item in sublist]
In [51]:
fake_observations = zeros([10000, len(observation)])
theories = zeros([10000, len(observation)])
for i in arange(10000):
obs = []
obs.append(fake_spt[i,:])# - dot(cloud.best_fit, spt.windows.T.tolist()))
obs.append(fake_wmap[i,:])# - dot(cloud.best_fit, wmap.windows.T.tolist()))
obs = [item for sublist in obs for item in sublist]
fake_observations[i,:] = obs
obs = []
obs.append(dot(draws[i,:], spt.windows.T))# - dot(cloud.best_fit, spt.windows.T.tolist()))
obs.append(dot(draws[i,:], wmap.windows.T))# - dot(cloud.best_fit, wmap.windows.T.tolist()))
obs = [item for sublist in obs for item in sublist]
theories[i,:] = obs
In [ ]:
measured_dof = dot(observation, real(lda.modes))
cloud_dof = dot(fake_observations, real(lda.modes))
theory_dof = dot(theories, real(lda.modes))
In [60]:
#junk, global_zscore = fit.test.global_statistic(cloud_dof, measured_dof)
#junk, max_location, max_zscore = fit.test.max_statistic(cloud_dof, measured_dof)
In [61]:
#print max_zscore
#print global_zscore
#print max_location
In [64]:
from numpy import std
numdof = 10
figure(figsize = (11,8.5), dpi = 400)
subplots_adjust(wspace = 0, hspace = 0)
for i in arange(numdof):
for j in arange(i+1):
if j == i:
ax = subplot(numdof, numdof, i + j*numdof + 1)
stdev = std(cloud_dof[:,i])
expected = mean(cloud_dof[:,i])
bins = linspace(-5 * stdev+expected, 5*stdev+expected, 50)
ax.hist(cloud_dof[:,i], bins, color = 'blue')
ax.plot([measured_dof[i], measured_dof[i]], [0, 500], color = 'red', linewidth = 1)
ax.autoscale_view('tight')
setp( ax.get_yticklabels(), visible=False)
setp( ax.get_xticklabels(), visible=False)
# for tick in ax.xaxis.get_major_ticks():
# tick.label.set_fontsize(12)
# tick.label.set_rotation('vertical')
else:
ax=subplot(numdof,numdof, i+j*numdof+1)
ax.scatter(cloud_dof[:,i], cloud_dof[:,j], s = 1, alpha = .2, linewidths = 0, color = 'blue')
ax.scatter(theory_dof[i], theory_dof[j], s = 1, alpha = 1, linewidths = 0, color = 'k')
ax.scatter(measured_dof[i], measured_dof[j], s = 10, color = 'red')
setp( ax.get_xticklabels(), visible=False)
setp( ax.get_yticklabels(), visible=False)
ax.autoscale_view('tight')
In [ ]:
nrun = np.load('nrun_response.npy')
In [ ]:
nrun.shape
spt.windows.shape
lda.modes.shape
In [ ]:
nrun_bandpowers = dot(nrun, spt.windows.T[:3001,:])
nrun_amps = dot(nrun_bandpowers, lda.modes[:47,:])
In [ ]:
nrun_amps
In [ ]:
bandpowers_reconstructed = dot(nrun_amps, lda.modes[:47,:].T)
In [ ]:
plot(nrun_bandpowers)
plot(bandpowers_reconstructed)
In [ ]:
plot(lda.modes[:47,:]);
In [ ]:
dot(lda.modes[:,1], lda.modes[:,1])
In [ ]:
In [ ]:
In [37]:
std(lcdm_dof, axis = 0 )
Out[37]:
In [36]:
std(cloud_dof, axis = 0)
Out[36]:
In [ ]: