Strong machine dependence of spey results (simplified likelihoods)

[WolfgangWaltenberger]

System Settings

In [2]: spey.about()
DEPRECATION: Loading egg at /home/walten/.venvs/311/lib/python3.11/site-packages/ml_likelihoods-0.0.0-py3.11.egg is deprecated. pip 24.3 will enforce this behaviour change. A possible replacement is to use pip for package installation.. Discussion can be found at https://github.com/pypa/pip/issues/12330
Name: spey
Version: 0.1.3
Summary: Smooth inference for reinterpretation studies
Home-page: https://github.com/SpeysideHEP/spey
Author: Jack Y. Araz
Author-email: [email protected]
License: MIT
Location: /home/walten/.venvs/311/lib/python3.11/site-packages
Requires: autograd, numpy, scipy, semantic-version, tqdm
Required-by: spey-pyhf

Platform info:            Linux-6.6.1-060601-generic-x86_64-with-glibc2.38
Python version:           3.11.6
Numpy version:            1.26.2
Scipy version:            1.10.0
DEPRECATION: Loading egg at /home/walten/.venvs/311/lib/python3.11/site-packages/ml_likelihoods-0.0.0-py3.11.egg is deprecated. pip 24.3 will enforce this behaviour change. A possible replacement is to use pip for package installation.. Discussion can be found at https://github.com/pypa/pip/issues/12330
Autograd version:         1.5
tqdm version:             4.66.1
semantic_version version: 2.10.0

Installed backend plug-ins:

- default_pdf.correlated_background (spey-0.1.3)
- default_pdf.effective_sigma (spey-0.1.3)
- default_pdf.poisson (spey-0.1.3)
- default_pdf.third_moment_expansion (spey-0.1.3)
- default_pdf.uncorrelated_background (spey-0.1.3)
- pyhf (spey-pyhf-0.1.1)
DEPRECATION: Loading egg at /home/walten/.venvs/311/lib/python3.11/site-packages/ml_likelihoods-0.0.0-py3.11.egg is deprecated. pip 24.3 will enforce this behaviour change. A possible replacement is to use pip for package installation.. Discussion can be found at https://github.com/pypa/pip/issues/12330
Name: spey-pyhf
Version: 0.1.1
Summary: pyhf plugin for spey interface
Home-page: https://github.com/SpeysideHEP/spey-pyhf
Author: Jack Y. Araz
Author-email: [email protected]
License: MIT
Location: /home/walten/.venvs/311/lib/python3.11/site-packages
Requires: pyhf, spey
Required-by: 

- pyhf.simplify (spey-pyhf-0.1.1)
- pyhf.uncorrelated_background (spey-pyhf-0.1.1)
- ml.likelihoods (ml-likelihoods-0.0.0)
DEPRECATION: Loading egg at /home/walten/.venvs/311/lib/python3.11/site-packages/ml_likelihoods-0.0.0-py3.11.egg is deprecated. pip 24.3 will enforce this behaviour change. A possible replacement is to use pip for package installation.. Discussion can be found at https://github.com/pypa/pip/issues/12330
Name: ml-likelihoods
Version: 0.0.0
Summary: 
Home-page: 
Author: 
Author-email: 
License: 
Location: /home/walten/.venvs/311/lib/python3.11/site-packages/ml_likelihoods-0.0.0-py3.11.egg
Requires: 
Required-by:

Describe the bug

The simplified likelihoods in spey show some strong machine dependence. E.g. the code snippet below sometimes produce the correct:

spey oUL(mu)=0.2842

on other machines -- same versions of spey, autograd, scipy, numpy -- it produces e.g.:

spey oUL(mu)=0.0000

I traced back the problem, it is imho due to hypotest_base.py around line 892:

sigma_mu = self.sigma_mu(muhat, expected=expected) if muhat != 0.0 else 1.0

on the machines with the right result muhat = 0.0, therefore sigma_mu becomes 1.0
on the problematic machines muhat = 4e-24 != 0.0, sigma_mu becomes 1.4e-21 and the low_init, hi_init bracket becomes super small, the ul is not found.
I think, muhat != 0.0 has to be replaced with a floating-point-correct (in-)equality comparison operation.

Source code

#!/usr/bin/env python3

obsN=[138.0, 91.0, 14.0, 3.0, 54.0, 38.0, 4.0, 0.0, 8.0, 2.0, 4.0, 0.0, 1.0, 3.0, 1.0, 1.0, 42.0, 6.0, 1.0, 4.0, 0.0, 0.0]
bg=[160.57, 90.439, 11.512, 2.7508, 53.529, 28.319, 2.5856, 2.6029, 5.0628, 2.1711, 0.063528, 0.88532, 2.6812, 1.2599, 0.41911, 0.67201, 33.567, 7.3285, 1.6546, 4.0289, 0.88189, 0.19966]
cov=[[174.5, 6.601, 0.2403, 0.163, 19.64, 6.506, 0.5148, 0.519, 1.754, 0.6346, 0.007154, 0.422, 1.397, 0.6174, 0.171, 0.286, -0.5543, -0.1305, 0.02306, -0.05027, -0.006686, 0.0109], [6.601, 89.09, 0.5931, -0.04077, 7.188, 7.695, 0.3474, 0.3852, 1.01, 0.4502, 0.009081, 0.2456, 0.8593, 0.4282, 0.1174, 0.235, -0.2041, -0.04318, 0.03952, -0.1545, -0.01002, 0.01648], [0.2403, 0.5931, 9.672, -0.01319, 0.5717, 0.4039, 0.4499, 0.03413, 0.08148, 0.03547, 0.001684, 0.07249, 0.07084, 0.04303, 0.01724, 0.02002, 0.03696, 0.0385, 0.02107, -0.01743, 0.001832, 0.003709], [0.163, -0.04077, -0.01319, 3.709, 0.2659, 0.1292, -0.001839, 0.5376, 0.03348, 0.00997, 0.001397, 0.05909, 0.01182, 0.005355, 0.006677, 0.03921, 0.01953, -0.01931, 0.02678, 0.01732, -0.002887, 0.004678], [19.64, 7.188, 0.5717, 0.2659, 69.77, 6.934, 0.5394, 0.8391, 1.98, 0.7593, 0.01766, 0.3011, 1.661, 0.7474, 0.216, 0.224, -0.2935, -0.01572, 0.04757, -0.001996, 0.01541, 0.01112], [6.506, 7.695, 0.4039, 0.1292, 6.934, 27.89, 0.3204, 0.3368, 1.036, 0.4105, 0.008852, 0.1561, 0.8648, 0.3931, 0.1036, 0.1509, -0.2699, -0.02976, 0.006113, -0.06488, 0.01266, 0.01277], [0.5148, 0.3474, 0.4499, -0.001839, 0.5394, 0.3204, 1.699, 0.01839, 0.08511, 0.02799, 0.0003195, 0.001698, 0.06698, 0.02765, 0.01019, 0.004009, -0.06498, -0.01117, 0.01824, 0.01453, 0.002514, 0.002796], [0.519, 0.3852, 0.03413, 0.5376, 0.8391, 0.3368, 0.01839, 5.76, 0.09795, 0.05031, 0.002012, 0.04714, 0.078, 0.03875, 0.01596, 0.02729, 0.0001083, -0.01461, 0.01978, 0.02572, 0.004099, 0.0067], [1.754, 1.01, 0.08148, 0.03348, 1.98, 1.036, 0.08511, 0.09795, 2.108, 0.1177, 0.002927, 0.03125, 0.5298, 0.1115, 0.02586, 0.03543, 0.03278, -0.006467, 0.005036, 0.002551, -0.002566, 0.0006324], [0.6346, 0.4502, 0.03547, 0.00997, 0.7593, 0.4105, 0.02799, 0.05031, 0.1177, 0.5019, 0.001326, 0.0148, 0.09351, 0.1455, 0.01469, 0.01873, 0.02184, 0.002662, 0.00445, -0.004747, 0.0002134, 0.0007696], [0.007154, 0.009081, 0.001684, 0.001397, 0.01766, 0.008852, 0.0003195, 0.002012, 0.002927, 0.001326, 0.0121, 0.0001844, 0.001437, 0.001037, 0.01401, 0.0009814, -0.0007443, -8.312e-05, 0.0009001, 0.0009795, 0.0001779, 0.0002314], [0.422, 0.2456, 0.07249, 0.05909, 0.3011, 0.1561, 0.001698, 0.04714, 0.03125, 0.0148, 0.0001844, 4.761, 0.03474, 0.01384, 0.007457, 1.983, 0.07728, 0.0002071, 0.009179, 0.07125, 0.007518, 0.004858], [1.397, 0.8593, 0.07084, 0.01182, 1.661, 0.8648, 0.06698, 0.078, 0.5298, 0.09351, 0.001437, 0.03474, 0.8703, 0.0914, 0.02294, 0.02833, 0.009643, -0.001688, 0.003651, -0.005076, -0.0001599, 0.001133], [0.6174, 0.4282, 0.04303, 0.005355, 0.7474, 0.3931, 0.02765, 0.03875, 0.1115, 0.1455, 0.001037, 0.01384, 0.0914, 0.2891, 0.01253, 0.01471, 0.00162, 0.002663, 0.00349, -0.002747, -0.001798, -0.0003085], [0.171, 0.1174, 0.01724, 0.006677, 0.216, 0.1036, 0.01019, 0.01596, 0.02586, 0.01469, 0.01401, 0.007457, 0.02294, 0.01253, 0.3721, 0.005811, 0.006178, 0.004303, 0.001909, 0.002492, 0.0006547, 0.0008085], [0.286, 0.235, 0.02002, 0.03921, 0.224, 0.1509, 0.004009, 0.02729, 0.03543, 0.01873, 0.0009814, 1.983, 0.02833, 0.01471, 0.005811, 1.776, 0.06913, 0.008452, 0.01022, 0.03104, 0.007245, 0.005361], [-0.5543, -0.2041, 0.03696, 0.01953, -0.2935, -0.2699, -0.06498, 0.0001083, 0.03278, 0.02184, -0.0007443, 0.07728, 0.009643, 0.00162, 0.006178, 0.06913, 33.32, 5.406, 1.608, 1.192, 0.0703, 0.07672], [-0.1305, -0.04318, 0.0385, -0.01931, -0.01572, -0.02976, -0.01117, -0.01461, -0.006467, 0.002662, -8.312e-05, 0.0002071, -0.001688, 0.002663, 0.004303, 0.008452, 5.406, 3.271, 0.5664, 0.07033, 0.262, 0.03766], [0.02306, 0.03952, 0.02107, 0.02678, 0.04757, 0.006113, 0.01824, 0.01978, 0.005036, 0.00445, 0.0009001, 0.009179, 0.003651, 0.00349, 0.001909, 0.01022, 1.608, 0.5664, 0.8431, 0.07743, 0.03581, 0.1028], [-0.05027, -0.1545, -0.01743, 0.01732, -0.001996, -0.06488, 0.01453, 0.02572, 0.002551, -0.004747, 0.0009795, 0.07125, -0.005076, -0.002747, 0.002492, 0.03104, 1.192, 0.07033, 0.07743, 1.903, 0.203, 0.01753], [-0.006686, -0.01002, 0.001832, -0.002887, 0.01541, 0.01266, 0.002514, 0.004099, -0.002566, 0.0002134, 0.0001779, 0.007518, -0.0001599, -0.001798, 0.0006547, 0.007245, 0.0703, 0.262, 0.03581, 0.203, 0.13, 0.04504], [0.0109, 0.01648, 0.003709, 0.004678, 0.01112, 0.01277, 0.002796, 0.0067, 0.0006324, 0.0007696, 0.0002314, 0.004858, 0.001133, -0.0003085, 0.0008085, 0.005361, 0.07672, 0.03766, 0.1028, 0.01753, 0.04504, 0.04159]]
nsig=[0.0, 0.0002892207, 0.0004011771, 0.005287607475, 0.000170267025, 0.000746376, 0.0007673678249999999, 0.0069832804500000005, 0.0007930245, 0.00267295905, 0.001835618475, 0.031107552225, 0.0008769918, 0.005873046149999999, 0.007456762725000001, 0.085919539725, 0.59161307046, 0.7652938326899998, 5.264143279499999, 0.5127043338, 0.63509786919, 4.421279522100001]
analysis='CMS-SUS-20-004'
lumi=137.0

import spey

                                                                               
stat_wrapper = spey.get_backend("default_pdf.correlated_background")           
speyModel = stat_wrapper( data = obsN, background_yields = bg,
    covariance_matrix = cov, signal_yields = nsig,
    xsection = [ x / lumi for x in nsig ], analysis = analysis )

print ( f"spey oUL(mu)={speyModel.poi_upper_limit( ):.4f}" ) 
print ( f"spey eUL(mu)={speyModel.poi_upper_limit( expected = spey.ExpectationType.aposteriori ):.4f}" ) 
# the upper limits

Tracebacks

No response

Expected behaviour

spey oUL(mu)=0.2842
spey eUL(mu)=0.3516

which is the result obtained on those machines that just happen to compute muhat to be == 0.0.

Additional information

This is for sure not the only part of spey that is machine dependent. However, one step at a time. Expect more
such trouble tickets in the future.

Existing GitHub issues

• I have searched existing GitHub issues to make sure the issue does not already exist.

[jackaraz]

Thanks @WolfgangWaltenberger; in the future, can you also report the machines that you observe the difference so that I can try to find the exact machine configuration? Note that if I don't have access to the machine I won't be able to replicate the issue :/ This particular issue seems to be related to machine precision and optimiser performance tagged to it. Is this solved for other Scipy versions?

Also, this might be due to Python 3.11, I haven't tested Spey with 3.11 yet.

The ultimate solution may be adding an iminuit optimiser, but that will take some time to implement.

[WolfgangWaltenberger]

Hey Jack, I can try and give you all that info. That said, statements "if floating_variable != specific_value:" are clearly, undoubtedly, horribly wrong. In this case, this is the culprit. It seems a given that it is not python 3.11, or a specific version of scipy or any of this. I can give you all the version numbers, but are you sure this is helping in any way? All digits of 64-bit floats after the 16th are a lottery. There will be other machine dependent results coming up that require all these versions, don't you worry ;)

[jackaraz]

Hi @WolfgangWaltenberger, if your solution solves it for all the machines I can definitely add it or maybe add a clip on the mu value? Let me know if there is a one-fits-all type of solution and I'll update the release ASAP.

[WolfgangWaltenberger]

All that said:
*) ubuntu 23.10, python3.11, scipy 1.10.0, numpy 1.26.2, it so happens that right now i get: oUL=0, eUL=inf. Could be different next time, on an ARM computer, etc etc etc.
*) on a centos 7, python3.10, scipy version 1.10.0, numpy 1.26.1, i currently get oUL=0.2842, eUL=0.3516, the right answers.
I do recall even different answers, but thats enough for now i guess.

[jackaraz]

How about

sigma_mu = self.sigma_mu(muhat, expected=expected) if not np.isclose(muhat, 0.0) else 1.0

[jackaraz]

Hi @WolfgangWaltenberger, I also implemented the update within the issue23 branch and attached a PR to it. On my machine, your example works perfectly without any interruption, and I get the following result:

spey oUL(mu)=0.0016
spey eUL(mu)=13.3741

Can you please confirm this on your side as well?

[WolfgangWaltenberger]

This is a realistic example, the data are from CMS-SUS-20-004, testing a very average simplified model point. So, eUL being 4 orders of magnitude bigger than oUL is absolutely impossible. The values I gave above were realistic. I will try to
illustrate with a plot.

[jackaraz]

This is a realistic example, the data are from CMS-SUS-20-004, testing a very average simplified model point. So, eUL being 4 orders of magnitude bigger than oUL is absolutely impossible. The values I gave above were realistic. I will try to illustrate with a plot.

Maybe the issue is originating from an earlier step? perhaps, $\hat{\mu}$, max llhd or asimov data?

[WolfgangWaltenberger]

Now I am the one who cannot reproduce. Using your patch I get on one machine
oUL(mu)=0.2842, eUL(mu)=0.3516 (thats sensible)
on another i get
oUL(mu)=0, eUL(mu)=0.9062.
pretty random. I tried to visualise nll(mu) in a sensible range, e.g. (-.5, 2.) but i get incompatible nonsense

[jackaraz]

HI @WolfgangWaltenberger, I need a more specific issue. I compared the tool with the previous application (including Bill's version) and got the same result. If you can point out what needs to be changed, I can implement it accordingly. Otherwise, I'm afraid I don't know how to help here.

[jackaraz]

@WolfgangWaltenberger is there any update on this? or is it solved?