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

Re-write Galaxy modelling code #157

Merged
merged 3 commits into from Aug 7, 2014

Conversation

Projects
None yet
2 participants
@adonath
Member

adonath commented Jul 26, 2014

Revision of the galaxy source modelling code, including the following changes:

  • Implementation of the spatial and velocity distribution models as astropy.modeling.Fittable1DModel
  • Use of astropy.units.Quantity throughout the code
  • Added tests
  • Added minimal docs
from astropy.table import Table
from astropy.utils.data import get_pkg_data_filename
filename = get_pkg_data_filename('data/atnf_sample.txt')

This comment has been minimized.

@cdeil

cdeil Jul 27, 2014

Member

This test fails because the data/atnf_sample.txt file is not installed:
https://travis-ci.org/gammapy/gammapy/jobs/30936682#L974

You probably need to add this file via get_package_data in setup_package.py as described here ... you can find a few examples in other parts of Gammapy or Astropy.

assert_allclose(pulsar.luminosity_spindown(t), reference)
def test_Pulsar_energy_integrated():

This comment has been minimized.

@cdeil

cdeil Jul 27, 2014

Member

There's a few tests on travis-ci like this one that fail because you didn't declare to skip them if scipy is not present:
https://travis-ci.org/gammapy/gammapy/jobs/30936680#L1636

plt.xlim(0, max_radius)
plt.ylim(0, 0.28)
plt.xlabel('Galactocentric Distance [kpc]')
plt.ylabel('Surface Density [a.u.]')

This comment has been minimized.

@cdeil

cdeil Jul 27, 2014

Member

You write 'Surface Density [a.u.]' here ... don't you have some well-defined units here?

plt.xlim(v_min, v_max)
plt.ylim(0, 0.004)
plt.xlabel('Velocity [km/s]')
plt.ylabel('PDF [a.u.]')

This comment has been minimized.

@cdeil

cdeil Jul 27, 2014

Member

Don't you have well-defined units here ... maybe PDF dP / dv (1 / (km / s)) or something?

@@ -9,14 +9,22 @@ Astrophysical source models (`gammapy.astro.source`)
Introduction
============
`gammapy.astro.source` implements some common astrophysical source and population models.
The `gammapy.astro.source` module contains classes of source models, which can be used for population synthesis of galactic gamma-ray sources.
Every source model class

This comment has been minimized.

@cdeil

cdeil Jul 27, 2014

Member

This sentence in the docs in unfinished.

TODO: describe
The `gammapy.astro.population` module provides a simple framework for population synthesis of
TeV sources.

This comment has been minimized.

@cdeil

cdeil Jul 27, 2014

Member

TeV -> gamma-ray

from gammapy.astro.source import SNR, SNRTrueloveMcKee
snr_models = [SNR, SNRTrueloveMcKee]
densities = [Quantity(1, 'cm^-3'), Quantity(0.1, 'cm^-3')]

This comment has been minimized.

@cdeil

cdeil Jul 27, 2014

Member

Can you write this instead here?

Quantity([1, 0.1], 'cm^-3')

If it works I would find it a bit simpler to read.

label=snr.__class__.__name__ + ' (n_ISM = {0})'.format(density.value),
linestyle=linestyle, color=color)
plt.xlabel('time [years]')
plt.ylabel('radius [pc]')#

This comment has been minimized.

@cdeil

cdeil Jul 27, 2014

Member

Remove # at the end of line.

ra, dec = coordinate.ra.deg, coordinate.dec.deg
# Add columns to table
table['distance'] = Column(distance, unit='pc', description='Distance observer to source center')
table['GLON'] = Column(glon, unit='deg', description='Galactic longitude')
table['GLAT'] = Column(glat, unit='deg', description='Galactic latitude')
table['VGLON'] = Column(v_glon.to('deg/Myr'), unit='deg/Myr', description='Velocity in Galactic longitude')

This comment has been minimized.

@cdeil

cdeil Jul 27, 2014

Member

Break line (put "description" on extra line)?

theta : `~astropy.units.Quantity`
Angle coordinate
amount: float
Amount of blurring of the position, given as a fraction of `radius`.

This comment has been minimized.

@cdeil

cdeil Jul 27, 2014

Member

If you're documenting the Parameters, maybe document the Returns as well?

This comment has been minimized.

@adonath

adonath Jul 31, 2014

Member

I started to document the function, than I realized that it is not visible to the end-user. What's the policy in this case? How should "private" functions be documented?

This comment has been minimized.

@cdeil

cdeil Aug 4, 2014

Member

There is no policy for documentation of private functions.

Writing code / docs such that you think it would be easiest to understand for someone else (or yourself in a year) is best. But there's also the 80:20 rule that says you should stop when it's in an OK state, which in this case, for well-structured code, means that it's OK to leave private functions without docstrings or just a one-line comment and instead move on to the next job.

dx, dy = cartesian(dr, dtheta)
return polar(x + dx, y + dy)
def _gc_correction(self, radius, theta, r_corr=Quantity(2.857, 'kpc')):

This comment has been minimized.

@cdeil

cdeil Jul 27, 2014

Member

If you're documenting the Parameters, maybe document the Returns as well?

max_age = 1e6
return simulate.make_cat_gal(nsources=nsources, rad_dis=YK04, vel_dis=H05, max_age=max_age)
return simulate.make_base_catalog_galactic(n_sources=n_sources, rad_dis=rad_dis, vel_dis=vel_dis, max_age=max_age)

This comment has been minimized.

@cdeil

cdeil Jul 27, 2014

Member

Break line.

'velocity_distributions',
]
# Simulation range used for random number drawing
v_range = 4000 # km/s
VMIN, VMAX = 0, 4000 # km/s

This comment has been minimized.

@cdeil

cdeil Jul 27, 2014

Member

Use Quantity here?

Density in velocity ``v``
amplitude : float
Value of the integral
sigma1 : float

This comment has been minimized.

@cdeil

cdeil Jul 27, 2014

Member

You can combine multiple parameters on one line ... I think here it makes sense:

sigma1, sigma2, w : float
    See model formula
SEC_TO_YEAR = Unit('second').to(Unit('year'))
YEAR_TO_SEC = 1. / SEC_TO_YEAR
from astropy.units import Unit, Quantity

This comment has been minimized.

@cdeil

cdeil Jul 27, 2014

Member

Remove Unit ... you don't use it.

-----
The spin-down luminosity is given by:
.. math::

This comment has been minimized.

@cdeil

cdeil Jul 27, 2014

Member

Should be OK to de-dent .. math:: by four spaces, no?

-----
The magnetic field is given by:
.. math::

This comment has been minimized.

@cdeil

cdeil Jul 27, 2014

Member

Again ... the .. math:: directive is unnecessarily indented by four spaces, I think.

def luminosity_spindown(self, t=None):
"""
Spin down luminosity (erg sec^-1) at age t (yr).

This comment has been minimized.

@cdeil

cdeil Jul 27, 2014

Member

No need to specify t (yr) ... it's a quantity now and doesn't need a specified input unit.

-----
The spin-down luminosity is given by:
.. math::

This comment has been minimized.

@cdeil

cdeil Jul 27, 2014

Member

Remove unneeded indentation of .. math::.

-----
The period is given by:
.. math::

This comment has been minimized.

@cdeil

cdeil Jul 27, 2014

Member

Dedent .. math::.

"""Total released energy (erg) at age t (yr).
def energy_integrated(self, t=None):
"""
Total released energy (erg) at age t (yr).

This comment has been minimized.

@cdeil

cdeil Jul 27, 2014

Member

Again ... no need to specify yr for t input because it's a quantity (I'll stop mentioning that now ... I see another case below where this can be removed).

__all__ = ['PWN']
from astropy.units import Unit, Quantity

This comment has been minimized.

@cdeil

cdeil Jul 27, 2014

Member

No need to import Unit (Eclipse shows that it's unused and thus a PEP8 violation)

class PWN(object):
"""
Simple pulsar wind nebula (PWN) evolution model.
Parameters

