Continuum BG Estimation

[ckarwin]

This is an initial method for estimating the background for continuum sources. It is still highly preliminary, and so I am opening this as a PR draft, in order to share the example notebook which has the initial algorithm hard-coded (available here). I plan to develop the ContinuumEstimation class leading up to the October meeting, in order to at least have a working version that can then be further develop. I think the ContinuumEstimation class should be part of a background_estimation base class, which will also contain a LineEstimation class and a BurstEstimation class. But we should discuss this and figure out what makes the most sense.

Here are slides summarizing the initial results from the method presented in the example tutorial:
continuum_bg_estimation.pdf. Below is a summary of the next steps:

Improving algorithm:

• Optimize cut for mask selection (so far only tested cumdist < 0.4)

• Use actual ARM measurement for mask selection

• Use more sophisticate inpainting method:

  • Inpainting Galactic Foreground Intensity and Polarization Maps Using Convolution Neural Networks (Puglisi+20, link)
    • Tested for small angular scales, but what about large scales?
  • SMILE inpainting algorithm: (Li 2014, link)
    • Based on solving the Laplace equation for each missing pixel by using the values of the surrounding pixels. Might not work so well for our purposes.
    • Also used in Macias+19 (link).
  • We used CNN for a somewhat similar problem: Deep learning models of the discrete component of the Galactic interstellar gamma-ray emission (link).
  • Develop our own, optimized for our specific problem (maybe employing some ML algorithm).
  • Other methods?

• Give more freedom to fit (first need to fix unrelated convergence problem).

• Vectorize code

More tests:

• Test method on extended/diffuse source.
• Test Crab fit with GDCE included in total BG.
• Test Crab fit with other point source(s) included in BG.
• Test Crab fit will all other astrophysical sources included in BG.

[codecov]

Codecov Report

Attention: Patch coverage is 64.28571% with 50 lines in your changes missing coverage. Please review.

Project coverage is 74.25%. Comparing base (0ea3e15) to head (b3877c3).
Report is 39 commits behind head on develop.

Files with missing lines	Patch %	Lines
...osipy/background_estimation/ContinuumEstimation.py	63.76%	50 Missing ⚠️

Files with missing lines	Coverage Δ
cosipy/__init__.py	100.00% <100.00%> (ø)
cosipy/background_estimation/__init__.py	100.00% <100.00%> (ø)
...osipy/background_estimation/ContinuumEstimation.py	63.76% <63.76%> (ø)

[hiyoneda]

load (typo?)

[ckarwin]

Yup, thanks! Fixed now.

[hiyoneda]

I have reviewed the code. I added some suggestions about the code.

[hiyoneda]

In the future (not in this PR), we can move this part outside this class because the point source calculation is also performed in several different classes, e.g., point source injector.

[ckarwin]

Yes, that makes sense.

[hiyoneda]

Can you use sc_orientation directly instead of the path to an orientation file?

[ckarwin]

Ok, done.

[hiyoneda]

Can you return psr here? I think it would be better than saving a PSR file outside this function: users receive psr from this function and save it by themselves if they want.

[ckarwin]

Yes, done.

[hiyoneda]

-> "...loading the pre-computed point source response ..." ?

[ckarwin]

Changed.

[hiyoneda]

if this_index+j >= self.psr.axes['PsiChi'].nbins-1:

[ckarwin]

More elegant, thanks! Changed.

[hiyoneda]

Is it possible to return arm_mask and sorted_indices? Assigning them as member parameters might be confusing because it will be hard to know from which sliced data self.arm_mask and sorted_indices are generated. My suggestion is that arm_mask and sorted_indices are explicitly generated for each loop in continuum_bg_estimation and they are used as input parameters of simple_inpainting.

[ckarwin]

I made your suggested changes.

[hiyoneda]

Following the above comment, can you add arm_mask and sorted_indices as input parameters?

[ckarwin]

Done.

[hiyoneda]

Can you move the line 92 inside continuum_bg_estimation because the background model is not yet generated at this point? And, maybe, estimated_bg can be just an output of the function continuum_bg_estimation.

[ckarwin]

Done.

[hiyoneda]

Can you return estimated_bg here? It would be better that estimated_bg is saved by users outside this function.

[ckarwin]

Done.

[ckarwin]

Thanks again for your review, @hiyoneda! It was very helpful. I implemented all your comments, and I also added unit tests. I think it can be merged if you think it's ok.

[hiyoneda]

@ckarwin Thanks. It looks good. I have two more comments on the notebook. After they are reflected, I think that it can be merged.

1. On the 3rd cell, the response file is downloaded, but the file is not consistent with the one we load later. Should it be SMEXv12.Continuum.HEALPixO3_10bins_log_flat.binnedimaging.imagingresponse.nonsparse_nside8.area.good_chunks_unzip.earthocc.zip?
2. At the last 2nd cell, the background file is generated, but only 1 Em bin and 2 Phi bins are considered. It is already mentioned, but can you add more information? Especially, how the values in other bins are filled (I suppose they are the same as the actually events), and how long it takes to run through all of bins roughly.

[ckarwin]

@hiyoneda Ok, it's been updated.

1. Yes, the correct file is SMEXv12.Continuum.HEALPixO3_10bins_log_flat.binnedimaging.imagingresponse.nonsparse_nside8.area.good_chunks_unzip.earthocc.zip
2. I added the following text:

Note that the current code has not yet been optimized for speed, as it uses a simple nested for loop. The time required to generate the estimated background using all bins is roughyly 4 hours. The option to use a subset of the Em bins and/or Phi bins may be useful for analyses that also use a given subset, but at this point the main motivation for this option is for demonstrational purposes, and when using this option, nothing is done with the other bins.

[hiyoneda]

Thanks! It can be merged.