Unexpected behavior from custom function SpecFromDat

[ckarwin]

I am analyzing the Crab from DC3 using the SpecFromDat custom function and noticed some odd behavior, as summarized below.

In the doc string, the normalization parameter K specifies: transformation : log10. This transformation leads to some odd behavior. First, when minimizing the likelihood, the normalization is taking negative values. The min value in the doc string is 1e-30, and so it appears as though negative numbers would not be allowed, but because of the log transformation, any input <1 maps to a negative value, which of course is non-physical. The other issue with the log transformation is that it's not clear for the user to know that this is happening, specifically when passing initial values for K. To fix this my suggestion is to simply remove the log transformation. I would also adjust the other parameters a bit: min, max, and delta.

The other issue with this function is that the normalization in the doc string is ph/cm2/s, however, the actual normalization that comes from the code is ph/cm2/s/keV. The former makes sense, since the spectrum is normalized, and K gives the integrated flux over the entire energy range. I think this can be fixed by simply removing lines 138-139 (and changing set_units to the correct units):

cosipy/cosipy/threeml/custom_functions.py

Line 138 in 0570719

if self._x_unit != None:

.

I think @krishnatejavedula worked on this function most recently. What do you think about these proposed changes? Do you think they make sense, or is there something else I'm not considering here?

[krishnatejavedula]

Hi @ckarwin,

I've updated the docstring for SpecFromDat in PR #241 to reflect the correct units for the normalization parameter K, based on our discussion a few months ago. At that time, we agreed that K should have units of ph/cm²/s, since it represents the total integrated flux.

However, while reviewing the GRB spectral fitting notebooks and other examples, I noticed that K is consistently defined with units of ph/cm²/s/keV (i.e., / u.cm / u.cm / u.s / u.keV). When I tried updating the units in the code to match the ph/cm²/s convention, it resulted in unit compatibility errors during execution. This suggests that 3ML internally expects K to be in ph/cm²/s/keV.

Because of this discrepancy, I had to revert the changes in the PR to maintain compatibility with existing notebooks. I'm not entirely sure where the discrepancy arises from and would like to discuss it further.

Regarding your earlier suggestion about redefining the parameter, would you prefer something along these lines?

K:
    desc: Normalization
    initial value: 1.0
    is_normalization: True
    min: 1e-5
    max: 1e3
    delta: 0.1
    units: ph/cm2/s/keV

[israelmcmc]

About the transformation:

• It's usually useful to let the algorithm explore unphysical regions, as long as the underlaying likelihood profile is well defined in those regions. It's up to the user to interpret the result later. For example, when fitting a background-only data, the fitted flux will be negative ~50% of the time just out of fluctuations, which the user can interpret later. If the source flux range is restricted to >0, the algorithm will have a hard time converging to 0.
• The main purpose of the transformation is helping the algorithm converge. The log transformation is usually useful when dealing with power laws and similar functions that change dramatically over the energy range. I believe that's why it is the default (which I think it's a good one) for the normalization in multiple functions in 3ML.
• This is just the default, but the users can change it if they see fit for the particular problem at hand. @ckarwin You can do e.g. spectrum.K.remove_transformation()

About the units:

• There are 2 "anchor" units that I think make sense for SpecFromDat: K should be ph/cm2/s (as @ckarwin pointed out) and the output of evaluate() should be ph/cm2/s/keV.
• Last time I checked Revised units for parameter K (ph/cm2/s) in extended source fit tutorial to match the units in the SpecFromDat file #241 that was not yet the case. I can check again if you've fixed that @krishnatejavedula, please let me know. It's been a while since I went through the PR, sorry.
• I don't think 3ML expects a particular value for K, as long as the units are worked out internally such that the output units make sense. These can be dealt with internally instead of changing the input units of K

[ckarwin]

@israelmcmc @krishnatejavedula

About the transformation:

It's usually useful to let the algorithm explore unphysical regions, as long as the underlaying likelihood profile is well defined in those regions. It's up to the user to interpret the result later. For example, when fitting a background-only data, the fitted flux will be negative ~50% of the time just out of fluctuations, which the user can interpret later. If the source flux range is restricted to >0, the algorithm will have a hard time converging to 0.
The main purpose of the transformation is helping the algorithm converge. The log transformation is usually useful when dealing with power laws and similar functions that change dramatically over the energy range. I believe that's why it is the default (which I think it's a good one) for the normalization in multiple functions in 3ML.

Ok, that's reasonable. I was initially having trouble getting the fit to converge with the log transformation. But I think I might have had some bugs in my pipeline (or maybe I wasn't patient enough), and now the fit works well. It gives the same results with and without the log transformation.

This is just the default, but the users can change it if they see fit for the particular problem at hand. @ckarwin You can do e.g. spectrum.K.remove_transformation()

Thanks! I was looking for something like that.

• There are 2 "anchor" units that I think make sense for SpecFromDat: K should be ph/cm2/s (as @ckarwin pointed out) and the output of evaluate() should be ph/cm2/s/keV.
• Last time I checked Revised units for parameter K (ph/cm2/s) in extended source fit tutorial to match the units in the SpecFromDat file #241 that was not yet the case. I can check again if you've fixed that @krishnatejavedula, please let me know. It's been a while since I went through the PR, sorry.
• I don't think 3ML expects a particular value for K, as long as the units are worked out internally such that the output units make sense. These can be dealt with internally instead of changing the input units of K

Yup, I agree. Reading @israelmcmc's comments in PR #241, I think this is a good description of how it should be done.

I also think we need to modify how the dat file is read. Currently it requires a very specific format (in terms of skipped lines in header, footer, etc), but this should be more general.

@krishnatejavedula are you able to implement these changes? I can do it if you prefer.

[ckarwin]

One more thing: It's probably better if the normalization uses an actual integral, instead of the approximation being used.

cosipy/cosipy/threeml/custom_functions.py

Line 134 in 0570719

dataFlux = dataFlux / np.sum(dataFlux * ewidths)

I think the approximation might not be as accurate as needed.

[israelmcmc]

@ckarwin About the format, there's this other pending PR by @krishnatejavedula: #257

[ckarwin]

Ah, thanks for pointing that out, @israelmcmc. It would be good to follow-up on these two PRs. @krishnatejavedula Let me know if I can help with anything.

[israelmcmc]

Indeed. Thanks @ckarwin!

[israelmcmc]

@krishnatejavedula Do you have an update on this issue? No pressure, I'm just cleaning up the issue list and trying to see where we stand.

[krishnatejavedula]

Hi @israelmcmc The issue is still open. I will try to close it by the end of this week.