Spectral Analysis of Extended Source in Galactic Coordinates

[ckarwin]

Building from the tools and examples from @israelmcmc, @eneights, and @fieldrog, I was able to successfully fit the diffuse 511 spectrum in the Galactic frame (using a simplified toy model). The example notebook is still a development version, but it should serve as a proof of concept. Some notes below:

Some of the things in the notebook still need to be moved to the cosipy source code.

In order to get the expected counts for an extended source, I am rotating the response to Galactic coordinates on the fly. In general, this can take some time depending on how many pixels are included in the calculation. It will be better to use the precomputed response in Galactic coordinates. I will be implementing this next.

I setup COSILike so that it should be able to take models consisting of multiple extended sources. I think it should be pretty straight forward to do this for the point sources as well, and to also allow for models to include both point sources and extended sources. Before doing this though, I think the point source analysis needs to be implemented in the Galactic frame, which I believe @eneights is already working on.

[israelmcmc]

Nice! I skimmed through the code changes, it seems that the only change needed to take both extended and point sources is changing = to += here https://github.com/cositools/cosipy/pull/114/files#diff-9c8499c20719aa0f7755d960aec888f802742b0caeee03533d6948978bfae32dR164

[ckarwin]

@israelmcmc Yes, that's right, essentially you should just need to change = to +=, and then of course make it compatible with the rest of the code. It still needs to be tested though. Some details:

I think we first need the point source analysis done in the Galactic frame. Am I correct that @eneights is already working on this?

The expectation for the point source seems to return a sparse histogram, whereas for the extended source it returns a dense histogram. This is any easy fix, but they need to be compatible before summing them.

The point source analysis in COSILike is setup to also allow for the source position to be optimized, but just for a single source. This will need to be modified as well -- at least eventually, but I'm not sure if we're really going to use that option for now. Tthat should also be a pretty easy change.

[eneights]

@ckarwin Yes, I’m working on point source spectral fitting in galactic coordinates. I’ll make a PR soon so you can see the changes in COSILike

[ckarwin]

Excellent, thanks @eneights!

[eneights]

Hey @ckarwin, I'm still working through the notebook, but just a couple of small changes while I'm waiting for it to run:

Add an option for the user to define a path to where the data is stored, so the files don't need to be in the same directory as the notebook Be more clear in the introduction about where the files are stored on wasabi. I had some trouble finding them Maybe add a cell to download the binned data from wasabi?

[ckarwin]

Thanks, @eneights. I'll get this implemented in the cleaned up notebook.

[eneights]

I ran through the notebook, and it works well! I just have a few minor comments:

The initial value of the background parameter is set to 1e-4, while the true value is ~1. I guess this came from the GRB spectral fit notebook, whereI had set it to 1e-4 since it's using the full 2 hours of background for the background model, but only 2 s of data. This isn’t causing any problems and it’s fitting correctly, so it doesn’t really matter in this case, but it’s something to keep in mind, especially when fitting more parameters later on.

You currently can’t run method 2 without first running some cells in method 1 (to define grid/grid_response).

Combining the unbinned files and binning the data is still using a lot of memory (up to ~150 GB), which @ckarwin is already working on fixing

[ckarwin]

Thanks, @eneights. I'm working on getting things updated to be compatible with your last PR.

[hiyoneda]

Thanks, Chris. I am running the notebook and posting comments. Sorry if some comments were already pointed out. In addition to the above 2 comments, I've got error when calculating the scatt map. I needed to modify get_scatt_map in SpacecraftFile.py at line390: from

x,y,z = attitudes[:-1].as_axes()

to

x,y,z = attitudes.as_axes()
x,y,z = x[:-1], y[:-1], z[:-1]

[ckarwin]

Hi @hiyoneda, I'll look into this, but @eneights and I both tested the notebook and didn't have this error. Are you using the most current version of cosipy in main?

[hiyoneda]

Yes, I think so.

[ckarwin]

@hiyoneda I'm in the middle of updating my PR to be compatible with Eliza's latest changes. Can you please wait till I get things updated, in the next day or two, before testing?

[eneights]

@hiyoneda Do you get the error when running scatt_map = ori.get_scatt_map(nside = response.nside * 2, coordsys = 'galactic')? I just ran that cell again and it works for me without having to change anything

[hiyoneda]

Oh, sorry. scatt_map = ori.get_scatt_map(inside = response.inside * 2, coordsys = 'galactic') works well without any changes now. The error was because I was using the main branch of scoords. After I reinstalled it using the develop branch, it worked.

@ckarwin I understood. I will wait until you finish updating.

[ckarwin]

@eneights @fieldrog @hiyoneda @israelmcmc I just updated my PR. The extended source fit should now be compatible with the point source fit (for galactic coordinates of course). I was able to successfully fit a model with two components. In principle, you should be able to fit the DC2 models following a similar approach. If you have a chance to test my example notebook that would be great! Or if you want to try to complete the DC2 challenge that would also be great! I also plan to try to do the DC2 challenge next.

Note: the only thing I didn't test directly are the cells for copying the data from wasabi. Hopefully it should be fine, but let me know if not.

[hiyoneda]

I could run through the notebook without errors.

[israelmcmc]

I marked it as a draft for now, just for tracking purposes. Feel free to remove the flag whenever it is ready.

[ckarwin]

I performed another fit variation with the notebook, where I perform the likelihood fit using the entire sky. In the fit I use both the extended source and the point source, and keep free just the normalizations, as well as the BG normalization. Here is a summary of the results:

computation time: 235.2 minutes

log(likelihood) = 1.527559e+07

best-fit norm (extended source): (4.6951 +/- 0.0025) x 10^-2 ph/cm2/s injected norm (extended source): 4.0 x 10^-2 ph/cm2/s

best-fit norm (point source): (0.0 +/- 2.0) x 10^-9 ph/cm2/s injected norm (point source): 1.0 x 10^-2 ph/cm2/s note: I have set the lower bound to be 0.

BG norm: (9.32 +/- 0.05) x 10^-1 injected: 1.0

The likelihood is significantly better than the other fit variation where I only use a subset of the pixels, which makes sense.

The point source goes to zero. I think that this actually is not very surprising. There is a high degree of degeneracy between the extended source and the point source. They both have the same spectral model, although the point source flux is 4x weaker. I think that in this case the point source is just not significant enough to be detected.

The normalization of the extended source is close to the injected value, although slightly higher. Correspondingly, the BG normalization is close to the expected value of 1.0, although slightly lower.

[hiyoneda]

I've also tried fitting the data using the entire sky. Since I was interested in whether we estimate a spatial model parameter, I used only the extended model by making gaussian.Gaussian_on_sphere.sigma free. The initial value was set to 10 degrees so as it is different from the true value (5 deg.) For the spectrum model, only F is free, and others are fixed.

The fitted spectral and background normalizations seem well. For the sigma of the 2d Gaussian model, it was derived as 3.786 +/- 0.007 deg. I think that the spatial model also fits well. It is smaller than the true value, but it may be acceptable because the simulation data also includes a point source at GC, which makes sigma small effectively when fitting the data with only the single component. Rather, I am curious why the error is so small.

In this test, I modified COSILike.py to make the calculation faster. I implemented a similar code that I am using in the image deconvolution module. As a result, fitting the entire sky took only 1 minute 11 seconds with my m1mac 64 GB. The essential parts are attached at the end of this comment (but sharing the code directly may be better...).

As Israel suggested in the meeting, in the near future, I think that it would be a good idea to prepare a kind of new class "ExtendedSourceResponse". The expected count calculation is better to be implemented as a part of the diffuse response handling class. Then, this computationally heavy part can be separated from the actual fitting class. (If the time still allows, I would like to work on that.)

# load the precomputed response in __init__
    def __init__(self, name, dr, data, bkg, sc_orientation,
...
...
...
    self.image_response = Histogram.open(self.precomputed_psr_file)

    def set_model(self, model):
...
        for name,source in extended_sources.items():

            # Set spectrum:
            # Note: the spectral parameters are updated internally by 3ML
            # during the likelihood scan.
            spectrum = source.spectrum.main.shape

            # just for Gaussian
            spectrum_unit = spectrum.F.unit / spectrum.sigma.unit

            eaxis = self.image_response.axes['Ei']

            flux = Quantity([integrate.quad(spectrum, lo_lim/lo_lim.unit, hi_lim/hi_lim.unit)[0] * spectrum_unit * lo_lim.unit for lo_lim,hi_lim in zip(eaxis.lower_bounds, eaxis.upper_bounds)])

            modelmap = ModelMap(nside = self.image_response.axes['NuLambda'].nside, energy_edges = eaxis.edges)
            modelmap.axes[0].label = 'NuLambda'

            npix = modelmap.axes[0].npix
            coords = modelmap.axes[0].pix2skycoord(np.arange(npix))
            normalized_map = source.spatial_shape(coords.l.deg, coords.b.deg) / u.sr

            modelmap[:] = np.tensordot(normalized_map, flux, axes = 0)

            # Get expectation using precomputed psr in Galactic coordinates:
            total_expectation = Histogram(edges = self.image_response.axes[2:])

            total_expectation += np.tensordot(modelmap.contents, self.image_response.contents, axes = ([0,1], [0,1])) * 4 * np.pi / npix * u.sr
...

[israelmcmc]

Hey @hiyoneda, it's very nice that you were able to leave the gaussian width free and it's working! Also the speed improvement is pretty cool.

Rather, I am curious why the error is so small. I'm not super surprised since I think the signal is quite strong. This causes the TS profile to be high and peaked. The key thing there is to remember that this is only statistical uncertainty, and doesn't include a systematic, which is likely large due to the coarse detector response. My guess is that if we run a goodness of fit or get the residual they would show up that the "real" error (stats+sys) should be larger.

The essential parts are attached at the end of this comment (but sharing the code directly may be better...). Another option is for you to submit a quick PR against @ckarwin's ckarwin:main branch, and once merge this PR will be updated automatically. Anyway is fine with me of course, whatever you two find it easier.

(If the time still allows, I would like to work on that ["ExtendedSourceResponse"].) 👍

[ckarwin]

This is great, @hiyoneda. I'm relieved to see that things are running much faster. Indeed, I figured the code could be much faster without looping over the entire sky each time. Some questions/comments:

1. Where is ModelMap coming from?
2. Can you're updated code fit multiple sources at once, including point sources and extended sources? If so, are you able to replicate the same fit that I did in my last comment, i.e. free just the normalizations of the two components, using the entire sky? It would be a good sanity check to make sure that we get the same results.
3. How do you want to proceed for merging our updates to get a working version for DC2? Based on what you sent, I think I can implement your changes in my PR fairly easily. The suggestion from @israelmcmc is also a good one. Or if you want to update things on your end that will also work. Keep in mind though that the code should be able to fit multiple sources at once (this is actually pretty straight forward). Please let me know what you think.

[hiyoneda]

I've submitted the PR to @ckarwin 's main branch just now. Cna you check it?

Where is ModelMap coming from?

This class was originally developed for the image deconvolution. It is now in cosipy.image_deconvolution.

Can you're updated code fit multiple sources at once, including point sources and extended sources?

I've checked the fitting test with one extended model and a point source as tested in the last part of Chris's notebook. The attached figure shows the results. While there is a slight difference, it is basically consistent. I haven't tested multiple point sources, which should be our next action item.

How do you want to proceed for merging our updates to get a working version for DC2?

First, please check my PR to @ckarwin 's main branch. Then, we can go forward to test multiple-point-source fitting. I feel like that it is better to modify the code to improve the calculation time, because it seems like a rotated response for each source is recalculated every time its expected counts are calculated. Also, we will merge @fieldrog 's 2d gaussian model and 511 keV models. I expect that this should be enough to do before the release of DC2. In the middle-term plan, we need to implement ExtendedSourceResponse to make things separated and modularized. Initially, I wanted to do this quickly, but I began to notice that we need to do more, e.g., determining the response format more consistently, etc. So maybe I rather will summarize task lists before that, including posting issues, and then start to work on that after DC2.

[ckarwin]

Excellent, thanks @hiyoneda! I'll get to this first thing next week. I'm happy to see we get essentially matching results for the extended src + point source fit over the entire sky.

[ckarwin]

Hi @hiyoneda, I checked your PR (against my main branch) and things are working well. Very nice work! I can merge it soon, and I will update my example notebook. Then, I think @fieldrog will just need to change the input models. First though, I have just a couple suggestions:

1. By default you always initialize the point source response and make a corresponding print statement: https://github.com/hiyoneda/cosipy/blob/2c6f42f5d4db8cf503aecfaee41ba35bd3dbcaf5/cosipy/threeml/COSILike.py#L181 But this isn't needed if the model doesn't have any point sources. The code is actually fine as it it, I would just change the print statement so that it only prints if there are actually point sources in the model.

2. It would be very convenient to have access to the total_expectation for each source (i.e. the histogram object), after the likelihood scan finishes. This would allow for easily making plots of the expected counts. Note that one can already access the total signal after the scan with cosi._signal, but this is a numpy array. I wouldn't change anything in the code, maybe just define an object of the class that gives you the expected counts for each source. Or maybe you know a better way?

[hiyoneda]

Thanks @ckarwin. Sounds good. I agree with your suggestions. As you said, defining an object to store the expected counts for each source seems the best way for the second point. Also, probably, it is better to redefine self._signal, expectation as histpy.Histogram s to make it consistent (but it might cause some conflict with other classes?)

[ckarwin]

@hiyoneda I think it is best to leave self._signal the way it is for now. The reason is that for point sources the expectation is a sparse histogram, and for extended sources it is a dense histogram. If your model contains both a point source and an extended source, then you'll get an error if you try to add the two matrices. I agree that it might be better to eventually make things consistent, but maybe it's better to implement this in the response class. If we have the expected counts for each source, the total signal (histogram object) can easily be obtained by summing all the sources.

[ckarwin]

I merged @hiyoneda PR into my main, and also implemented the two comments I mentioned above (https://github.com/cositools/cosipy/pull/114#issuecomment-1864826445). Additionally, I updated the extended source example notebook. I also tested Eliza's two notebooks (SpectralFit_Crab.ipynb and SpectralFit_Line.ipynb), and verified that they are still working ok after all the changes. I think this PR is ready to be merged into main, after being checked by @eneights, who has agreed to review it by Jan 10th. If anybody else can also check it, that would be great.

The only thing remaining is to add a specific 511 model from DC2. My suggestion is to add this at the end of the current example notebook, i.e. don't change anything in the first two examples, but just add a third example at the bottom.