Add Cosmoglobe DR2 model

[angela-bonato]

In the latest version of the Cosmoglobe Data Release 2, the zodiacal light model uses the parametrization of Kelsall et al. (1998) with the phase function from Hong et al. (1985). For this reason, I added a new model in zodipy called "cosmoglobe" which features the parametrization and the best-fit parameters derived from the DR2 analysis. I also set the parameters of the existing zodipy "dirbe" model to those listed in Kelsall et al (1998).
I was able to successfully run the new models and produce plots and maps with them, but I have not yet updated the test functions or the documentation. This still needs to be done.

[MetinSa]

Happy to see that the cosmoglobe model is finally being added here! Well done 😊
I left a comment regarding the difference between the DIRBE values seen here and in the paper.

You should also run the ruff linter on the project and the tests (See the readme for how to do that).

Additionally, you should update the docs/ and include the cosmoglobe model as one of the available models in zodipy.

And finally, once the PR is merged we need to publish a new version of cosmoglobe to PIP.

You should request a review from me again when the tests run green and your are happy with the PR. Feel free to ask if there is anything you need 😄

[MetinSa]

Im aware that these values slightly differ from the values presented in the K98 paper. I found these values by reading the DIRBE source code which was sent to us at some point and decided to use these more spesific values instead. But given the uncertianties in the model I dont think it matters much. But if using these "rounded" values significant change the end results we should consider adding a new DIRBE model instead with these values so for backward compatibility.

[MetinSa]

@angela-bonato any updates here? 👀

[angela-bonato]

Thank you for the inputs @MetinSa and I'm sorry for not answering before. I'm on it now and I think in a couple of days I will be able to fix it and let you revise it.

[codecov]

Codecov Report

❌ Patch coverage is 98.52941% with 1 line in your changes missing coverage. Please review.
✅ Project coverage is 99.48%. Comparing base (f8e6bad) to head (1f1fd55).
⚠️ Report is 16 commits behind head on main.

Files with missing lines	Patch %	Lines
zodipy/unpack_model.py	95.23%	1 Missing ⚠️

Additional details and impacted files

@@             Coverage Diff             @@
##              main      #56      +/-   ##
===========================================
- Coverage   100.00%   99.48%   -0.52%     
===========================================
  Files           14       14              
  Lines          708      782      +74     
===========================================
+ Hits           708      778      +70     
- Misses           0        4       +4     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[angela-bonato]

I ran the tests listed in the README and fixed those. Then I changed the documentation to add in the list of available models the new one, I don't think there is something else that needs to be updated.

For the DIRBE model I changed the parameters values to match those in the main branch. The difference between using them or the values listed in the original K98 paper is about 0.5% so I think we can just keep what goes with the current tests.

There are coverage issues that I am not sure how to fix. My understanding is they might be caused by the fact that maybe additional tests specifically written for the new model are expected, or maybe I just got something wrong when using the DIRBE model as a basis to write the Cosmoglobe model, but in any case I need help with those.

[MetinSa]

The coverage should probably not be a blocker so I will disable it when I get to work. Looks good otherwise! Have you tested and seen if it reproduces commander over the dirbe bandpasses for instance?

[angela-bonato]

No, I haven't. Do you mean like taking the difference between the zodi moldel as fitted by Commander and the one simulated by Zodipy in all the bands or is there something more sophisticated I could do?

[MetinSa]

No, I haven't. Do you mean like taking the difference between the zodi moldel as fitted by Commander and the one simulated by Zodipy in all the bands or is there something more sophisticated I could do?

Yeah. Produce a zodi map in zodipy with the new cosmoglobe model using one of the dirbe bandpasses and compare it to the corresponding zodi map from the commander run (for instance band 6). If its "close enough" we know that the model in zodipy is doing what we want it to do.

[MetinSa]

Regarding the code coverage, we should add a new coverage.yml to the root of the repository. Its documented here https://docs.codecov.com/docs/commit-status#informational, but its contents should look something like this:

coverage:
  status:
    project:
      default:
        informational: true
    patch:
      default:
        informational: true

Could you try that and see if it stops failing the checks?

[angela-bonato]

No, I haven't. Do you mean like taking the difference between the zodi moldel as fitted by Commander and the one simulated by Zodipy in all the bands or is there something more sophisticated I could do?

Yeah. Produce a zodi map in zodipy with the new cosmoglobe model using one of the dirbe bandpasses and compare it to the corresponding zodi map from the commander run (for instance band 6). If its "close enough" we know that the model in zodipy is doing what we want it to do.

This is what I am getting for bands 5 and 6:

I am afraid this does not look good. Some of the structures are probably due to the difference between having to deal with simulated data and working with TODs, but to me it seems also zodi stuff is involved. I'll try to debug the code again and look also at the difference per-component to see if there's one particularly bad. I am talking about this with Duncan as well so I hope we will sort this out, if you have any suggestions please be my guest :)

[MetinSa]

How are you specifying the spacecraft coordinates when performing the zodipy sim? In commander i believe there are some optimizations where a cache is used to reuse certain zodi evaluations if the spacecraft has moved less than some fraction on the sky. Zodipy will perform the actual evaluation per specified sky coordinate.

If you could post the zodipy code that would be useful

[angela-bonato]

Ths is a script that Duncan gave me to produce mission averaged maps and I just adapted it to my needs.

comm_tod = TODLoader("/mn/stornext/d5/data/duncanwa/DIRBE/hdf_files", "")

NSIDE = 512
band_num = 6
scan_id  = 1

band = f'{band_num:02}'
scan = f'{scan_id:06}'

comm_tod.init_file(f"DIRBE_{band}_nside512_V22", "")
# pixel values in galactic coordinates
pix = comm_tod.load_field(f'{scan}/{band}_A/pix')
# boresight angle in radians
psi = comm_tod.load_field(f'{scan}/{band}_A/psi')
# flags used for processing. Bitmask is currently 26615
flags = comm_tod.load_field(f'{scan}/{band}_A/flag')
# calibrated instrument observations
cios = comm_tod.load_field(f'{scan}/{band}_A/tod')

# Time in MJD
t_i = comm_tod.load_field(f'{scan}/common/time')[0]
t_f = comm_tod.load_field(f'{scan}/common/time_end')[0]

model_c = zodipy.Model(25*u.micron, name='cosmoglobe')

FLAG0 = 26615

inds = np.bitwise_and(flags, FLAG0) == 0

b_map_c = np.zeros(hp.nside2npix(NSIDE))
hits_map = np.zeros(hp.nside2npix(NSIDE))

for scan_id in tqdm(np.arange(1, 285+1)):
    scan = f'{scan_id:06}'
    pix = comm_tod.load_field(f'{scan}/{band}_A/pix')
    flags = comm_tod.load_field(f'{scan}/{band}_A/flag')
    t_i = comm_tod.load_field(f'{scan}/common/time')[0]
    inds = np.bitwise_and(flags, FLAG0) == 0

    lon, lat = hp.pix2ang(NSIDE, pix,lonlat=True)
    skycoord = SkyCoord(lon*u.deg, lat*u.deg, obstime=Time(t_i, format='mjd'), frame="galactic")

    emission_c = model_c.evaluate(skycoord)
    
    for i in range(len(inds)):
        if inds[i]:
            b_map_c[pix[i]] += emission_c[i].value
            hits_map[pix[i]] += 1

zero_mask = (hits_map==0)
hits_map[zero_mask] = hp.UNSEEN

mhat_c = b_map_c/hits_map
mhat_c[zero_mask] = hp.UNSEEN
hp.fitsfunc.write_map("tot_mission_c_main.fits", mhat_c, dtype = np.float32, overwrite = True)

[MetinSa]

I suppose you don't have the zodi timestream from commander? If you do I would have tried to bin those in the same loop to validate that the binning isn't the issue. And maybe look at them in the time domain for comparison aswell.

Also, if you want this script to be a bit faster and have enough ram, you could reconstruct the full Lon lat array for the whole 285 days and do a single zodi evaluation. Alternatively do each scan like you do now but multiprocessed.

Zodipy also supports an obstime per pointing (the obstime can be an array of times) so you wouldn't have to pretend that the spacecraft sits still for a full day before jumping a full degree like you do here. That would give you a more accurate zodi map when comparing with Tod's (if that's relevant). I'm not sure how the obstime/obspos is solved in commander but it probably doesn't account for the full diff anyway.

[angela-bonato]

I don't really know the details of how the scanning strategy was handled in both cases, so I think I need more time to address that. However, I made some maps that maybe can help us better understand the sources of error.

First row, left column is the total K98 model which has been tested on Zodipy so I think this map should give us an idea about the error given by the simulated scanning strategy, next to it there is the total Cosmoglobe model and then the components we fit (cloud, band1, band2, band3). All maps are in band 6 (25 microns).

I think there is something going on with the zero level as well as the Ecliptic plane, so it is not a mere question of how mission averaged maps are simulated. I'll talk about this with the Cosmoglobe group and will write here as soon as I have news.