Create spanning templates for DC2

[jbkalmbach]

The goal here is to produce a spanning template set using the 30 band SEDs of DC2. @sschmidt23 has provided a set of the Brown et al. templates with a linear interpolation over the strongest emission line regions. I will use these to create a template set that better matches the DC2 color space.

[sschmidt23]

After turning in my review, I forgot something, @jbkalmbach do the Magtrue[band]_lsstz0 colors include internal reddenning? i.e. are they the equivalent of the "magnitude:LSST[band]:observed:dustAtlas" magnitudes? If the z=0 colors do not include internal reddening then this might explain why the colors are still so blue.

Also, I was thinking about SEDs that we could add to the base 129 Brown SEDs and remembered that I have two single stellar population models with low/no dust and low metallicity that were used with BPZ for the HUDF to mimic super blue galaxies at z~5-7. I interpolated those two templates on to the same grid as the Brown SEDs and re-ran Bryce's notebook, and the blue end color matching is much better (see attached plot).

I can add these two templates to the data/brown_templates/ directory on the issue/2/ branch for people to play around with as well, if others are curious

[sschmidt23]

Here's a plot of the two new templates in blue compared to the two bluest Brown SEDs in red/orange, normalized to have the same flux at 8000 angstroms. The steeper slope of the ssp5Myr_z008 template is bluer than even most of the DC2 data, it makes sense that this would help the PCA capture blue galaxies. And the ssp25Myr_z008 template looks similar in blue slope to the bluest Brown SED, but maybe the slight diversity in the SED blueward of the Lyman limit helps as well?

[sschmidt23]

I created another branch u/sschmidt23/add_blue_templates (note: edited branch name for clarity) so that these aren't included by default. The two new templates are in data/brown_templates, and Bryce's notebook should run fine with no changes in this branch.

[jbkalmbach]

From the gcr-catalogs reader for the cosmoDC2 catalogs it looks like the Mag_true_{}_lsst_z0 bands are equivalent to LSST_filters/magnitude:LSST_{}:rest:dustAtlas.

[evevkovacs]

That's correct. Mag_true* are supposed to be the un-lensed rest-frame magnitudes, including host extinction. Please let me know if you have any questions.

[jbkalmbach]

Ok, @sschmidt23 @morriscb I have updated the notebook based upon the comments from the PR. I think I incorporated all of them but let me know if there's still something you guys would like to see or that I overlooked.

In terms of final results here are the results from running just with LSST colors and then with 26 of the tophat bands. The plots show the rest frame LSST colors of the new templates created from each starting color space. It seems to me that they are pretty close, but the ones that ran the Gaussian Process in the LSST space might be slightly preferable based upon their performance covering the bluer colors. Let me know what you all think and we can move the preferred set of new templates into a common folder for others.

w/ Tophat bands:

w/ LSST colors:

[sschmidt23]

I finally got back to looking in to this problem, and I think I figured out at least one of the problems: it looks like the DC2 galaxies have strange dust extinction for high values of Av. I added a notebook showing an example here: https://github.com/LSSTDESC/pz_pdf/blob/issue/2/templates_for_DC2/doc/notebooks/Test_ext_vs_unext_vs_dustmodels.ipynb For high values of Av, the dust extinction does not follow a smooth Cardelli/Fitzpatrick/Calzetti-like curve in terms of how much the flux is extincted, rather there is a plateau at ~8000-12000 angstroms and a small jump at ~7500 angstroms that induces the too red in i-z and too blue in z-y effect that we're seeing in the figure above.
Here is an example SED showing the 30 top-hat bands both unextincted (green) and extincted (black)

The ratio of unreddened to reddened flux is also off in a few other wavelengths (notably the blue end), which may be causing some of the other color discrepancies that we're seeing between DC2 colors and typical galaxies such as the Brown SEDs or CWW models. Here is a plot of the ratio of the unextincted to extincted flux compared to a couple common dust laws (normalized to 1 at reddest bin with an eyeballed Av value).
You can see that it does not look much like a standard dust law. Other galaxies with high Av had similar (but not exactly the same) patterns for the ratio.

We can talk about maybe trying to figure out why Galacticus is doing this, but if we want to match these colors, I think we can back out a version of this wonky dust law and apply it to some of the Brown SEDs, which should at least allow us to get some templates to feed into ESP to generate a spanning set.

[janewman-pitt-edu]

It would be good to communicate the issue to I guess Eve and Andrew H.? I agree that if we know the dust law we potentially could apply it.

[evevkovacs]

@sschmidt23 @janewman-pitt-edu Which version of cosmoDC2 are you using?

[sschmidt23]

Yep, I was just getting ready to compose an email to Andrew and Eve before going into too much more detail here.
@evevkovacs on Friday I switched to "cosmoDC2_v1.1.4_small", which may have eliminated some even stranger looking cases that I saw when we were looking at the earlier version of the catalog "cosmoDC2_v1.0_9556". Here is one such example from the v1.0 catalog where there were a few galaxies with a big chunk of flux removed by dust at ~7000 angstroms that had even stranger colors as a result.

But the other strange looking examples are still present in the new 1.1.4 catalogs

[sschmidt23]

Actually, this issue is probably a better place to discuss than email so we have all the discussion in one place. Much of this is just me poking at the data, so comments from others would be super helpful and quite welcome.

I have comments sprinkled into the notebook mentioned above (and pushed to the issue/2/templates_for_DC2 branch): https://github.com/LSSTDESC/pz_pdf/blob/issue/2/templates_for_DC2/doc/notebooks/Test_ext_vs_unext_vs_dustmodels.ipynb so people can look through that for more background info.

In addition to @aphearin I'll mention @abensonca here to see if he can comment on whether this behavior looks familiar for Galacticus at high Av values. Note that I haven't confirmed that this happens for all high Av galaxies, I am selecting from a region of rest frame z-y vs i-z that looks quite different from what we expect from typical galaxies, so it may only happen some of the time.

Again, I think we can pretty easily make some dust models that fit these galaxies so we can generate spanning templates. One of my concerns, though, is that having a feature at ~7500 angstroms in some galaxies may mean that you can now distinguish some low-z galaxies unambiguously that otherwise might have been degenerate colors with some high redshift templates, so our photo-z estimates might be a bit optimistic for these galaxies.

And as I mentioned above, I've been focused on z-y vs i-z rest frame, the 1000-2500 Angstrom portion of the extinction curve also looks very odd for these galaxies, which is most likely resulting in some odd colors for the blue end of some galaxies (at redshifts above 0.5 or so when these wavelengths shift into the LSST bands).

[abensonca]

I suspect that this behavior is related to the limited scope of the dust models that we used here. Unlike typical dust modeling in which an A_V is determined and then a single extinction curve is applied to each model galaxy, we've used a library of 3D radiative transfer models to compute the attenuation based on the physical properties of each galaxy. While in principle a better approach (i.e. should correctly capture the effects of inclination, morphology, etc.) the model we used had some significant limitations.

Primarily it was based on some old calculations (from 1999) which meant that the grid of models was very coarsely sampled - requiring us to do a lot of interpolation and extrapolation to get the attenuation for each galaxy at each wavelength.

This probably explains the "plateau" feature you're seeing: The table of dust models had tabulation points in wavelength at 8059Å and 12369Å for example. This probably causes the plateau over that wavelength range.

Additionally, while building an improved model (see below) I discovered some clearly incorrect values in the tabulated dust attenuations. These are only a small number of points (and i haven't been able to track back through 20 years of history of the tabulated data to figure out if these were bugs in the original models, or just transcription errors when somebody built the files) but could possibly explain some of the sharp features seen in the curves.

I recently built a new dust tabulation (see here) which is: built using a modern radiative transfer code, tabulates attenuation over a much wider range of wavelengths and galactic properties, and tabulates on a much finer grid in all of these properties. These new models should resolve the problems you're finding with the current implementation.

It would be useful if we could identify one or two of the problem galaxies in the original Galacticus outputs - I'd like to take these and confirm that I can reproduce the bad dust behavior with the original dust models, and then see how the new dust model performs in these cases.

[sschmidt23]

Ok, thanks @abensonca , a grid of models with some tabulated points is exactly what I suspected based on what I'm seeing. As I said, I think we can pretty easily factor this into our SED models for photo-z's, it was just puzzling since we couldn't reproduce these colors with the standard dust laws, glad we now understand the cause.

I pushed an update to the notebook and added 'galaxy_id' and printed a few out so you can query for them in the GCRCatalog for cosmoDC2_v1.1.4_small catalog. In the notebook I was selecting from a very narrow range of rest frame colors, so most of the galaxies are one of two types, one blue galaxy and one much redder. The blue ones have a derived A_v >~2 while the redder galaxies have a derived A_v more similar to 1.5-1.6. As I said, I think you can probably find examples just by selecting on high values of derived A_v in the catalog. I would just do a filters = ['A_v' > 1.5] or something in the get_quantities GCR method.

You can see the effects of the discretization of the dust law in many of the galaxies at all colors, I think the effects are just much more noticeable in the high A_v ones because the flux is so much more attenuated, and since many galaxies have very little dust, it's not a problem there. Though, we may want to look at maybe using more points to interpolate the dust law in the future.

Actually, more striking than the plateau is the blue end differences from typical Cardell/Calzetti-like curves. Is there any evidence for such differences in observed galaxy samples at higher extinctions, higher redshifts, or other areas of parameter space? In photo-z algorithms we usually just pick one or two of the dust laws, compute a quick grid and call it a day. If we have a few other options that could actually be helpful in getting better redshift estimates.

[abensonca]

Thanks - I'll pull some example galaxies from cosmoDC2_v1.1.4_small and run some tests with the new dust grids.

I don't know if the blue end behavior is realistic for high redshift, high extinction galaxies. If we could figure out if there are any samples that could answer that question it would be very useful - I'm always concerned that the lack of constraints on the dust modeling means that it could easily have unrealistic behavior.

[janewman-pitt-edu]

@abensonca : are you looking for something different from Naveen Reddy's work (e.g. https://arxiv.org/abs/1705.09302 )?

[abensonca]

No, that looks like a great place to start. I'll takea look at this.

[janewman-pitt-edu]

I think he's done other work on this too but that was the first from him to pop up on google. Naveen is almost local for you of course.

[abensonca]

Yes, and I'm actually meeting with him tomorrow so will get his input on this.

[abensonca]

I've been unable to log in to jupyter-dev so I'm not yet able to pull properties for the galaxies in your list above. But, just to show an example of a comparison between the old and new dust models, here's a case for a galaxy with high extinction. The black line is the old dust model, the magenta line is the new dust model. You can see that the discreteness of the tabulation of attenuations in the old model (black curve) very clearly. But in the new model (magenta curve) the attenuation is smoothly varying as a function of wavelength. (The dust properties differ between the two models, so the attenuations differ.)

[sschmidt23]

I've been working on this a bit more, by pulling out two of the dust laws as shown above and creating some extra templates to feed into Bryce's notebook to give the Gaussian Process a better starting place, I'm able to generate much better (but still not perfect) templates. The notebook showing this has been added to a new branch: u/sschmidt23/galacticus_dust, and the notebook computing some templates is here: https://github.com/LSSTDESC/pz_pdf/blob/u/sschmidt23/galacticus_dust/doc/notebooks/extrapolated_dc2_templates_galacticus_SEDs.ipynb

This plot from the notebook shows the much better coverage in the rest frame colors, a few areas not covered, but much better than previous efforts:

I output 100 k-means SEDs from this set to test it as a template set with BPZ, results are ok but not great with sigma_z/1+z~0.045 for a simple point-like z_peak estimate (note: I did not add uncertainties to the input photometry, and set the supposed error to a very small number). Plotting the observed frame colors of the 100 templates it looks like we are missing the colors for on the outskirts of the color blobs, so 100 k-means SEDs may not be enough to summarize this set.

I am re-running the notebook now with a few tweaks and will save 200 k-means SEDs, hopefully that will be a little better in terms of how well it spans.

[janewman-pitt-edu]

There are methods which can select points at the edges of distributions -- unfortunately I don't remember exactly off-hand (might be methods that make the points repel each other and evolve the distribution to stasis, or that choose the point most distant from any current point). Adding some points selected that way might be helpful.

[janewman-pitt-edu]

Maybe could use scipy.spatial.convexhull ? or concave hull (https://github.com/sebastianbeyer/concavehull )? or https://github.com/tomilov/quickhull ?

[sschmidt23]

I tried something like that for DC1 and it over-emphasized the outliers. I was looking through the notebook and noticed that there was a cut down to the the 5000 "best fit" Gauss Process SEDs. It turns out that a lot of those were in the central parts of the color distribution. I switched that out to 50,000 and it looked better in terms of covering the rest frame colors, though 200 still looked better than 100. Here's the rest frame colors for the 200 k-means templates for 50k "best fit" SEDs. I will run BPZ on these tomorrow and see if the results look better, I expect they will.

[jbkalmbach]

I looked through the notebook. It looks like you haven't committed the latest version since it was still only using the 5000 best models and 100 templates instead of 50000/200 as you posted but I'm assuming everything else is the same. My main question was on the templates you used in the training. It says there are ~1400 of them. Did you just generate additional versions of the existing templates with new dust models and include them all, is that where they all came from?

Anyway, everything looks really good!

[sschmidt23]

Yes, the 1400 templates include the 129 Brown SEDs, plus the Brown SEDs with several values of the extinction for the two empirical dust models that I pulled out of the DC2 data (5 Av values x 2 extinctions x 129 templates = 1290 + original 129 is where the 1400 comes from). As you can see from the blue points in the plot above, many of those are overkill, and too red compared to the DC2 data, but the GP properly ignores those when coming up with the final data set.

I have not committed the new notebook, which does simply switch the "best" models to 50,000 and k=200 for the k-means. I just did a run with BPZ and the sigma_z, bias, and outlier rate do not improve much, looking at the observed colors, they did not improve a great deal. While they were better, they still missed a couple of areas of observed color space.

I instead did what it mentioned at the top of your original notebook: used the original 100 k-means SEDs to supplement the 129 Brown SEDs, and results did improve. With that total 229 SED set the sigma_z is 0.041, bias is -0.0097 and outler rate is 5.4%, and I expect the outlier rate might actually improve with realistic magnitude errors, as right now with tiny errors in my quick runs the chi^2 values are all more extreme, and thus the prior does not have much of an effect because the likelihoods are all very peaky.

Here is a plot with both the k-means colors (black) and the Brown colors (blue).

I think we can call this good enough for our first round of BPZ based photo-z's, I'll train up a prior with this template set and start working on getting a photo-z for the whole cosmoDC2v1.1.4 data and run 1.2i data sets for people to play with.

[janewman-pitt-edu]

Are the red points with errors or without? If without, we clearly don't have a spanning/appropriate template set.
Possibly we should be just making the templates from the same spectral synthesis galacticus uses? Or perhaps if we choose ~1000 representative objects Eve could send us more detailed SEDs for them?

[sschmidt23]

The red points are without errors, so the template set is not completely spanning. I wasn't too concerned about this, actually I was rather glad that the PCA/GP didn't completely cover the input colors, since DC2 is supposed to have some level of incompleteness. I figured this incompleteness, while not realistic, would be a proxy for our eventual realistic incompleteness model, while a sigma_z of 0.04 was not too scary for placeholder p(z) for people to have access to through the GCR initially.

Are you thinking that we want a complete set of SEDs that cover the cosmoDC2 dataset before we then work on removing certain portions of the training data as an incompleteness model? We can try adding some of the SEDs from either/both the SPHEREx set that we received or ask Eve et al for some representative SEDs to see if that fills in the remaining uncovered areas of parameter space, but I was worried that it might take a while to get that settled.

What do you think, is this "good enough" for an initial photo-z estimate, or should we go for a spanning template set? Here is a z_peak vs spec-z plot for 100k cosmodc2 galaxies with at i<25.3 for the 229 templates

[janewman-pitt-edu]

I think this is mostly OK except for z~0.2/2.8 degeneracy which seems too strong. Is there a prior here? Might be helpful to use a low alpha for visualization. Sorry for the slow response, I've been tied up with DOE proposal.

[janewman-pitt-edu]

The completely falling apart at ~2.5 also doesn't seem right. Is there any Madau-type reddening model involved?

[sschmidt23]

I am almost positive that there is mo Madau reddening, both because I am fairly sure I asked about this when we started work on the photo-z's, and also because the u-band colors don't start getting redder as redshift increases, i.e. there is no wholesale reddening of the u-g colors at z>2 that we would expect with Madau extinction.

In the pz vs sz plot above, I was running with very tiny errors (0.005 in all bands since the data did not come with errors), even with these tiny errors there is still the large degeneracy in colors between low and high redshift, as a lot of the galaxies are very very blue and basically look like a power law. With tiny error bars the chi^2 are high and the likelihoods very peaky for faint, high redshift galaxies. To make this more realistic, I simulated 10-year depth uncertainties and plugged those into BPZ (I didn't actually scatter the magnitudes, I kept the true mags but told BPZ the expected errors). This should let a prior do a better job at influencing the redshift for faint galaxies, which should help for DC2, since there is a somewhat unrealistic lack of faint blue galaxies at low redshifts.

But, things still didn't look great, so I added a small grid of Bruzual & Charlot SEDs to the GP/PCA notebook and re-ran, which generated a template set that spans a little better.

I retrained the prior and ran with and without the prior on the data where the 10-year depth errors were fed in, results are better but still not perfect, for point estimates the sample statistics are: with prior: sigma_z/(1+z) = 0.033 bias = 0.002 cat. outlier rate= 5.02% (>3sigma) without prior: sigma_z/(1+z) = 0.036 bias = 0.0005 cat. outlier rate = 6.66% (since sigma is bigger the drop is more than the 1.6% difference)

Point estimate plots are above, no prior applied on left, with prior on the right.