API changes

[drvdputt]

Here is an almost working version of the new API. Good enough for other people to start testing.

The old code is still there. In fact, the new Model class internally works with PAHFITBase and param_info.

I have also written a first test for this model. It makes some fit results, converts the changes back to the features table, and then converts the features table back to the fit model, and checks if the original and reconstructed model are the same.

Anyone testing this is encouraged to share their testing code, as inspiration for including more tests.

[drvdputt]

Attaching my notebook with some things I tested during development. Use it as a way to get started. dev-testing.ipynb.zip

[drvdputt]

Getting close in terms of functionality. See discussion at #231 . Next step is removing obsolete code from helpers.py.

[drvdputt]

initialize_model and fit_spectrum were removed from helpers. I have edited the run_pahfit and plot_pahfit command line scripts to use the new api. Doing so was also a good test for the save/load functionality. Here is an example of the new model, saved in ecsv format (with extension changed so github lets me upload it)

M101_Nucleus_irs_output.txt

Note the metadata, which currently stores the instrument name and redshift, so that the model can be loaded in again.

Edit: plot_pahfit has been confirmed to work too. Will now write a test for the save/load functionality.

[drvdputt]

Important detail: when features get trimmed away (== ignored during the setup of the astropy model), what should we do with the features table once we get the fit result with the reduced model?

Right now, any features that were not included in the model, just have their rows in the features table unchanged. Removing them would be an option too.

[drvdputt]

Writing down problems revealed by testing here:

Wavelength check

I have been getting the 1% warning, or even the 3% error, with MIRI data sometimes. @BethanyRS was also getting the warning with ISO speed 2 data. So it seems quite common that the data is a little bit outside the specification.

Redshift

I have a very rudimentary implementation for the redshift, which just shifts the data to the rest frame (so that data wavelength matches features table wavelength).

The problem with this is that the FWHM determination from the instrument specification, uses the feature wavelengths, and therefore the rest frame wavelength, instead of the observed wavelength. We should think about how redshift should interact with

• trimming the lines based on the wavelength coverage of the instrument
• determining the fwhm based on the instrument model
• plotting observational data + model

Storing fitted FWHM of lines in segment overlap region.

Right now, the lines that fall in instrument segment overlap regions are given variable widths. But these widths are not backported to the feature table, because they are considered an instrumental effect. Therefore, we can not inspect if that part of the fitting is working well (aside from parsing the resulting astropy model).

[BethanyRS]

Just for show, here is an example of fitting a spectrum with a range of wavelengths missing. To create this graph, I took some ISO data in the repository and deleted sections of it. The fit to the data still present seems reasonable.

[BethanyRS]

Wavelength check

I have been getting the 1% warning, or even the 3% error, with MIRI data sometimes. @BethanyRS was also getting the warning with ISO speed 2 data. So it seems quite comment that data is a little bit outside the specification.

To add to this, I was getting the 1% error for "Orion_D5_ISO-SWS_merged.ipac". I tried the other merged ISO-SWS files and they seem to all give the 3% error.

[jdtsmith]

Important detail: when features get trimmed away (== ignored during the setup of the astropy model), what should we do with the features table once we get the fit result with the reduced model?

Right now, any features that were not included in the model, just have their rows in the features table unchanged. Removing them would be an option too.

Hmm, good point. Perhaps we should simply mask their values in the masked array so no-one gets confused. Presumably the "guess" function will also do this.

[jdtsmith]

Writing down problems revealed by testing here:

Wavelength check

I have been getting the 1% warning, or even the 3% error, with MIRI data sometimes. @BethanyRS was also getting the warning with ISO speed 2 data. So it seems quite common that the data is a little bit outside the specification.

This probably means our wave-bounds are inaccurate in the instrument packs, which isn't terribly surprising. I based these on some tables I found on JDox I think, but it's possible those were pre-launch or otherwise out of date. If you or someone could ask around on the MIRI team for any updates to that, I'd appreciate. It is important to pass spectra in the observed frame, btw (I presume people are).

Redshift

I have a very rudimentary implementation for the redshift, which just shifts the data to the rest frame (so that data wavelength matches features table wavelength).

You must also multiply the flux by (1+z) to avoid power non-conservation.

The problem with this is that the FWHM determination from the instrument specification, uses the feature wavelengths, and therefore the rest frame wavelength, instead of the observed wavelength. We should think about how redshift should interact with

The way I do this on the new model side is as follows:

• trimming the lines based on the wavelength coverage of the instrument

I red-shift the features' line wavelengths on input as I check within_segment.

• determining the fwhm based on the instrument model

Same.

• plotting observational data + model

Haven't solved this, but I'd be OK always plotting in the rest-frame (especially during the fit), but I suppose it could be an option of plot.

We can also come back to redshift once the API PR is merged.

Storing fitted FWHM of lines in segment overlap region.

Right now, the lines that fall in instrument segment overlap regions are given variable widths. But these widths are not backported to the feature table, because they are considered an instrumental effect. Therefore, we can not inspect if that part of the fitting is working well (aside from parsing the resulting astropy model).

I think a good idea is to update the features table for any non-fixed line FWHM, regardless of whether that was passed in by the user, or implicitly set due to the presence of stitched spectra.

[jdtsmith]

Wavelength check

I have been getting the 1% warning, or even the 3% error, with MIRI data sometimes. @BethanyRS was also getting the warning with ISO speed 2 data. So it seems quite comment that data is a little bit outside the specification.

To add to this, I was getting the 1% error for "Orion_D5_ISO-SWS_merged.ipac". I tried the other merged ISO-SWS files and they seem to all give the 3% error.

Can you check with @jancami ? He made the ISO instrument packs. Again some ambiguity on "where a spectrum ends" is not uncommon.

[jdtsmith]

Just for show, here is an example of fitting a spectrum with a range of wavelengths missing. To create this graph, I took some ISO data in the repository and deleted sections of it. The fit to the data still present seems reasonable.

This example brings up a subtlety: for any dust features not well covered, they are pretty much unconstrained, but PAHFIT will go right ahead and fit them. Ideally the covariance matrix will show you just how uncertain they are, but we may need a heuristic to alert the user here.

[jdtsmith]

BTW, @drvdputt, I added features.util, which provides a few convenience functions for asking a parameter if it is fixed, missing and its min/max value. Thought you might like to use them in your Features plumbing so pushed that up.

[drvdputt]

BTW, @drvdputt, I added features.util, which provides a few convenience functions for asking a parameter if it is fixed, missing and its min/max value. Thought you might like to use them in your Features plumbing so pushed that up.

Rebase on master was trivial. Although I don't have immediate use for them, I do have a prototype model builder that skips the param_info step entirely (would remove MANY lines of code), where they might be useful. But I have not yet added those changes, because some things like the initial guessing also work with param_info, so those parts would have to be rewritten too.

I have now removed the save() function from PAHFITBase, because the files created are the old output format, and can't be read in anymore. It made use of the stuff in feature_strengths though, to add a few extra columns. So we might want to add some line strenght, equivalent width, and pah feature strength stuff to the new model class.

[jdtsmith]

I have now removed the save() function from PAHFITBase, because the files created are the old output format, and can't be read in anymore.

What did we decide about a model object "saving" itself? For it to store a Features table, or be a Features table?

It made use of the stuff in feature_strengths though, to add a few extra columns. So we might want to add some line strenght, equivalent width, and pah feature strength stuff to the new model class.

Yes, these should be methods not table columns. In particular, since covariance matrix will be in meta, to get the uncertainty of a given line, you'll have to consult that matrix. Will also provide a Features capability to co-add rows with proper covariance error propagation. I'd also like to move model and geometry out of the table itself and into model.meta['attenuation']. Since we control how Features tables print we can fix that up to make it look nicer. Only issue to consider is when saving the table, it's maybe inconvenient to have values and their uncertainties in different locations.

[drvdputt]

"Ignored" features are now masked out in the features table. The mask is applied in _parse_astropy_results, i.e. when the fit (or guess) results are being written to the table. The masking and unmasking happens via a new member function of Features, which uses some internal information to make sure only the relevant columns are masked out. The mask of the bounds are left untouched, as those have a special meaning (fixed / not fixed)

This mask is ignored during "guess" or "fit", as the features to be included might need to be reconsidered at that point.

Where it becomes useful is

1. When the user looks at the table contents by eye
2. For "plot", to make sure no irrelevant features are plotted
3. For the analysis functions, e.g. line strengths.

With this, the better our "feature trimming" algorithm works, the easier it will be to make good quality plots (no weird features like the plot shown in this thread by Bethany), and make sure that no badly fitted features make it into the analysis functions.

[drvdputt]

There's another subtlety with redshift:

If we shift the data to the rest frame, then all the features in that data will get narrower (wavelength points get closer together!). We should therefore adjust the FWHM in the model, so that it better represents the shifted data. In terms of equations, the correct calculation is

x_0_obs = x_0 * (1+z)
fwhm_obs= instrument.fwhm(x_0_obs)
fwhm = fwhm_obs / (1+z)

# or all at once
fwhm = instrument.fwhm(x_0 * (1+z)) / (1+z)

Of course, the "right" way to do things, would be to leave the data as-is, and include the redshift in the model (even if it's just fixed). I.e. a purely "forward modeling" approach. But that might not be straightforward with the attenuation etc (unless astropy model has some kind of wavelength shift operator).

[jdtsmith]

"Ignored" features are now masked out in the features table. The mask is applied in _parse_astropy_results, i.e. when the fit (or guess) results are being written to the table. The masking and unmasking happens via a new member function of Features, which uses some internal information to make sure only the relevant columns are masked out. The mask of the bounds are left untouched, as those have a special meaning (fixed / not fixed)

This mask is ignored during "guess" or "fit", as the features to be included might need to be reconsidered at that point.

Where it becomes useful is

1. When the user looks at the table contents by eye
2. For "plot", to make sure no irrelevant features are plotted
3. For the analysis functions, e.g. line strengths.

With this, the better our "feature trimming" algorithm works, the easier it will be to make good quality plots (no weird features like the plot shown in this thread by Bethany), and make sure that no badly fitted features make it into the analysis functions.

This is good insight here. And if a user wants to omit features from a given pack, they can just remove them from the table altogether. Especially good to "reconsider" each feature on each new fit, just in case.

If we shift the data to the rest frame, then all the features in that data will get narrower (wavelength points get closer together!). We should therefore adjust the FWHM in the model, so that it better represents the shifted data. In terms of equations, the correct calculation is...

Yeah good point. Maybe we should just perform the fit in the observed frame (and perhaps provide another axis with rest-frame value in the plotting).

[jdtsmith]

But that might not be straightforward with the attenuation etc (unless astropy model has some kind of wavelength shift operator).

Hmmm, ok. This is a real issue, since our attenuation laws have finite wavelength range of applicability. That's perhaps a another good reason to eventually parameterize our attenuation curves going forward, so they work at any wavelength (even if with questionable scientific validity). So for now your shift-forward/shift-back approach seems sensible.

I once experimented with freeing up redshift as a bounded parameter; it got ugly.

[karllark]

(unless astropy model has some kind of wavelength shift operator).

You mean like https://docs.astropy.org/en/stable/api/astropy.modeling.functional_models.RedshiftScaleFactor.html

[codecov-commenter]

Codecov Report

:exclamation: No coverage uploaded for pull request base (dev_yaml-scipack+api@4b04c95). Click here to learn what that means. The diff coverage is n/a.

@@                   Coverage Diff                   @@
##             dev_yaml-scipack+api     #233   +/-   ##
=======================================================
  Coverage                        ?   63.02%           
=======================================================
  Files                           ?       14           
  Lines                           ?     1025           
  Branches                        ?        0           
=======================================================
  Hits                            ?      646           
  Misses                          ?      379           
  Partials                        ?        0           

Help us with your feedback. Take ten seconds to tell us how you rate us. Have a feature suggestion? Share it here.

[drvdputt]

There are a bunch of conflicts now, because I merged in the master branch (so I could get the recent instrument fix and documentation additions). Rebasing the target branch PAHFIT:dev_yaml-scipack+api onto master should fix this.

[jdtsmith]

Quick summary of TBD's based on today's call, to get this merged soon:

• [ ] Merge master into the target devbranch, to eliminate conflicts.
• [ ] Omit 'instrumentname' and 'redshift' from pahfit.model.Model.
• [ ] Accept only Spectrum1D's with sp.meta['instrument'] (I prefer that because it's a bit shorter to type). Hopefully other tooling will eventually do this for the user.
• [ ] Also respect the redshift of the Spectrum1D object [see 1].
• [ ] Give guess and fit optional redshift-override flags.
• [ ] Pass Model's __repr__ through to its Features' __repr__, so a model is "self-printing".
• [ ] Update the example notebook using the slightly simplified API.
• [ ] Update docs (such as they are) to reflect all this.

Then once merged into dev, everyone should give the dev branch a good test with the notebook, make more doc/etc. PRs. Soon after this, we can "flip the switch" on the new API and replace main with it, making a new release.

[1] Here's how I'm doing it (untested):

        z = sp.spectral_axis.redshift.value if redshift is None else redshift

Longer term TBD (after merge):

• More YAML packs (@ThomasSYLai what's the state of the 2.5-5µm pack components?).
• override Features row flag to override the instrument pack (for FWHM).
• expose uncertainties (via feat.meta['covariance'], once available, directly in the printout).
• Tiny convenience subclass of Spectrum1D to make setting instrument more sensible (it's already got redshift). I.e. hide meta as an implementation detail.

[jdtsmith]

BTW @drvdputt, I was wrong about the Row's meta, it's not distinct, but just aliases to the full table's meta (see below). We could create our own Row class which has its own row metadata, or just create a dictionary keyed by name in the full table meta (probably easier). I'm thinking to use feature_table.meta['fit_details'] as the name space. Let me know what you think.

Update: or, as you are doing, keep Model-specific attributes like fit_info in the model object, and don't bother with the table. So Features.meta contains just fit-result related info (like covariance matrix).

    @property
    def meta(self):
        return self._table.meta

[drvdputt]

Update with latest push:

• [x] Merge master into the target devbranch, to eliminate conflicts.

• [x] Omit 'instrumentname' and 'redshift' from pahfit.model.Model.

• [X] Accept only Spectrum1D's with sp.meta['instrument'] (I prefer that because it's a bit shorter to type). Hopefully other tooling will eventually do this for the user.

• [X] Also respect the redshift of the Spectrum1D object [see 1].

• [X] Give guess and fit optional redshift-override flags.

• [X] Pass Model's repr through to its Features' repr, so a model is "self-printing".

• [X] Update the example notebook using the slightly simplified API.

• [ ] Update docs (such as they are) to reflect all this.

• Redshift can be overridden everywhere it is used.

• The print function also writes out a string showing if the model was fit or not (could contain more info in the future). I've also added a _reprhtml function (also mostly a wrapper), for prettier printing in notebooks.

• An error will be raised if no instrument is set.

• I adjusted the demonstration notebook to these changes. Will post a share link on the slack. Anyone testing this can go ahead and try it. Documentation can be changed after some feedback.

[jdtsmith]

I merged master into the dev branch. Not sure if that was the issue given the remaining conflicts. You might need to rebase/refactor your changes on top of the dev branch.

[drvdputt]

Strange that there were these conflicts, as I could merge into master without problems. The merge was not easy, but here it is. Good thing I wrote/reworked those tests. The demo notebook also still works.

[jdtsmith]

Thanks. It would be good if we understood what happened. Sounds like you pulled some commits from master into your PR that came after the dev branch was created into your PR. I don't see why that would cause such serious conflicts.

Anyway, high time we get this PR in (to the dev branch at least). Anything else left to do? A few linter complaints:

pahfit/base.py:51:38: E712 comparison to False should be 'if cond is False:' or 'if not cond:' pahfit/base.py:52:38: E712 comparison to False should be 'if cond is False:' or 'if not cond:' pahfit/base.py:90:42: E712 comparison to False should be 'if cond is False:' or 'if not cond:' pahfit/base.py:266:25: E126 continuation line over-indented for hanging indent pahfit/base.py:270:21: E126 continuation line over-indented for hanging indent

[drvdputt]

I can go ahead and fix those. By the way, the windows test isn't working because it doesn't allow the creation of files on /tmp (or the windows-equivalent thereof). @karllark You told me you know a solution for this?

[drvdputt]

Alternative solution: just disable the windows tests?

[karllark]

I can go ahead and fix those. By the way, the windows test isn't working because it doesn't allow the creation of files on /tmp (or the windows-equivalent thereof). @karllark You told me you know a solution for this?

I found this search to provide a few possible (simple) solutions. Worth a try. I'm sure we'll hear it doesn't work on windows fairly quickly - so we should put some effort in keeping the windows tests.

https://www.google.com/search?q=tempfile.NamedTemporaryFile+windows+fails&oq=tempfile.NamedTemporaryFile+windows+fails&aqs=chrome..69i57j33i160l2.3848j0j7&sourceid=chrome&ie=UTF-8

[jdtsmith]

One of my students and my postdoc use Windows! Gah.

[drvdputt]

I found the cause: tempfile opens the file and returns a file handle. I used file_handle.name, and then the astropy write function also tried to do open(). Windows doesn't like it when you do open() on a file that is already open.

Solution: use a temporary directory, and put a file with a fixed name in that. Directory and its contents are deleted when going out of the with tempfile.TemporaryDirectory() as d block.

Finally, I changed something in the documentation to make linkcheck work, and cleaned up the last bit of formatting.