This comment has been minimized.

@cdeil

cdeil Jul 27, 2014

Member

Document all parameters here?

Reference: http://adsabs.harvard.edu/abs/2006ARA%26A..44...17G (Formula 8).
TODO: shouldn't we assume constant fraction of SNR radius?

This comment has been minimized.

@cdeil

cdeil Jul 27, 2014

Member

What about this "TODO"?
Do you have a suggestion what would be best here?

-----
During the free expansion phase the radius of the PWN evolves like:
.. math::

This comment has been minimized.

@cdeil

cdeil Jul 27, 2014

Member

De-dent .. math:: (I'll stop mentioning that now ... just check once in the html output that it really works before you go through and change it).

t = self.age
else:
raise ValueError('Need time variable or age attribute.')
# 4e3 years is a typical value that works for fsolve

This comment has been minimized.

@cdeil

cdeil Jul 27, 2014

Member

This comment is out of place here ... the fsolve call is further up.

# TODO: The following PWN model should be adapted to use gammafit classes.
class PWNElectronSpectrum(PWN):

This comment has been minimized.

@cdeil

cdeil Jul 27, 2014

Member

Can you open a new issue that this class needs to be finished or removed so that we don't forget about it?

__all__ = ['SNR', 'SNR_Truelove']
from astropy.units import Unit, Quantity

This comment has been minimized.

@cdeil

cdeil Jul 27, 2014

Member

Remove Unit.

Ejecta mass (g)
T_stop : float
t_stop : `~astropy.units.Quantity`
Post-shock temperature where gamma-ray emission stops (K)

This comment has been minimized.

@cdeil

cdeil Jul 27, 2014

Member

Remove (K) ... default input unit no longer needed because it's a Quantity now.

raise ValueError('Need time variable or age attribute.')
r = np.where(t > self.sedov_taylor_begin, self._radius_sedov_taylor(t).to('cm'),
self._radius_free_expansion(t).to('cm'))
return Quantity(r, Unit('cm'))

This comment has been minimized.

@cdeil

cdeil Jul 27, 2014

Member

No need to write Unit here ... this will work the same:

Quantity(r, 'cm')
t = self.age
else:
raise ValueError('Need time variable or age attribute.')
r = np.where(t > self.sedov_taylor_begin, self._radius_sedov_taylor(t).to('cm'),

This comment has been minimized.

@cdeil

cdeil Jul 27, 2014

Member

Add line break after comma? (I find this a bit more readable.)

tf = 200 * term1 * term2 * term3
def radius(self, t=None):
"""
Outer shell radius (pc) at age t (yr).

This comment has been minimized.

@cdeil

cdeil Jul 27, 2014

Member

pc -> cm (seems to have been changed to cm, no?)

return CM_TO_PC * np.select([t < t_core, t >= t_core], [R_1, R_RS])
def luminosity_tev(self, t=None, energy_min=Quantity(1, 'TeV')):
"""
Gamma-ray luminosity above 1 TeV (ph s^-1 cm^-2) at age t (yr).

This comment has been minimized.

@cdeil

cdeil Jul 27, 2014

Member

Document energy_min parameter and correct description here (It's not always "luminosity above 1 TeV" any more).

term_0 = energy_min / Quantity(1, 'TeV')
term_1 = self.e_sn / Quantity(1e51, 'erg')
term_2 = self.rho_ISM / (Quantity(1, 'cm^-3') * const.m_p)
L = self.theta * term_0 ** (-1.1) * term_1 * term_2

This comment has been minimized.

@cdeil

cdeil Jul 27, 2014

Member

Assign 1.1 to a variable, or maybe even write this?

spectral_index = 2.1

and then use it as term_0 ** (1 - spectral_index)?

Notes
-----
The beginning of the ST phase of the SNR is defined by the condition,

This comment has been minimized.

@cdeil

cdeil Jul 27, 2014

Member

ST -> Sedov-Taylor (abbrev is not defined and gain to introduce it is too low).

Notes
-----
The end of the ST phase of the SNR is defined by the condition, that the

This comment has been minimized.

@cdeil

cdeil Jul 27, 2014

Member

OK ... ST is used again here ... still I think it's better to be explicit and write Sedov-Taylor.

self.r_c = self.m_ejecta ** (1. / 3) * self.rho_ISM ** (-1. / 3)
self.t_c = self.e_sn ** (-1. / 2) * self.m_ejecta ** (5. / 6) * self.rho_ISM ** (-1. / 3)
def radius(self, t=None):

This comment has been minimized.

@cdeil

cdeil Jul 27, 2014

Member

Same comments apply here that I made for the other SNR class.

__all__ = ['coordinate_iau_format',
'ra_iau_format',
'dec_iau_format']
'dec_iau_format',
'as_quantity']

This comment has been minimized.

@cdeil

cdeil Jul 27, 2014

Member

This doesn't belong in gammapy.catalog ... should go somewhere in gammapy.utils ... maybe gammapy.utils.array is better than starting a new sub-module?
(This should go away anyways with Astropy 1.0 when Table / Quantity integration happens (hopefully)).

This comment has been minimized.

@cdeil

cdeil Jul 27, 2014

Member

Or maybe simply not add it to __all__ so that it doesn't show up in the docs, i.e. it's a "hidden function" and it can stay here?

r = np.sqrt(x ** 2 + y_prime ** 2)
v_glon = (-y_prime * vx + x * vy) / r ** 2
v_glat = (vz / (np.sqrt(1 - (z / d) ** 2) * d) - np.sqrt(vx ** 2 + vy ** 2 + vz ** 2) * z /

This comment has been minimized.

@cdeil

cdeil Jul 27, 2014

Member

I find such very long formulas very hard to read.
Python makes temp variables anyways, so without performance hit you can use temp variables (or split the expression into more lines if you don't like temp variables)

@@ -12,8 +12,14 @@ def normalize(func, x_min, x_max):
return normalized_func
def pdf(func):

This comment has been minimized.

@cdeil

cdeil Jul 27, 2014

Member

The name is not specific enough ... maybe this?

radial_density_to_pdf
radial_pdf_to_density
@cdeil

This comment has been minimized.

Member

cdeil commented Jul 27, 2014

Thanks for cleaning all this up!

I did a review just based on reading the code.

Feel free to merge after making travis-ci pass and addressing the comment, or ping me again if you want me to have a second look or actually try this stuff out.

@cdeil

This comment has been minimized.

Member

cdeil commented Jul 27, 2014

One thing that I don't like is that the radial distribution classes like FaucherKaspi2006 are simply listed in gammapy.astro.source and their name doesn't contain GalacticRadialDistribution or something that describes what they actually are.

The names get super-long if you write GalacticRadialDistributionFaucherKaspi2006, so I can see why you didn't do it, but still, it's not nice to have classes where the name doesn't describe at all what they are.
@adonath What do you think?

@cdeil

This comment has been minimized.

Member

cdeil commented Aug 4, 2014

@adonath - I think Ellis wants to produce the main results for the SciNeGhe proceeding this week ... maybe it's possible to merge this on Thursday or Friday latest? (you can always make new PRs or smaller changes directly in master)

Axel Donath added some commits Jul 17, 2014

adonath added a commit that referenced this pull request Aug 7, 2014

@adonath adonath merged commit f30420f into gammapy:master Aug 7, 2014

1 check passed

continuous-integration/travis-ci The Travis CI build passed
Details
@cdeil

This comment has been minimized.

Member

cdeil commented Aug 7, 2014

Thanks!

@cdeil cdeil changed the title from Galaxy modelling code revision to Re-write Galaxy modelling code Apr 8, 2015

@cdeil cdeil added the feature label Apr 8, 2015

@cdeil cdeil added this to the 0.1 milestone Apr 8, 2015

@adonath adonath deleted the adonath:astropy_population_models branch Nov 20, 2018

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