threeML / hawc_hal

HAWC Accelerated Likelihood - python-only framework for HAWC data analysis
BSD 3-Clause "New" or "Revised" License
11 stars 21 forks source link

Plotting residuals (seemingly) incompatible with non-int bins #60

Closed cbrisboi closed 3 years ago

cbrisboi commented 3 years ago

Describe the bug When using hawc_hal with the 2D bins, if one wants to plot the binned count residuals one get an error

  File "minimal_crab.py", line 71, in <module>
    fig = hawc.display_spectrum()
  File "/data/disk01/home/cbrisboi/hawc_software/miniconda3/envs/hal_prod/lib/python3.7/site-packages/hawc_hal/HAL.py", line 410, in display_spectrum
    return self._plot_spectrum(net_counts, yerr, model_only, residuals, residuals_err)
  File "/data/disk01/home/cbrisboi/hawc_software/miniconda3/envs/hal_prod/lib/python3.7/site-packages/hawc_hal/HAL.py", line 415, in _plot_spectrum
    planes = np.array(self._active_planes, dtype=int)
ValueError: invalid literal for int() with base 10: '1c'

To Reproduce

Here is a minimal example produced from the Markarian example in the README.

from hawc_hal import HAL, HealpixConeROI
import matplotlib.pyplot as plt
from threeML import *

# Define the ROI                                                                                                                             
ra_crab, dec_crab = 83.63, 22.01
data_radius = 3.0
model_radius = 8.0

roi = HealpixConeROI(data_radius=data_radius,
                     model_radius=model_radius,
                     ra=ra_crab,
                     dec=dec_crab)

# Instance the plugin                                                                                                                        

response ="/data/scratch/userspace/kmalone/sweets-for-systematics/stages/bobs-ednas-thresh0p2/detRes-GP-allDecs.root"
maptree ="/data/scratch/userspace/kmalone/datasets/energy/maps/combined/dataset-ch103-ch603/rejiggered-files/maptree-ch103-ch603.root"

hawc = HAL("HAWC",
           maptree,
           response,
           roi)

# Use from bin 1 to bin 9                                                                                                                                                                                                                       
hawc.set_active_measurements(bin_list=['1c', '1d', '1e', '2c', '2d', '2e', '3c', '3d', '3e', '3f', \
                                      '4c', '4d', '4e', '4f', '4g', '5d', '5e', '5f', '5g', '5h', \
                                      '6d', '6e', '6f', '6g', '6h', '7f', '7g', '7h', '7i', '8f', \
                                      '8g', '8h', '8i', '8j', '9g', '9h', '9i', '9j', '9k', '9l'])

# Define model as usual                                                                                                                      
spectrum = Log_parabola()
source = PointSource("crab", ra=ra_crab, dec=dec_crab, spectral_shape=spectrum)

spectrum.piv = 1 * u.TeV
spectrum.piv.fix = True

spectrum.K = 1e-14 / (u.TeV * u.cm ** 2 * u.s)  # norm (in 1/(keV cm2 s))                                                                    
spectrum.K.bounds = (1e-25, 1e-19)  # without units energies are in keV                                                                      

spectrum.beta = 0  # log parabolic beta                                                                                                      
spectrum.beta.bounds = (-4., 2.)

spectrum.alpha = -2.5  # log parabolic alpha (index)                                                                                         
spectrum.alpha.bounds = (-4., 2.)

model = Model(source)

data = DataList(hawc)

jl = JointLikelihood(model, data, verbose=False)
jl.set_minimizer("minuit")
param_df, like_df = jl.fit()

# See the model in counts space and the residuals                                                                                            
fig = hawc.display_spectrum()
# Save it to file                                                                                                                            
fig.savefig("hal_crab_residuals.png")

Expected behavior I expect a plot to be output, and saved. A clear and concise description of what you expected to happen.

Log files

[WARNING ] The GSL library or the pygsl wrapper cannot be loaded. Models that depend on it will not be available.
[WARNING ] The ebltable package is not available. Models that depend on it will not be available
[INFO    ] Starting 3ML!

WARNING RuntimeWarning: numpy.ndarray size changed, may indicate binary incompatibility. Expected 16 from C header, got 88 from PyObject

[WARNING ] The cthreeML package is not installed. You will not be able to use plugins which require the C/C++ interface (currently HAWC)
[WARNING ] Could not import plugin HAWCLike.py. Do you have the relative instrument software installed and configured?
[WARNING ] Could not import plugin FermiLATLike.py. Do you have the relative instrument software installed and configured?

WARNING DeprecationWarning: `np.object` is a deprecated alias for the builtin `object`. To silence this warning, use `object` by itself. Doing this will not modify any behavior and is safe. 
Deprecated in NumPy 1.20; for more details and guidance: https://numpy.org/devdocs/release/1.20.0-notes.html#deprecations

[INFO    ] Creating singleton for /data/scratch/userspace/kmalone/sweets-for-systematics/stages/bobs-ednas-thresh0p2/detRes-GP-allDecs.root

WARNING DeprecationWarning: `np.object` is a deprecated alias for the builtin `object`. To silence this warning, use `object` by itself. Doing this will not modify any behavior and is safe. 
Deprecated in NumPy 1.20; for more details and guidance: https://numpy.org/devdocs/release/1.20.0-notes.html#deprecations

[WARNING ] We have set the min_value of crab.spectrum.main.Log_parabola.K to 1e-99 because there was a postive transform
[INFO    ] set the minimizer to minuit
[INFO    ] set the minimizer to MINUIT
[WARNING ] The current value of the parameter K (1.0) was above the new maximum 1e-19.
Best fit values:

                                                         result  \
parameter                                                         
crab.spectrum.main.Log_parabola.K      (3.76 +/- 0.07) x 10^-20   
crab.spectrum.main.Log_parabola.alpha          -2.363 +/- 0.028   
crab.spectrum.main.Log_parabola.beta    (1.08 +/- 0.09) x 10^-1   

                                                  unit  
parameter                                               
crab.spectrum.main.Log_parabola.K      1 / (cm2 keV s)  
crab.spectrum.main.Log_parabola.alpha                   
crab.spectrum.main.Log_parabola.beta                    

Correlation matrix:

 col0  col1  col2
----- ----- -----
 1.00 -0.72 -0.48
-0.72  1.00  0.92
-0.48  0.92  1.00

Values of -log(likelihood) at the minimum:

       -log(likelihood)
HAWC       66270.311832
total      66270.311832

Values of statistical measures:

     statistical measures
AIC         132546.623714
BIC         132579.890109
Traceback (most recent call last):
  File "minimal_crab.py", line 67, in <module>
    fig = hawc.display_spectrum()
  File "/data/disk01/home/cbrisboi/hawc_software/miniconda3/envs/hal_prod/lib/python3.7/site-packages/hawc_hal/HAL.py", line 410, in display_spectrum
    return self._plot_spectrum(net_counts, yerr, model_only, residuals, residuals_err)
  File "/data/disk01/home/cbrisboi/hawc_software/miniconda3/envs/hal_prod/lib/python3.7/site-packages/hawc_hal/HAL.py", line 415, in _plot_spectrum
    planes = np.array(self._active_planes, dtype=int)
ValueError: invalid literal for int() with base 10: '1c'

Desktop (please complete the following information): (not a desktop, but its the machine I'm running on)

Additional context I should have some time to investigate this next week.

henrikef commented 3 years ago

Hm, I see. The tests only run on fHit maps. Do we have any maps/IRFs with 2D binning that are publicly available? Then we could add this to the unit tests.

In any case, try just removing the dtype=int. I think that current versions of matplotlib should be able to handle alphanumeric data points.

henrikef commented 3 years ago

Fixed in #64