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

Fix and improvements to TS map tool #456

Merged
merged 7 commits into from Feb 23, 2016

Conversation

Projects
None yet
3 participants
@OlgaVorokh
Member

OlgaVorokh commented Feb 17, 2016

Working on my first task: fixing documentation of https://gammapy.readthedocs.org/en/latest/detect/index.html

The python code works well, but the command line tool
$ cd gammapy-extra/datasets/fermi_survey
$ gammapy-ts-image all.fits.gz ts_map_0.00.fits --scale 0
Needs more changes than just renaming.

I can attach other fixes for this documentation page to this pull request later or make a new one.

@OlgaVorokh

This comment has been minimized.

Owner

OlgaVorokh commented on gammapy/scripts/image_ts.py in 7165503 Feb 17, 2016

it conflicted with the signature of image_ts

def image_ts(input_file, output_file, psf, model, scales, downsample, residual,
morphology, width, overwrite):

I did a small investigation: the image_ts used to have a threshold parameter, but then it got deleted(not sure when since file was renamed and github.com history only works if you don't move file around)
I found an old commit gammapy@67c0bad where compute_ts_map have been replaced by compute_ts_map_multiscale, which doesn't have "threshold" parameter.
Based on this I believe threshold is not needed anymore.

@OlgaVorokh

This comment has been minimized.

Owner

OlgaVorokh commented on gammapy/scripts/image_ts.py in 7165503 Feb 17, 2016

Path() doesn't have split() method. Moreover, I am not sure why we even create folder with exist_ok=False, because in the tutorial https://gammapy.readthedocs.org/en/latest/detect/index.html

gammapy-ts-image all.fits.gz ts_map_0.00.fits --scale 0

the folder is always going to be our current working directory, thus failing

@adonath

This comment has been minimized.

Member

adonath commented Feb 17, 2016

Thanks for working on this! compute_ts_map was not replaced, but rather compute_ts_map_multiscale was added (it is still part of the code base, if you take a look...) and is the main function called in the command line tool. The threshold parameter was obviously forgotten to pass to the image_ts function. Could you fix that instead of removing it? (I must admit that it was probably forgotten, because it wasn't used at all... but it's still good to have it)

The command line tool outputs either one file or several files, depending on the number of scales, that are given. The way the code used to work was, that if several scales are given, the parent folder of the output file was taken and all output files were written in there, for convenience the mkdir command was added to create the parent folder, if it doesn't exist yet. This way one could specify the output file like gaussian/result.fits, without that the gaussian folder existed yet. Maybe you can try to get that behaviour again using the Path() class. Here's the pathlib documentation as a starting point: https://docs.python.org/3/library/pathlib.html#

@OlgaVorokh

This comment has been minimized.

Member

OlgaVorokh commented Feb 17, 2016

Thanks for the quick reply!

By "replaced" I meant that in the main function image_ts(), located in https://github.com/gammapy/gammapy/blob/master/gammapy/scripts/image_ts.py we used to call compute_ts_map and now call compute_ts_map_multiscale instead after your old pull request :)

But according to code and documentation, unlike compute_ts_map http://gammapy.readthedocs.org/en/latest/api/gammapy.detect.compute_ts_map_multiscale.html doesn't accept threshold parameter

@adonath

This comment has been minimized.

Member

adonath commented Feb 17, 2016

Ok, I see what you mean. If you take a look again at the code, compute_ts_map_multiscale takes *args and **kwargs that are directly passed to compute_ts_map, so there's no need to repeat the arguments like method, threshold or optimizer. But I just noticed that this behaviour is not documented... should maybe be added.

@cdeil cdeil added this to the 0.4 milestone Feb 18, 2016

@cdeil

This comment has been minimized.

Member

cdeil commented Feb 18, 2016

@OlgaVorokh - Thanks for the pull request!

@adonath - Since you wrote this code started reviewing this, I'll stay out of it ... assigning this PR to you. FYI: @OlgaVorokh is considering applying for the analysis project for GSoC with Gammapy, and I suggested to her trying out the example from the gammapy.detect docs page and fixing some small things as a first pull request, like here the command line tool is called gammapy-ts-image, which should be changed to the current name gammapy-image-ts.

@OlgaVorokh and @adonath - Up to you how far you want to take this PR (just the docs page fix, or further code or test changes).

@adonath

This comment has been minimized.

Member

adonath commented Feb 18, 2016

@OlgaVorokh

If you want to add a unit test, I would suggest the following:

Currently there's only one test for the compute_ts_map function. It uses a test dataset, that is provided here and asserts that the correct numbers are obtained. Now it would make sense to write a test, that calls the command line tool's main function (image_ts_main), with the corresponding arguments and write the result to a temporary file. Then open the file and check that the numbers from test_compute_ts_map, mentioned above, are reproduced. You can find an example how to do that here. The main task will be to get the test data into the right format i.e. a single fits file. To achieve this you can either write a helper function in the test file or extend this script and push the file to the gammapy-extra repository. Let us know, when have further questions or got stuck.

@adonath adonath referenced this pull request Feb 18, 2016

Closed

Further TS map improvements #223

6 of 10 tasks complete
@OlgaVorokh

This comment has been minimized.

Member

OlgaVorokh commented on gammapy/detect/test_statistics.py in c07509a Feb 19, 2016

it used to fail on the example dataset because positions were empty -- usually kernel was too big in one dimension, like (50, 50) for image (25, 200)

This comment has been minimized.

Member

adonath replied Feb 19, 2016

OK, that shouldn't happen. Can you add the code that fails as a test to this PR? And we'll see how to fix it.

@adonath

This comment has been minimized.

Member

adonath commented on gammapy/scripts/image_ts.py in c07509a Feb 19, 2016

You still need to pass the threshold parameter here as a keyword argument.

@adonath

This comment has been minimized.

Member

adonath commented on gammapy/detect/test_statistics.py in c07509a Feb 19, 2016

except statements should be always used with the corresponding error. Like e.g. except ValueError:, otherwise the code can silently fail and doesn't give a traceback to the user. In this case the try and except statement even covers ~50-60 line of code, where in principle in any line could occur an error. By catching any exception the information on what went wrong is hidden from the user, which is bad. Can you remove that overall try and except statement? We might add it later to the corresponding line, where the error occurs, if necessary. Please add the test first.

@cdeil

This comment has been minimized.

Member

cdeil commented Feb 19, 2016

@OlgaVorokh - FYI: I'll make the Gammapy 0.4 release this Sunday.
Having this in would be nice, but if you don't have time in the coming days, no worries, it'll go in 0.5 in a few months.

@OlgaVorokh

This comment has been minimized.

Member

OlgaVorokh commented Feb 19, 2016

@cdeil I don't have much time this weekend, but I'll try :)

@adonath you are right, I used this mainly for debug reasons in case if there are more issues(and I am logging the exception into INFO).

Unfortunately I don't have time right now(I'll add test in ~14 hours), but you can reproduce the error by cloning this branch and running
cd gammapy/gammapy-extra/datasets/fermi_survey/
gammapy-image-ts all.fits.gz output/ts_map_0.00.fits --scales 0 1 2 3

--scales 0 and 1 work, 2 and 3 fail and you will see the exception in the log
INFO Exception during computation of scale 3.0, perhaps the kernel size is bigger than the imageand positions are empty because of that [gammapy.detect.test_statistics]
Traceback (most recent call last):
File "/home/ubuntu/gammapy/gammapy/detect/test_statistics.py", line 215, in compute_ts_map_multiscale
kernel, _args, *_kwargs)
File "/home/ubuntu/gammapy/gammapy/detect/test_statistics.py", line 381, in compute_ts_map
j, i = zip(*positions)

@adonath

This comment has been minimized.

Member

adonath commented Feb 19, 2016

@OlgaVorokh Yes, most probably the kernel is to large. The test data is set up with 200x200 pixel and a pixel size of 0.02 deg / pixel. I.e. the total image covers a region of 4x4 deg. The scale parameter is given in deg. So rather choose something smaller like 0, 0.05 0.1, 0.2. For the test a single scale of 0 is already sufficient at first.

@OlgaVorokh

This comment has been minimized.

Member

OlgaVorokh commented Feb 19, 2016

@adonath Yes, this makes sense. I've actually made a debug output to print image and kernel sizes to be absolutely sure.

So what should behaviour be when user specifies several scales and some of them are too large? Right now it fails with cryptic message on line j,i = zip(*positions), so I think we could put call compute_ts_map into try/except block and catch this kind of exceptions, like I am doing now, but with more specific catch?

@adonath

This comment has been minimized.

Member

adonath commented Feb 20, 2016

I don't think it's necessary to catch this error at all. The normal and intended use case is, for sure, that the kernel is (much) smaller than the input data. If necessary, one could require that the data shape is larger by a factor of 2 in every dimension then the kernel shape. So I'd rather prefer to check the inputs at the beginning of compute_ts_map, and give a specific error message, if the kernel is too large, and describe this requirement in the docstring. I definitely want to avoid the the try and except statement around the whole compute_ts_map function, because the algorithm is rather complex and could fail at different places with the same error.

@OlgaVorokh

This comment has been minimized.

Owner

OlgaVorokh commented on docs/detect/index.rst in c37b02f Feb 22, 2016

fails with ubuntu@ip-172-31-15-21:~/gammapy/gammapy-extra/datasets/source_diffuse_separation/galactic_simulations$ gammapy-detect-iterative --counts fermi_counts.fits --background fermi_diffuse.fits --exposure fermi_exposure_gal.fits output_fits output_regions

^[[1;2DINFO:gammapy.scripts.detect_iterative:Reading counts map: fermi_counts.fits

INFO:gammapy.scripts.detect_iterative:Reading background map: fermi_diffuse.fits

INFO:gammapy.scripts.detect_iterative:Reading exposure map: fermi_exposure_gal.fits

INFO:gammapy.scripts.detect_iterative:Number of scales: 3

INFO:gammapy.scripts.detect_iterative:DEG_PER_PIX: 0.1

INFO:gammapy.scripts.detect_iterative:Scales in deg: [0.1, 0.2, 0.4]

INFO:gammapy.scripts.detect_iterative:Scales in pix: [ 1. 2. 4.]

Traceback (most recent call last):

File "/home/ubuntu/anaconda2/envs/py27/bin/gammapy-detect-iterative", line 9, in

load_entry_point('gammapy', 'console_scripts', 'gammapy-detect-iterative')()

File "/home/ubuntu/gammapy/gammapy/scripts/detect_iterative.py", line 30, in detect_iterative_main

detect_iterative(**vars(args))

File "/home/ubuntu/gammapy/gammapy/scripts/detect_iterative.py", line 71, in detect_iterative

detector.run()

File "/home/ubuntu/gammapy/gammapy/detect/iterfind.py", line 113, in run

self.compute_iter_maps()

File "/home/ubuntu/gammapy/gammapy/detect/iterfind.py", line 151, in compute_iter_maps

background += self.model_excess(self.sources)

ValueError: operands could not be broadcast together with shapes (103,2001) (2,101,2001) (103,2001)

@OlgaVorokh

This comment has been minimized.

Owner

OlgaVorokh commented on docs/detect/index.rst in c37b02f Feb 22, 2016

Also fails with ubuntu@ip-172-31-15-21:~/gammapy/gammapy-extra/test_datasets/unbundled/poisson_stats_image$ gammapy-detect-iterative --background background.fits.gz --counts counts.fits.gz --exposure exposure.fits.gz output_fits/ output_regions/

INFO:gammapy.scripts.detect_iterative:Reading counts map: counts.fits.gz

INFO:gammapy.scripts.detect_iterative:Reading background map: background.fits.gz

INFO:gammapy.scripts.detect_iterative:Reading exposure map: exposure.fits.gz

INFO:gammapy.scripts.detect_iterative:Number of scales: 3

INFO:gammapy.scripts.detect_iterative:DEG_PER_PIX: 0.02

INFO:gammapy.scripts.detect_iterative:Scales in deg: [0.1, 0.2, 0.4]

INFO:gammapy.scripts.detect_iterative:Scales in pix: [ 5. 10. 20.]

Traceback (most recent call last):

File "/home/ubuntu/anaconda2/envs/py27/bin/gammapy-detect-iterative", line 9, in

load_entry_point('gammapy', 'console_scripts', 'gammapy-detect-iterative')()

File "/home/ubuntu/gammapy/gammapy/scripts/detect_iterative.py", line 30, in detect_iterative_main

detect_iterative(**vars(args))

File "/home/ubuntu/gammapy/gammapy/scripts/detect_iterative.py", line 71, in detect_iterative

detector.run()

File "/home/ubuntu/gammapy/gammapy/detect/iterfind.py", line 113, in run

self.compute_iter_maps()

File "/home/ubuntu/gammapy/gammapy/detect/iterfind.py", line 151, in compute_iter_maps

background += self.model_excess(self.sources)

TypeError: Cannot cast ufunc add output from dtype('float64') to dtype('>i2') with casting rule 'same_kind'

@OlgaVorokh

This comment has been minimized.

Owner

OlgaVorokh commented on docs/detect/index.rst in c37b02f Feb 22, 2016

I think I can fix the second one, but I just wantet to ask @cdeli :)

@OlgaVorokh

This comment has been minimized.

Member

OlgaVorokh commented Feb 22, 2016

Oh, made a typo in username.

@cdeil, I copied documentation for gammapy-detect-iterative, but it doesn't work on samples with different errors(mentioned above), so I wanted to ask how to proceed with this -- is there are testing data for gammapy-detect-iterative or should I attempt to make it work on theses testcases?

@adonath Please review whenever you have a free minute :)

@cdeil

This comment has been minimized.

Member

cdeil commented Feb 23, 2016

This PR looks good to me.

@adonath - Ready to merge?

The iterative source detect needs more work and should be done in a separate PR.
But I think it's OK to keep the description in the docs for now, no need to remove it again.

As I wrote here
https://github.com/gammapy/gammapy/blob/master/gammapy/detect/iterfind.py#L13
the iterfind code was an afternoon hack to try it out.
I don't even remember if I finished it to a state where it ever worked.

But such algorithms are things we want in gammapy.detect and that we'd like to work on with a GSoC student this summer. So @OlgaVorokh, feel free to try to make it work if you have time. But please also email me offline when you're available for a chat about GSoC in general.

Concerning a test dataset, the Fermi image is pretty realistic for what we'd like to analysis.
It's not ideal for debugging and a test here:
https://github.com/gammapy/gammapy/blob/master/gammapy/detect/tests/test_iterfind.py#L8
because the correct answer isn't known.

So for debugging / testing I'd suggest to simulate an image with a few Gaussian sources (start with one) with known location, flux, extension on flat background and see if you can find those and recover their source parameters correctly.
There's some nice utility functions to make test catalogs and images of Gaussian sources here.

@adonath

This comment has been minimized.

Member

adonath commented Feb 23, 2016

@OlgaVorokh Thanks! Looks good to me as well. Are you still interested in implementing a test for the command line tool? If so, I would suggest to do this in a separate PR so we can merge this PR.

cdeil added a commit that referenced this pull request Feb 23, 2016

Merge pull request #456 from OlgaVorokh/documentation_fixes
Fix missing parameter and call of non-existing method in image_ts.py

@cdeil cdeil merged commit 192d6ab into gammapy:master Feb 23, 2016

2 checks passed

continuous-integration/appveyor/pr AppVeyor build succeeded
Details
continuous-integration/travis-ci/pr The Travis CI build passed
Details

@cdeil cdeil changed the title from Fix missing parameter and call of non-existing method in image_ts.py to Fix and improvements to TS map tool Feb 23, 2016

@cdeil

This comment has been minimized.

Member

cdeil commented Feb 23, 2016

@OlgaVorokh - Thanks!

I've added a changelog entry and added you as Gammapy contributor in c30bf84

If I run the example I get this warning:

eduroam-3-163:fermi_survey deil$ time gammapy-image-ts all.fits.gz ts_map_0.00.fits --scale 0 
INFO     Reading all.fits.gz [gammapy.scripts.image_ts]
INFO     Reading psf.json [gammapy.scripts.image_ts]
INFO     Computing TS map for scale 0.000 deg and Gaussian2D morphology. [gammapy.detect.test_statistics]
INFO     No down sampling used. [gammapy.detect.test_statistics]
INFO     Using 8 cores to compute TS map. [gammapy.detect.test_statistics]
/Users/deil/Library/Python/3.5/lib/python/site-packages/gammapy-0.4.dev2240-py3.5-macosx-10.11-x86_64.egg/gammapy/detect/test_statistics.py:380: RuntimeWarning: invalid value encountered in multiply
  TS = np.empty(counts.shape) * np.nan
/Users/deil/Library/Python/3.5/lib/python/site-packages/gammapy-0.4.dev2240-py3.5-macosx-10.11-x86_64.egg/gammapy/detect/test_statistics.py:381: RuntimeWarning: invalid value encountered in multiply
  amplitudes = np.empty(counts.shape) * np.nan
/Users/deil/Library/Python/3.5/lib/python/site-packages/gammapy-0.4.dev2240-py3.5-macosx-10.11-x86_64.egg/gammapy/detect/test_statistics.py:382: RuntimeWarning: invalid value encountered in multiply
  niter = np.empty(counts.shape) * np.nan
INFO     TS map computation took 11.3 s 
 [gammapy.detect.test_statistics]
INFO     Writing ts_map_0.00.fits [gammapy.scripts.image_ts]

real    0m13.265s
user    1m21.671s
sys 0m1.942s

I don't see the warning if I just run this:

$ python -c 'import numpy as np; np.empty((3,2)) * np.nan'

@adonath @OlgaVorokh - if you can reproduce the warning or have an idea how to get rid of it, please do. If not, no worries ... it's not important.

@cdeil

This comment has been minimized.

Member

cdeil commented Feb 23, 2016

I now do see the warning if I have more pixels in the empty array:

$ python -c 'import numpy as np; np.empty((100, 300)) * np.nan'
-c:1: RuntimeWarning: invalid value encountered in multiply

I'm changing this from np.empty to np.ones directly in master now to get rid of the warning.

@adonath

This comment has been minimized.

Member

adonath commented Feb 23, 2016

I don't see the warning with Python 3.5 and Numpy 1.10.4. With Numpy > 1.8 I think you can use np.full(shape, np.nan), which is a bit more elegant. Or np.tile(np.nan, shape) should work as well for older numpy versions.

@cdeil

This comment has been minimized.

Member

cdeil commented Feb 23, 2016

I think the warning can appear if there's strange floats in the array allocated with np.empty.
So it's normal that the issue only appears sometimes, maybe even only on some machines.

I've decided to use np.ones() * np.nan -- see 5e03244

@OlgaVorokh OlgaVorokh deleted the OlgaVorokh:documentation_fixes branch Feb 23, 2016

@OlgaVorokh

This comment has been minimized.

Member

OlgaVorokh commented Feb 24, 2016

@adonath I just realized from reading your comment that I forgot to actually add my test to the commit. Oh well, I'll open a new PR

@cdeil

This comment has been minimized.

Member

cdeil commented Feb 24, 2016

If you do open a new PR, could you do some further cleanup?

I saw some code for TS computation using iminuit as an option, but this should probably be removed.
(If it's working, OK to keep as one of the options, but it's probably broken and not worth having as another option, right?)

@adonath

This comment has been minimized.

Member

adonath commented Feb 24, 2016

@cdeil, @OlgaVorokh As we still see some artefacts at the boundaries of some of our maps, I'd like to keep iminuit for the moment as a debug option (that's how it used it in the past). Later it can be most probably removed, yes.

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