search

Starfish-develop / Starfish

Tools for Flexible Spectroscopic Inference
https://starfish.readthedocs.io/en/latest/
BSD 3-Clause "New" or "Revised" License
68 stars 22 forks source link

Calibrated Spectra #139

Closed mileslucas closed 3 years ago

mileslucas commented 3 years ago

fixes #138

TODO list:

This PR incorporates methods for absolute flux calibration throughout the emulator and spectral model. The interpolator is added to the Emulator following @gully's code https://github.com/gully/Starfish/blob/fe9879b67a12ddd1a101cf07c60de25bdb46a705/Starfish/emulator.py#L367

This has changed the data structure for the HDF5 files of the emulator, so I've noted a breaking change in the versioning (0.3 -> 0.4).

The only interface changes to the emulator are the new method emu.norm_factor(params) which just returns the interpolated scaling factor, and I've added a norm keyword argument to load_flux--

flux = emu.load_flux(params)
absolute_flux = emu.load_flux(params, norm=True)

Inside the spectrum model, I've added a keyword and property called normalize, which controls whether absolute flux calibration is applied to the flux--

model = SpectrumModel(
    "F_SPEX_emu.hdf5",
    data,
    grid_params=[6800, 4.2, 0],
    Av=0,
    log_scale=-13.477,
    normalize=True,
    global_cov=dict(log_amp=38, log_ls=2),
)
# or after instantiation
model.normalize = True

In addition, for convenience I've made it a little easier to get a starting point for log_scale. If you don't have any prior info, you can instantiate it without the log_scale parameter like 

model = SpectrumModel(
    "F_SPEX_emu.hdf5",
    data,
    grid_params=[6800, 4.2, 0]
)

but now you'll see in the representation the value that is fit by renorm internally (within the model.__call__ method)

SpectrumModel
-------------
Data: Example Spectrum
Emulator: F_SPEX_emu
Log Likelihood: None

Parameters
  T: 6800
  logg: 4.2
  Z: 0
  log_scale: -0.020519147372786654 (fit)

the log_scale is not in the parameter dictionary, so it won't be fit by the user, but this way you can still see it. If you want to be really cheeky and steal the value, just access _log_scale

model = SpectrumModel(...)
# call once so _log_scale gets set
model()
model["log_scale"] = model._log_scale

@iyera22 can you give this a try and let me know how it works for you? I'd like to find any bugs or pain points before we push and release a new version. You can update to this version within your virtualenv with

$ pip install -e git+https://github.com/iancze/Starfish@ml/scaling#egg=astrostarfish

like I mentioned above, the saved emulator structure is different, so you'll have to retrain your emulator. If that is a serious show-stopper, let me know and I can walk through exactly what data to edit/add in your saved HDF5 file to make it compatible with the new emulator structure.

@gully @iancze can you take a look through the Emulator and model code and see if anything seems egregiously incorrect? Thankfully not much has changed, but I'm also eager for comments on the documentation and any tests.

aishaiyer commented 3 years ago

Thank you @mileslucas !! Yes, this version works fine on my end. The emulator reconstructed absolute flux vs raw spectrum (of my stellar model) has fractional errors less ~1ppt over the bandpass I'm focusing on (comparable to your plot in #138 but also since I'm working with lower resolution). Additionally, I agree with @iancze, the log_scale parameter defined this way (to be used for physical (R/D)^2 scaling) is definitely super sensitive to Teff, I think this also answers my question on why the Teff posteriors are badly constrained from my plot in #138 when the log_scale is a large negative number, particularly a problem when working with H- or V-band calibrated cool star spectra in my case. Perhaps with ZJ's+Gully's upcoming version release one could constrain both the log_scale and the Teff to be well sampled?