In [1]:
%load_ext watermark
%watermark -u -d -v -p numpy,matplotlib,scipy,pandas,sklearn,mlxtend
In [2]:
%matplotlib inline
from __future__ import division, print_function
from collections import defaultdict
import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib.ticker import FormatStrFormatter
import matplotlib.gridspec as gridspec
import healpy as hp
import dask
from dask import delayed, multiprocessing
from dask.diagnostics import ProgressBar
import pyprind
from scipy.special import erfcinv
import comptools as comp
import comptools.analysis.plotting as plotting
import comptools.anisotropy.anisotropy as anisotropy
import comptools.anisotropy.teststatistic as ts
color_dict = comp.analysis.get_color_dict()
In [29]:
config = ['IC86.2011', 'IC86.2012', 'IC86.2013', 'IC86.2014', 'IC86.2015']
years_str = '2011-2015'
composition='all'
n_side = 64
scale = 3
smooth = 0.0
n_bins = 60
decmax = -55
low_energy = True
In [30]:
figures_dir = os.path.join(comp.paths.figures_dir, 'anisotropy', 'random-trials',
'smooth-{}_RAbins-{}_decmax-{}'.format(int(smooth), int(n_bins), int(decmax)))
if not os.path.isdir(figures_dir):
print('Making the directory {}'.format(figures_dir))
os.makedirs(figures_dir)
In [31]:
n_total = anisotropy.get_num_events(config=config, composition='all', decmax=decmax, low_energy=low_energy)
n_light = anisotropy.get_num_events(config=config, composition='light', decmax=decmax, low_energy=low_energy)
n_heavy = anisotropy.get_num_events(config=config, composition='heavy', decmax=decmax, low_energy=low_energy)
In [32]:
n_light/n_total, n_heavy/n_total
Out[32]:
In [33]:
kwargs_relint = {'config': config, 'low_energy': low_energy, 'smooth': smooth, 'scale': scale, 'decmax': decmax}
In [34]:
relint_all = anisotropy.get_map(name='relint', composition='all', **kwargs_relint)
relint_light = anisotropy.get_map(name='relint', composition='light', **kwargs_relint)
relint_heavy = anisotropy.get_map(name='relint', composition='heavy', **kwargs_relint)
In [35]:
kwargs_plot_relint = {'smooth': smooth, 'scale': scale, 'decmax': decmax}
In [36]:
print(relint_all.max())
print(relint_all[relint_all != hp.UNSEEN].min())
title = 'Relative Intensity [$\\times \ 10^{}$]'.format('{'+str(-scale)+'}')
fig, ax = anisotropy.plot_skymap(relint_all, color_palette='RdBu_r', symmetric=True, llabel=years_str,
cbar_title=title, cbar_min=-2.5, cbar_max=2.5, polar=True, **kwargs_plot_relint)
outfile = 'IC86-{}_relint_{}_nside-{}_smooth-{:0.1f}_decmax-{}.png'.format(years_str, 'all', n_side, smooth, decmax)
plt.savefig(os.path.join(figures_dir, outfile))
In [11]:
print(relint_light.max())
print(relint_light[relint_light != hp.UNSEEN].min())
title = 'Relative Intensity [$\\times \ 10^{}$]'.format('{'+str(-scale)+'}')
fig, ax = anisotropy.plot_skymap(relint_light, color_palette='RdBu_r', symmetric=True, llabel=years_str,
cbar_title=title, cbar_min=-2.5, cbar_max=2.5, polar=True, **kwargs_plot_relint)
outfile = 'IC86-{}_relint_{}_nside-{}_smooth-{:0.1f}_decmax-{}.png'.format(years_str, 'light', n_side, smooth, decmax)
plt.savefig(os.path.join(figures_dir, outfile))
In [12]:
print(relint_heavy.max())
print(relint_heavy[relint_heavy != hp.UNSEEN].min())
title = 'Relative Intensity [$\\times \ 10^{}$]'.format('{'+str(-scale)+'}')
fig, ax = anisotropy.plot_skymap(relint_heavy, color_palette='RdBu_r', symmetric=True, llabel=years_str,
cbar_title=title, cbar_min=-2.5, cbar_max=2.5, polar=True, **kwargs_plot_relint)
outfile = 'IC86-{}_relint_{}_nside-{}_smooth-{:0.1f}_decmax-{}.png'.format(years_str, 'heavy', n_side, smooth, decmax)
plt.savefig(os.path.join(figures_dir, outfile))
In [13]:
kwargs_sig = {'config': config, 'low_energy': low_energy, 'smooth': smooth, 'scale': None, 'decmax': decmax}
In [14]:
sig_all = anisotropy.get_map(name='sig', composition='all', **kwargs_sig)
sig_light = anisotropy.get_map(name='sig', composition='light', **kwargs_sig)
sig_heavy = anisotropy.get_map(name='sig', composition='heavy', **kwargs_sig)
In [15]:
kwargs_plot_sig = {'smooth': smooth, 'scale': None, 'decmax': decmax}
In [16]:
print(sig_all.max())
print(sig_all[sig_all != hp.UNSEEN].min())
title = 'Significance [$\mathrm{\sigma_{LM}}$]'
fig, ax = anisotropy.plot_skymap(sig_all, color_palette='RdBu_r', symmetric=True, llabel=years_str,
cbar_title=title, cbar_min=-4.1, cbar_max=4.1, polar=True, **kwargs_plot_sig)
outfile = 'IC86-{}_sig_{}_nside-{}_smooth-{:0.1f}_decmax-{}.png'.format(years_str, 'all', n_side, smooth, decmax)
plt.savefig(os.path.join(figures_dir, outfile))
In [17]:
print(sig_light.max())
print(sig_light[sig_light != hp.UNSEEN].min())
title = 'Significance [$\mathrm{\sigma_{LM}}$]'
fig, ax = anisotropy.plot_skymap(sig_light, color_palette='RdBu_r', symmetric=True, llabel=years_str,
cbar_title=title, cbar_min=-4.1, cbar_max=4.1, polar=True, **kwargs_plot_sig)
outfile = 'IC86-{}_sig_{}_nside-{}_smooth-{:0.1f}_decmax-{}.png'.format(years_str, 'light', n_side, smooth, decmax)
plt.savefig(os.path.join(figures_dir, outfile))
In [18]:
print(sig_heavy.max())
print(sig_heavy[sig_heavy != hp.UNSEEN].min())
title = 'Significance [$\mathrm{\sigma_{LM}}$]'
fig, ax = anisotropy.plot_skymap(sig_heavy, color_palette='RdBu_r', symmetric=True, llabel=years_str,
cbar_title=title, cbar_min=-4.1, cbar_max=4.1, polar=True, **kwargs_plot_sig)
outfile = 'IC86-{}_sig_{}_nside-{}_smooth-{:0.1f}_decmax-{}.png'.format(years_str, 'heavy', n_side, smooth, decmax)
plt.savefig(os.path.join(figures_dir, outfile))
[ back to top ]
In [37]:
kwargs_relint_proj = dict(kwargs_relint)
kwargs_relint_proj['scale'] = None
Get relative intensity (statistical) error plot
In [38]:
relint_all = anisotropy.get_map(name='relint', composition='all', **kwargs_relint_proj)
relint_all_err = anisotropy.get_map(name='relerr', composition='all', **kwargs_relint_proj)
Get projected relative intensity
In [39]:
# ri_all, ri_all_err, ra, ra_err = anisotropy.get_proj_relint(relint_all, relint_all_err,
# n_bins=n_bins, decmax=decmax)
data_all = anisotropy.get_map(name='data', composition='all', **kwargs_relint_proj)
ref_err = anisotropy.get_map(name='ref', composition='all', **kwargs_relint_proj)
ri_all, ri_all_err, ra, ra_err = anisotropy.get_binned_relint(data_all, ref_err,
n_bins=n_bins, decmax=decmax)
In [40]:
n_dof = ri_all.shape[0]
chi2_all = np.sum(ri_all**2 / ri_all_err**2)
chi2_all_red = chi2_all / n_dof
chi2_all
Out[40]:
In [41]:
from scipy.stats import chi2
In [42]:
p_all = chi2.sf(chi2_all, n_dof, loc=0, scale=1)
s_all = erfcinv(2*p_all)*np.sqrt(2)
print(p_all, s_all)
In [43]:
chi2_samples = chi2.rvs(n_dof, loc=0, scale=1, size=int(2e4), random_state=2)
fig, ax = plt.subplots()
counts, bins, _ = ax.hist(chi2_samples, bins=75, alpha=0.7, label='$\mathrm{\chi^2}$ distribution '+'({} d.o.f)'.format(n_dof))
ax.axvline(chi2_all, marker='None', ls='-.', color='C1', lw=1.5, label='Observed $\mathrm{\chi^2}$ ' + '({:0.2f}$\\sigma$)'.format(s_all))
ax.set_ylabel('Counts')
ax.set_xlabel('$\mathrm{\chi^2}$')
ax.grid()
# ax.legend()
ax.legend(bbox_to_anchor=(0., 1.02, 1., .102), loc=3,
ncol=2, mode="expand", borderaxespad=0., frameon=False)
outfile = 'chi2-dist_proj-RI-all.png'
plt.savefig(os.path.join(figures_dir, outfile))
plt.show()
In [44]:
chi2.sf(chi2_all, ri_all.shape[0]-1, loc=0, scale=1)
Out[44]:
In [45]:
popt_all, perr_all, _ = anisotropy.get_proj_fit_params(np.deg2rad(ra), ri_all, sigmay=ri_all_err, l=3)
amp_all = popt_all[1]
amp_err_all = perr_all[1]
phase_all = np.rad2deg(popt_all[2])
phase_err_all = np.rad2deg(perr_all[2])
In [46]:
fig, ax = plt.subplots()
# ax.errorbar(ra, ri_all, yerr=ri_all_err, marker='.', ls='None', c='C2', label='all')
ra_bins = np.linspace(0, 360, n_bins + 1)
plotting.plot_steps(ra_bins, ri_all, yerr=ri_all_err, color='C2', label='all', fillalpha=0.2, ax=ax)
# ax.plot(ra, anisotropy.cos_fit_func(np.deg2rad(ra), *popt_all[:3]), color='C2', marker='None')
ax.axhline(0, marker='None', ls='-.', c='k')
ax.set_ylabel('$\mathrm{\langle RI \\rangle}$')
ax.set_xlabel('RA [ $^{\circ}$]')
ax.grid()
# ax.set_ylim(-1.5e-3, 1.5e-3)
ax.set_xlim(0, 360)
ax.invert_xaxis()
ax.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
# ax.legend()
# all_amp_str = 'Amp = {:0.2e} +/- {:0.1e}'.format(amp_all, amp_err_all)
# all_phase_str = 'Phase = {:0.2f} {} +/- {:0.2f} {}'.format(phase_all, '$^{\circ}$', phase_err_all, '$^{\circ}$')
# ax.text(250, 5.0e-3, all_amp_str + '\n' + all_phase_str)
outfile = 'IC86-{}_proj-RI-all.png'.format(years_str)
plt.savefig(os.path.join(figures_dir, outfile))
plt.show()
In [47]:
# relint_light = anisotropy.get_map(name='relint', composition='light', **kwargs_relint_proj)
# relint_light_err = anisotropy.get_map(name='relerr', composition='light', **kwargs_relint_proj)
# ri_light, ri_light_err, ra, ra_err = anisotropy.get_proj_relint(relint_light, relint_light_err,
# n_bins=n_bins, decmax=decmax)
In [48]:
data_light = anisotropy.get_map(name='data', composition='light', **kwargs_relint_proj)
ref_light = anisotropy.get_map(name='ref', composition='light', **kwargs_relint_proj)
ri_light, ri_light_err, ra, ra_err = anisotropy.get_binned_relint(data_light, ref_light,
n_bins=n_bins, decmax=decmax)
In [49]:
# relint_heavy = anisotropy.get_map(name='relint', composition='heavy', **kwargs_relint_proj)
# relint_heavy_err = anisotropy.get_map(name='relerr', composition='heavy', **kwargs_relint_proj)
# ri_heavy, ri_heavy_err, ra, ra_err = anisotropy.get_proj_relint(relint_heavy, relint_heavy_err,
# n_bins=n_bins, decmax=decmax)
In [50]:
data_heavy = anisotropy.get_map(name='data', composition='heavy', **kwargs_relint_proj)
ref_heavy = anisotropy.get_map(name='ref', composition='heavy', **kwargs_relint_proj)
ri_heavy, ri_heavy_err, ra, ra_err = anisotropy.get_binned_relint(data_heavy, ref_heavy,
n_bins=n_bins, decmax=decmax)
In [51]:
chi2_light = np.sum(ri_light**2 / ri_light_err**2)
chi2_heavy = np.sum(ri_heavy**2 / ri_heavy_err**2)
p_light = chi2.sf(chi2_light, n_dof, loc=0, scale=1)
s_light = erfcinv(2*p_light)*np.sqrt(2)
print(p_light, s_light)
p_heavy = chi2.sf(chi2_heavy, n_dof, loc=0, scale=1)
s_heavy = erfcinv(2*p_heavy)*np.sqrt(2)
print(p_heavy, s_heavy)
In [52]:
def calc_chi2(ri_0, ri_1, ri_err_0, ri_err_1):
chi2 = np.sum((ri_0 - ri_1)**2/(ri_err_0**2 + ri_err_1**2)) / ri_0.shape[0]
return chi2
In [53]:
calc_chi2(ri_light, ri_heavy, ri_light_err, ri_heavy_err)
Out[53]:
In [54]:
light_skymap_file = anisotropy.get_skymap_files(config=config, n_side=n_side, composition='light',
low_energy=low_energy)
heavy_skymap_file = anisotropy.get_skymap_files(config=config, n_side=n_side, composition='heavy',
low_energy=low_energy)
chi2_data = ts.get_proj_RI_red_chi2(light_skymap_file, heavy_skymap_file, smooth=smooth, n_bins=n_bins, decmax=decmax)
print(chi2_data)
In [55]:
fig, ax = plt.subplots()
plotting.plot_steps(ra_bins, ri_light, yerr=ri_light_err, color='C0', label='light', fillalpha=0.2, ax=ax)
plotting.plot_steps(ra_bins, ri_heavy, yerr=ri_heavy_err, color='C1', label='heavy', fillalpha=0.2, ax=ax)
ax.axhline(0, marker='None', ls='-.', c='k')
ax.set_ylabel('$\mathrm{\langle RI \\rangle}$')
ax.set_xlabel('RA [ $^{\circ}$]')
ax.set_title('$\mathrm{\chi^2_{red}}$' + ' = {:0.1f}'.format(chi2_data))
ax.legend()
ax.grid()
ax.set_xlim(ra_bins.min(), ra_bins.max())
ax.invert_xaxis()
ax.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
outfile = 'IC86-{}_proj-RI-comps.png'.format(years_str)
plt.savefig(os.path.join(figures_dir, outfile))
plt.show()
In [56]:
fig, ax = plt.subplots()
plotting.plot_steps(ra_bins, ri_light, yerr=ri_light_err, color='C0',
label='light ({:0.2f}$\sigma$)'.format(s_light), fillalpha=0.2, ax=ax)
plotting.plot_steps(ra_bins, ri_heavy, yerr=ri_heavy_err, color='C1',
label='heavy ({:0.2f}$\sigma$)'.format(s_heavy), fillalpha=0.2, ax=ax)
plotting.plot_steps(ra_bins, ri_all, yerr=ri_all_err, color='C2',
label='all ({:0.2f}$\sigma$)'.format(s_all), fillalpha=0.2, ax=ax)
ax.axhline(0, marker='None', ls='-.', c='k')
ax.set_ylabel('$\mathrm{\langle RI \\rangle }$')
ax.set_xlabel('RA [ $^{\circ}$]')
ax.grid()
ax.set_xlim(0, 360)
ax.invert_xaxis()
ax.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
ax.legend(bbox_to_anchor=(0., 1.02, 1., .102), loc=3,
ncol=3, mode="expand", borderaxespad=0., frameon=False)
outfile = 'IC86-{}_proj-RI.png'.format(years_str)
plt.savefig(os.path.join(figures_dir, outfile))
plt.show()
In [63]:
df_test_stat = ts.load_test_stats(config=config, low_energy=low_energy, smooth=smooth, n_bins=n_bins, decmax=decmax)
df_test_stat.head()
Out[63]:
In [64]:
def get_proj_trial_dict(config, trial_num):
proj_dict = {}
proj_dict['trial'] = trial_num
trial_files_0 = []
trial_files_1 = []
for c in config:
map_dir = os.path.join(comp.paths.comp_data_dir, c + '_data',
'anisotropy/random_trials')
sample_0_file = os.path.join(map_dir, 'random_split_0_trial-{}.fits'.format(trial_num))
trial_files_0.append(sample_0_file)
sample_1_file = os.path.join(map_dir, 'random_split_1_trial-{}.fits'.format(trial_num))
trial_files_1.append(sample_1_file)
relint_0 = anisotropy.get_map(files=trial_files_0, name='relint', **kwargs_relint_proj)
relint_1 = anisotropy.get_map(files=trial_files_1, name='relint', **kwargs_relint_proj)
relerr_0 = anisotropy.get_map(files=trial_files_0, name='relerr', **kwargs_relint_proj)
relerr_1 = anisotropy.get_map(files=trial_files_1, name='relerr', **kwargs_relint_proj)
ri_0, ri_0_err, ra, _ = anisotropy.get_proj_relint(relint_0, relerr_0, n_bins=n_bins, decmax=decmax)
ri_1, ri_1_err, ra, _ = anisotropy.get_proj_relint(relint_1, relerr_1, n_bins=n_bins, decmax=decmax)
proj_dict['ri_0'] = ri_0
proj_dict['ri_0_err'] = ri_0_err
proj_dict['ri_1'] = ri_1
proj_dict['ri_1_err'] = ri_1_err
proj_dict['chi2'] = ts.get_proj_RI_red_chi2(trial_files_0, trial_files_1, smooth=smooth,
n_bins=n_bins, decmax=decmax)
return proj_dict
In [65]:
n_trials = 16
proj_dicts = [delayed(get_proj_trial_dict)(config, trial_num) for trial_num in range(n_trials)]
proj_df = delayed(pd.DataFrame.from_records)(proj_dicts, index='trial')
In [66]:
with ProgressBar() as bar:
proj_df = proj_df.compute(get=multiprocessing.get, num_workers=10)
proj_df.head()
Out[66]:
In [67]:
n_side = np.sqrt(n_trials).astype(int)
fig, axarr = plt.subplots(n_side, n_side, figsize=(10, 8), sharex=True, sharey=True)
for trial, ax in zip(proj_df.index, axarr.flatten()):
proj_trial = proj_df.loc[trial]
ra_bins = np.linspace(0, 360, n_bins + 1)
ax = plotting.plot_steps(ra_bins, proj_trial['ri_0'], yerr=proj_trial['ri_0_err'], color='C3', ax=ax)
ax = plotting.plot_steps(ra_bins, proj_trial['ri_1'], yerr=proj_trial['ri_1_err'], color='C4', ax=ax)
ax.set_title('$\mathrm{\chi^2_{red}}$' + ' = {:0.1f}'.format(proj_trial['chi2']))
ax.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
ax.set_xlim(ra_bins.min(), ra_bins.max())
# ax.set_ylim(-1e-2, 1e-2)
ax.invert_xaxis()
ax.grid()
for i, row in enumerate(axarr):
for j, cell in enumerate(row):
if i == len(axarr) - 1:
cell.set_xlabel('RA [ $^{\circ}$]')
if j == 0:
cell.set_ylabel('$\mathrm{\langle RI \\rangle }$')
plt.tight_layout()
plt.savefig(os.path.join(figures_dir, 'projRI-random-trials-grid.png'))
plt.show()
In [68]:
test_stat = df_test_stat['proj_RI_red_chi2']
pval = np.sum(test_stat > chi2_data)/len(test_stat)
print(pval)
significance = erfcinv(2*pval)*np.sqrt(2)
print(significance)
In [69]:
fig, ax = plt.subplots()
# chi2_max = 1.0
chi2_max = 5
chi2_bins = np.linspace(0, chi2_max, 100)
counts = np.histogram(df_test_stat['proj_RI_red_chi2'], bins=chi2_bins)[0]
# ax = plotting.plot_steps(chi2_bins, counts, yerr=np.sqrt(counts))
ax = plotting.plot_steps(chi2_bins, counts, yerr=np.sqrt(counts), label='Random trials')
ax.axvline(chi2_data, marker='None', ls='-.', color='C2', lw=1.5,
label='Light-heavy split \n ($\mathrm{\chi^2_{red}} = '+'${:0.2f})'.format(chi2_data))
ax.set_xlabel('$\mathrm{\chi^2_{red}}$')
ax.set_ylabel('Counts')
ax.set_title('IC86 {} random trials'.format(years_str))
# ax.text(0.3, 4, 'p-value $= {:g}$'.format(pval) + \
# '\n significance $= {:0.2f}\sigma$'.format(significance))
ax.set_ylim(1e-1)
ax.set_xlim(0, chi2_max)
ax.set_yscale('log', nonposy='clip')
ax.legend()
ax.grid()
plt.savefig(os.path.join(figures_dir, 'random-trials-chi2_{}.png'.format(years_str)))
plt.show()
In [70]:
fig, ax = plt.subplots()
# chi2_max = 1.0
chi2_max = 5*n_bins
chi2_bins = np.linspace(0, chi2_max, 100)
counts = np.histogram(df_test_stat['proj_RI_red_chi2']*n_bins, bins=chi2_bins, density=True)[0]
# ax = plotting.plot_steps(chi2_bins, counts, yerr=np.sqrt(counts))
ax = plotting.plot_steps(chi2_bins, counts, label='Random trials')
x = np.linspace(0, 5*n_bins, 100)
ax.plot(x, chi2.pdf(x, df=36))
# ax.axvline(chi2_data, marker='None', ls='-.', color='C2', lw=1.5,
# label='Light-heavy split \n ($\mathrm{\chi^2_{red}} = '+'${:0.2f})'.format(chi2_data))
ax.set_xlabel('$\mathrm{\chi^2_{red}}$')
ax.set_ylabel('Counts')
ax.set_title('IC86 {} random trials'.format(years_str))
# ax.text(0.3, 4, 'p-value $= {:g}$'.format(pval) + \
# '\n significance $= {:0.2f}\sigma$'.format(significance))
# ax.set_ylim(1e-1)
ax.set_xlim(0, chi2_max)
# ax.set_yscale('log', nonposy='clip')
ax.legend()
ax.grid()
# plt.savefig(os.path.join(figures_dir, 'random-trials-chi2_{}.png'.format(years_str)))
plt.show()
In [ ]:
chi2.pdf(np.linspace(0, 5, 100), df=36)
In [ ]: