Power Drude and Gaussian

[drvdputt]

Here is the next pull request. I have just rebased this branch on dev.

Replaces #225 , and this should be the final part of #280 .

Changes with respect to the user experience:

• the power column of the Features table now has a unit
• power now actually means power (no longer amplitude)
• fitted power values will be reasonable numbers
• only spectra compatible with MJy / sr are supported

[drvdputt]

I changed a use of trapz to trapezoid to avoid a deprecation warning from numpy. But apparently that's a numpy 2.0 thing. If require numpy >=2.0 as a dependency, that should fix the tests (do we want this?).

[jdtsmith]

Thanks for working on this. RE: numpy 2.0, we'll have to go there at some point, so might as well make ourselves compliant now. Will add review separately.

[jdtsmith]

Looks good overall. I do want to make sure we "save room" for easily adding support for flux units (and not just intensity units). I'm also not sure I understand the "assume all units are the internal PAHFIT units". Rather than assume, can't we enforce this by directly converting? I guess this is at the level of Fitter so it's a safe assumption. But if/when we do accept flux units, we'll need Fitter's to be versant in units.flux_power.

[jdtsmith]

Hi @drvdputt, would be great to get this merged to dev and then to master soon.

[drvdputt]

I pushed some changes.

For the numpy trapz (<2.0.0) vs trapezoid (>=2.0.0) issue, I have found a workaround: use scipy.integrate.trapezoid instead.

To deal with the flux/intensity split, there will have to be some logic in certain places.

• In PowerDrude and PowerGaussian, a different value of intensity_amplitude_factor will have to be chosen.
• In the initial guess, see the comment above.
• In Model._convert_spec_data, I currently validate the user input units, and then convert the quantities. This function could be made to return a flag that tells us which version of the above switches to use

I should note that I find that most of the differences come down to the differences in ratios between flux_density / flux_power and intensity_power / intensity. If those two quantities were the same, things would be more trivial. But at the moment, the chosen ratios are as follows

In [14]: (u.flux_density / u.flux_power).decompose()
Out[14]: Unit("1e-07 s")

In [15]: (u.intensity / u.intensity_power).decompose()
Out[15]: Unit("1e-10 s")

[jdtsmith]

For the numpy trapz (<2.0.0) vs trapezoid (>=2.0.0) issue, I have found a workaround: use scipy.integrate.trapezoid instead.

Did they really remove one name and add another? Sigh.

I should note that I find that most of the differences come down to the differences in ratios between flux_density / flux_power and intensity_power / intensity. If those two quantities were the same, things would be more trivial. But at the moment, the chosen ratios are as follows

Can you expand on this? Why is the ratio significant? If it simplifies treating the two flavors the same, we could make both of these ratios 1e-9, for example.

[jdtsmith]

In PowerDrude and PowerGaussian, a different value of intensity_amplitude_factor will have to be chosen.

BTW, with my suggestion above of normalizing on intensity for Fitter, this wouldn't be the case. It just means Model would have to be able to do/undo the conversion.

[drvdputt]

The ratio matters because it affect the value of intensity_amplitude_factor, e.g. in this code

 intensity_amplitude_factor = (
        (2 * units.intensity_power * units.wavelength / (constants.c * np.pi))
        .to(units.intensity)
        .value
    )

Multiplying with units.intensity_power and applying .to(units.intensity) effectively constitutes a multiplication with units.intensity_power / units.intensity. So if this value is the same as units.flux_power / units.flux_density, then the intensity_amplitude_factor will have the same value in both unit cases.

If we go with the canonical solid angle route, then we will need to be careful that the tabulated and plotted values still make sense for the user. Although we have something like that almost in place. Model.tabulate() already converts its output to the units of the spectrum provided by the user. Putting in an conditional conversion that applies this solid angle should be straightforward.

Not sure what would be expected for the power values in a model that was written to disk. Perhaps the same solid angle factor can be applied when writing/reading.

[jdtsmith]

The ratio matters because it affect the value of intensity_amplitude_factor, e.g. in this code

I see, so if we are careful with our units, the factor needed to go between power and amplitude for Drude's/Gaussian's would be agnostic as to whether we have intensity of flux units. I can see an advantage there for Fitter's: they don't have to worry about it. The same would of course hold if we use a canonical solid angle to force everything to be intensity in Fitter space. I suppose that would be cast as units.canon_solid_angle. Model should hide all of these details from the user.

Why don't we punt for now on the question of a canonical solid angle or "special units". A Model saved to disk will presumably also include the "desired output units" so will always "do the right thing". Ideally a given Model could be fitted to an intensity spectrum, then re-fitted to a flux spectrum (or over-plotted, etc.) without issue. I think a canonical solid angle would work for that, but need to think it through more. Not a blocker for this power PR though. I'll open an issue.

See #297 (for discussion of another PR).

[jdtsmith]

Hi @drvdputt, I'd like to get this PR buttoned up and onto dev so we can do a bit of testing before merging to master. Can you let me know what remains to be done?

[drvdputt]

I reread the comments and looked over my code changes again, and I think no more changes are needed at this point. I think we can merge to dev, and the other minor issues can be addressed in smaller/quicker pull requests.

With this merged, that concludes all of my major changes, and the dev branch should be open for new additions by other contributors. And any other changes I develop in my experimental branches will likely be small enough to copy paste and adapt.

[jdtsmith]

Great, thanks for your work on this, Dries. I've merged to dev. @drvdputt & @sarduval, can you give dev some good testing, in particular checking the power column (which is now no longer a misrepresentation)?

In terms of other changes, I'd like to restore our old practice of smaller PRs doing ~one thing, vetted and merged to master, so if you have any small changes let's wait for master to settle. Once dev gets a thorough testing, I plan to merge it to master, and we can then close a bunch of old issues that have been addressed. Then from there, I'll build out the much-discussed SPFitter (which is laying in disarray and in the form of old notes on my disk). I'll probably have questions about our Model⇒Fitter API at that point; will open an issue for that for posterity.