In [1]:
%load_ext watermark
%watermark -u -d -v -p numpy,matplotlib,scipy,pandas,sklearn,mlxtend


last updated: 2017-07-20 

CPython 2.7.13
IPython 5.3.0

numpy 1.12.1
matplotlib 2.0.2
scipy 0.19.0
pandas 0.20.1
sklearn 0.18.1
mlxtend 0.6.0

In [47]:
%matplotlib inline
from __future__ import division, print_function
from collections import defaultdict
import os
import glob
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import healpy as hp
from scipy.stats import ks_2samp
import multiprocessing as mp
from scipy.special import erfcinv
import pyprind
import dask
from dask import delayed, multiprocessing
from dask.diagnostics import ProgressBar

from icecube import astro

import comptools as comp
import comptools.analysis.plotting as plotting
import comptools.anisotropy.anisotropy as anisotropy

color_dict = comp.analysis.get_color_dict()

In [70]:
config = ['IC86.2011', 'IC86.2012', 'IC86.2013', 'IC86.2014','IC86.2015']
n_side = 64
smooth_rad = 5.0
scale = 3
decmax = -55
years_str = '2011-2015'
composition='all'
low_energy = True
n_bins = 36

In [71]:
df = anisotropy.get_test_stats(config=config, low_energy=low_energy, smooth=smooth_rad, n_bins=n_bins)


Loading /data/user/jbourbeau/composition/anisotropy_random_trials/teststat-2011-2012-2013-2014-2015_smooth-5.0_RAbins-36_lowenergy.hdf...

In [72]:
df['proj_RI_red_chi2'][df['proj_RI_red_chi2'].isnull()].index


Out[72]:
Int64Index([], dtype='int64')

In [73]:
df['proj_RI_red_chi2'].describe()


Out[73]:
count    1000.000000
mean        1.103543
std         0.334419
min         0.311463
25%         0.866288
50%         1.062819
75%         1.310826
max         2.564191
Name: proj_RI_red_chi2, dtype: float64

In [74]:
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 [75]:
def get_chi2_data(config):
    kwargs_relint = {'smooth': smooth_rad, 'scale': None, 'decmax': -55, 'low_energy': True}

    relint_light = anisotropy.get_map(config=config, name='relint', composition='light', **kwargs_relint)
    relint_heavy = anisotropy.get_map(config=config, name='relint', composition='heavy', **kwargs_relint)

    relerr_light = anisotropy.get_map(config=config, name='relerr', composition='light', **kwargs_relint)
    relerr_heavy = anisotropy.get_map(config=config, name='relerr', composition='heavy', **kwargs_relint)

    ri_light, ri_err_light, ra, _ = anisotropy.get_proj_relint(relint_light, relerr_light, n_bins=n_bins)
    ri_heavy, ri_err_heavy, ra, _ = anisotropy.get_proj_relint(relint_heavy, relerr_heavy, n_bins=n_bins)

    chi2_data = calc_chi2(ri_light, ri_heavy, ri_err_light, ri_err_heavy)
    
    return chi2_data

In [76]:
def calc_sig_ks_pval(config):

    kwargs_relint = {'smooth': smooth_rad, 'scale': None, 'decmax': -55, 'low_energy': True}

    # Get 2D Li-Ma significance maps
    sig_light = anisotropy.get_map(config=config, name='sig', composition='light', **kwargs_relint)
    sig_heavy = anisotropy.get_map(config=config, name='sig', composition='heavy', **kwargs_relint)

    # Construct masks
    is_good_mask_light = (sig_light != hp.UNSEEN) & ~np.isnan(sig_light)
    is_good_mask_heavy = (sig_heavy != hp.UNSEEN) & ~np.isnan(sig_heavy)

    # Calculate ks test statistic and corresponding p-value
    ks_statistic, pval = ks_2samp(sig_light[is_good_mask_light], sig_heavy[is_good_mask_heavy])
    
    # Calculate ~area between cumulative distributions
    bins = np.linspace(0, 5, 50)
#     counts_light = np.histogram(sig_light[is_good_mask_light], bins=bins)[0]
#     counts_heavy = np.histogram(sig_heavy[is_good_mask_heavy], bins=bins)[0]

    counts_light = np.histogram(np.abs(sig_light[is_good_mask_light]), bins=bins)[0]
    counts_heavy = np.histogram(np.abs(sig_heavy[is_good_mask_heavy]), bins=bins)[0]

    sig_cumsum_diff_area = np.abs(np.cumsum(counts_light) - np.cumsum(counts_heavy)).sum()

    return ks_statistic, pval, sig_cumsum_diff_area

In [77]:
plt.errorbar(ra, ri_light, yerr=ri_err_light)
plt.errorbar(ra, ri_heavy, yerr=ri_err_heavy)
plt.grid()
plt.show()



NameErrorTraceback (most recent call last)
<ipython-input-77-c14c2754d326> in <module>()
----> 1 plt.errorbar(ra, ri_light, yerr=ri_err_light)
      2 plt.errorbar(ra, ri_heavy, yerr=ri_err_heavy)
      3 plt.grid()
      4 plt.show()

NameError: name 'ra' is not defined

In [78]:
chi2_data = get_chi2_data(config)
sig_ks_dist_data, sig_pval_data, sig_cumsum_diff_area = calc_sig_ks_pval(config)

In [79]:
sig_ks_dist_data, sig_pval_data, sig_cumsum_diff_area


Out[79]:
(0.03701241134751776, 0.0040066424150784503, 5307)

In [80]:
test_stat = df['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)


0.0
inf

In [81]:
fig, ax = plt.subplots()
chi2_max = 4
chi2_bins = np.linspace(0, chi2_max, 100)
counts = np.histogram(df['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 2011-2015 random trials')
# 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(comp.paths.figures_dir, 'random-trials-chi2_IC86.2011-2015.png'))
plt.show()



In [16]:
df['sig_ks_dist'].plot(kind='hist', bins=50)


Out[16]:
<matplotlib.axes._subplots.AxesSubplot at 0x7fc8b8128410>

In [17]:
df['sig_cumsum_diff_area'].plot(kind='hist', bins=50)


Out[17]:
<matplotlib.axes._subplots.AxesSubplot at 0x7fc8b5f6f4d0>

In [18]:
np.log10(df['sig_pval']).plot(kind='hist', bins=50, logy=True)


Out[18]:
<matplotlib.axes._subplots.AxesSubplot at 0x7fc8b5f1fc50>

In [113]:
fig, ax = plt.subplots()
sig_pval_bins = np.linspace(-300, 0, 100)
counts = np.histogram(np.log10(df['sig_pval']), bins=sig_pval_bins)[0]
# ax = plotting.plot_steps(chi2_bins, counts, yerr=np.sqrt(counts))
ax = plotting.plot_steps(sig_pval_bins, counts, yerr=np.sqrt(counts), label='Random trials')
ax.axvline(np.log10(sig_pval_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('KS p-value')
ax.set_ylabel('Counts')
ax.set_title('IC86 2011-2015 random trials')
# ax.text(1, 5, 'p-value $= {:g}$'.format(pval) + \
#              '\n significance $= {:0.2f}\sigma$'.format(significance))
# ax.set_ylim((1e-1, 1e2))
# ax.set_xlim(0, chi2_max)
# ax.set_yscale("log", nonposy='clip')
ax.legend()
ax.grid()
# plt.savefig(os.path.join(comp.paths.figures_dir, 'random-trials-chi2_IC86.2012-2014.png'))
plt.show()



In [35]:
df['sig_ks_dist'].describe()


Out[35]:
count    1000.000000
mean        0.044340
std         0.018044
min         0.009304
25%         0.031028
50%         0.042115
75%         0.055242
max         0.116578
Name: sig_ks_dist, dtype: float64

In [47]:
fig, ax = plt.subplots()
sig_ks_bins = np.linspace(0, 0.4, 75)
counts = np.histogram(df['sig_ks_dist'], bins=sig_ks_bins)[0]
# ax = plotting.plot_steps(chi2_bins, counts, yerr=np.sqrt(counts))
ax = plotting.plot_steps(sig_ks_bins, counts, yerr=np.sqrt(counts), label='Random trials')
ax.axvline(sig_ks_dist_data, marker='None', ls='-.', color='C2', lw=1.5,
           label='Light-heavy split')
#            label='Light-heavy split \n ($\mathrm{\chi^2_{red}} = '+'${:0.2f})'.format(chi2_data))
ax.set_xlabel('KS distance')
ax.set_ylabel('Counts')
ax.set_title('IC86 2011-2015 random trials')

test_stat = df['sig_ks_dist']
pval = np.sum(test_stat > sig_ks_dist_data)/len(test_stat)
significance = erfcinv(2*pval)*np.sqrt(2)

ax.text(0.06, 2, 'p-value $= {:g}$'.format(pval) + \
             '\n significance $= {:0.2f}\sigma$'.format(significance))
ax.set_ylim(0)
# ax.set_xlim(0, chi2_max)
ax.set_yscale("log", nonposy='clip')
ax.set_ylim(1e-1)
ax.legend()
ax.grid()
plt.savefig(os.path.join(comp.paths.figures_dir, 'random-trials-ks-dist_IC86.2011-2015.png'))
plt.show()



In [115]:
test_stat = df['sig_ks_dist']
pval = np.sum(test_stat > sig_ks_dist_data)/len(test_stat)
print(pval)
significance = erfcinv(2*pval)*np.sqrt(2)
print(significance)


0.184
0.900225985701

In [125]:
fig, ax = plt.subplots()
sig_cumsum_diff_bins = np.linspace(0, 70000, 100)
counts = np.histogram(df['sig_cumsum_diff_area'], bins=sig_cumsum_diff_bins)[0]
# ax = plotting.plot_steps(chi2_bins, counts, yerr=np.sqrt(counts))
ax = plotting.plot_steps(sig_cumsum_diff_bins, counts, yerr=np.sqrt(counts), label='Random trials')
ax.axvline(sig_cumsum_diff_area, 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('Cumulative sum area proxy')
ax.set_ylabel('Counts')
ax.set_title('IC86 2011-2015 random trials')
# ax.text(1, 5, 'p-value $= {:g}$'.format(pval) + \
#              '\n significance $= {:0.2f}\sigma$'.format(significance))
# ax.set_ylim((1e-1, 1e2))
# ax.set_xlim(0, chi2_max)
# ax.set_yscale("log", nonposy='clip')
ax.legend()
ax.grid()
# plt.savefig(os.path.join(comp.paths.figures_dir, 'random-trials-chi2_IC86.2012-2014.png'))
plt.show()



In [126]:
test_stat = df['sig_cumsum_diff_area']
pval = np.sum(test_stat > sig_cumsum_diff_area)/len(test_stat)
print(pval)
significance = erfcinv(2*pval)*np.sqrt(2)
print(significance)


0.155
1.01522203322

In [40]:
fig, ax = plt.subplots()
sig_ks_bins = np.linspace(0, 0.5, 100)
configs = [['IC86.2012'], ['IC86.2012', 'IC86.2013', 'IC86.2014']]
# configs = [['IC86.2012'], ['IC86.2012', 'IC86.2013'], ['IC86.2012', 'IC86.2013', 'IC86.2014']]
for config, color in zip(configs, ['C0', 'C1', 'C2']):
    print('On {}...'.format(config))
    df = anisotropy.get_test_stats(config=config, low_energy=True)
    counts = np.histogram(df['sig_ks_dist'], bins=sig_ks_bins)[0]
    ax = plotting.plot_steps(sig_ks_bins, counts, yerr=np.sqrt(counts), label='{} years'.format(len(config)), color=color, ax=ax)
#     chi2_data = get_chi2_data(config)
#     ax.axvline(chi2_data, marker='None', ls='-.', color=color, lw=1.5)
# #                label='Light-heavy split \n ($\mathrm{\chi^2_{red}} = '+'${:0.2f})'.format(chi2_data))
# ax.axvline(6.1, marker='None', ls='-.', color=color, lw=1.5)
ax.set_xlabel('KS distance')
ax.set_ylabel('Counts')
ax.set_title('1000 random trials')
#     ax.text(1, 5, 'p-value $= {:g}$'.format(pval) + \
#                  '\n significance $= {:0.2f}\sigma$'.format(significance))
# ax.set_ylim((1e-1, 1e2))
# ax.set_xlim(0, chi2_max)
# ax.set_yscale("log", nonposy='clip')
ax.legend()
ax.grid()
# fig.text(0.5, 0.0, 'RA [ $^{\circ}$]', ha='center', fontsize=16)
# fig.text(0.0, 0.5, '$\mathrm{\langle RI \\rangle }$', va='center', rotation='vertical', fontsize=16)
plt.tight_layout()
# plt.savefig(os.path.join(comp.paths.figures_dir, 'trials_compare-running-sum.png'))
plt.show()


On ['IC86.2012']...
Loading /data/user/jbourbeau/composition/anisotropy_random_trials/teststat_2012_lowenergy.hdf...
On ['IC86.2012', 'IC86.2013', 'IC86.2014']...
Loading /data/user/jbourbeau/composition/anisotropy_random_trials/teststat_2012_2013_2014_lowenergy.hdf...

In [41]:
fig, ax = plt.subplots()
sig_pval_bins = np.linspace(-300, 0, 100)
configs = [['IC86.2012'], ['IC86.2012', 'IC86.2013', 'IC86.2014']]
# configs = [['IC86.2012'], ['IC86.2012', 'IC86.2013'], ['IC86.2012', 'IC86.2013', 'IC86.2014']]
for config, color in zip(configs, ['C0', 'C1', 'C2']):
    print('On {}...'.format(config))
    df = anisotropy.get_test_stats(config=config, low_energy=True)
    counts = np.histogram(np.log10(df['sig_pval']), bins=sig_pval_bins)[0]
    ax = plotting.plot_steps(sig_pval_bins, counts, yerr=np.sqrt(counts), label='{} years'.format(len(config)), color=color, ax=ax)
#     chi2_data = get_chi2_data(config)
#     ax.axvline(chi2_data, marker='None', ls='-.', color=color, lw=1.5)
# #                label='Light-heavy split \n ($\mathrm{\chi^2_{red}} = '+'${:0.2f})'.format(chi2_data))
# ax.axvline(6.1, marker='None', ls='-.', color=color, lw=1.5)
ax.set_xlabel('KS p-value')
ax.set_ylabel('Counts')
ax.set_title('1000 random trials')
#     ax.text(1, 5, 'p-value $= {:g}$'.format(pval) + \
#                  '\n significance $= {:0.2f}\sigma$'.format(significance))
# ax.set_ylim((1e-1, 1e2))
# ax.set_xlim(0, chi2_max)
# ax.set_yscale("log", nonposy='clip')
ax.legend()
ax.grid()
# fig.text(0.5, 0.0, 'RA [ $^{\circ}$]', ha='center', fontsize=16)
# fig.text(0.0, 0.5, '$\mathrm{\langle RI \\rangle }$', va='center', rotation='vertical', fontsize=16)
plt.tight_layout()
# plt.savefig(os.path.join(comp.paths.figures_dir, 'trials_compare-running-sum.png'))
plt.show()


On ['IC86.2012']...
Loading /data/user/jbourbeau/composition/anisotropy_random_trials/teststat_2012_lowenergy.hdf...
On ['IC86.2012', 'IC86.2013', 'IC86.2014']...
Loading /data/user/jbourbeau/composition/anisotropy_random_trials/teststat_2012_2013_2014_lowenergy.hdf...

In [4]:
fig, ax = plt.subplots()
chi2_max = 25
chi2_bins = np.linspace(0, 30, 100)
for config, color in zip(['IC86.2012', 'IC86.2013', 'IC86.2014'],
                         ['C0', 'C1', 'C2']):
    df = anisotropy.get_test_stats(config=config, low_energy=True)
    counts = np.histogram(df.chi2, bins=chi2_bins)[0]
    ax = plotting.plot_steps(chi2_bins, counts, yerr=np.sqrt(counts), label=config, color=color)
#     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('1000 random trials')
#     ax.text(1, 5, 'p-value $= {:g}$'.format(pval) + \
#                  '\n significance $= {:0.2f}\sigma$'.format(significance))
    ax.set_ylim((1e-1, 1e2))
    ax.set_xlim(0, chi2_max)
    ax.set_yscale("log", nonposy='clip')
    ax.legend()
    ax.grid()
plt.savefig(os.path.join(comp.paths.figures_dir, 'random-trials-chi2_compare-years.png'))
plt.show()


/home/jbourbeau/.virtualenvs/composition/lib/python2.7/site-packages/matplotlib/figure.py:1743: UserWarning: This figure includes Axes that are not compatible with tight_layout, so its results might be incorrect.
  warnings.warn("This figure includes Axes that are not "

In [11]:
fig, axarr = plt.subplots(2, 3, figsize=(10, 8), sharex=True, sharey=True)
for config, ax in zip(['IC86.2011', 'IC86.2012', 'IC86.2013', 'IC86.2014', 'IC86.2015'],
                            axarr.flatten()):
    print('On {}...'.format(config))
    kwargs_relint = {'smooth': 20, 'scale': None, 'decmax': -55, 'low_energy': True}

    relint_light = anisotropy.get_map(config=config, name='relint', composition='light', **kwargs_relint)
    relint_heavy = anisotropy.get_map(config=config, name='relint', composition='heavy', **kwargs_relint)
    relerr_light = anisotropy.get_map(config=config, name='relerr', composition='light', **kwargs_relint)
    relerr_heavy = anisotropy.get_map(config=config, name='relerr', composition='heavy', **kwargs_relint)

    ri_light, ri_err_light, ra, _ = anisotropy.get_proj_relint(relint_light, relerr_light, n_bins=24)
    ri_heavy, ri_err_heavy, ra, _ = anisotropy.get_proj_relint(relint_heavy, relerr_heavy, n_bins=24)
    
    ax.errorbar(ra, ri_light, yerr=ri_err_light, label='light', marker='.', ls=':')
    ax.errorbar(ra, ri_heavy, yerr=ri_err_heavy, label='heavy', marker='.', ls=':')
    ax.axhline(0, marker='None', ls='-.', c='k')
    
    chi2_data = calc_chi2(ri_light, ri_heavy, ri_err_light, ri_err_heavy)
    ax.set_title('{}: '.format(config)+'$\mathrm{\chi^2_{red}} =$'+ ' {:0.1f}'.format(chi2_data))
    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 }$')
    ax.grid()
    ax.legend()
# fig.text(0.5, 0.0, 'RA [ $^{\circ}$]', ha='center', fontsize=16)
# fig.text(0.0, 0.5, '$\mathrm{\langle RI \\rangle }$', va='center', rotation='vertical', fontsize=16)
plt.tight_layout()
plt.savefig(os.path.join(comp.paths.figures_dir, 'light-heavy-split_compare-years.png'))
plt.show()


On IC86.2011...
On IC86.2012...
On IC86.2013...
On IC86.2014...
On IC86.2015...

In [24]:
fig, axarr = plt.subplots(2, 3, figsize=(10, 8), sharex=True, sharey=True)
configs = [['IC86.2011'], ['IC86.2011', 'IC86.2012'], ['IC86.2011', 'IC86.2012', 'IC86.2013'], 
           ['IC86.2011', 'IC86.2012', 'IC86.2013', 'IC86.2014'],
           ['IC86.2011', 'IC86.2012', 'IC86.2013', 'IC86.2014', 'IC86.2015']]
for config, ax in zip(configs, axarr.flatten()):
    print('On {}...'.format(config))
    kwargs_relint = {'smooth': 20, 'scale': None, 'decmax': -55, 'low_energy': True}

    relint_light = anisotropy.get_map(config=config, name='relint', composition='light', **kwargs_relint)
    relint_heavy = anisotropy.get_map(config=config, name='relint', composition='heavy', **kwargs_relint)
    relerr_light = anisotropy.get_map(config=config, name='relerr', composition='light', **kwargs_relint)
    relerr_heavy = anisotropy.get_map(config=config, name='relerr', composition='heavy', **kwargs_relint)

    ri_light, ri_err_light, ra, _ = anisotropy.get_proj_relint(relint_light, relerr_light, n_bins=24)
    ri_heavy, ri_err_heavy, ra, _ = anisotropy.get_proj_relint(relint_heavy, relerr_heavy, n_bins=24)
    
    ax.errorbar(ra, ri_light, yerr=ri_err_light, label='light', marker='.', ls=':')
    ax.errorbar(ra, ri_heavy, yerr=ri_err_heavy, label='heavy', marker='.', ls=':')
    ax.axhline(0, marker='None', ls='-.', c='k')
    
    popt_light, perr, _ = anisotropy.get_proj_fit_params(np.deg2rad(ra), ri_light, sigmay=ri_err_light)
    popt_heavy, perr, _ = anisotropy.get_proj_fit_params(np.deg2rad(ra), ri_heavy, sigmay=ri_err_heavy)
    ax.plot(ra, anisotropy.cos_fit_func(np.deg2rad(ra), *popt_light[:3]), marker='None', color=color_dict['light'])
    ax.plot(ra, anisotropy.cos_fit_func(np.deg2rad(ra), *popt_heavy[:3]), marker='None', color=color_dict['heavy'])
    
    chi2_data = calc_chi2(ri_light, ri_heavy, ri_err_light, ri_err_heavy)
    ax.set_title('{} years: '.format(len(config))+'$\mathrm{\chi^2_{red}} =$'+ ' {:0.1f}'.format(chi2_data))
    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 }$')
    ax.grid()
    ax.legend()
# fig.text(0.5, 0.0, 'RA [ $^{\circ}$]', ha='center', fontsize=16)
# fig.text(0.0, 0.5, '$\mathrm{\langle RI \\rangle }$', va='center', rotation='vertical', fontsize=16)
plt.tight_layout()
plt.savefig(os.path.join(comp.paths.figures_dir, 'light-heavy-split_compare-running-sum.png'))
plt.show()


On ['IC86.2011']...
On ['IC86.2011', 'IC86.2012']...
On ['IC86.2011', 'IC86.2012', 'IC86.2013']...
On ['IC86.2011', 'IC86.2012', 'IC86.2013', 'IC86.2014']...
/home/jbourbeau/cr-composition/comptools/anisotropy/anisotropy.py:188: RuntimeWarning: invalid value encountered in sqrt
  perr = np.sqrt(np.diag(pcov))
On ['IC86.2011', 'IC86.2012', 'IC86.2013', 'IC86.2014', 'IC86.2015']...

In [14]:
fig, ax = plt.subplots()
chi2_max = 25
chi2_bins = np.linspace(0, 30, 100)
configs = [['IC86.2012'], ['IC86.2012', 'IC86.2013'], ['IC86.2012', 'IC86.2013', 'IC86.2014']]
for config, color in zip(configs, ['C0', 'C1', 'C2']):
    print('On {}...'.format(config))
    df = anisotropy.get_test_stats(config=config, low_energy=True)
    counts = np.histogram(df.chi2, bins=chi2_bins)[0]
    ax = plotting.plot_steps(chi2_bins, counts, yerr=np.sqrt(counts), label='{} years'.format(len(config)), color=color, ax=ax)
    chi2_data = get_chi2_data(config)
    ax.axvline(chi2_data, marker='None', ls='-.', color=color, lw=1.5)
#                label='Light-heavy split \n ($\mathrm{\chi^2_{red}} = '+'${:0.2f})'.format(chi2_data))
ax.axvline(6.1, marker='None', ls='-.', color=color, lw=1.5)
ax.set_xlabel('$\mathrm{\chi^2_{red}}$')
ax.set_ylabel('Counts')
ax.set_title('1000 random trials')
#     ax.text(1, 5, 'p-value $= {:g}$'.format(pval) + \
#                  '\n significance $= {:0.2f}\sigma$'.format(significance))
ax.set_ylim((1e-1, 1e2))
ax.set_xlim(0, chi2_max)
ax.set_yscale("log", nonposy='clip')
ax.legend()
ax.grid()
# fig.text(0.5, 0.0, 'RA [ $^{\circ}$]', ha='center', fontsize=16)
# fig.text(0.0, 0.5, '$\mathrm{\langle RI \\rangle }$', va='center', rotation='vertical', fontsize=16)
plt.tight_layout()
plt.savefig(os.path.join(comp.paths.figures_dir, 'trials_compare-running-sum.png'))
plt.show()


On ['IC86.2012']...
Loading /data/user/jbourbeau/composition/anisotropy_random_trials/teststat_2012_lowenergy.hdf...
On ['IC86.2012', 'IC86.2013']...
Loading /data/user/jbourbeau/composition/anisotropy_random_trials/teststat_2012_2013_lowenergy.hdf...
On ['IC86.2012', 'IC86.2013', 'IC86.2014']...
Loading /data/user/jbourbeau/composition/anisotropy_random_trials/teststat_2012_2013_2014_lowenergy.hdf...

In [31]:
plt.plot(counts.cumsum())
plt.xlabel('')
plt.grid()



In [31]:
n_bins


Out[31]:
36

In [33]:
np.sum(test_stat > chi2_data)


Out[33]:
19

In [8]:
map_dir = '/data/user/jbourbeau/composition/IC86.2013_data/anisotropy/random_trials/'
sample_0_file_pattern = os.path.join(map_dir, 'random_split_0_trial-*.fits')
sample_1_file_pattern = os.path.join(map_dir, 'random_split_1_trial-*.fits')
infiles_sample_0 = sorted(glob.glob(sample_0_file_pattern))
infiles_sample_1 = sorted(glob.glob(sample_1_file_pattern))

# kwargs_relint = {'smooth': 20, 'scale': None, 'decmax': -55}
n_trials = 2

In [66]:
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)

    kwargs_relint = {'smooth': smooth_rad, 'scale': None, 'decmax': -55}

    relint_0 = anisotropy.get_map(files=trial_files_0, name='relint', **kwargs_relint)
    relint_1 = anisotropy.get_map(files=trial_files_1, name='relint', **kwargs_relint)

    relerr_0 = anisotropy.get_map(files=trial_files_0, name='relerr', **kwargs_relint)
    relerr_1 = anisotropy.get_map(files=trial_files_1, name='relerr', **kwargs_relint)

    ri_0, ri_0_err, ra, _ = anisotropy.get_proj_relint(relint_0, relerr_0, n_bins=n_bins)
    ri_1, ri_1_err, ra, _ = anisotropy.get_proj_relint(relint_1, relerr_1, n_bins=n_bins)

    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

    chi2 = calc_chi2(ri_0, ri_1, ri_0_err, ri_1_err)
    proj_dict['chi2'] = chi2
        
    return proj_dict

In [67]:
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 [68]:
with ProgressBar() as bar:
# with dask.set_options(pool=mp.Pool(10)) as _, ProgressBar() as bar:
    proj_df = proj_df.compute(get=multiprocessing.get, num_workers=10)

proj_df.head()


[########################################] | 100% Completed | 13.2s
Out[68]:
chi2 ri_0 ri_0_err ri_1 ri_1_err
trial
0 0.969489 [0.00321165542178, 0.000115269095227, 0.002368... [0.00173124362307, 0.00173141300553, 0.0017505... [0.00199736147034, 0.00289491003121, 0.0006607... [0.00210537594418, 0.00215605149437, 0.0021451...
1 1.122090 [0.00176563989164, 0.00219092747269, 0.0031199... [0.00172009168706, 0.00174824950129, 0.0017579... [0.00510266249962, 0.00116202822101, 0.0008925... [0.00214209904845, 0.00213806157986, 0.0021452...
2 1.061886 [0.00432630999139, 0.00424003554803, 0.0032409... [0.00172987657282, 0.00175660827852, 0.0017526... [0.00263422364449, -0.000793998166798, 0.00049... [0.00212915307852, 0.00212418977155, 0.0021524...
3 0.920603 [0.00314843476466, 0.00354549898345, 0.0042647... [0.0017287056899, 0.00175922691732, 0.00176760... [0.00461875620905, 0.00115157374417, -0.000173... [0.00213539922731, 0.00214522605925, 0.0021455...
4 1.293273 [0.00211656755968, 0.000759990585587, 0.001169... [0.00172794015787, 0.00174516104541, 0.0017496... [0.00623097487189, 0.0039378548647, 0.00177354... [0.00213115860417, 0.0021501422825, 0.00215605...

In [69]:
# n_side = np.sqrt(n_trials).astype(int)
n_side = 4
# n_bins = 72
# fig, ax = plt.subplots(n, n, figsize=(10, 8), sharex=True, sharey=True)
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(comp.paths.figures_dir, 'random-trials-grid.png'))
plt.show()



In [33]:
n_trials = 16
sig_trials_0 = []
sig_trials_1 = []
for trial_num in pyprind.prog_bar(range(n_trials)):
    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)

    kwargs_relint = {'smooth': smooth_rad, 'scale': None, 'decmax': -55}
#     kwargs_relint = {'smooth': 20, 'scale': None, 'decmax': -55}

    # Get 2D Li-Ma significance maps
    sig_0 = anisotropy.get_map(files=trial_files_0, name='sig', **kwargs_relint)
    sig_1 = anisotropy.get_map(files=trial_files_1, name='sig', **kwargs_relint)
    
    sig_trials_0.append(sig_0)
    sig_trials_1.append(sig_1)


0% [################] 100% | ETA: 00:00:00
Total time elapsed: 00:00:46

In [34]:
n_side = 4
# fig, ax = plt.subplots(n, n, figsize=(10, 8), sharex=True, sharey=True)
fig, axarr = plt.subplots(n_side, n_side, figsize=(10, 8), sharex=True, sharey=False)
for trial, ax in enumerate(axarr.flatten()):
    trial = int(trial)
    print('On trial {}...'.format(trial))
    
    # Get 2D Li-Ma significance maps
    sig_0 = sig_trials_0[trial]
    sig_1 = sig_trials_1[trial]
    
#     sig_0 = np.abs(sig_trials_0[trial])
#     sig_1 = np.abs(sig_trials_1[trial])
    
    # Construct masks
    is_good_mask_0 = (sig_0 != hp.UNSEEN) & ~np.isnan(sig_0)
    is_good_mask_1 = (sig_1 != hp.UNSEEN) & ~np.isnan(sig_1)

    bins = np.linspace(-5, 5, 50)
#     bins = np.linspace(0, 5, 50)
    counts_0, _, _ = ax.hist(sig_0[is_good_mask_0], bins=bins, alpha=0.6, color='C3', label='Split 0')
    counts_1, _, _ = ax.hist(sig_1[is_good_mask_1], bins=bins, alpha=0.6, color='C4', label='Split 1')
    ks_statistic, pval = ks_2samp(sig_0[is_good_mask_0], sig_1[is_good_mask_1])
    ax.set_title('KS p-value = {:0.1e}'.format(pval))

#     ax.set_title('$\mathrm{\chi^2_{red}}$'+' = {:0.2f}'.format(chi2_trials[trial]))
    ax.grid()
#     ax.set_xlim(0, 5)
    #     ax.legend()

    for i, row in enumerate(axarr):
        for j, cell in enumerate(row):
            if i == len(axarr) - 1:
                cell.set_xlabel('$\mathrm{\sigma_{LM}}$')
            if j == 0:
                cell.set_ylabel('Counts')

plt.tight_layout()
    # fig.text(0.5, 0.00, 'RA [ $^{\circ}$]', ha='center')
    # fig.text(0.00, 0.5, '$\mathrm{\langle RI \\rangle }$', va='center', rotation='vertical')
plt.savefig(os.path.join(comp.paths.figures_dir, 'random-trials-grid.png'))
plt.show()


On trial 0...
On trial 1...
On trial 2...
On trial 3...
On trial 4...
On trial 5...
On trial 6...
On trial 7...
On trial 8...
On trial 9...
On trial 10...
On trial 11...
On trial 12...
On trial 13...
On trial 14...
On trial 15...

In [141]:
n_side = 4
# fig, ax = plt.subplots(n, n, figsize=(10, 8), sharex=True, sharey=True)
fig, axarr = plt.subplots(n_side, n_side, figsize=(10, 8), sharex=True, sharey=False)
for trial, ax in enumerate(axarr.flatten()):
    trial = int(trial)
    print('On trial {}...'.format(trial))
    
    # Get 2D Li-Ma significance maps
    sig_0 = sig_trials_0[trial]
    sig_1 = sig_trials_1[trial]
    
#     sig_0 = np.abs(sig_trials_0[trial])
#     sig_1 = np.abs(sig_trials_1[trial])
    
    # Construct masks
    is_good_mask_0 = (sig_0 != hp.UNSEEN) & ~np.isnan(sig_0)
    is_good_mask_1 = (sig_1 != hp.UNSEEN) & ~np.isnan(sig_1)

#     bins = np.linspace(0, 5, 50)
    bins = np.linspace(-5, 5, 50)
    counts_0 = np.histogram(sig_0[is_good_mask_0], bins=bins)[0]
    counts_1 = np.histogram(sig_1[is_good_mask_1], bins=bins)[0]
    
    ax.plot(np.cumsum(counts_0), alpha=0.7, color='C2', label='Split 0')
    ax.plot(np.cumsum(counts_1), alpha=0.7, color='C3', label='Split 1')
    
    ax.set_title(np.abs(np.cumsum(counts_0) - np.cumsum(counts_1)).sum())


#     ax.set_title('$\mathrm{\chi^2_{red}}$'+' = {:0.2f}'.format(chi2_trials[trial]))
    ax.grid()
#     ax.set_xlim(0, 5)
    #     ax.legend()

    for i, row in enumerate(axarr):
        for j, cell in enumerate(row):
            if i == len(axarr) - 1:
                cell.set_xlabel('$\mathrm{\sigma_{LM}}$')
            if j == 0:
                cell.set_ylabel('Counts')

plt.tight_layout()
    # fig.text(0.5, 0.00, 'RA [ $^{\circ}$]', ha='center')
    # fig.text(0.00, 0.5, '$\mathrm{\langle RI \\rangle }$', va='center', rotation='vertical')
#     plt.savefig(os.path.join(comp.paths.figures_dir, 'random-trials-grid.png'))
plt.show()


On trial 0...
On trial 1...
On trial 2...
On trial 3...
On trial 4...
On trial 5...
On trial 6...
On trial 7...
On trial 8...
On trial 9...
On trial 10...
On trial 11...
On trial 12...
On trial 13...
On trial 14...
On trial 15...

In [77]:
np.cumsum(counts_0)


Out[77]:
array([  190.,   399.,   627.,   829.,  1038.,  1263.,  1474.,  1725.,
        1959.,  2192.,  2444.,  2667.,  2887.,  3087.,  3254.,  3412.,
        3562.,  3697.,  3822.,  3950.,  4074.,  4184.,  4289.,  4371.,
        4428.,  4471.,  4501.,  4510.,  4512.,  4512.,  4512.,  4512.,
        4512.,  4512.,  4512.,  4512.,  4512.,  4512.,  4512.,  4512.,
        4512.,  4512.,  4512.,  4512.,  4512.,  4512.,  4512.,  4512.,
        4512.])

In [9]:
kwargs_relint = {'smooth': 20, 'scale': None, 'decmax': -55}

file_0 = infiles_sample_0[501]
print(file_0)
file_1 = infiles_sample_1[501]
print(file_1)

relint_0 = anisotropy.get_map(files=file_0, name='relint', **kwargs_relint)
relint_1 = anisotropy.get_map(files=file_1, name='relint', **kwargs_relint)

relerr_0 = anisotropy.get_map(files=file_0, name='relerr', **kwargs_relint)
relerr_1 = anisotropy.get_map(files=file_1, name='relerr', **kwargs_relint)

ri_0, ri_err_0, ra, ra_err = anisotropy.get_proj_relint(relint_0, relerr_0, n_bins=24)
ri_1, ri_err_1, ra, ra_err = anisotropy.get_proj_relint(relint_1, relerr_1, n_bins=24)

plt.errorbar(ra, ri_0, yerr=ri_err_0)
plt.errorbar(ra, ri_1, yerr=ri_err_1)
plt.grid()
plt.show()


/data/user/jbourbeau/composition/IC86.2013_data/anisotropy/random_trials/random_split_0_trial-55.fits
/data/user/jbourbeau/composition/IC86.2013_data/anisotropy/random_trials/random_split_1_trial-55.fits

In [4]:
kwargs_data = {'config': config, 'low_energy': low_energy, 'smooth': False, 'scale': False, 'decmax': decmax}

In [14]:
n_total = anisotropy.get_num_events(config=config, composition='all', decmax=decmax)
n_light = anisotropy.get_num_events(config=config, composition='light', decmax=decmax)
n_heavy = anisotropy.get_num_events(config=config, composition='heavy', decmax=decmax)

In [15]:
n_light/n_total, n_heavy/n_total


Out[15]:
(0.61598432615372278, 0.38401543147863343)

In [69]:
data_heavy = anisotropy.get_map('data', composition='heavy', **kwargs_data)
data_heavy[ data_heavy < 0 ] = 0
'{:g}'.format(data_heavy.sum())


Out[69]:
'1.1091e+07'

In [36]:
kwargs_relint = {'config': config, 'low_energy': low_energy, 'smooth': smooth_rad, 'scale': scale, 'decmax': decmax}

In [37]:
relint_all = anisotropy.get_map('relint', composition='all', **kwargs_relint)
relint_light = anisotropy.get_map('relint', composition='light', **kwargs_relint)
relint_heavy = anisotropy.get_map('relint', composition='heavy', **kwargs_relint)

In [38]:
print(rel_int_all.max())
print(rel_int_all[rel_int_all != hp.UNSEEN].min())
title = 'Relative Intensity [$\\times \ 10^{}$]'.format('{'+str(-scale)+'}')
fig, ax = anisotropy.plot_skymap(relint_all, smooth=smooth_rad, decmax=decmax, scale=scale,
                                 color_palette='RdBu_r', symmetric=True, llabel='IC86 2012-2015', 
                                 cbar_title=title, cbar_min=-2.5, cbar_max=2.5, polar=True)
outfile = 'IC86-{}_relint_{}_nside-{}_smooth-{:0.1f}.png'.format(years_str, 'all', n_side, smooth_rad)
plt.savefig(os.path.join(comp.paths.figures_dir, outfile))


2.8661449384
-3.43987296025

In [39]:
print(rel_int_light.max())
print(rel_int_light[rel_int_light != hp.UNSEEN].min())
title = 'Relative Intensity [$\\times \ 10^{}$]'.format('{'+str(-scale)+'}')
fig, ax = anisotropy.plot_skymap(relint_light, smooth=smooth_rad, decmax=decmax, scale=scale,
                                 color_palette='RdBu_r', symmetric=True, llabel='IC86 2012-2015', 
                                 cbar_title=title, cbar_min=-2.5, cbar_max=2.5, polar=True)
outfile = 'IC86-{}_relint_{}_nside-{}_smooth-{:0.1f}.png'.format(years_str, 'light', n_side, smooth_rad)
plt.savefig(os.path.join(comp.paths.figures_dir, outfile))


4.28225073662
-3.7855295012

In [40]:
print(rel_int_heavy.max())
print(rel_int_heavy[rel_int_heavy != hp.UNSEEN].min())
title = 'Relative Intensity [$\\times \ 10^{}$]'.format('{'+str(-scale)+'}')
fig, ax = anisotropy.plot_skymap(relint_heavy, smooth=smooth_rad, decmax=decmax, scale=scale,
                                 color_palette='RdBu_r', symmetric=True, llabel='IC86 2012-2015', 
                                 cbar_title=title, cbar_min=-2.5, cbar_max=2.5, polar=True)
outfile = 'IC86-{}_relint_{}_nside-{}_smooth-{:0.1f}.png'.format(years_str, 'heavy', n_side, smooth_rad)
plt.savefig(os.path.join(comp.paths.figures_dir, outfile))


2.59085836656
-3.28462140504

In [46]:
kwargs_sig = {'config': config, 'low_energy': low_energy, 'smooth': smooth_rad, 'scale': None, 'decmax': decmax}

In [47]:
sig_all = anisotropy.get_map('sig', composition='all', **kwargs_sig)
sig_light = anisotropy.get_map('sig', composition='light', **kwargs_sig)
sig_heavy = anisotropy.get_map('sig', composition='heavy', **kwargs_sig)

In [48]:
print(sig_all.max())
print(sig_all[sig_all != hp.UNSEEN].min())
title = 'Significance [$\mathrm{\sigma_{LM}}$]'
fig, ax = anisotropy.plot_skymap(sig_all, smooth=smooth_rad, decmax=decmax, scale=scale,
                                 color_palette='RdBu_r', symmetric=True, llabel='IC86 2012-2015', 
                                 cbar_title=title, cbar_min=-4.1, cbar_max=4.1, polar=True)
outfile = 'IC86-{}_sig_{}_nside-{}_smooth-{:0.1f}.png'.format(years_str, 'all', n_side, smooth_rad)
plt.savefig(os.path.join(comp.paths.figures_dir, outfile))


4.74007594155
-5.26517211233

In [49]:
print(sig_light.max())
print(sig_light[sig_light != hp.UNSEEN].min())
title = 'Significance [$\mathrm{\sigma_{LM}}$]'
fig, ax = anisotropy.plot_skymap(sig_light, smooth=smooth_rad, decmax=decmax, scale=scale,
                                 color_palette='RdBu_r', symmetric=True, llabel='IC86 2012-2015', 
                                 cbar_title=title, cbar_min=-4.1, cbar_max=4.1, polar=True)
outfile = 'IC86-{}_sig_{}_nside-{}_smooth-{:0.1f}.png'.format(years_str, 'light', n_side, smooth_rad)
plt.savefig(os.path.join(comp.paths.figures_dir, outfile))


5.55383911119
-4.48361748687

In [50]:
print(sig_heavy.max())
print(sig_heavy[sig_heavy != hp.UNSEEN].min())
title = 'Significance [$\mathrm{\sigma_{LM}}$]'
fig, ax = anisotropy.plot_skymap(sig_heavy, smooth=smooth_rad, decmax=decmax, scale=scale,
                                 color_palette='RdBu_r', symmetric=True, llabel='IC86 2012-2015', 
                                 cbar_title=title, cbar_min=-4.1, cbar_max=4.1, polar=True)
outfile = 'IC86-{}_sig_{}_nside-{}_smooth-{:0.1f}.png'.format(years_str, 'heavy', n_side, smooth_rad)
plt.savefig(os.path.join(comp.paths.figures_dir, outfile))


3.15578385989
-2.85770816159

In [30]:
rel_int_diff = anisotropy.get_relint_diff(**kwargs_relint)
print(rel_int_diff.max())
print(rel_int_diff[rel_int_diff != hp.UNSEEN].min())


2.38421078863
-2.57748060774

In [31]:
title = '$\mathrm{\Delta \ RI}$'
fig, ax = anisotropy.plot_skymap(rel_int_diff, smooth=smooth_rad, decmax=decmax, scale=scale,
                                 color_palette='BrBG', symmetric=True, llabel='IC86 2012-2013', 
                                 cbar_title=title, cbar_min=-2.5, cbar_max=2.5, polar=True)
outfile = 'IC86-{}_RI-diff_nside-{}_smooth-{:0.1f}.png'.format(years_str, n_side, smooth_rad)
plt.savefig(os.path.join(comp.paths.figures_dir, outfile))


Projected relative intensity

[ back to top ]


In [4]:
kwargs_relint = {'config': config, 'low_energy': low_energy, 'smooth': smooth_rad, 'scale': None, 'decmax': decmax}

In [5]:
relint_all = anisotropy.get_map('relint', composition='all', **kwargs_relint)
relint_all_err = anisotropy.get_map('relerr', composition='all', **kwargs_relint)

In [80]:
ri_all, ri_all_err, ra, ra_err = anisotropy.get_proj_relint(relint_all, relint_all_err, n_bins=30)

In [81]:
fig, ax = plt.subplots()
ax.errorbar(ra, ri_all, yerr=ri_all_err, marker='.', ls='None', c='C2', label='all')
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(-2e-3, 2e-3)
ax.invert_xaxis()
ax.legend()
outfile = 'IC86-{}_proj-RI-all_nside-{}_smooth-{:0.1f}.png'.format(years_str, n_side, smooth_rad)
plt.savefig(os.path.join(comp.paths.figures_dir, outfile))
plt.show()



In [82]:
relint_light = anisotropy.get_map('relint', composition='light', **kwargs_relint)
relint_light_err = anisotropy.get_map('relerr', composition='light', **kwargs_relint)
ri_light, ri_light_err, ra, ra_err = anisotropy.get_proj_relint(relint_light, relint_light_err, n_bins=30)

In [83]:
relint_heavy = anisotropy.get_map('relint', composition='heavy', **kwargs_relint)
relint_heavy_err = anisotropy.get_map('relerr', composition='heavy', **kwargs_relint)
ri_heavy, ri_heavy_err, ra, ra_err = anisotropy.get_proj_relint(relint_heavy, relint_heavy_err, n_bins=30)

In [84]:
popt_light, perr_light, chi2 = anisotropy.get_proj_fit_params(np.deg2rad(ra), ri_light, sigmay=ri_light_err, l=3)
popt_heavy, perr_heavy, chi2 = anisotropy.get_proj_fit_params(np.deg2rad(ra), ri_heavy, sigmay=ri_heavy_err, l=3)

In [85]:
fig, ax = plt.subplots()
ax.errorbar(ra, ri_light, yerr=ri_light_err, marker='.', ls='None', label='light')
ax.errorbar(ra, ri_heavy, yerr=ri_heavy_err, marker='.', ls='None', label='heavy')
ax.axhline(0, marker='None', ls='-.', c='k')
ax.set_ylabel('$\mathrm{\langle RI \\rangle }$')
ax.set_xlabel('RA [ $^{\circ}$]')
ax.legend()
ax.grid()
ax.invert_xaxis()
ax.set_ylim(-2e-3, 2e-3)
outfile = 'IC86-{}_proj-RI-comps_nside-{}_smooth-{:0.1f}.png'.format(years_str, n_side, smooth_rad)
plt.savefig(os.path.join(comp.paths.figures_dir, outfile))
plt.show()



In [86]:
fig, ax = plt.subplots()
ax.errorbar(ra, ri_light, yerr=ri_light_err, marker='.', ls='None', label='light')
ax.errorbar(ra, ri_heavy, yerr=ri_heavy_err, marker='.', ls='None', label='heavy')
ax.errorbar(ra, ri_all, yerr=ri_all_err, marker='.', ls='None', label='all')
ax.axhline(0, marker='None', ls='-.', c='k')
ax.set_ylabel('$\mathrm{\langle RI \\rangle }$')
ax.set_xlabel('RA [ $^{\circ}$]')
ax.legend()
ax.grid()
ax.invert_xaxis()
ax.set_ylim(-2e-3, 2e-3)
outfile = 'IC86-{}_proj-RI_nside-{}_smooth-{:0.1f}.png'.format(years_str, n_side, smooth_rad)
plt.savefig(os.path.join(comp.paths.figures_dir, outfile))
plt.show()



In [87]:
fig, ax = plt.subplots()
ax.errorbar(ra, ri_light, yerr=ri_light_err, marker='.', ls='None', label='light')
ax.errorbar(ra, ri_heavy, yerr=ri_heavy_err, marker='.', ls='None', label='heavy')
ax.errorbar(ra, ri_all, yerr=ri_all_err, marker='.', ls='None', label='all')
ax.axhline(0, marker='None', ls='-.', c='k')
ax.plot(ra, anisotropy.cos_fit_func(np.deg2rad(ra), *popt_light), color='C0', marker='None')
ax.plot(ra, anisotropy.cos_fit_func(np.deg2rad(ra), *popt_heavy), color='C1', marker='None')
ax.set_ylabel('$\mathrm{\langle RI \\rangle }$')
ax.set_xlabel('RA [ $^{\circ}$]')
ax.legend()
ax.grid()
ax.invert_xaxis()
ax.set_ylim(-2e-3, 2e-3)
outfile = 'IC86-{}_proj-RI-fit_nside-{}_smooth-{:0.1f}.png'.format(years_str, n_side, smooth_rad)
plt.savefig(os.path.join(comp.paths.figures_dir, outfile))
plt.show()



In [88]:
amp_light = popt_light[1]
amp_heavy = popt_heavy[1]
amp_err_light = perr_light[1]
amp_err_heavy = perr_heavy[1]

phase_light = np.rad2deg(popt_light[2])
phase_heavy = np.rad2deg(popt_heavy[2])
phase_err_light = np.rad2deg(perr_light[2])
phase_err_heavy = np.rad2deg(perr_heavy[2])

In [89]:
fig, ax = plt.subplots()
ax.errorbar(ra, ri_light, yerr=ri_light_err, marker='.', ls='None', label='light')
ax.errorbar(ra, ri_heavy, yerr=ri_heavy_err, marker='.', ls='None', label='heavy')
ax.errorbar(ra, ri_all, yerr=ri_all_err, marker='.', ls='None', label='all')
ax.axhline(0, marker='None', ls='-.', c='k')
ax.plot(ra, anisotropy.cos_fit_func(np.deg2rad(ra), *popt_light[:3]), color='C0', marker='None')
ax.plot(ra, anisotropy.cos_fit_func(np.deg2rad(ra), *popt_heavy[:3]), color='C1', marker='None')
ax.set_ylabel('$\mathrm{\langle RI \\rangle }$')
ax.set_xlabel('RA [ $^{\circ}$]')
ax.legend()
ax.grid()
ax.invert_xaxis()
ax.set_ylim(-2e-3, 2e-3)

light_amp_str = 'light amp = {:0.1g} +/- {:0.1g}'.format(amp_light, amp_err_light)
light_phase_str = 'light phase = {:0.2f} {} +/- {:0.2f} {}'.format(phase_light, '$^{\circ}$', phase_err_light, '$^{\circ}$')
ax.text(350, 0.00225, light_amp_str + '\n' + light_phase_str)

heavy_amp_str = 'heavy amp = {:0.1g} +/- {:0.1g}'.format(amp_heavy, amp_err_heavy)
heavy_phase_str = 'heavy phase = {:0.2f} {} +/- {:0.2f} {}'.format(phase_heavy, '$^{\circ}$', phase_err_heavy, '$^{\circ}$')
ax.text(175, 0.00225, heavy_amp_str + '\n' + heavy_phase_str)

outfile = 'IC86-{}_proj-RI-dipolefit_nside-{}_smooth-{:0.1f}.png'.format(years_str, n_side, smooth_rad)
plt.savefig(os.path.join(comp.paths.figures_dir, outfile))
plt.show()


Cross-check: random splitting of events

[ back to top ]


In [66]:
config = ['IC86.2012', 'IC86.2013', 'IC86.2014', 'IC86.2015']
n_side = 64
smooth_rad = 20.0
scale = 3
decmax = -55
figures_dir = '/home/jbourbeau/public_html/figures/composition-anisotropy/cross-check-random-split'
years_str = '2012-2015'

In [67]:
rel_int_0 = anisotropy.get_map('relint', config=config, composition='random_0',
                             smooth=smooth_rad, scale=scale, decmax=decmax)
rel_int_1 = anisotropy.get_map('relint', config=config, composition='random_1',
                             smooth=smooth_rad, scale=scale, decmax=decmax)
print(rel_int_0.max())
print(rel_int_0[rel_int_0 != hp.UNSEEN].min())

print(rel_int_1.max())
print(rel_int_1[rel_int_1 != hp.UNSEEN].min())


2.65429851244
-3.36021713829
2.58918508496
-3.14294109658

In [68]:
title = 'Relative Intensity [$\\times \ 10^{}$]'.format('{'+str(-scale)+'}')
fig, ax = anisotropy.plot_skymap(rel_int_0, smooth=smooth_rad, decmax=decmax, scale=scale,
                                 color_palette='RdBu_r', symmetric=True, llabel='IC86 2012-2013', 
                                 cbar_title=title, cbar_min=-2.5, cbar_max=2.5, polar=True)
outfile = 'IC86-{}_relint_{}_nside-{}_smooth-{:0.1f}.png'.format(years_str, 'random_0', n_side, smooth_rad)
plt.savefig(os.path.join(figures_dir, outfile))



In [69]:
title = 'Relative Intensity [$\\times \ 10^{}$]'.format('{'+str(-scale)+'}')
fig, ax = anisotropy.plot_skymap(rel_int_1, smooth=smooth_rad, decmax=decmax, scale=scale,
                                 color_palette='RdBu_r', symmetric=True, llabel='IC86 2012-2013', 
                                 cbar_title=title, cbar_min=-2.5, cbar_max=2.5, polar=True)
outfile = 'IC86-{}_relint_{}_nside-{}_smooth-{:0.1f}.png'.format(years_str, 'random_1', n_side, smooth_rad)
plt.savefig(os.path.join(figures_dir, outfile))



In [70]:
sig_0 = anisotropy.get_map('sig', composition='random_0', config=config, smooth=smooth_rad, decmax=decmax)
sig_1 = anisotropy.get_map('sig', composition='random_1', config=config, smooth=smooth_rad, decmax=decmax)

In [71]:
print(sig_0.max())
print(sig_0[sig_0 != hp.UNSEEN].min())

print(sig_1.max())
print(sig_1[sig_1 != hp.UNSEEN].min())


3.10673742569
-3.66635747678
3.68417830318
-3.60931721794

In [72]:
title = 'Significance [$\mathrm{\sigma_{LM}}$]'
fig, ax = anisotropy.plot_skymap(sig_0, smooth=smooth_rad, decmax=decmax, scale=scale,
                                 color_palette='RdBu_r', symmetric=True, llabel='IC86 2012-2013', 
                                 cbar_title=title, cbar_min=-3.5, cbar_max=3.5, polar=True)
outfile = 'IC86-{}_sig_{}_nside-{}_smooth-{:0.1f}.png'.format(years_str, 'random_0', n_side, smooth_rad)
plt.savefig(os.path.join(figures_dir, outfile))



In [73]:
title = 'Significance [$\mathrm{\sigma_{LM}}$]'
fig, ax = anisotropy.plot_skymap(sig_1, smooth=smooth_rad, decmax=decmax, scale=scale,
                                 color_palette='RdBu_r', symmetric=True, llabel='IC86 2012-2013', 
                                 cbar_title=title, cbar_min=-3.5, cbar_max=3.5, polar=True)
outfile = 'IC86-{}_sig_{}_nside-{}_smooth-{:0.1f}.png'.format(years_str, 'random_1', n_side, smooth_rad)
plt.savefig(os.path.join(figures_dir, outfile))



In [74]:
relint = anisotropy.get_map('relint', config=config, composition='all',
                             smooth=smooth_rad, scale=None, decmax=decmax)
relint_err = anisotropy.get_map('relerr', config=config, composition='all',
                             smooth=smooth_rad, scale=None, decmax=decmax)

In [75]:
ri_all, ri_all_err, ra, ra_err = anisotropy.get_proj_relint(relint, relint_err)

In [76]:
fig, ax = plt.subplots()
ax.errorbar(ra, ri_all, yerr=ri_all_err, marker='.', ls='None')
ax.axhline(0, marker='None', ls='-.', c='k')
ax.set_ylabel('$\mathrm{\langle RI \\rangle }$')
ax.set_xlabel('RA [ $^{\circ}$]')
ax.grid()
ax.invert_xaxis()
plt.show()



In [91]:
relint_0 = anisotropy.get_map('relint', config=config, composition='random_0',
                             smooth=smooth_rad, scale=None, decmax=decmax)
relint_err_0 = anisotropy.get_map('relerr', config=config, composition='random_0',
                             smooth=smooth_rad, scale=None, decmax=decmax)
ri_0, ri_err_0, ra, ra_err = anisotropy.get_proj_relint(relint_0, relint_err_0)

In [92]:
relint_1 = anisotropy.get_map('relint', config=config, composition='random_1',
                             smooth=smooth_rad, scale=None, decmax=decmax)
relint_err_1 = anisotropy.get_map('relerr', config=config, composition='random_1',
                             smooth=smooth_rad, scale=None, decmax=decmax)
ri_1, ri_err_1, ra, ra_err = anisotropy.get_proj_relint(relint_1, relint_err_1)

In [93]:
popt_0, perr, chi2 = anisotropy.get_proj_fit_params(np.deg2rad(ra), ri_0, sigmay=ri_err_0)
popt_1, perr, chi2 = anisotropy.get_proj_fit_params(np.deg2rad(ra), ri_1, sigmay=ri_err_1)

In [94]:
fig, ax = plt.subplots()
ax.errorbar(ra, ri_0, yerr=ri_err_0, marker='.', ls='None', label='Subset 0')
ax.errorbar(ra, ri_1, yerr=ri_err_1, marker='.', ls='None', label='Subset 1')
ax.axhline(0, marker='None', ls='-.', c='k')
ax.set_ylabel('$\mathrm{\langle RI \\rangle }$')
ax.set_xlabel('RA [ $^{\circ}$]')
ax.legend()
ax.grid()
ax.invert_xaxis()
plt.show()



In [95]:
fig, ax = plt.subplots()
ax.errorbar(ra, ri_0, yerr=ri_err_0, marker='.', ls='None', label='Subset 0')
ax.errorbar(ra, ri_1, yerr=ri_err_1, marker='.', ls='None', label='Subset 1')
ax.errorbar(ra, ri_all, yerr=ri_all_err, marker='.', ls='None', label='all')
ax.axhline(0, marker='None', ls='-.', c='k')
ax.set_ylabel('$\mathrm{\langle RI \\rangle }$')
ax.set_xlabel('RA [ $^{\circ}$]')
ax.legend()
ax.grid()
ax.invert_xaxis()
plt.show()



In [96]:
fig, ax = plt.subplots()
ax.errorbar(ra, ri_0, yerr=ri_err_0, marker='.', ls='None', label='Subset 0')
ax.errorbar(ra, ri_1, yerr=ri_err_1, marker='.', ls='None', label='Subset 1')
ax.errorbar(ra, ri_all, yerr=ri_all_err, marker='.', ls='None', label='all')
ax.axhline(0, marker='None', ls='-.', c='k')
ax.plot(ra, anisotropy.cos_fit_func(np.deg2rad(ra), *popt_0), color='C0', marker='None')
ax.plot(ra, anisotropy.cos_fit_func(np.deg2rad(ra), *popt_1), color='C1', marker='None')
ax.set_ylabel('$\mathrm{\langle RI \\rangle }$')
ax.set_xlabel('RA [ $^{\circ}$]')
ax.legend()
ax.grid()
ax.invert_xaxis()
plt.show()



In [97]:
fig, ax = plt.subplots()
ax.errorbar(ra, ri_0, yerr=ri_err_0, marker='.', ls='None', label='Subset 0')
ax.errorbar(ra, ri_1, yerr=ri_err_1, marker='.', ls='None', label='Subset 1')
ax.errorbar(ra, ri_all, yerr=ri_all_err, marker='.', ls='None', label='all')
ax.axhline(0, marker='None', ls='-.', c='k')
ax.plot(ra, anisotropy.cos_fit_func(np.deg2rad(ra), *popt_0[:3]), color='C0', marker='None')
ax.plot(ra, anisotropy.cos_fit_func(np.deg2rad(ra), *popt_1[:3]), color='C1', marker='None')
ax.set_ylabel('$\mathrm{\langle RI \\rangle }$')
ax.set_xlabel('RA [ $^{\circ}$]')
ax.legend()
ax.grid()
ax.invert_xaxis()
ax.set_ylim(-2e-3, 2e-3)
outfile = 'IC86-{}_proj-RI-dipolefit_random_nside-{}_smooth-{:0.1f}.png'.format(years_str, n_side, smooth_rad)
plt.savefig(os.path.join(figures_dir, outfile))
plt.show()



In [98]:
'{:g}'.format(popt_0[1]), '{:g}'.format(popt_1[1])


Out[98]:
('0.0010178', '0.00113404')

In [99]:
np.rad2deg(popt_0[2]), np.rad2deg(popt_1[2])


Out[99]:
(-71.862836300904505, -61.870525912281074)

In [ ]:


In [3]:
# Load full DataFrame for config
mypaths = comp.get_paths()
df_file = os.path.join(mypaths.comp_data_dir, 'IC86.2013' + '_data',
                       'anisotropy_dataframe.hdf')
with pd.HDFStore(df_file) as store:
    data_df = store['dataframe']

In [6]:
data_df.head()


Out[6]:
lap_zenith lap_azimuth start_time_mjd pred_comp
0 0.206773 5.480972 56414.452644 light
1 0.181282 2.019701 56414.452706 heavy
2 0.220391 3.160242 56414.452744 heavy
3 0.408630 1.879472 56414.452762 light
4 0.061005 1.460350 56414.452762 light

In [5]:
n_side = 64
n_pix = hp.nside2npix(n_side)
bdt_score_skymap = np.zeros(n_pix)

for idx, row in data_df.iterrows():
    
    local_zenith = row['lap_zenith']
    local_azimuth = row['lap_azimuth']
    local_time = row['start_time_mjd']

    ra, dec = astro.dir_to_equa(local_zenith, local_azimuth, local_time)
    hp_theta, hp_phi = anisotropy.equatorial_to_healpy(ra, dec)
    pix = hp.ang2pix(n_side, hp_theta, hp_phi)
    bdt_score_skymap[pix] += row['score']



KeyErrorTraceback (most recent call last)
<ipython-input-5-bbcae491435d> in <module>()
----> 1 data_df.apply(lambda row: astro.dir_to_equa(row['lap_zenith'], row['lap_azimuth'], row['start_time_mjd']))

/home/jbourbeau/.virtualenvs/composition/lib/python2.7/site-packages/pandas/core/frame.pyc in apply(self, func, axis, broadcast, raw, reduce, args, **kwds)
   4358                         f, axis,
   4359                         reduce=reduce,
-> 4360                         ignore_failures=ignore_failures)
   4361             else:
   4362                 return self._apply_broadcast(f, axis)

/home/jbourbeau/.virtualenvs/composition/lib/python2.7/site-packages/pandas/core/frame.pyc in _apply_standard(self, func, axis, ignore_failures, reduce)
   4454             try:
   4455                 for i, v in enumerate(series_gen):
-> 4456                     results[i] = func(v)
   4457                     keys.append(v.name)
   4458             except Exception as e:

<ipython-input-5-bbcae491435d> in <lambda>(row)
----> 1 data_df.apply(lambda row: astro.dir_to_equa(row['lap_zenith'], row['lap_azimuth'], row['start_time_mjd']))

/home/jbourbeau/.virtualenvs/composition/lib/python2.7/site-packages/pandas/core/series.pyc in __getitem__(self, key)
    599         key = com._apply_if_callable(key, self)
    600         try:
--> 601             result = self.index.get_value(self, key)
    602 
    603             if not is_scalar(result):

/home/jbourbeau/.virtualenvs/composition/lib/python2.7/site-packages/pandas/core/indexes/base.pyc in get_value(self, series, key)
   2426         try:
   2427             return self._engine.get_value(s, k,
-> 2428                                           tz=getattr(series.dtype, 'tz', None))
   2429         except KeyError as e1:
   2430             if len(self) > 0 and self.inferred_type in ['integer', 'boolean']:

pandas/_libs/index.pyx in pandas._libs.index.IndexEngine.get_value (pandas/_libs/index.c:4363)()

pandas/_libs/index.pyx in pandas._libs.index.IndexEngine.get_value (pandas/_libs/index.c:4046)()

pandas/_libs/index.pyx in pandas._libs.index.IndexEngine.get_loc (pandas/_libs/index.c:4912)()

KeyError: ('lap_zenith', u'occurred at index lap_zenith')

In [80]:
fig, axarr = plt.subplots(1, 2)



In [82]:
ax = axarr[0]

In [ ]:
ax.s

In [ ]: