Thick Target (and thermal) Update - #291
Conversation
…d_equivalent_units()
…e. Similar to other functions.
|
|
||
| __all__ = ["ThickTarget", "ThinTarget"] | ||
|
|
||
| FLOAT_EPSILON_FOR_POWER_LAW = 1 + FLOAT_EPSILON * 1e30 |
There was a problem hiding this comment.
Where does the FLOAT_EPSILON * 1e30 come from, trail and error or something else?
There was a problem hiding this comment.
Trial and error unfortunately. I needed a value as close to 1 as I could get but any of the Python epsilons added to 1 just returned 1 annoyingly. The FLOAT_EPSILON * 1e30 was the lowest value I could find that would return 1+(a little something).
I.e., 1 - (1 + FLOAT_EPSILON)==0 when evaluated in the model but 1 - (1 + FLOAT_EPSILON * 1e30)==FLOAT_EPSILON * 1e30.
I'm more than open to suggestions on using something better but I couldn't find anything :/
| scaled_thick_eflux_units = u.def_unit("scaled_thick_eflux_units", 1e35 * (u.electron * u.s**-1)) | ||
| scaled_thin_eflux_units = u.def_unit("scaled_thin_eflux_units", 1e55 * (u.electron * u.cm ** (-2) * u.s**-1)) | ||
| scaled_em_units = u.def_unit("scaled_em_units", 1e49 * (u.cm ** (-3))) |
There was a problem hiding this comment.
On the naming these aren't scaled they are the thick_eflux units etc, no?
There was a problem hiding this comment.
I did go between using the word "scaled" and "normalised". Do you think I should have stuck to "normalised" instead?
The thick_eflux units should be electrons/second but the models, for convenience and fitting, also want to be able to accept something more reasonable like these scaled units. Since I made the edits to make the models more unit aware, a user can just use the models normally with the correct values and units, i.e. 1e35 electrons/second, but could also pass (as the fitting will) 1 scaled_thick_eflux_units.
I suppose it's more of a philosophical question. Do we want the models to natively use these units, then scaled_thick_eflux_units becomes thick_eflux_units, etc, or do we want the docs to just say these models use normal units like electrons/second then only use scaled_thick_eflux_units internally?
My initial thoughts went with the latter, but I don't mind working more on it.
|
@samaloney , your work in #295 looks like it applies some of the things I did in here (like the I think the only other things I've done here is
I think it would be easier to merge #295 first then deal with the merge conflicts in this PR rather than the other way. This is ultimately a simple PR in comparison to be able to just make the edits again. |
KriSun95
added
Solar Physics Knowledge
2 participants
This PR got a bit away from me but it was all in the investigation of trying to figure out why the Thick target model would return NaNs sometimes.
So, what I've done is:
norm_thick_target_eflux_units,norm_thin_target_eflux_units, andnorm_thermal_emission_measure_units.5e49 cm^-3, if we like but internally the value will be converted to the scaled version of5 new_unitswherenew_untisis1e49 cm^-3. This is unit aware and so conversions can easily happen that the code knows about rather than it being in some places in the doc string that5 cm^-3really means5e49 cm^-3.5e49 cm^-3and5 new_unitswork the same and as intended.RunTimeWarningsas they should.<<when the model evaluates. This means we don't need to haveif/elsestatements checking all the input units then assuming the units if they're not there.The above were fixes I made to keep investigating the weird NaN behaviour I've seen before, removing cases where this has happened. I then ran a fit, producing the plot below (data=thermal+thick source, model=thermal+thick model, thick=just the thick target model contribution):
From here, the thick target is returning one bin around 20 keV with a NaN.
Investigating this further, I think there are two important points:
energy_edges = np.arange(2, 15, 0.1) << u.keVinstead does not return a NaN in that bin.Still to do:
Following the Astropy Parameter docs, add a getter/setter to the flux and emission measure parameters (those that are scaled) so that if the parameter is edited outside initialising the model or calling evaluate (e.g., just doing something like.emission_measure = 3e50<<u.cm**-3) the units will be converted to the new scaled units for use internally to the model.The code I used to make the plot is here:
With output (the arrays being printed are those that contain a NaN during the parameter search during fitting with the parameters producing the array printed on the line before):