threeML / hawc_hal

HAWC Accelerated Likelihood - python-only framework for HAWC data analysis
BSD 3-Clause "New" or "Revised" License
11 stars 22 forks source link

Model residual maps #17

Closed cbrisboi closed 5 years ago

cbrisboi commented 5 years ago

This branch should now properly implement making model and residuals maps. I also have added a script to convert the hdf5 files to fits for compatibility with older hawc tools (this could probably use some better info from the origin maptrees, but it approximately works so far.

codecov[bot] commented 5 years ago

Codecov Report

Merging #17 into master will increase coverage by 0.08%. The diff coverage is 93.93%.

@@            Coverage Diff            @@
##           master     #17      +/-   ##
=========================================
+ Coverage   92.81%   92.9%   +0.08%     
=========================================
  Files          39      39              
  Lines        1532    1592      +60     
=========================================
+ Hits         1422    1479      +57     
- Misses        110     113       +3
cbrisboi commented 5 years ago

I made a commit to reflect your feedback, before realizing that its much easier to read on github then via email. After discussing with @henrikef I also made a model map(but did not commit it) which might be used to produce some unit tests using the already provided detector response file in hal.

colasri commented 5 years ago

There is a map already in pa-pub.umd.edu:/data/disk01/home/giacomov/test_data, which is used in the HAL test:

# Test the HAWC plugin
# First copy the data in /data/disk01/home/giacomov/test_data at UMD to
# the directory test_data
mkdir hawc_test_data
rsync --progress 'pa-pub.umd.edu:/data/disk01/home/giacomov/test_data/*' hawc_test_data
export HAWC_3ML_TEST_DATA_DIR=$SOFTWARE_BASE/3ML_tests/hawc_test_data
pytest -vv --pyargs threeML -k "hawc"

Could you use the same map for your test, or do you really need a new file?

cbrisboi commented 5 years ago

Colas - I don't as yet have a test. It was something @henrikef told me I should perhaps look into for the future.

I was talking to @henrikef earlier today, she suggested I generate the hd5 file for the test directly, and so that would only require adding the script, rather than a whole new map. The idea would be to check that the model/residual maps were generated properly, since I haven't looked at those maps I cant say for sure. I don't have a preference per se, but generating the needed maptree on the fly might be a good plan for bandwidth and things. I'm not sure what the best plan would be.

giacomov commented 5 years ago

In general a test should accomplish two things:

While for the first point any data is good, for the second one you need a dataset that yields a predictable answer, so you can check that the answer is what you expect. I have no preference on whether to use or not a new map. In the latter case, please do make it hd5 and with a fairly small ROI (< 20 deg) so that the size of the file is small.

I will wait for the implementation of the tests before merging this in.

cbrisboi commented 5 years ago

My idea for a test is to generate a model map using the geminga_maptree for background and the detector response file we have in hal already. Thats slightly more portable than adding a whole map. I'm thinking to just read in the file and compare that the model that was saved is the same as what is read back in. Does that sound reasonable?

giacomov commented 5 years ago

Yes, it does! Make sure to cover in the test also corner cases and exceptions. Have a look at the other tests for inspiration, and feel free to ask questions either here or on Slack.

giacomov commented 5 years ago

@cbrisboi I am looking at your test so far. I understand that it is work in progress, but I thought I can give you some more pointers to avoid too much work on your side.

Now you can understand that all the code you have up until line 33 can be avoided.

You can get exactly the same by doing:

def test_model_residual_maps(geminga_roi, geminga_maptree, geminga_response):

    hawc = HAL("HAWC", geminga_maptree, geminga_response, geminga_roi)

    # Now you can continue as before
    hawc.set_active_measurements(1, 9)

    # Here you can put the rest of your code
    ...
cbrisboi commented 5 years ago

@giacomov I'm not sure that I fully addressed the your statement that the code up to line 33 can be avoided, but I hope what I did commit is clear. The test is in a function now, and uses a fixture as input.