Adding an arbitrary offset and scaling to atmospheric models

[gabrielastro]

When fitting photometry or spectra (as in https://species.readthedocs.io/en/latest/tutorials/fitting_model_spectra.html), would it be possible, for exploratory purposes, to have the possibility of fitting for an abitrary linear flux term and for a scaling of the theoretical spectrum (e.g. PHOENIX, BT-Settl, ATMO, Exo-REM, etc.)? This can be currently partly done by adding a blackbody component when fitting (using disk_teff and disk_radius in bounds= for FitModel()) but this component cannot yield a negative flux; if the model has the same spectral features (as strong peak-to-peak) as the data but is too bright, adding a blackbody will not work.

Thanks, and I hope this makes sense :)!

[tomasstolker]

That should be possible! Not sure if I fully understand the blackbody part. Would you fit in this case an arbitrary scaling instead of using the radius and parallax/distance parameters for the scaling? And then as second parameter an offset applied to the model spectrum?

[gabrielastro]

Using a blackbody is just a partial workaround but I would like to fit a regular atmospheric model but which can be modified by an offset and a scaling. Currently, one can only scale the model spectrum (through the radius), but not offset it arbitarily. So, yes, I would fit for an arbitary scaling and offset. It would be sensible for the offset to be a linear function of wavelength (and not only a constant).

[tomasstolker]

Ah I see, so you used the blackbody for the offset. Should the scaling and offset parameters be sampled linearly? Or logarithmic?

[gabrielastro]

Exactly! I guess the scaling of the atmospheric model spectrum should be sampled logarithmically since it should not be negative. Actually, for this one can use the radius (but sampled logarithmically for that reason); it will not have physical meaning but this is ok for these purposes. The new part is the linear offset function, whose slope and intercept should be sampled linearly since they can be positive or negative. I guess their needed range is not clear in general and cannot be derived physically (and will depend on the atmospheric model). Maybe useful parameters scales are Delta_offset = (mean of starting theoretical spectrum) - (mean of data spectrum) and Delta_slope = (mean slope of theoretical spectrum = quick, rough least-squares fit) - (mean slope of data), i.e., starting at offset = 0 but stepping around with Delta_offset, and similarly for slope)?

[tomasstolker]

I have added the flux_scaling, log_flux_scaling, and flux_offset (in W m-2 um-1) parameters in commit 4804241eba3e8bf16babf21660cbb38bd02eb050. These parameters can be added to bounds in FitModel. They are also supported by ReadModel. For now I have not included the offset slope, perhaps it is not needed.

[gabrielastro]

Thanks a lot! Sorry for the delay. I wanted first to check "quickly" simple, classical fitting and then report but I have not quite succeeded yet, so I want to come back to this once the rest works…

[gabrielastro]

Finally! I am able to try it out. Thanks a lot again. However, I get:

KeyError: 'radius'
Exception ignored on calling ctypes callback function: <function run.<locals>.loglike at 0x7f88492a4160>
Traceback (most recent call last):
  File "[…]/.local/lib/python3.9/site-packages/pymultinest/run.py", line 228, in loglike
    return LogLikelihood(cube, ndim, nparams)
  File "[…]/species/fit/fit_model.py", line 2018, in _lnlike_multinest
    return self._lnlike_func(params)
  File "[…]/species/fit/fit_model.py", line 1375, in _lnlike_func
    params[self.cube_index["radius"]],

i.e., the radius does not exist as a free parameter anymore yet it would be needed for logg. (By the way, it would be better if the Exception of KeyError were not ignored; can this be done from species or would it require the user to recompile pymultinest?)

I guess a solution is instead to let the distance (1/parallax) be a completely free parameter, or be the real distance times a free parameter. Would this break something else, though?

[tomasstolker]

The radius is not needed for logg, so should be removed from bounds. Then it will hopefully work!

[gabrielastro]

The problem is that the 'radius' is getting set in species, not by the user. I do:

Spekkalpar = {'teff': (1000, 1900), 'log_flux_scaling': (-24.0, -18), 'flux_offset': (-1e-15, 1e-15)}
fit = FitModel(…, bounds=Spekkalpar)

but FitMode prints:

Fitting 5 parameters:
   - teff
   - logg
   - feh
   - log_flux_scaling
   - flux_offset

Uniform priors (min, max):
   - teff = (1000, 1900)
   - log_flux_scaling = (-24.0, -18)
   - flux_offset = (-1e-15, 1e-15)
   - logg = (3.0, 5.5)
   - feh = (-0.6, 0.3)

and indeed, after the call to FitModel, Spekkalpar contains additionally 'logg': (3.0, 5.5), 'feh': (-0.6, 0.3).

[tomasstolker]

Do you use a mass prior? That should not be used if the radius is not a free parameter.

[gabrielastro]

Oops! Thank you, indeed I was. I thought I had deactivated that but apparently it was still there. (Or was it because the mass prior remained in an internal variable in FitModel between a first call with normal_prior=… and a second one with normal_prior=None in the same session?) Maybe warn the user if the input is non-consistent in this way… In any case, thanks a lot!

[gabrielastro]

Sorry, two additions: 1.) It works :partying_face:! Thanks! 2.) It would be good/needed to add automatically the a_flux and b_flux parameters to the plot_spectrum legend, as a reminder that they were used.

[tomasstolker]

I changed the implementation a bit. Any parameter that was used by a model should be included by default now (apart from parallax and distance).

[gabrielastro]

Hm, thanks but for the legend, I would prefer if the default did not include the mass and luminosity because they are only derived and the legend entries get too long, forcing hand tweaking (that option is very good but up to now I was very happy with the default, which is practical) or tempting users to shrink the font, which is usually a bad idea

[tomasstolker]

Sure! I will remove the mass and luminosity by default.

[gabrielastro]

Small further suggestions while you are at it (thanks so much!):

• Have the b_flux be in "pretty scientific notation" (not with 9.7e-17 but rather 9.7×10⁻¹⁷)
• Let the negative numbers (e.g. a_flux typically) be in math mode, which gives a minus sign (longer than a hyphen) :nerd_face:
• At the end of the fit, print to screen not only the median but also the ±1 sigma values (this is more general than this thread)?