In [ ]:
import os
import pickle
import cPickle as cpkl
import csv
import sys
import matplotlib as mpl
mpl.use('PS')
import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse
mpl.rcParams['mathtext.rm'] = 'serif'
plt.rcParams["mathtext.fontset"] = "dejavuserif"
mpl.rcParams['font.family'] = 'serif'
mpl.rcParams['font.serif'] = 'Times New Roman'
title = 18
label = 16
mpl.rcParams['text.usetex'] = False
mpl.rcParams['axes.titlesize'] = title
mpl.rcParams['axes.labelsize'] = label
mpl.pyplot.rcParams['xtick.labelsize'] = label
mpl.pyplot.rcParams['ytick.labelsize'] = label
mpl.rcParams['figure.subplot.left'] = 0.2
mpl.rcParams['figure.subplot.right'] = 0.9
mpl.rcParams['figure.subplot.bottom'] = 0.2
mpl.rcParams['figure.subplot.top'] = 0.9
mpl.rcParams['figure.subplot.wspace'] = 0.5
mpl.rcParams['figure.subplot.hspace'] = 0.5
mpl.rcParams['savefig.dpi'] = 250
mpl.rcParams['savefig.format'] = 'pdf'
mpl.rcParams['savefig.bbox'] = 'tight'
import numpy as np
import scipy.stats as sps
import chippr
import chippr.plot_utils as pu
%matplotlib inline
In [ ]:
import matplotlib.hatch
import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse, Polygon
house_path = Polygon(
[[-0.3, -0.4], [0.3, -0.4], [0.3, 0.1], [0., 0.4], [-0.3, 0.1]],
closed=True, fill=False).get_path()
class CustomHatch(mpl.hatch.Shapes):
"""
Custom hatches defined by a path drawn inside [-0.5, 0.5] square.
Identifier 'c'.
"""
filled = True
size = 1.0
path = house_path
def __init__(self, hatch, density):
self.num_rows = (hatch.count('c')) * density
self.shape_vertices = self.path.vertices
self.shape_codes = self.path.codes
mpl.hatch.Shapes.__init__(self, hatch, density)
mpl.hatch._hatch_types.append(CustomHatch)
In [ ]:
# print(np.shape(cl_input))
# n_tomobins = np.shape(cl_input)[1] - 1
# 1 column of redshifts, 4 tomographic bins, 350 redshift bins
# too many redshift bins, I'll set the truth to that binned down
In [ ]:
chippr_bins = 35
colors = ['k', pu.colors[2], pu.colors[-3], pu.colors[0]]#'#332288', '#117733', '#999933', '#FF6700']
In [ ]:
cl_input = np.genfromtxt('../results/CosmoLike/nz_histo.txt')
n_tomobins = np.shape(cl_input)[1] - 1
fine_dif = np.mean(cl_input.T[0])
cl_data = np.empty((chippr_bins, 5))
factor = len(cl_input.T[0]) / chippr_bins
In [ ]:
for i in range(chippr_bins):
cl_data[i][0] = cl_input[i * factor][0]
cl_data[i][1:] = np.sum(cl_input[i * factor:(i + 1) * factor, 1:], axis=0) / (factor)
np.savetxt('../results/CosmoLike/nz_format_test.txt', cl_data)
In [ ]:
# this is the truth!
cl_data = np.empty((chippr_bins, 5))
for i in range(chippr_bins):
cl_data[i][1:] = np.sum(cl_input[i * 10:(i + 1) * 10, 1:], axis=0)
cl_data[i][0] = cl_input[i * 10][0]
for i in range(n_tomobins):
plt.plot(cl_data.T[0], cl_data.T[i+1])
# np.savetxt('../results/CosmoLike/nz_format_test.txt', cl_data)
In [ ]:
bin_ends = np.genfromtxt('../results/CosmoLike/0single_lsst-in-thesis/data/metadata.txt')
bin_difs = bin_ends[1:] - bin_ends[:-1]
bin_mids = (bin_ends[1:] + bin_ends[:-1]) / 2.
In [ ]:
print(bin_ends)
In [ ]:
all_true_nzs = [bin_mids]
for i in range(n_tomobins):
true_zs = np.genfromtxt('../results/CosmoLike/' + str(i) + 'single_lsst-in-thesis/data/true_vals.txt').T[0]
true_nz = np.histogram(true_zs, bins=bin_ends)[0] / float(len(true_zs)) / bin_difs
all_true_nzs.append(true_nz)
all_true_nzs = np.array(all_true_nzs)
np.savetxt('../results/CosmoLike/true_nz.txt', all_true_nzs.T)
In [ ]:
for i in range(n_tomobins):
plt.plot(all_true_nzs[0], all_true_nzs[i+1])
In [ ]:
cl_input.T[0]
In [ ]:
f = plt.figure(figsize=(10, 5))
ax = f.add_subplot(1, 1, 1)
ax.set_title(r'')
ax.set_xlabel(r'$z$', fontsize=18)
ax.set_ylabel(r'$n(z)$', fontsize=18)
ax.plot([-1., -2], [0., 1.], color='k', lw=1., alpha=1., label=r'underlying tomographic $n(z)$')
ax.plot([-1., -2], [0., 1.], color='k', lw=2., alpha=0.5, linestyle='-', dashes=[3, 2], label=r'binned tomographic $n(z)$')
ax.plot([-1., -2], [0., 1.], color='k', lw=2., alpha=1., linestyle='-', dashes=[1, 1], label=r'drawn tomographic $n(z)$')
for i in range(n_tomobins):
ax.plot(cl_input.T[0], cl_input.T[i+1], color=colors[i], lw=1., alpha=1.)
norm_factor = np.sum(bin_difs * cl_data.T[i+1])
print(norm_factor)
pu.plot_step(ax, cl_data.T[0], cl_data.T[i+1] / norm_factor, c=colors[i], w=2., s='--', a=0.5, d=[(0,(3, 2))])
pu.plot_step(ax, bin_ends, all_true_nzs[i+1], c=colors[i], w=2., s='--', a=0.5, d=[(0,(1, 1))])
ax.set_xlim(-0.1, 3.6)
ax.legend(loc='upper right', fontsize=20)
f.savefig('cosmolike_inputs.png', bbox_inches='tight', pad_inches = 0, dpi=250)
In [ ]:
all_stk_nzs, all_mle_nzs, all_map_nzs = [bin_mids], [bin_mids], [bin_mids]
for i in range(n_tomobins):
with open('../results/' + str(i) + 'single_lsst-in-thesis/results/nz.p', 'rb') as test_file:
test_info = cpkl.load(test_file)
all_stk_nzs.append(np.exp(test_info['estimators']['log_stacked_nz']))
all_mle_nzs.append(np.exp(test_info['estimators']['log_mmle_nz']))
all_map_nzs.append(np.exp(test_info['estimators']['log_mmap_nz']))
all_stk_nzs = np.array(all_stk_nzs)
all_mle_nzs = np.array(all_mle_nzs)
all_map_nzs = np.array(all_map_nzs)
np.savetxt('../results/CosmoLike/stack_nz.txt', all_stk_nzs.T)
np.savetxt('../results/CosmoLike/mmle_nz.txt', all_mle_nzs.T)
np.savetxt('../results/CosmoLike/mmap_nz.txt', all_map_nzs.T)
In [ ]:
for i in range(n_tomobins):
plt.plot(all_true_nzs[0], all_true_nzs[i+1], c='k')
plt.plot(all_true_nzs[0], all_stk_nzs[i+1], c=pu.colors[0])
plt.plot(all_true_nzs[0], all_mle_nzs[i+1], c='b')
plt.plot(all_true_nzs[0], all_map_nzs[i+1], c=pu.colors[-3])
In [ ]:
testnames = ['true_nz', 'stack_nz', 'mmle_nz', 'mmap_nz']
In [ ]:
all_covariances = []
all_invcovariances = []
magfactors = []
for testname in testnames:
covariances = np.zeros((200, 200)) + sys.float_info.min
# with open(os.path.join('../results/CosmoLike', 'Cl_cov.nz'+testname+'.txt'), 'rb') as cosmolike_file:
# cosmolike_reader = csv.reader(cosmolike_file, delimiter=' ')
# cosmolike_reader.next()
# for row in cosmolike_reader:
# # covariance matrix is filled with positive values << floating point precision.
# # I'm going to inflate them and deflate before reporting results
# covariances[int(row[0])][int(row[1])] = float(row[2])
# covariances[int(row[1])][int(row[0])] = float(row[2])
covariance_table = np.genfromtxt('../results/CosmoLike/Cl_cov.'+testname+'.txt')
magfactor = np.asarray(np.min(covariance_table.T[-1]))
print(magfactor)
magfactors.append(magfactor)
for row in covariance_table:
# print(row)
covariances[int(row[0])][int(row[1])] = row[2]# / magfactor
covariances[int(row[1])][int(row[0])] = row[2]# / magfactor
# covariances = covariances - sys.float_info.min
# assert(np.all(covariances >= 0.))
# covariances += sys.float_info.min
# covariances = covariances / magfactor
assert(np.all(covariances > 0.))
all_covariances.append(covariances)
invcovariances = np.linalg.inv(covariances)
mask = np.where((invcovariances < sys.float_info.min))
invcovariances[mask] = sys.float_info.min
print(np.shape(invcovariances))
assert(np.all(invcovariances > 0.))
all_invcovariances.append(invcovariances)
In [ ]:
# # plt.imshow(np.log(covariances))
# # print(covariances)
# invcovariances = np.linalg.pinv(covariances)# / sys.float_info.epsilon) * sys.float_info.epsilon
# # print(invcovariances)
# # # invcovariances = np.linalg.pinv(1.e15 * covariances) * 1.e15
# # # print(invcovariances)
# # plt.imshow(np.log(invcovariances * sys.float_info.epsilon))
# # print(invcovariances * sys.float_info.epsilon)
In [ ]:
# Cl, dCldOm, dClds8, dCldns, dCldw0, dCldwa, dCldOb, dCldH0 = np.zeros(200), np.zeros(200), np.zeros(200), np.zeros(200), np.zeros(200), np.zeros(200), np.zeros(200)
all_deriv_info = []
for testname in testnames:
# deriv_info = np.zeros((200, 7)) + sys.float_info.min
# with open(os.path.join('../results/CosmoLike', 'Cl_derivs.nz'+testname+'.txt'), 'rb') as cosmolike_file:
# cosmolike_reader = csv.reader(cosmolike_file, delimiter=' ')
# cosmolike_reader.next()
# i = 0
# while i < 200:
# for row in cosmolike_reader:
# for j in range(7):
# deriv_info[i][j] = float(row[j])
# i += 1
deriv_info = np.genfromtxt('../results/CosmoLike/Cl_derivs.'+testname+'.txt')#, skip_header=0)
# print(np.shape(deriv_info))
deriv_info = np.vstack((np.zeros(7), deriv_info))
deriv_info = deriv_info.T
all_deriv_info.append(deriv_info)
In [ ]:
# # deriv_info = deriv_info.T
# dCldOm, dClds8, dCldns, dCldw0, dCldwa, dCldOb, dCldH0 = deriv_info[0], deriv_info[1], deriv_info[2], deriv_info[3], deriv_info[4], deriv_info[5], deriv_info[6]
In [ ]:
all_fisher = []
all_invfisher = []
for k in range(len(testnames)):
ls = np.arange(200)
fisher = np.eye(7)
for i in range(7):
for j in range(i+1):
fisher[i][j] = np.sum((2. * ls[1:] + 1.) * all_deriv_info[k][i][1:] *
(all_invcovariances[k][1:, 1:]) * all_deriv_info[k][j][1:])
fisher[j][i] = fisher[i][j]
all_fisher.append(fisher)
invfisher = np.linalg.pinv(fisher)#np.linalg.inv(fisher/sys.float_info.epsilon) * sys.float_info.epsilon
all_invfisher.append(invfisher)
In [ ]:
# # print(fisher/sys.float_info.epsilon)
# # plt.imshow(fisher/sys.float_info.epsilon)
# invfisher = np.linalg.inv(fisher/sys.float_info.epsilon) * sys.float_info.epsilon
# plt.imshow(np.log(all_fisher[0]))
# plt.colorbar()
# print(invfisher)
In [ ]:
def make_ellipse_params(fisher):
diag_elems = np.diag(fisher)
term1 = (diag_elems[:, np.newaxis] + diag_elems[np.newaxis, :]) / 2.
term2 = np.sqrt((diag_elems[:, np.newaxis] - diag_elems[np.newaxis, :])**2 / 4. + fisher * fisher.T)
a = np.sqrt(term1 + term2)
b = np.sqrt(term1 - term2)
t = np.arctan((fisher + fisher.T) / np.abs(diag_elems[:, np.newaxis] - diag_elems[np.newaxis, :])) / 2.
try:
assert(np.all(a >= 0.))
except:
print(a)
try:
assert(np.all(b >= 0.))
except:
print(b)
return(a, b, t)
In [ ]:
fish_params = []
for k in range(len(testnames)):
fisher = all_invfisher[k]
(a, b, t) = make_ellipse_params(fisher)
fish_params.append((a, b, t))
fish_params = np.array(fish_params)
Calculate the KLD between ellipses
In [ ]:
def kld(p, q):
term1 = np.trace(np.matmul(np.linalg.pinv(q), p))
# print(term1)
term2 = float(len(q))
term3 = np.log(np.linalg.det(q) / np.linalg.det(p))
# print(term3)
return (term1 - term2 + term3) / 2.
In [ ]:
def area(a, b):
return np.pi * a * b
In [ ]:
for k in range(1, len(testnames)):
print(testnames[k])
the_kld = kld(all_invfisher[0], all_invfisher[k])
print(the_kld)
In [ ]:
colors = ['k', pu.colors[0], pu.colors[2], pu.colors[4]]
styles = [[(0, (1, 0.001))], [(0, (2, 2))], [(0, (1, 2))], [(0, (2, 1))]]
names = ['true', 'stacked', 'CHIPPR', 'modes']
hatches = [None, '--', None, '||']
fills = [False, False, True, False]
alphas = [1., 1., 0.5, 1.]
keys = ['dCl/dOm', 'dCl/ds8', 'dCl/dns', 'dCl/dw0', 'dCl/dwa', 'dCl/dOb', 'dCl/dH0']
fancy = [r'$dC_{l}/d\Omega_{m}$', r'$dC_{l}/ds_{8}$', r'$dC_{l}/dn_{s}$', r'$dC_{l}/dw_{0}$', r'$dC_{l}/dw_{a}$', r'$dC_{l}/d\Omega_{b}$',r'$dC_{l}/dH_{0}$']
mini_inds = [0, 3, 4, -1]
mini_keys = [keys[mini_ind] for mini_ind in mini_inds]
mini_fancy = [fancy[mini_ind] for mini_ind in mini_inds]
# keys = [r'$\Omega_{m}$', r'$S_{8}$', r'$n_{s}$', r'$w_{0}$', r'$w_{a}$', r'$\Omega_{b}$', r'$H_{0}$']
In [ ]:
def mycorner(fish_params, keys, names, hatches, fills, alphas, labels=None, extra=''):
if labels==None:
labels = keys
ncol = len(keys)
fig = plt.figure(figsize=(ncol*(5), ncol*5))
ax = [[fig.add_subplot(ncol, ncol, ncol * i + j + 1) for j in range(i+1)] for i in range(ncol)]
to_keep = range(ncol)#[0, 1, 2, 4]
extrema = np.zeros(ncol)
for k in range(len(testnames)):
if k != 0:
the_kld = '\n KLD='+"{0:.3e}".format(kld(all_invfisher[0], all_invfisher[k]))
else:
the_kld = ''
axl = fig.add_subplot(ncol, ncol, ncol+ncol-1)
axl.text(0.1, 0.9 - 0.25 * k, names[k]+r' $n(z)$'+the_kld, fontsize=30, color=colors[k])
axl.set_xticks([])
axl.set_yticks([])
# fisher = fishers_info[k]
a, b, t = fish_params[k]
extrema = np.max(np.sqrt(a), axis=0) / 2.
for i in range(ncol):
lim_val_i = np.sqrt(all_invfisher[k][i][i])
for j in range(i + 1):#to_keep[:k]:#range(i+1):
lim_val_j = np.sqrt(all_invfisher[k][j][j])
# ax[i][j].set_xlim(-1. * extrema[i], extrema[i])
if i == j:
lim_val = np.sqrt(all_invfisher[k][i][i])
x_grid = np.linspace(-3. * lim_val_i, 3. * lim_val_i, 100)
func = sps.norm(0., lim_val)
ax[i][j].plot(x_grid, func.pdf(x_grid),
color=colors[k], label=names[k], alpha=1., linestyle=styles[k][0], linewidth=3.)
ax[i][j].set_xlabel(labels[i], fontsize=24)
ax[i][j].set_yticks([])
else:
ellipse = Ellipse(xy=(0., 0.), height=2.*a[i][j], width=2.*b[i][j], angle=(90. - t[i][j]*180./np.pi),
alpha=alphas[k], color=colors[k], linestyle=styles[k][0], lw=2., hatch=hatches[k], fill=fills[k])
ax[i][j].add_artist(ellipse)
if j == 0:
ax[i][j].set_ylabel(labels[i], fontsize=24)
if i == ncol-1:
ax[i][j].set_xlabel(labels[j], fontsize=24)
# note x and y swapping here. . .
extreme_y = a[i][j]*abs(np.cos(t[i][j]))+b[i][j]*abs(np.sin(t[i][j]))
extreme_x = a[i][j]*abs(np.sin(t[i][j]))+b[i][j]*abs(np.cos(t[i][j]))
ax[i][j].set_xlim(-1.5*extreme_y, 1.5*extreme_y)
ax[i][j].set_ylim(-1.5*extreme_x, 1.5*extreme_x)
plt.subplots_adjust(hspace=0., wspace=0.)
plt.savefig('../results/CosmoLike/final_plot'+extra+'.png', dpi=100, bbox_inches='tight', pad_inches = 0)
return
In [ ]:
# mycorner(fish_params, mini_keys, names, mini_fancy, 'mini')
mycorner(fish_params, keys, names, hatches, fills, alphas, fancy)
In [ ]:
In [ ]:
test_dir = '../results/single'
In [ ]:
simulated_posteriors = chippr.catalog(params='single.txt', loc=test_dir)
In [ ]:
data = simulated_posteriors.read(loc='data', style='.txt')
In [ ]:
params = chippr.utils.ingest('single.txt')
def check_prob_params(params):
"""
Sets parameter values pertaining to components of probability
Parameters
----------
params: dict
dictionary containing key/value pairs for probability
Returns
-------
params: dict
dictionary containing key/value pairs for probability
"""
if 'prior_mean' not in params:
params['prior_mean'] = 'interim'
else:
params['prior_mean'] = params['prior_mean'][0]
if 'no_prior' not in params:
params['no_prior'] = 0
else:
params['no_prior'] = int(params['no_prior'][0])
if 'no_data' not in params:
params['no_data'] = 0
else:
params['no_data'] = int(params['no_data'][0])
return params
params = check_prob_params(params)
def set_up_prior(data, params):
"""
Function to create prior distribution from data
Parameters
----------
data: dict
catalog dictionary containing bin endpoints, log interim prior, and log
interim posteriors
params: dict
dictionary of parameter values for creation of prior
Returns
-------
prior: chippr.mvn object
prior distribution as multivariate normal
"""
zs = data['bin_ends']
log_nz_intp = data['log_interim_prior']
log_z_posts = data['log_interim_posteriors']
z_difs = zs[1:]-zs[:-1]
z_mids = (zs[1:]+zs[:-1])/2.
n_bins = len(z_mids)
n_pdfs = len(log_z_posts)
a = 1.# / n_bins
b = 20.#1. / z_difs ** 2
c = a / n_pdfs
prior_var = np.eye(n_bins)
for k in range(n_bins):
prior_var[k] = a * np.exp(-0.5 * b * (z_mids[k] - z_mids) ** 2)
prior_var += c * np.identity(n_bins)
prior_mean = log_nz_intp
prior = chippr.mvn(prior_mean, prior_var)
if params['prior_mean'] == 'sample':
new_mean = prior.sample_one()
prior = chippr.mvn(new_mean, prior_var)
print(params['prior_mean'], prior_mean, new_mean)
else:
print(params['prior_mean'], prior_mean)
return (prior, prior_var)
(prior, cov) = set_up_prior(data, params)
In [ ]:
with open(os.path.join('../results/single/data', 'true_params.p'), 'r') as true_file:
true_nz_info = pickle.load(true_file)
In [ ]:
true_nz_info['amps']
In [ ]:
true_funcs = [chippr.discrete(np.array([true_nz_info['bins'][0], true_nz_info['bins'][-1]]), true_nz_info['amps'])]
true_amps = true_nz_info['amps']
# # print(true_amps)
true_nz = true_funcs[0]#chippr.gmix(true_amps, true_funcs, limits=(min(true_nz_info['bins']), max(true_nz_info['bins'])))
grid_mids = (true_nz_info['bins'][1:] + true_nz_info['bins'][:-1])/2.
plt.plot(grid_mids, true_nz.evaluate(grid_mids))
In [ ]:
nz = chippr.log_z_dens(data, prior, truth=true_nz, loc='../results/single', vb=True)
In [ ]:
nz.read('nz.p')
In [ ]:
bins_to_write = np.linspace(0.0101, 3.5001, 350)
empty_bins = np.random.random((350, 4))/350.
In [ ]:
bin_mids = (nz.info['bin_ends'][1:] + nz.info['bin_ends'][:-1])/2.
with open(os.path.join('../results/single/results', 'nz_mmle_test.txt'), 'wb') as cosmolike_file:
# cosmolike_file.write(zip(bin_mids, np.exp(nz.info['estimators']['log_mmle_nz']))
cosmolike_writer = csv.writer(cosmolike_file, delimiter=' ')
cosmolike_writer.writerows(zip(bin_mids, np.exp(nz.info['estimators']['log_mmle_nz'])))
How to run cosmolike. . .
In [ ]:
z_0 = 0.3
def smooth_func(z):
return 1/(2 * z_0) * (z/z_0)**2 * np.exp(-z/z_0)
zs = np.linspace(0., 1., 100)
nz = smooth_func(zs[:-1])
nz /= np.dot(nz, zs[1:]-zs[:-1])
plt.plot(zs[:-1], nz)
In [ ]:
print(np.dot(true_funcs[0].evaluate(grid_mids), true_nz_info['grid'][1:]-true_nz_info['grid'][:-1]))
In [ ]:
print(np.dot(true_funcs[0].distweights, true_funcs[0].dbins))
In [ ]:
func = sps.norm(0.25, 0.05)
print(func.std())
x = np.linspace(0., 1., 100)
plt.plot(x, func.pdf(x))
In [ ]:
n_bins = 100
z_mids = np.linspace(0., 1., n_bins)
a = 1.
b = 20.#mid-scale variations, larger means more peaks
c = 1.e-6#longest-scale variation, lower increases amplitude relative to small-scale
prior_var = np.eye(n_bins)
for k in range(n_bins):
prior_var[k] = a * np.exp(-0.5 * b * (z_mids[k] - z_mids) ** 2)
prior_var += c * np.identity(n_bins)
prior_mean = np.zeros(n_bins)
prior = chippr.mvn(prior_mean, prior_var)
In [ ]:
samples = prior.sample(7)
for each in samples:
plt.plot(z_mids, each)
In [ ]:
In [ ]: