Skip to content
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 multi-scale TS image computation #234

Merged
merged 3 commits into from Feb 5, 2015

Conversation

@adonath
Copy link
Member

@adonath adonath commented Jan 24, 2015

I now there is still a bit of work to be done here (put the script in a main function etc., change to back astropy.nddata.utils.extract_array, docs missing...), however I think its ready for a first review.

return C_0 - C_1, amplitude, niter


def _amplitude_bounds(counts, background, model):
"""

This comment has been minimized.

@cdeil

cdeil Jan 24, 2015
Member

Add one-line description what this function does.

if scale > 0:
from astropy.convolution import convolve
sigma = scale / (BINSZ * factor)
if morphology == 'Gaussian2D':

This comment has been minimized.

@cdeil

cdeil Jan 24, 2015
Member

Always end such if -elif ... with an else and throw a ValueError.
If you don't the user will mistype the string and then your function will fail in weird hard-to-understand ways later.

def compute_ts_map(counts, background, exposure, kernel, flux=None,
method='root', optimizer='Brent', parallel=True,
def compute_ts_map_multiscale(maps, psf_parameters, scales=[0], downsample='auto',
residual=False, morphology='Gaussian2D', width=0.3,

This comment has been minimized.

@cdeil

cdeil Jan 24, 2015
Member

The width parameter is not documented.
Probably 0.3 is something specific to your application ... if so this parameter should not have a default.

You should also say what happens to args and kwargs in a sentence ... if those are really needed ... it usually makes for complex code to be flexible in that way.

for scale, factor in zip(scales, downsample):
logging.info('Computing (residual) TS map for scale {0} deg'.format(scale))

# Sample down and require that scale parameters is at least 5 pix

This comment has been minimized.

@cdeil

cdeil Jan 24, 2015
Member

This function is quite long ... ideally every function should fit on a screen.

Can you split it up, e.g. extract this part into a _downsample helper function?

Number of iterations map
runtime : float
Time needed to compute TS map.
"""
def read(self, filename):

This comment has been minimized.

@cdeil

cdeil Jan 24, 2015
Member

If possible I'd like read to always be a @staticmethod like Table.read in astropy or EnergyDependentTablePSF.read in Gammapy.

shape = maps[0].data.shape
logging.info('Reading {0}'.format(args.input))
maps = fits.open(args.input)
psf_parameters = json.load(open(args.psf))

This comment has been minimized.

@cdeil

cdeil Jan 24, 2015
Member

Add logging.info that you're reading the PSF file.

@@ -36,89 +45,48 @@
help='Input excess model FITS file name')
parser.add_argument('--threshold', type=float, default=None,
help="Minimal required initial (before fitting) TS value, that the fit is done at all.")
parser.add_argument('--overwrite', action='store_true',
help='Overwrite existing output file?')
parser.add_argument('--clobber', action='store_true',

This comment has been minimized.

@cdeil

cdeil Jan 24, 2015
Member

I've opted for consistent use of overwrite in Gammapy instead of clobber.
I was 50:50 on this, but astropy.table.Table uses overwrite and I think the plan is to use it consistently in Astropy, even if now astropy.io.fits uses clobber.
OK?

parser.add_argument('--scale', type=float, default=0,
help='Scale of the kernel in deg.')
parser.add_argument('--morphology', type=str, default='Gaussian2D',
help="Which source morphology to use for TS calculation. Either 'Gaussian2D' or 'Shell2D'.")

This comment has been minimized.

@cdeil

cdeil Jan 24, 2015
Member

Wrap line

@cdeil
Copy link
Member

@cdeil cdeil commented Jan 24, 2015

I've left some stylistic comments inline.

You need to fix the tests:
https://travis-ci.org/gammapy/gammapy/jobs/48174547#L2459

Once that passes and these little things are fixed it's up to you when you want to merge.
I'm working on spectra and gammapy.data classes, so no rush from my side.

It would be nice if you moved the script to gammapy/scripts and put it in functions so that it's re-usable and shows up in the docs ....

@@ -196,7 +319,7 @@ def _ts_value(position, counts, exposure, background, kernel, flux,
Background map.
exposure : array
Exposure map.
kernel : `astropy.convolution.Kernel2D`
kernel : astropy.convolution.core.Kernel2D

This comment has been minimized.

@cdeil

cdeil Jan 25, 2015
Member

Why did you change this?
I think the old version is the correct one that creates a HTML link and the new one doesn't?
(see https://gammapy.readthedocs.org/en/latest/development/index.html#where-should-i-import-from)

This comment has been minimized.

@adonath

adonath Jan 26, 2015
Author Member

I had merge conflicts and just copied the files, I guess that's where this "backward" changes come from. Thanks for the comments so far! I will address them soon and also add an example to gammapy-extra.

import os

from astropy.utils.exceptions import AstropyDeprecationWarning
warnings.filterwarnings('ignore', category=AstropyDeprecationWarning)

This comment has been minimized.

@hamogu

hamogu Jan 26, 2015

If this is brand-new code, is it not dangerous to ignore astropy depreciation warnings? If you use a feature that's deprecated already, there is a good chance it will be removed in one of the next astropy versions.

This comment has been minimized.

@cdeil

cdeil Jan 26, 2015
Member

Agreed ... since we only support Astropy 1.0 (master at the moment) in Gammapy it should be possible to write code that doesn't raise any deprecation warnings.
@adonath If you get deprecation warnings from unrelated parts of Gammapy you don't have to fix them in this PR, but please don't silence them either.

This comment has been minimized.

@adonath

adonath Jan 26, 2015
Author Member

Yes, there were unrelated deprecation warnings and I included this local fix to have a clean output of the command line tool, when I was testing. Just forgot to remove it...Thanks!

exposure_slice = extract_array(exposure, kernel.shape, position).astype(float)
counts_slice = extract_array_2d(counts, kernel.shape, position).astype(float)
background_slice = extract_array_2d(background, kernel.shape, position).astype(float)
exposure_slice = extract_array_2d(exposure, kernel.shape, position).astype(float)

This comment has been minimized.

@hamogu

hamogu Jan 26, 2015

Can you tell me why you switched back from astropy.nddata.utils.extract_array to imageutils.extract_array_2d? imageutils has been merged with astropy and will not be developed actively any longer. I'm currently working on a PR to improve astropy.nddata.utils.extract_array so please let me know if you found any differences between the two that matter for you.

This comment has been minimized.

@adonath

adonath Jan 26, 2015
Author Member

This change was not intended. As I had some merge conflicts, I just copied the changed files, but obviously they included some old stuff. Thanks for the hint!

@adonath adonath force-pushed the adonath:ts_computation_multiscale branch 4 times, most recently from 1a0298a to 2bab6fb Jan 27, 2015
+ kernel.shape[0] / 2)
source_kernel = Model2DKernel(model, x_size=x_size, mode='oversample')
else:
raise ValueError('Unknown morphology model.')

This comment has been minimized.

@cdeil

cdeil Jan 30, 2015
Member

Change to:

ValueError('Unknown morphology: {}'.format(morphology))
index = np.where(ts[:, :, i] == ts_max)
scale_max[index] = scale
niter_max[index] = niter[:, :, i][index]
amplitude_max[index] = amplitude[:, :, i][index]

This comment has been minimized.

@cdeil

cdeil Jan 30, 2015
Member

Add empty line.

for name, order in zip(['ts', 'amplitude', 'niter'], [1, 1, 0]):
ts_results[name] = upsample_2N(ts_results[name], factor,
order=order, shape=shape)
multiscale_result.append(ts_results)

This comment has been minimized.

@cdeil

cdeil Jan 30, 2015
Member

Add empty line.

if niter > MAX_NITER:
logging.warning('Exceeded maximum number of function evaluations!')
return np.nan, amplitude, niter
raise ValueError('Invalid fitting method.')

This comment has been minimized.

@cdeil

cdeil Jan 30, 2015
Member

Add empty line.


def _amplitude_bounds(counts, background, model):
"""
Compute bounds for the root of `f_cash_root`.

This comment has been minimized.

@cdeil

cdeil Jan 30, 2015
Member

for the -> for the

@@ -375,6 +376,47 @@ def coordinates(image, world=True, lon_sym=True, radians=False, system=None):
return lon, lat


def detect_peaks(image, threshold, mask=None, world=True):

This comment has been minimized.

@cdeil

cdeil Jan 30, 2015
Member

Maybe wait with this PR until the photutils version is merged in photutils and then use it here?

from __future__ import (absolute_import, division, print_function,
unicode_literals)
from ..utils.scripts import get_parser
import os

This comment has been minimized.

@cdeil

cdeil Jan 30, 2015
Member

Always sort imports:

  • Python stdlib
  • numpy
  • scipy
  • Astropy
  • Gammapy
ts_image(**vars(args))


def ts_image(input, psf, model, scales, downsample, residual, morphology,

This comment has been minimized.

@cdeil

cdeil Jan 30, 2015
Member

input is a Python keyword ... you shouldn't use it for variable names.
Maybe use data here or data_filename?

@cdeil
Copy link
Member

@cdeil cdeil commented Jan 30, 2015

I've left some cosmetic inline comments ... this looks good to me know, but tests should be added and if possible the peak finder from photutils should be used.

Feel free to merge this yourself next week!

@adonath adonath force-pushed the adonath:ts_computation_multiscale branch from 2bab6fb to b2778c8 Feb 5, 2015
adonath added a commit that referenced this pull request Feb 5, 2015
Multiscale ts image computation
@adonath adonath merged commit 67c0bad into gammapy:master Feb 5, 2015
1 check passed
1 check passed
continuous-integration/travis-ci The Travis CI build passed
Details
@cdeil cdeil mentioned this pull request Feb 5, 2015
6 of 10 tasks complete
@cdeil cdeil changed the title Multiscale ts image computation Add multi-scale TS image computation Apr 8, 2015
@cdeil cdeil added the feature label Apr 8, 2015
@cdeil cdeil added this to the 0.2 milestone Apr 8, 2015
@adonath adonath deleted the adonath:ts_computation_multiscale branch Nov 20, 2018
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Linked issues

Successfully merging this pull request may close these issues.

None yet

3 participants