In [1]:

    

import model


    

In [2]:

    

primordial = model.model('primordial')
lcdm = model.model('lcdm')


    

In [3]:

    

primordial.get_dcldphi()
lcdm.get_dcldphi()


    

In [4]:

    

primordial.dcldphi = primordial.dcldphi[:2501,:]
lcdm.dcldphi = lcdm.dcldphi[:2501,:]


    

In [5]:

    

from scipy.ndimage.filters import gaussian_filter1d as gaussian_filter


    

In [6]:

    

lcdm.dcldphi_smooth = gaussian_filter(lcdm.dcldphi, sigma = 20, axis = 0)
primordial.dcldphi_smooth = gaussian_filter(primordial.dcldphi, sigma = 20, axis = 0)


    

In [7]:

    

plot(primordial.dcldphi_smooth);
primordial.save_dcldphi()
lcdm.save_dcldphi() 
#.dcldphi / dot(lcdm.dcldphi, lcdm.dcldphi).reshape(-1,1)


    

In [8]:

    

row_sums = sqrt((lcdm.dcldphi**2).sum(axis = 0))
lcdm.dcldphi = lcdm.dcldphi /row_sums
plot(lcdm.dcldphi);


    

In [9]:

    

primordial.get_covariance()
lcdm.get_covariance()


    

In [10]:

    

imshow(lcdm.covariance)
colorbar()
import data


    

In [11]:

    

experiment = data.data('planck_like')


    

In [12]:

    

experiment.get_covariance()


    

In [13]:

    

import compression_scheme


    

In [14]:

    

lda = compression_scheme.compression_scheme(
                                            'LDA', 
                                            dcldphi = primordial.dcldphi_smooth, 
                                            noise_cov = experiment.covariance, 
                                            supermodel_cov = primordial.covariance 
                                            #kappa = 1.0e10, 
                                            #SM_cov = lcdm.covariance
                                            )


    

In [15]:

    

lda.get_compression_modes()


    

In [16]:

    

plot(lda.modes[5:,19]);


    

In [17]:

    

import simulate
from numpy.random import multivariate_normal as normal
#observation = experiment.fiducial_cl[:2501] + normal(np.zeros(2501), experiment.covariance[:2501,:2501])


    

In [18]:

    

clobs = loadtxt('../clobs.csv', delimiter = ',')
#plot(clobs[:,2:])
observation = experiment.fiducial_cl[:2501].copy()
observation += normal(np.zeros(2501), experiment.covariance[:2501,:2501])
observation2 = experiment.fiducial_cl[:2501].copy()
observation2[clobs[:,0].astype('int')] = clobs[:,4]
plot(observation)


    

    
Out[18]:





[<matplotlib.lines.Line2D at 0x3ba3890>]






    

In [19]:

    

cloud = simulate.gauss_approx(observation = observation, 
                              covariance = experiment.covariance, 
                              deltaCl = lcdm.dcldphi_smooth)


    

In [20]:

    

cloud.get_bestfit()
plot(cloud.best_fit - experiment.fiducial_cl[:2501])
ylim(-100,100)


    

    
Out[20]:





(-100, 100)






    

In [21]:

    

plot(cloud.best_fit)
plot(experiment.fiducial_cl[:2501], alpha = 1)
plot(cloud.observation, alpha = 0.3, color = 'k')


    

    
Out[21]:





[<matplotlib.lines.Line2D at 0x7d55ed0>]






    

In [22]:

    

x = array([cloud.sample() for i in arange(10000)])
x.shape


    

    
Out[22]:





(10000, 2501)

In [23]:

    

#plot(x.T, color = 'b', alpha = .3);
#plot(cloud.observation, color = 'k', alpha = 0.3);
#ylim(1,1e4);


    

In [24]:

    

#plot(array([x[i,:]/cloud.best_fit for i in arange(10000)]).T, color = 'b', alpha = 0.3);
#plot(cloud.observation/cloud.best_fit, color = 'k', alpha = 0.3);
#ylim(0,2);


    

In [25]:

    

#lda.modes.shape
#x.shape
sigma2 = diag(cloud.covariance)
sigma2[0] = 1
xwithnoise = array([x[j,:] + array([random.normal(scale = sqrt(sigma2[i])) for i in arange(0,2501)]) for j in arange(10000)])


    

    




-c:4: FutureWarning: Numpy has detected that you (may be) writing to an array returned
by numpy.diagonal or by selecting multiple fields in a record
array. This code will likely break in the next numpy release --
see numpy.diagonal or arrays.indexing reference docs for details.
The quick fix is to make an explicit copy (e.g., do
arr.diagonal().copy() or arr[['f0','f1']].copy()).

In [26]:

    

cloud.amps = dot(lda.modes.T, xwithnoise.T)
#cloud.amps.shape


    

In [27]:

    

observed_amp = dot(lda.modes.T, cloud.observation)
#observed_amp.shape


    

In [28]:

    

from scipy.linalg import sqrtm
means = mean(cloud.amps, axis = 1)
covariance = cov(cloud.amps)
covariance.shape


    

    
Out[28]:





(30, 30)

In [29]:

    

z = dot(sqrtm(covariance), observed_amp)
zcloud = dot(sqrtm((covariance)), cloud.amps)
zcloud.shape


    

    
Out[29]:





(30, 10000)

In [29]:

    




    

In [30]:

    

plot(cloud.amps, color = 'k', alpha = .01);
plot(observed_amp, color = 'b', alpha = 1);
#test = dot(x[:10,:], lda.modes);
#plot(abs(test.T), color = 'b');
#plot(abs(observed_amp));


    

In [46]:

    

distance = dot((observed_amp - means).T, dot(inv(covariance), (observed_amp - means).T))
test_distance = dot((test_amp - means).T, dot(inv(covariance), (test_amp - means).T))


    

    




---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-46-04ceb450bcb8> in <module>()
----> 1 distance = dot((observed_amp - means).T, dot(inv(covariance), (observed_amp - means).T))
      2 test_distance = dot((test_amp - means).T, dot(inv(covariance), (test_amp - means).T))

ValueError: operands could not be broadcast together with shapes (30) (100) 

In [40]:

    

from numpy.random import randint
test_amp = dot(lda.modes.T, observation2)
k = randint(0,10000)
figure(figsize = (100,100), dpi = 100)
subplots_adjust(wspace=0,
                    hspace = 0)
for i in arange(30):
    subplot(30,30, 31*i+1)
    hist, edges = histogram(cloud.amps[i,:], bins = 100)
    hist = hist/float(max(hist))
    window = ones(int(2))/float(2)
    xvals = convolve(edges, window, 'valid')
    plot(xvals, hist)
    plot([cloud.amps[i,k], cloud.amps[i,k]], [0,1], color = 'yellow', linewidth = 4)
    plot([observed_amp[i], observed_amp[i]], [0,1], color = 'red', linewidth = 4)
    gca().axes.get_xaxis().set_ticks([])
    gca().axes.get_yaxis().set_ticks([])
    xlabel('mode %i'%i)
    ylabel('mode %i'%i)
    
    for j in arange(i):
        subplot(30,30,30*j + i+1)
        scatter(cloud.amps[i,:], cloud.amps[j,:], s = 1, alpha = .1);
        scatter(observed_amp[i], observed_amp[j], color = 'red');
        scatter(cloud.amps[i,k], cloud.amps[j,k], color = 'yellow');
        #scatter(test_amp[i], test_amp[j], color = 'purple')
        #if ((30*j + i + 1) % 31 == 2):
        #    
        #    ylabel('mode %i'%j)
        gca().axes.get_xaxis().set_ticks([])
        gca().axes.get_yaxis().set_ticks([])
savefig('/home/btfollin/triangle.png')


    

In [41]:

    

plot(array([dot(lda.modes, cloud.amps[:,i]) - dot(lda.modes, observed_amp) for i in arange(10000)]).T, color = 'k', alpha = .01);
#plot(dot(lda.modes, observed_amp))
ylabel(r'$C_l^{cloud}  - C_l^{observed}$--with noise', fontsize = 20);
title('S_N projected power; no exotic perturbations');


    

In [45]:

    

semilogy(lda.eigspectrum[:30]/max(lda.eigspectrum))


    

    
Out[45]:





[<matplotlib.lines.Line2D at 0x1761d950>]






    

In [ ]: