New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add flux point computation using Lafferty & Wyatt (1995) #128

Merged
merged 3 commits into from Jul 10, 2014

Conversation

Projects
None yet
2 participants
@ellisowen
Contributor

ellisowen commented Jun 11, 2014

Implements #50 .

@ellisowen

This comment has been minimized.

Contributor

ellisowen commented Jun 11, 2014

Pull request as you asked - It's not clear to me what is complete and what is work in progress as it's been a few weeks since I worked on this. I'll make this more complete over the next few days.

@cdeil cdeil added this to the 0.1 milestone Jun 11, 2014

@cdeil cdeil added the feature label Jun 11, 2014

@cdeil

This comment has been minimized.

Member

cdeil commented Jun 11, 2014

Please rebase this again against the master branch.
I just merged e45ca57 with a ton of whitespace changes which will probably create extra rebase conflicts.

@ellisowen

This comment has been minimized.

Contributor

ellisowen commented Jun 13, 2014

Still in progress - adding unit tests & docstrings. The actual Lafferty flux point implementation looks a bit fishy, so I will investigate if I actually did this correctly or if it was just a quick hack while I worked on the API.

@@ -33,25 +456,22 @@ def __init__(self, n_on, n_off, a_on, a_off):
@property
def alpha(self):
r"""Background exposure ratio.

This comment has been minimized.

@cdeil

cdeil Jun 13, 2014

Member

Looks like you introduced an indentation problem here and in a bunch of other places ... please review the diff on Github and fix those cases.

@@ -1,8 +1,8 @@
# Licensed under a 3-clause BSD style license - see LICENSE.rst
from __future__ import print_function, division
from numpy.testing import assert_allclose
from ..data import Stats
from ..data import *

This comment has been minimized.

@cdeil

cdeil Jun 13, 2014

Member

import * is forbidden in Astropy / Gammapy (and you shouldn't use it in your scripts either).
Either list directly the things you want to import, or use from .. import data and then use data.thing if there are very many things you want to import.

return values
if __name__ == '__main__':

This comment has been minimized.

@cdeil

cdeil Jun 13, 2014

Member

Did you commit test_script.py by accident?
If you actually want it to become part of Gammapy you need a better name ... usually the convention is to have for each file something.py a corresponding test file tests/test_something.py so that it's clear where to find the unit tests for stuff from something.py.

class ApperturePhotometry(object):
def __init__(self, lat, lon, radius, psf=None, counts=None, background=None, excess=None, exposure=None, bg_model=None, flux=None):

This comment has been minimized.

@cdeil

cdeil Jun 13, 2014

Member

Please roughly respect the 80-character limit and PEP8 guidelines.
If you activate the PEP8 plugin in the Eclipse Pydev settings as warnings your editor will point things out as you type.

@staticmethod
def from_ascii(filename):
"""File should be an ascii table in two columns; energy and flux of the event
Headers 'Energy' and 'Flux' should be stated

This comment has been minimized.

@cdeil

cdeil Jun 13, 2014

Member

Another case of screwed up indentation.

'compute_total_stats']
class ImportsExports(object):

This comment has been minimized.

@cdeil

cdeil Jun 13, 2014

Member

Having a class with name ImportsExports or LoadInputCubes is a bad idea ... there should be classes representing the objects your code deals with (e.g. PowerLawSpectrumModel or FluxPoints or ...) and they can have static or class methods to load various formats ... look e.g. at the TablePSF or GammaSpectralCube I recently wrote for guidance.

__all__ = ['Stats', 'make_stats', 'combine_stats',
__all__ = ['ImportsExports', 'CalcSED', 'ApperturePhotometry',

This comment has been minimized.

@cdeil

cdeil Jun 13, 2014

Member

Aperture with on p only.

from astropy.wcs import WCS
radius = radius / image.header['CDELT2']
if cube == 'Yes':
from gammapy.image import cube_to_image

This comment has been minimized.

@cdeil

cdeil Jun 13, 2014

Member

Gammapy and astropy imports always at the top of the file.

@ellisowen

This comment has been minimized.

Contributor

ellisowen commented Jun 17, 2014

Just spotted some other things to change - will commit again shortly.

@ellisowen

This comment has been minimized.

Contributor

ellisowen commented Jun 17, 2014

The data file 'read' function I think can be improved (and tests added), but I will do this later in a separate pull request.

@ellisowen

This comment has been minimized.

Contributor

ellisowen commented Jun 17, 2014

I think this is OK now - Please can you review?

@cdeil

This comment has been minimized.

Member

cdeil commented Jun 17, 2014

Here's the API we discussed just now:
http://nbviewer.ipython.org/gist/cdeil/c49b786a9e6ceca19c95

The documentation of the format should be in gammapy/docs/dataformats/spectra.rst ... thinking about it some more we don't even need our own read and write methods ... users can just use the Table methods and rename columns to be consistent if needed themselves. It's the most flexible approach for the user and the lest work for us.

@ellisowen

This comment has been minimized.

Contributor

ellisowen commented Jun 18, 2014

Implements as we discussed yesterday - seems to give consistent conversion between integral and differential fluxes and back again. Still need to check whether the actual values are correct though.

I still need to:

  • Add an example
  • Add unit tests
  • Update docstrings
  • Change syntax to be consistent with that used in gammapy.spectrum

Just wondered if the implementation seems roughly OK to you/if there is anything missing or that should obviously be changed?

@ellisowen

This comment has been minimized.

Contributor

ellisowen commented Jun 19, 2014

Still working on example & some unit tests.

raise ValueError('Unknown method {0}'.format(y_method))
energy_col = Column(name='ENERGY', data=energy)
table.add_column(energy_col)

This comment has been minimized.

@cdeil

cdeil Jun 19, 2014

Member

You can use this way to add columns to tables:

table['column_name'] = numpy_array
return XYMETHODS[y_method](yint, xmin, xmax, x_method, model)
def _x_log_center(xmin, xmax, function=None, table=None):

This comment has been minimized.

@cdeil

cdeil Jun 19, 2014

Member

Use scipy.stats.gmean to compute log center positions

@ellisowen

This comment has been minimized.

Contributor

ellisowen commented Jun 24, 2014

There are still a couple of minor things I need to sort out. Will do this later today (& will also squash commits then too)

@ellisowen

This comment has been minimized.

Contributor

ellisowen commented Jun 24, 2014

The unit test test_compute_differential_flux_points() currently fails due to insufficient info about published data points it compares to i.e powerlaw index/spectral model/how energies are calculated in bins (although the values aren't wildly off). I guess comparison against a better understood dataset is the way forward here? Any suggestions of any?

Otherwise, can you please review?

@cdeil

This comment has been minimized.

Member

cdeil commented Jun 24, 2014

I'll have a look tomorrow. Can you please rebase this via git rebase master and then after resolving the conflicts force push this branch again to Github.

from .flux_point import *
from . import fitting_utils
from .inverse_compton import *

This comment has been minimized.

@cdeil

cdeil Jun 25, 2014

Member

inverse_compton.py has been removed in the meantime in master ... you need to rebase again and resolve this merge conflict (or remove this line).

'energy_bounds_equal_log_spacing',
'np_to_pha',
]
def linear_extrapolator(x_vals, y_vals):

This comment has been minimized.

@cdeil

cdeil Jun 25, 2014

Member

I'm not sure if this is useful ... can't the same thing be done e.g. with
numpy.polyfit or numpy.linalg.leastsq?

If this actually is something that can't easily be done with numpy or scipy, then OK to put it in gammapy, but please somewhere in gammapy/utils instead of gammapy/spectrum/utils.py because it's not spectrum-specific.

@cdeil

This comment has been minimized.

Member

cdeil commented Jun 25, 2014

Please have a look at the coverage of flux_points.py and add unit tests for the red lines (most importantly y_method == 'Model' is not tested at the moment):

python setup.py test --coverage
open htmlcov/index.html 
Useful for automatic processing.
"""
XYMETHODS['LogCenter'] = _x_log_center

This comment has been minimized.

@cdeil

cdeil Jun 25, 2014

Member

Please change to snake_case everywhere instead of CamelCase for the option strings, i.e. log_center, lafferty, ... I think this is done consistently in numpy, scipy, astropy and doing it differently here will just lead to users mis-typing it all the time.

# Licensed under a 3-clause BSD style license - see LICENSE.rst
"""Demonstrates integral flux to differential with Lafferty method)
"""
from gammapy.spectrum.flux_point import compute_differential_flux_points

This comment has been minimized.

@cdeil

cdeil Jun 25, 2014

Member

I'm trying to get a flatter namespace hierarchy in gammapy ... in this case it should be

from gammapy.spectrum import compute_differential_flux_points
"""
def from_dict(self, data):
emins = np.asanyarray(table['ENERGY_MIN'])

This comment has been minimized.

@cdeil

cdeil Jun 30, 2014

Member

how about using the same names as the table columns here?
I.e. energy_min, energy_max and int_flux?
IMO makes code a little bit easier to read and debug.

pass
table.meta['spectral_index'] = spectral_index

This comment has been minimized.

@cdeil

cdeil Jun 30, 2014

Member

Is it possible to add a description like "Spectral index assumed in the DIFF_FLUX computation"?

from flux_point_demo import get_flux_tables, plot_flux_points
from astropy.modeling.models import ExponentialCutoffPowerLaw1D
from astropy.table import Table
import gc

This comment has been minimized.

@cdeil

cdeil Jun 30, 2014

Member

Importing gc and calling gc_collect should not be necessary.
Why did you add it?
(If there are cyclic references that currently make this necessary, we have to restructure the code to avoid them, but I don't think there are!?)

y = ExponentialCutoffPowerLaw1D.eval(x, 10, 5, 1, 0.5)
function = lambda x: ExponentialCutoffPowerLaw1D.eval(x, 10, 5, 1, 0.5)
# Set the x-bins
xmin = [0.1, 0.3, 1, 3, 10, 30]

This comment has been minimized.

@cdeil

cdeil Jun 30, 2014

Member

Use the same names in the example as are required in the table, i.e. energy_min, energy_max

# Define the function
x = np.arange(0.01, 100, 0.001)
y = ExponentialCutoffPowerLaw1D.eval(x, 10, 5, 1, 0.5)
function = lambda x: ExponentialCutoffPowerLaw1D.eval(x, 10, 5, 1, 0.5)

This comment has been minimized.

@cdeil

cdeil Jun 30, 2014

Member

Do this instead:

spectral_model = ExponentialCutoffPowerLaw1D(10, 5, 1, 0.5)
y = spectral_model(x)

and then later pass the spectral_model callable directly instead of creating a function that wraps it.

class FluxPoints(object):
"""List of flux points.
def compute_differential_flux_points(table, spectral_index=None, x_method='lafferty',
y_method='PowerLaw', model=None):

This comment has been minimized.

@cdeil

cdeil Jun 30, 2014

Member

Please use power_law consistently instead of plaw as you now do in the file names.

lafferty_flux, log_flux = get_flux_tables(table, y_method, function)
axarr[0].loglog(x, ((x ** 2) * y))
axarr[0].loglog(lafferty_flux['ENERGY'], ((lafferty_flux['ENERGY'] ** 2) * lafferty_flux['DIFF_FLUX']), marker='D', color='k', linewidth=0, ms=5, label='Lafferty Method')

This comment has been minimized.

@cdeil

cdeil Jun 30, 2014

Member

Please reformat to roughly respect the 80 character line length limit (and from now on have PEP8 warnings on in your editor to realise it directly while you are coding).

plt = plot_flux_points(table, x, y, function, xmin, xmax, 'model')
plt.legend()
plt.savefig('model_example.png')

This comment has been minimized.

@cdeil

cdeil Jun 30, 2014

Member

I think the call to plt.savefig should be removed ... the Sphinx plot directive takes care of saving the plot and creating a link? See example here:
http://astropy.readthedocs.org/en/latest/modeling/index.html

Or is there a reason you need to handle this yourself here?

@cdeil

This comment has been minimized.

Member

cdeil commented Jul 1, 2014

Here's the analytical check where the point should go for power-laws we just coded together:
http://nbviewer.ipython.org/gist/cdeil/ff29110c1257fde3e06a

Begins implementation of Lafferty and Wyatt flux point method and API…
… changes to flux_point.py

Adds xfail test command to added unit test
@ellisowen

This comment has been minimized.

Contributor

ellisowen commented Jul 2, 2014

How does this look?

Where to stick your Spectral Points?
====================================
The ``gammapy.spectrum`` module offers a number of options for positioning data points within an energy band. This example offers a comparison between

This comment has been minimized.

@cdeil

cdeil Jul 2, 2014

Member

Use single-ticks and check that this generated a working link in the html docs.

====================================
The ``gammapy.spectrum`` module offers a number of options for positioning data points within an energy band. This example offers a comparison between
the log center and Lafferty & Wyatt (described in [Lafferty1994]_) methods. See ``gammapy.spectrum.flux_point`` for documentation on usage.

This comment has been minimized.

@cdeil

cdeil Jul 2, 2014

Member

This should be

`gammapy.spectrum.compute_differential_flux_points`

so that a link to that function in the html docs is generated.

# Compute x point
if x_method == 'table':
energy = np.array(table['ENERGY'])
elif x_method == 'log_center':

This comment has been minimized.

@cdeil

cdeil Jul 2, 2014

Member

Please add a unit test for the x_method == 'log_center' case.

This comment has been minimized.

@ellisowen

ellisowen Jul 2, 2014

Contributor

We removed this before as the code just calls into scipy - is it still worth including this?

e2 = np.max(energy_max)
# TODO: we shouldn't have to re-code the power-law function here.
# Try to find a way to re-use the existing powerlaw model.
model = lambda x: (x / (1 - g)) * (((e2 / x) ** (1 - g)) - ((e1 / x) ** (1 - g)))

This comment has been minimized.

@cdeil

cdeil Jul 2, 2014

Member

Now that you have an analytical expression for _energy_lafferty_power_law this power-law model here is never used and can be removed?

This comment has been minimized.

@ellisowen

ellisowen Jul 2, 2014

Contributor

No, this is still required to determine the y-position, line 82

energy = _energy_lafferty_power_law(energy_min, energy_max,
spectral_index)
else:
energy = np.array(_x_lafferty(table['ENERGY_MIN'],

This comment has been minimized.

@cdeil

cdeil Jul 2, 2014

Member

Please add a unit test that covers this case ... it's the most important one, passing a general model spectrum, no?

This comment has been minimized.

@ellisowen

ellisowen Jul 2, 2014

Contributor

Can I add this as a separate unit test, or as a second part to the existing test_compute_differential_flux_points() unit test?

assert_allclose(points.eyl[13], 2.724e-14)
assert_allclose(points.eyh[13], 2.512e-14)
from astropy.table import Table
from ..flux_point import _x_lafferty, _integrate, _ydiff_excess_equals_expected, compute_differential_flux_points, _energy_lafferty_power_law

This comment has been minimized.

@cdeil

cdeil Jul 2, 2014

Member

Either do

from .. import flux_point

and use flux_point._x_lafferty etc everywhere,
or break the line like this:

from ..flux_point import (_x_lafferty, _integrate, _ydiff_excess_equals_expected, 
                                      compute_differential_flux_points, _energy_lafferty_power_law)
temp = e1 ** (1 - g) - e2 ** (1 - g)
self.y = self.yint / self.x / temp / (g - 1)
return self.y
from .powerlaw import _conversion_factor

This comment has been minimized.

@cdeil

cdeil Jul 2, 2014

Member

Put imports from astropy and gammapy at the top.

Power Law Assumption
--------------------
In many "real world" examples, the nature of the true underlying frequency distribution is unknown and can only be estimated, either by using a

This comment has been minimized.

@cdeil

cdeil Jul 2, 2014

Member

I would suggest to always talk about "spectrum" instead of "frequency distribution" ... it's the more common term.

In many "real world" examples, the nature of the true underlying frequency distribution is unknown and can only be estimated, either by using a
fitting algorithm with data points or by assuming a certain spectral form. In this example, the true underlying frequency distribution
(being a piecewise power law function of indices -1 through to -5 for bins of increasing energy) is shown in blue for illustration. This would be

This comment has been minimized.

@cdeil

cdeil Jul 2, 2014

Member

always say "spectral indices" instead of only "indices" to be more clear and remove the minus here ... the gamma parameter in the formulas is positive because the minus is explicitly put in the model formula.

@ellisowen

This comment has been minimized.

Contributor

ellisowen commented Jul 10, 2014

Please can you review?

@cdeil

This comment has been minimized.

Member

cdeil commented Jul 10, 2014

Thanks!

cdeil added a commit that referenced this pull request Jul 10, 2014

@cdeil cdeil merged commit 6110eca into gammapy:master Jul 10, 2014

1 check failed

continuous-integration/travis-ci The Travis CI build could not complete due to an error
Details

@ellisowen ellisowen deleted the ellisowen:flux_point_method branch Jul 10, 2014

@cdeil

This comment has been minimized.

Member

cdeil commented Jul 19, 2014

@cdeil cdeil changed the title from Flux point method to Add flux point computation using Lafferty & Wyatt (1995) Apr 8, 2015

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment