Plotting

[henrikef]

I propose some updates to plotting routines...

1. Added a function to joint_likelihood.py that will return 1D/2D likelihood profiles/contours for all parameters/combinations of parameters for convenience.

2. Added an option to plot_point_source_spectra() to also plot extended source spectra (this will just plot whatever was set as the spectral shape, i.e. the 'spatially integrated' spectrum.

3a. Moved convergence_plots() and plot_chains() from bayesian_analysis.py to analysis_results.py, so that the plots can be produced from analysis results loaded from disk.
3b. Changed convergence_plots() to use the same number of samples for the bootstrapping average as for the sliding window, I hope that was the original intent.

[grburgess]

The tests are failing here. Can you show an example in pseudo-code of what you want to achieve? My original goal with this code was to have a separate handler for extended sources. Maybe we can coordinate and start working on that. I just have not had access to extended source analysis.

[giacomov]

@henrikef thanks for these nice contributions.

Maybe to make things easier on the testing/debugging you could divide this into two pull requests: one with the plotting things that do not use extended source (which I could merge right now since it is passing the tests), and another one with the addition that is failing now (the extended source part). This is just a suggestion, if you prefer to deal with everything at the same time that's fine too.

In any case, please take care of adding tests for all functionalities introduced.

Thanks!

[codecov-io]

Codecov Report

Merging #241 into master will decrease coverage by 0.29%.
The diff coverage is 57.74%.

@@            Coverage Diff            @@
##           master     #241     +/-   ##
=========================================
- Coverage   75.43%   75.14%   -0.3%     
=========================================
  Files          94       94             
  Lines       10979    11058     +79     
=========================================
+ Hits         8282     8309     +27     
- Misses       2697     2749     +52

[giacomov]

@grburgess has been working on the plotting interface up to now, I'll leave to him to review and merge when ready.

In the meantime thanks to @henrikef for this contribution!

[giacomov]

@henrikef @grburgess what is the status on this?

[henrikef]

Hi everyone!

I'd still like to get this pull request merged if possible. The travis-ci build is failing but as far as I can tell, the only test that fails is threeML/test/test_joint_likelihood_set.py::test_joint_likelihood_set_parallel for MacOS. I don't think it has anything to do with the changes in my branch as I also saw the same failure for a different pull request where I literally changed only one line having to do with reading from fits files. The tests run fine on my laptop and for the linux build.

Regarding plotting extended source spectra: @grburgess , what was the motivation for having separate handlers for extended sources? Wouldn't it just mean a lot of code would have to be duplicated? Maybe I'm missing something, but I didn't see anything in the current code that was particular to point sources. Me and several colleagues are using 3ML to fit extended sources to HAWC data, and we need to be able to plot the spectra. So far, some people are using my fork/branch, but it's annoying to have to keep it up to date whenever there are other changes to 3ML. But if there's something wrong with my approach of using the exact same methods for extended sources as for point sources, please let me know!

[grburgess]

My main concern was having this be general enough so that things can be expanded upon and that plotting of point sources was different from extended sources. That being said, is it possible for you to provide an example of the usage so that we can see how it works?

I trust that you and @giacomov have a much better handle on how extended sources should be plotted, so I really have to rely on your judgement to see if this fits into a bigger framework. Does this make sense?

[henrikef]

Maybe I'm being dense, but I don't understand what you mean by

that plotting of point sources was different from extended sources.

Could you explain? I'm talking about plotting the 'total' flux here, not the surface brightness (flux per solid angle).

I can't really upload HAWC datasets, and I think the fermipy plugin cannot handle extended sources yet, but I include an example with made-up results below:

import threeML
import astromodels
import numpy as np
import astropy.units as u
import matplotlib as mpl
import matplotlib.pyplot as plt

mpl.rc("font", family="serif", size=14)
from threeML.io.progress_bar import progress_bar

sedUnit = u.TeV  / ( u.cm**2 * u.s)
fluxUnit = 1  / (u.TeV* u.cm**2 * u.s)
E = np.logspace(np.log10(1), np.log10(100), 100) * u.TeV


specP = astromodels.Powerlaw()
sourceP = astromodels.PointSource( "pointsource", ra=20., dec=20., spectral_shape = specP )

specP.K = 2e-13*fluxUnit
specP.K.fix = False
specP.index = -2.5
specP.index.fix = False
specP.piv = 7*u.TeV
specP.piv.fix=True

specE = astromodels.Powerlaw()
shape = astromodels.Disk_on_sphere()

sourceE = astromodels.ExtendedSource( "extended_source", spatial_shape=shape, spectral_shape = specE )

specE.K = 3e-13*fluxUnit
specE.K.fix=False
specE.index = -2.7
specE.index.fix =  False
specE.piv = 7*u.TeV
specE.piv.fix=True
shape.radius = 1*u.degree
shape.radius.fix =False
shape.lat0  = 22*u.degree
shape.lat0.fix = True
shape.lon0 = 22*u.degree
shape.lon0.fix = True

model = astromodels.Model(sourceP, sourceE)

#print model.free_parameters.keys()
#['pointsource.spectrum.main.Powerlaw.K', 'pointsource.spectrum.main.Powerlaw.index', 'extended_source.Disk_on_sphere.radius', 'extended_source.spectrum.main.Powerlaw.K', 'extended_source.spectrum.main.Powerlaw.index']


uncertainties=[5e-14*fluxUnit, 0.3, 0.2*u.degree, 7e-14*fluxUnit, 0.4]

cov_list = []
for unc, (k, p)  in zip(uncertainties, model.free_parameters.items()):
  if p.unit != u.dimensionless_unscaled:
    unc = unc.to(p.unit).value
  if not p.has_transformation():
    cov_list.append(unc**2)
  else:  
    trafo=p.transformation.forward
    cov_list.append( (trafo( p.value + unc ) - trafo( p.value ))*(trafo( p.value )-trafo( p.value - unc )) )
      
cov_matrix=np.diag( np.array(cov_list) )
results = threeML.analysis_results.MLEResults(model, cov_matrix, {})  

fig, ax = plt.subplots(1, 1, figsize=(10, 5))

threeML.plot_point_source_spectra( results, ene_min=1*u.TeV, ene_max=100*u.TeV, flux_unit=sedUnit, grid=True, subplot=ax, include_extended=True, num_ene=20 )

fig.savefig("Spectra_including_extended.png" )

fig, ax = plt.subplots(1, 1, figsize=(10, 5))

threeML.plot_point_source_spectra( results, ene_min=1*u.TeV, ene_max=100*u.TeV, flux_unit=sedUnit, grid=True, subplot=ax, include_extended=False, num_ene=20 )

fig.savefig("Spectra_without_extended.png" )

[giacomov]

Reading the discussion above I suggest these modifications:

• plot_point_source_spectra -> plot_spectra
• change the plot_spectra calling sequence to include two optional parameters ra and dec. Within the code we should check the type of the requested source, if it is a point source these parameters are ignored, if it is an extended source these parameters are mandatory. This would allow to plot the spectrum of the extended source at any position, so it will work both for extended source with the same spectrum everywhere (like in @henrikef example) and for sources where the spectrum changes point to point
• after that, the only difference between a point source and an extended source would be in the calling sequence of the source. You would use spectrum = source(energies) for point sources, while spectrum = source(ra, dec, energies) for extended sources. This can be handled with a simple if. After this, all else is the same.

I don't know all the details of the code, but I think these suggestions would avoid code duplication and at the same time allow for the possibility of plotting extended source spectra.

[grburgess]

Aha. Now I understand. I agree /J. Michael

…

On 25 Apr 2018, at 21:33, Giacomo Vianello ***@***.***> wrote: Reading the discussion above I suggest these modifications: plot_point_source_spectra -> plot_spectra change the plot_spectra calling sequence to include two optional parameters ra and dec. Within the code we should check the type of the requested source, if it is a point source these parameters are ignored, if it is an extended source these parameters are mandatory. This would allow to plot the spectrum of the extended source at any position, so it will work both for extended source with the same spectrum everywhere (like in @henrikef example) and for sources where the spectrum changes point to point after that, the only difference between a point source and an extended source would be in the calling sequence of the source. You would use spectrum = source(energies) for point sources, while spectrum = source(ra, dec, energies) for extended sources. This can be handled with a simple if. After this, all else is the same. I don't know all the details of the code, but I think these suggestions would avoid code duplication and at the same time allow for the possibility of plotting extended source spectra. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or mute the thread.

[henrikef]

Hi @giacomov and @grburgess !

So the issue is with extended sources where the spectrum depends on the RA and Dec? Could you maybe provide an example of that for me, as I've never really worked with them.

Anyway, changing the name of the calling function is easy, I can do that :)

Adding RA/Dec parameters less so. By default, the plotting function loops through all or a given subset of the sources in the model. So we have to think of a user-friendly way to supply the requested coordinates for extended sources then. After that, I really have no idea for how to get at the flux at a given point. I'll look into it but I'd appreciate some help.

[colasri]

That would be great, what's the status here?

[grburgess]

Hi all, I'm sorry for falling behind on this PR. I am a bit behind on where things are. I think it would be very useful to have a skype between all of us to sort out what is the best approach as this will be very important for many instruments.

[colasri]

While I understand the concern about more general sources where the spectral shape could depend on the location (is there implementation of this somewhere in astromodels?), I think what Henrike implemented here addresses 99% (or 100%) of the current need.
I didn't know of this branch until recently, I'll probably start using it too if not merged.

[henrikef]

So I looked into adding an option to plot the spectrum at a given direction as @giacomov had asked. I honestly couldn't figure out a good way to add it. But if someone else wants to help me, that'd be great! If not, and not having that option is a deal breaker, I'll try to keep this branch up to date with the 3ML master. As it is, it works great for plotting the total spectrum for sources with separable spatial and spectral components.

[giacomov]

Let me fix the issue with HAL and I'll get to this. There is only so much I can do in a given day ;-)

The change to the plotting function that would allow plotting any source is trivial as explained above, I just need the time to sit down, do it and test it. Hopefully tomorrow...