In [1]:

    

%load_ext autoreload
%autoreload 2


    

In [63]:

    

from __future__ import division, print_function
import os
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
from matplotlib.colors import LogNorm
import seaborn as sns

import comptools as comp

color_dict = comp.get_color_dict()

sns.set_context(context='paper', font_scale=1.5)
# sns.set_context(context='paper', font_scale=1.75)

%matplotlib inline


    

In [3]:

    

config = 'IC86.2012'
num_groups = 2
comp_list = comp.get_comp_list(num_groups=num_groups)


    

In [4]:

    

df_sim = comp.load_sim(config=config,
                       energy_cut_key='MC_log_energy', 
                       energy_reco=False,
                       log_energy_min=None,
                       log_energy_max=None,
                       test_size=0,
                       verbose=True)


    

    




[                                        ] | 0% Completed |  0.1s





    




/home/jbourbeau/cr-composition/.virtualenv/lib/python2.7/site-packages/dask/base.py:835: UserWarning: The get= keyword has been deprecated. Please use the scheduler= keyword instead with the name of the desired scheduler like 'threads' or 'processes'
  warnings.warn("The get= keyword has been deprecated. "






    




[########################################] | 100% Completed |  8.0s
[########################################] | 100% Completed |  0.1s

In [5]:

    

energy_mask = df_sim['MC_log_energy'] >= 6.4


    

In [50]:

    

h, xedges, yedges = np.histogram2d(df_sim.loc[:, 'log_s125'].values,
                                   df_sim.loc[:, 'MC_log_energy'].values,
#                                    bins=[log_125_bins, log_MC_energy_bins])
                                   bins=(50, 50))


    

In [51]:

    

h = h / h.sum(axis=0)


    

In [52]:

    

h.sum(axis=0)


    

    
Out[52]:





array([1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.])

In [53]:

    

h = np.rot90(h)
h = np.flipud(h)
h = np.ma.masked_where(h == 0, h)


    

In [54]:

    

import matplotlib.colors as colors


    

In [55]:

    

extent = [xedges[0], xedges[-1], yedges[0], yedges[-1]]
extent


    

    
Out[55]:





[-0.7863706656992693, 2.1775389075026235, 5.316134256754441, 7.999821416140386]

In [64]:

    

fig, ax = plt.subplots()
im = ax.imshow(h,
               extent=extent,
               origin='lower',
               interpolation='none',
               cmap='viridis',
               aspect='auto')
plt.colorbar(im, label='$\mathrm{P(E_{true} | S_{125})}$')

ax.set_ylabel('$\mathrm{\log_{10}(E_{true}/GeV)}$')
ax.set_xlabel('$\mathrm{\log_{10}(S_{125})}$')

ax.grid()

outfile = os.path.join(comp.paths.figures_dir, 'laputop_performance', 's125_vs_MC_energy.png')
comp.check_output_dir(outfile)
plt.savefig(outfile)

plt.show()


    

In [25]:

    

fig, ax = plt.subplots()
im = ax.imshow(h_norm, origin='lower',
#                extent=extent,
              )
ax.set_ylabel('$\mathrm{\log_{10}(E_{true}/GeV)}$')
ax.set_xlabel('$\mathrm{\log_{10}(S_{125})}$')
# plt.colorbar(im, label='$\mathrm{P(E_{true} | S_{125})}$')

plt.xticks(log_125_bins)

# outfile = os.path.join(comp.paths.figures_dir, 'laputop_performance', 's125_vs_MC_energy.png')
# comp.check_output_dir(outfile)
# plt.savefig(outfile)

plt.show()


    

    





KeyboardInterruptTraceback (most recent call last)
<ipython-input-25-7c3bcb5e202b> in <module>()
     13 # plt.savefig(outfile)
     14 
---> 15 plt.show()

/home/jbourbeau/cr-composition/.virtualenv/lib/python2.7/site-packages/matplotlib/pyplot.pyc in show(*args, **kw)
    251     """
    252     global _show
--> 253     return _show(*args, **kw)
    254 
    255 

/home/jbourbeau/cr-composition/.virtualenv/lib/python2.7/site-packages/ipykernel/pylab/backend_inline.pyc in show(close, block)
     34     try:
     35         for figure_manager in Gcf.get_all_fig_managers():
---> 36             display(figure_manager.canvas.figure)
     37     finally:
     38         show._to_draw = []

/home/jbourbeau/cr-composition/.virtualenv/lib/python2.7/site-packages/IPython/core/display.pyc in display(*objs, **kwargs)
    309             publish_display_data(data=obj, metadata=metadata, **kwargs)
    310         else:
--> 311             format_dict, md_dict = format(obj, include=include, exclude=exclude)
    312             if not format_dict:
    313                 # nothing to display (e.g. _ipython_display_ took over)

/home/jbourbeau/cr-composition/.virtualenv/lib/python2.7/site-packages/IPython/core/formatters.pyc in format(self, obj, include, exclude)
    171             md = None
    172             try:
--> 173                 data = formatter(obj)
    174             except:
    175                 # FIXME: log the exception

<decorator-gen-9> in __call__(self, obj)

/home/jbourbeau/cr-composition/.virtualenv/lib/python2.7/site-packages/IPython/core/formatters.pyc in catch_format_error(method, self, *args, **kwargs)
    215     """show traceback on failed format call"""
    216     try:
--> 217         r = method(self, *args, **kwargs)
    218     except NotImplementedError:
    219         # don't warn on NotImplementedErrors

/home/jbourbeau/cr-composition/.virtualenv/lib/python2.7/site-packages/IPython/core/formatters.pyc in __call__(self, obj)
    332                 pass
    333             else:
--> 334                 return printer(obj)
    335             # Finally look for special method names
    336             method = get_real_method(obj, self.print_method)

/home/jbourbeau/cr-composition/.virtualenv/lib/python2.7/site-packages/IPython/core/pylabtools.pyc in <lambda>(fig)
    245 
    246     if 'png' in formats:
--> 247         png_formatter.for_type(Figure, lambda fig: print_figure(fig, 'png', **kwargs))
    248     if 'retina' in formats or 'png2x' in formats:
    249         png_formatter.for_type(Figure, lambda fig: retina_figure(fig, **kwargs))

/home/jbourbeau/cr-composition/.virtualenv/lib/python2.7/site-packages/IPython/core/pylabtools.pyc in print_figure(fig, fmt, bbox_inches, **kwargs)
    129 
    130     bytes_io = BytesIO()
--> 131     fig.canvas.print_figure(bytes_io, **kw)
    132     data = bytes_io.getvalue()
    133     if fmt == 'svg':

/home/jbourbeau/cr-composition/.virtualenv/lib/python2.7/site-packages/matplotlib/backend_bases.pyc in print_figure(self, filename, dpi, facecolor, edgecolor, orientation, format, **kwargs)
   2210                     orientation=orientation,
   2211                     dryrun=True,
-> 2212                     **kwargs)
   2213                 renderer = self.figure._cachedRenderer
   2214                 bbox_inches = self.figure.get_tightbbox(renderer)

/home/jbourbeau/cr-composition/.virtualenv/lib/python2.7/site-packages/matplotlib/backends/backend_agg.pyc in print_png(self, filename_or_obj, *args, **kwargs)
    511 
    512     def print_png(self, filename_or_obj, *args, **kwargs):
--> 513         FigureCanvasAgg.draw(self)
    514         renderer = self.get_renderer()
    515         original_dpi = renderer.dpi

/home/jbourbeau/cr-composition/.virtualenv/lib/python2.7/site-packages/matplotlib/backends/backend_agg.pyc in draw(self)
    431             # if toolbar:
    432             #     toolbar.set_cursor(cursors.WAIT)
--> 433             self.figure.draw(self.renderer)
    434             # A GUI class may be need to update a window using this draw, so
    435             # don't forget to call the superclass.

/home/jbourbeau/cr-composition/.virtualenv/lib/python2.7/site-packages/matplotlib/artist.pyc in draw_wrapper(artist, renderer, *args, **kwargs)
     53                 renderer.start_filter()
     54 
---> 55             return draw(artist, renderer, *args, **kwargs)
     56         finally:
     57             if artist.get_agg_filter() is not None:

/home/jbourbeau/cr-composition/.virtualenv/lib/python2.7/site-packages/matplotlib/figure.pyc in draw(self, renderer)
   1464                 try:
   1465                     self.tight_layout(renderer,
-> 1466                                       **self._tight_parameters)
   1467                 except ValueError:
   1468                     pass

/home/jbourbeau/cr-composition/.virtualenv/lib/python2.7/site-packages/matplotlib/figure.pyc in tight_layout(self, renderer, pad, h_pad, w_pad, rect)
   2273         kwargs = get_tight_layout_figure(
   2274             self, self.axes, subplotspec_list, renderer,
-> 2275             pad=pad, h_pad=h_pad, w_pad=w_pad, rect=rect)
   2276         self.subplots_adjust(**kwargs)
   2277 

/home/jbourbeau/cr-composition/.virtualenv/lib/python2.7/site-packages/matplotlib/tight_layout.pyc in get_tight_layout_figure(fig, axes_list, subplotspec_list, renderer, pad, h_pad, w_pad, rect)
    326                                      subplot_list=subplot_list,
    327                                      ax_bbox_list=ax_bbox_list,
--> 328                                      pad=pad, h_pad=h_pad, w_pad=w_pad)
    329 
    330     if rect is not None:

/home/jbourbeau/cr-composition/.virtualenv/lib/python2.7/site-packages/matplotlib/tight_layout.pyc in auto_adjust_subplotpars(fig, renderer, nrows_ncols, num1num2_list, subplot_list, ax_bbox_list, pad, h_pad, w_pad, rect)
    113 
    114         tight_bbox_raw = union([ax.get_tightbbox(renderer) for ax in subplots
--> 115                                 if ax.get_visible()])
    116         tight_bbox = TransformedBbox(tight_bbox_raw,
    117                                      fig.transFigure.inverted())

/home/jbourbeau/cr-composition/.virtualenv/lib/python2.7/site-packages/matplotlib/axes/_base.pyc in get_tightbbox(self, renderer, call_axes_locator)
   4168             bb.append(self._right_title.get_window_extent(renderer))
   4169 
-> 4170         bb_xaxis = self.xaxis.get_tightbbox(renderer)
   4171         if bb_xaxis:
   4172             bb.append(bb_xaxis)

/home/jbourbeau/cr-composition/.virtualenv/lib/python2.7/site-packages/matplotlib/axis.pyc in get_tightbbox(self, renderer)
   1145         ticks_to_draw = self._update_ticks(renderer)
   1146 
-> 1147         self._update_label_position(renderer)
   1148 
   1149         # go back to just this axis's tick labels

/home/jbourbeau/cr-composition/.virtualenv/lib/python2.7/site-packages/matplotlib/axis.pyc in _update_label_position(self, renderer)
   1904         # get bounding boxes for this axis and any siblings
   1905         # that have been set by `fig.align_xlabels()`
-> 1906         bboxes, bboxes2 = self._get_tick_boxes_siblings(renderer=renderer)
   1907 
   1908         x, y = self.label.get_position()

/home/jbourbeau/cr-composition/.virtualenv/lib/python2.7/site-packages/matplotlib/axis.pyc in _get_tick_boxes_siblings(self, renderer)
   1889         for nn, axx in enumerate(grp.get_siblings(self.axes)):
   1890             ticks_to_draw = axx.xaxis._update_ticks(renderer)
-> 1891             tlb, tlb2 = axx.xaxis._get_tick_bboxes(ticks_to_draw, renderer)
   1892             bboxes.extend(tlb)
   1893             bboxes2.extend(tlb2)

/home/jbourbeau/cr-composition/.virtualenv/lib/python2.7/site-packages/matplotlib/axis.pyc in _get_tick_bboxes(self, ticks, renderer)
   1128         for tick in ticks:
   1129             if tick.label1On and tick.label1.get_visible():
-> 1130                 extent = tick.label1.get_window_extent(renderer)
   1131                 ticklabelBoxes.append(extent)
   1132             if tick.label2On and tick.label2.get_visible():

/home/jbourbeau/cr-composition/.virtualenv/lib/python2.7/site-packages/matplotlib/text.pyc in get_window_extent(self, renderer, dpi)
    920             raise RuntimeError('Cannot get window extent w/o renderer')
    921 
--> 922         bbox, info, descent = self._get_layout(self._renderer)
    923         x, y = self.get_unitless_position()
    924         x, y = self.get_transform().transform_point((x, y))

/home/jbourbeau/cr-composition/.virtualenv/lib/python2.7/site-packages/matplotlib/text.pyc in _get_layout(self, renderer)
    307                 w, h, d = renderer.get_text_width_height_descent(clean_line,
    308                                                         self._fontproperties,
--> 309                                                         ismath=ismath)
    310             else:
    311                 w, h, d = 0, 0, 0

/home/jbourbeau/cr-composition/.virtualenv/lib/python2.7/site-packages/matplotlib/backends/backend_agg.pyc in get_text_width_height_descent(self, s, prop, ismath)
    230             fontsize = prop.get_size_in_points()
    231             w, h, d = texmanager.get_text_width_height_descent(
--> 232                 s, fontsize, renderer=self)
    233             return w, h, d
    234 

/home/jbourbeau/cr-composition/.virtualenv/lib/python2.7/site-packages/matplotlib/texmanager.pyc in get_text_width_height_descent(self, tex, fontsize, renderer)
    499         else:
    500             # use dviread. It sometimes returns a wrong descent.
--> 501             dvifile = self.make_dvi(tex, fontsize)
    502             with dviread.Dvi(dvifile, 72 * dpi_fraction) as dvi:
    503                 page = next(iter(dvi))

/home/jbourbeau/cr-composition/.virtualenv/lib/python2.7/site-packages/matplotlib/texmanager.pyc in make_dvi(self, tex, fontsize)
    363                 self._run_checked_subprocess(
    364                     ["latex", "-interaction=nonstopmode", "--halt-on-error",
--> 365                      texfile], tex)
    366             for fname in glob.glob(basefile + '*'):
    367                 if not fname.endswith(('dvi', 'tex')):

/home/jbourbeau/cr-composition/.virtualenv/lib/python2.7/site-packages/matplotlib/texmanager.pyc in _run_checked_subprocess(self, command, tex)
    333             report = subprocess.check_output(command,
    334                                              cwd=self.texcache,
--> 335                                              stderr=subprocess.STDOUT)
    336         except subprocess.CalledProcessError as exc:
    337             raise RuntimeError(

/home/jbourbeau/cr-composition/.virtualenv/lib/python2.7/site-packages/subprocess32.pyc in check_output(*popenargs, **kwargs)
    332     if 'stdout' in kwargs:
    333         raise ValueError('stdout argument not allowed, it will be overridden.')
--> 334     process = Popen(stdout=PIPE, *popenargs, **kwargs)
    335     try:
    336         output, unused_err = process.communicate(timeout=timeout)

/home/jbourbeau/cr-composition/.virtualenv/lib/python2.7/site-packages/subprocess32.pyc in __init__(self, args, bufsize, executable, stdin, stdout, stderr, preexec_fn, close_fds, shell, cwd, env, universal_newlines, startupinfo, creationflags, restore_signals, start_new_session, pass_fds)
    612                                     c2pread, c2pwrite,
    613                                     errread, errwrite,
--> 614                                     restore_signals, start_new_session)
    615             except:
    616                 # The cleanup is performed within the finally block rather

/home/jbourbeau/cr-composition/.virtualenv/lib/python2.7/site-packages/subprocess32.pyc in _execute_child(self, args, executable, preexec_fn, close_fds, pass_fds, cwd, env, universal_newlines, startupinfo, creationflags, shell, p2cread, p2cwrite, c2pread, c2pwrite, errread, errwrite, restore_signals, start_new_session)
   1209                                 errread, errwrite,
   1210                                 errpipe_read, errpipe_write,
-> 1211                                 restore_signals, start_new_session, preexec_fn)
   1212                         self._child_created = True
   1213                     else:

KeyboardInterrupt: 

In [13]:

    

df_sim.log_s125.plot(kind='hist', bins=100)
plt.show()


    

In [20]:

    

log_125_bins = np.linspace(0, 2.5, 100)


    

In [15]:

    

df_sim.MC_log_energy.plot(kind='hist', bins=100)
plt.show()


    

In [16]:

    

log_MC_energy_bins = np.linspace(6.4, 8.0, 100)


    

In [21]:

    

plt.hist2d(df_sim.loc[energy_mask, 'log_s125'].values,
           df_sim.loc[energy_mask, 'MC_log_energy'].values,
           bins=[log_125_bins, log_MC_energy_bins])
plt.show()


    

In [ ]: