gistfile1.txt

from __future__ import division
import numpy as np
from matplotlib import pyplot as plt
import cPickle as pickle
import scipy.linalg as la
import scipy

from commander.sphere.mmajor import lm_to_idx

nside = 16
npix = 12 * nside**2
lmax = 3 * nside - 1
nsh = (lmax + 1)**2

with file('simple_lowres_system.pickle') as f:
    Y, YhW, beam, S_diag, rms, mask = pickle.load(f)



P = Y * beam[None, :]# * np.sqrt(S_diag)[None, :]
P = P[mask == 1, :]
P = P[:, S_diag != 0]

Sinv = np.diag(1 / S_diag[S_diag != 0])

noise = (rms**2)[mask == 1]
N = np.diag(noise)

npix, nsh = P.shape

from joblib import Memory
mem = Memory(cachedir='joblib')

@mem.cache
def plot_A(N, Sinv, P):
    print 1
    A = np.vstack([
        np.hstack([N, P]),
        np.hstack([P.T, -Sinv])])
    print 2
    i = np.arange(A.shape[0])
    m = np.hstack([np.diag(N), np.diag(Sinv)])
    eig = la.eigvalsh(A, np.diag(m))
    print 3
    return A, eig


#P *= np.sqrt(S_diag)[S_diag != 0][None, :]
#P *= (1 / np.sqrt(noise))[:, None]

#A = plot_A(np.eye(npix), np.eye(nsh), P)

A, eig = plot_A(N, Sinv, P)
plt.semilogy(np.abs(sorted(eig)))


U = np.sqrt(1 / noise)[:, None] * P * np.sqrt(S_diag[S_diag != 0])
Q = np.dot(U.T, U)
i = np.arange(Q.shape[0])
Q[i, i] += 1
v = la.eigvalsh(Q, np.diag(Q[i, i]))
print v[-1] / v[0]
plt.semilogy(np.abs(sorted(v)))