Error when interpolating 0 to Sn for nld when energy is outside the bounds of nld.E[-1]

[hannahberg]

With my data, I do not have data for the NLD at E=0 MeV. When I try to normalize my gSF, the extrapolation fails as it is out of bounds. OMpy version '1.1.0.dev0+ce29c32'.

----------------------------------------------------------
ValueError               Traceback (most recent call last)
<ipython-input-119-f68f11161851> in <module>
----> 1 gsfnorm.normalize()
      2 gsfnorm.plot()

~/.local/lib/python3.8/site-packages/ompy/normalizer_gsf.py in normalize(self, gsf, normalizer_nld, alpha, nld, nld_model, norm_pars, num)
    201 
    202         # experimental Gg and calc. both in meV
--> 203         B_norm = self.norm_pars.Gg[0] / self.Gg_before_norm()
    204         # propagate uncertainty of D0
    205         B_norm_unc = B_norm * self.norm_pars.D0[1] / self.norm_pars.D0[0]

~/.local/lib/python3.8/site-packages/ompy/normalizer_gsf.py in Gg_before_norm(self)
    259         self.LOG.debug("Method to compute Gg: %s", self.method_Gg)
    260         if self.method_Gg == "standard":
--> 261             return self.Gg_standard()
    262         else:
    263             # placeholder if normalizations with

~/.local/lib/python3.8/site-packages/ompy/normalizer_gsf.py in Gg_standard(self)
    338             return integral
    339 
--> 340         integral = integrate()
    341 
    342         # factor of 2 because of equi-parity `spinsum` instead of `spinsum(π)`,

~/.local/lib/python3.8/site-packages/ompy/normalizer_gsf.py in integrate()
    333             Eg, stepsize = self.norm_pars.E_grid()
    334             Ex = self.norm_pars.Sn[0] - Eg
--> 335             integral = (np.power(Eg, 3) * wgsf(Eg) * wnld(Ex)
    336                         * self.SpinSum(Ex, self.norm_pars.Jtarget))
    337             integral = np.sum(integral) * stepsize

~/.local/lib/python3.8/site-packages/ompy/normalizer_gsf.py in wnld(E)
    325 
    326         def wnld(E) -> ndarray:
--> 327             return fnld(E, self.nld, self.nld_model)
    328 
    329         def wgsf(E) -> ndarray:

~/.local/lib/python3.8/site-packages/ompy/normalizer_gsf.py in fnld(E, nld, nld_model)
    592 
    593     conds = [E <= nld.E[-1], E > nld.E[-1]]
--> 594     return np.piecewise(E, conds, [fexp, nld_model(E[conds[-1]])])
    595 
    596 

<__array_function__ internals> in piecewise(*args, **kwargs)

~/.local/lib/python3.8/site-packages/numpy/lib/function_base.py in piecewise(x, condlist, funclist, *args, **kw)
    612             vals = x[condlist[k]]
    613             if vals.size > 0:
--> 614                 y[condlist[k]] = item(vals, *args, **kw)
    615 
    616     return y

~/.local/lib/python3.8/site-packages/ompy/library.py in <lambda>(zz)
    296     logy = np.log(yy)
    297     lin_interp = interp1d(xx, logy, kind='linear', **kwargs)
--> 298     log_interp = lambda zz: np.exp(lin_interp(zz))  # noqa
    299     return log_interp
    300 

~/.local/lib/python3.8/site-packages/scipy/interpolate/polyint.py in __call__(self, x)
     72         """
     73         x, x_shape = self._prepare_x(x)
---> 74         y = self._evaluate(x)
     75         return self._finish_y(y, x_shape)
     76 

~/.local/lib/python3.8/site-packages/scipy/interpolate/interpolate.py in _evaluate(self, x_new)
    657         y_new = self._call(self, x_new)
    658         if not self._extrapolate:
--> 659             below_bounds, above_bounds = self._check_bounds(x_new)
    660             if len(y_new) > 0:
    661                 # Note fill_value must be broadcast up to the proper size

~/.local/lib/python3.8/site-packages/scipy/interpolate/interpolate.py in _check_bounds(self, x_new)
    686         # !! Could provide more information about which values are out of bounds
    687         if self.bounds_error and below_bounds.any():
--> 688             raise ValueError("A value in x_new is below the interpolation "
    689                              "range.")
    690         if self.bounds_error and above_bounds.any():

ValueError: A value in x_new is below the interpolation range.

[vetlewi]

To be more specific, if the lowest energy bin for the experimental NLD is larger than zero (ie. nld.E[0] > 0) then the interpolation made with https://github.com/oslocyclotronlab/ompy/blob/5166b6dc4156e8041c7716ed2a73d827e814bc27/ompy/normalizer_gsf.py#L591

will fail since the integration goes from 0 to Sn.

A quick workaround could be:

fexp = None
    if nld.E[0] > 0:
        fexp = log_interp1d(np.concatenate(([0], nld.E)),
                            np.concatenate(([0], nld.values)))
    else:
        fexp = log_interp1d(nld.E, nld.values)

but is not very elegant

[fzeiser]

I'd say that this is not a bu, but a feature :), somehow at least. I think the question/problem is: You need nld(Ex=0), otherwise, how do you normalize gsf with Gg?

So you somehow need to make up values for nld(Ex=0).

[fzeiser]

A quick workaround could be:

fexp = None
    if nld.E[0] > 0:
        fexp = log_interp1d(np.concatenate(([0], nld.E)),
                            np.concatenate(([0], nld.values)))
    else:
        fexp = log_interp1d(nld.E, nld.values)

but is not very elegant

I agree, this is a quick fix, but I'm not sure I would want to make it a feature, instead of maybe having a few open issue and a "galleri" of different ways to fix this. Don't know what you think about it.

[vetlewi]

I'd say that this is not a bu, but a feature :), somehow at least. I think the question/problem is: You need nld(Ex=0), otherwise, how do you normalize gsf with Gg?

So you somehow need to make up values for nld(Ex=0).

Both yes and no. I'm thinking about the case where your first bin are at say 0.1 MeV, then you will get this error. Although I guess it is kinda bad choice of binning if the first bin is at 0.1 MeV?

[fzeiser]

To be more specific: I think that it is more clear if you do something along the lines of what Vetle suggested above ((np.concatenate(([0], nld.E, np.concatenate(([0], nld.values)=) -- but in the main script. So one does it before sending it to NormalizerNLD(nld=...) or so.

[vetlewi]

A quick workaround could be:

fexp = None
    if nld.E[0] > 0:
        fexp = log_interp1d(np.concatenate(([0], nld.E)),
                            np.concatenate(([0], nld.values)))
    else:
        fexp = log_interp1d(nld.E, nld.values)

but is not very elegant

I agree, this is a quick fix, but I'm not sure I would want to make it a feature, instead of maybe having a few open issue and a "galleri" of different ways to fix this. Don't know what you think about it.

Might be smart, but the error message should be more descriptive. Maybe give a link to the issue on GitHub?

[fzeiser]

I'd say that this is not a bu, but a feature :), somehow at least. I think the question/problem is: You need nld(Ex=0), otherwise, how do you normalize gsf with Gg? So you somehow need to make up values for nld(Ex=0).

Both yes and no. I'm thinking about the case where your first bin are at say 0.1 MeV, then you will get this error. Although I guess it is kinda bad choice of binning if the first bin is at 0.1 MeV?

So anyhow it will depend both on what is you central value and the bin width.

[fzeiser]

Might be smart, but the error message should be more descriptive. Maybe give a link to the issue on GitHub?

Yep, thought about this, too, when I wrote it :+1: !

[vetlewi]

Maybe something like

ValueError: nld does not extend to 0 MeV. Please see https://github.com/oslocyclotronlab/ompy/issues/170 for more info.

[fzeiser]

I can make a suggestion and write up a small example of how I would include it in the main script.

[fzeiser]

So one solution if you miss e.g. the lowest bin is to come up with a value for it (interpolate / find from known levels / ...).

fill_value = np.array([0, 0])  # (Ex, value) 
for nld in extractor.nld:
    nld.E = np.concatenate(([fill_value[0]], nld.E))
    nld.values = np.concatenate(([fill_value[1]], nld.values))

# diff = np.diff(nld.E)
# assert len(set(diff)) == 1, "nld does not have equal energy spacing (anymore)"

You may want to make sure that you preserve a linear energy spacing, so maybe you have to all several bins, and not just one. For example, you may (have to) append [0, 10] and [200, 20] to the nld (with the given value, assumes that the nld is 10 at 0 and 20 at 200.)

[fzeiser]

Alternatively, one can create a child class from NormalizerGSF and simply change the functions fnld and/or fgsf. There one can e.g. pass arguments for loginterpolate, which bases on scipy.interpolate.interp1D file fill_value= 0 or fill_value=extrapolate -- or create more complex functions.

[vetlewi]

Maybe reopen so that it's visible?

[vetlewi]

Alternatively, one can create a child class from NormalizerGSF and simply change the functions fnld and/or fgsf. There one can e.g. pass arguments for loginterpolate, which bases on scipy.interpolate.interp1D file fill_value= 0 or fill_value=extrapolate -- or create more complex functions.

Example:

import ompy as om
class MyNormalizerGSF(om.NormalizerGSF):
    @staticmethod
    def fnld(E, nld, nld_model):
        fexp = log_interp1d(nld.E, nld.values, fill_value=(0, np.max(nld.values), bounds_error=False)
        conds = [E <= nld.E[-1], E > nld.E[-1]]
        return np.piecewise(E, conds, [fexp, nld_model(E[conds[-1]])])

nldnorm = MyNormalizerGSF(...)

[fzeiser]

Maybe reopen so that it's visible?

People get an explicit link to this issue now, so I don't think that this is necessary.

[vetlewi]

One reason this may happen could be that when running the extractor the extracted nld/gsf may become nan for energies above 0 MeV. Maybe add a check at the end of the Extractor.decompose method that issues a warning whenever the extracted NLD has nan for excitation energies at or above 0 MeV.

[fzeiser]

Good idea. I will add it.