`HashPRNG` (`stdnormal`) gives different values in different architectures

[ikrommyd]

This correction which is a random number generator correction

gives different values on different architectures. It appears that it's only stdnormal while stdflat is fine.
Arm64:

In [49]: x.layout
Out[49]:
<RecordArray length="25">
    <field index="0" key="pt">
        <NumpyArray format="f" shape="25" data="29.56 60.0685 51.705 10.6553 8.59726 ... 26.9614 28.5527 36.3018 15.607 7.67689" at="0x600000cc7950"/>
    </field>
    <field index="1" key="r9">
        <NumpyArray format="f" shape="25" data="0.881836 0.879883 0.889648 0.997559 0.95752 ... 0.82373 0.672852 0.939453 0.981934 0.874512" at="0x600000cc79c0"/>
    </field>
    <field index="2" key="eta">
        <NumpyArray format="f" shape="25" data="1.83179 -0.292542 -0.904053 -2.18555 1.64502 ... -1.76221 0.784546 2.02441 -0.0594559 0.38092" at="0x600000cc72c0"/>
    </field>
    <field index="3" key="event">
        <NumpyArray format="L" shape="25" data="3749762 3749777 3749777 3749768 3749768 ... 3749856 3749862 3749861 3749875 3749883" at="0x600001ce5ba0"/>
    </field>
</RecordArray>

In [50]: rng.evaluate("stdnormal", x.pt, x.r9, x.eta, x.event)
Out[50]:
array([-0.26883039,  0.28672353, -0.57431689,  2.63626883, -0.83464608,
        0.90847055,  0.45283389, -0.41724338,  0.30672475,  0.60342825,
       -1.5795927 , -0.89867654, -1.61526027, -0.54997331, -0.00784177,
        2.00220383,  1.21941401, -0.36870525,  0.50568698,  0.28755395,
       -0.94606972,  2.23859349,  1.80704412, -1.85895929,  0.69082279])

In [51]:

In [51]: rng.evaluate("stdflat", x.pt, x.r9, x.eta, x.event)
Out[51]:
array([0.93352256, 0.87892313, 0.78439512, 0.12686755, 0.46446161,
       0.05117117, 0.46069956, 0.99490182, 0.49302054, 0.25524578,
       0.00626383, 0.50793706, 0.18885202, 0.93374629, 0.95984817,
       0.21819859, 0.77767419, 0.90941815, 0.82405071, 0.66864401,
       0.30018088, 0.69552463, 0.88240465, 0.85859066, 0.41470936])

In [52]: exit
(higgs-dna) ➜  ~ uname -m
arm64

versus x86:

In [14]: x.layout
Out[14]:
<RecordArray length="25">
    <field index="0" key="pt">
        <NumpyArray format="f" shape="25" data="29.56 60.0685 51.705 10.6553 8.59726 ... 26.9614 28.5527 36.3018 15.607 7.67689" at="0x55c8524d07c0"/>
    </field>
    <field index="1" key="r9">
        <NumpyArray format="f" shape="25" data="0.881836 0.879883 0.889648 0.997559 0.95752 ... 0.82373 0.672852 0.939453 0.981934 0.874512" at="0x55c8524d08d0"/>
    </field>
    <field index="2" key="eta">
        <NumpyArray format="f" shape="25" data="1.83179 -0.292542 -0.904053 -2.18555 1.64502 ... -1.76221 0.784546 2.02441 -0.0594559 0.38092" at="0x55c8521ee4a0"/>
    </field>
    <field index="3" key="event">
        <NumpyArray format="L" shape="25" data="3749762 3749777 3749777 3749768 3749768 ... 3749856 3749862 3749861 3749875 3749883" at="0x55c8522e6190"/>
    </field>
</RecordArray>

In [15]: rng.evaluate("stdnormal", x.pt, x.r9, x.eta, x.event)
Out[15]:
array([ 0.49944807,  1.39058204, -0.02216484, -0.08730098,  0.63380507,
       -0.74493107,  1.17528214, -1.06762684, -0.4320127 , -0.93337605,
       -0.17058287,  0.63699673,  0.98020086, -0.44526196,  0.73744931,
       -0.09207296, -0.3074458 , -0.88816113,  1.80998704, -1.1880895 ,
       -2.18910314, -1.64921238,  0.60019342, -2.07503974,  0.15367366])

In [16]: rng.evaluate("stdflat", x.pt, x.r9, x.eta, x.event)
Out[16]:
array([0.93352256, 0.87892313, 0.78439512, 0.12686755, 0.46446161,
       0.05117117, 0.46069956, 0.99490182, 0.49302054, 0.25524578,
       0.00626383, 0.50793706, 0.18885202, 0.93374629, 0.95984817,
       0.21819859, 0.77767419, 0.90941815, 0.82405071, 0.66864401,
       0.30018088, 0.69552463, 0.88240465, 0.85859066, 0.41470936])

In [17]: exit
u%                                                                                                                                                                                                                                           (coffea_0.7.x) ➜  ~ uname -m
x86_64

I'm attaching the correction file
EGMScalesSmearing_Pho_2022postEE2G.v1.json.gz

and also here's how I created this pseudo array to test this from our favorite nano_dy.root file

import awkward as ak
import uproot
import correctionlib

events = uproot.open("work/coffea_dev/coffea/tests/samples/nano_dy.root")["Events"].arrays()
x = ak.flatten(ak.zip({"pt": events.Electron_pt, "r9": events.Electron_r9, "eta": events.Electron_eta, "event": events.event, "run": events.run}))

cset = correctionlib.CorrectionSet.from_file("EGMScalesSmearing_Pho_2022postEE2G.v1.json.gz")
rng = cset["EGMRandomGenerator"]
print(rng.evaluate("stdnormal", x.pt, x.r9, x.eta, x.event))
print(rng.evaluate("stdflat", x.pt, x.r9, x.eta, x.event))

Other corrections like

In [9]: cset["EGMScaleVsRun_2022postEE"].evaluate(x.run, x.eta)
Out[9]:
array([1.00814164, 1.01591229, 1.03152823, 1.01894987, 0.99801081,
       1.03152823, 1.01894987, 1.01591229, 1.02502155, 1.02363706,
       1.0188297 , 1.01591229, 1.02372885, 1.01894987, 1.01514328,
       1.01811612, 1.01514328, 1.02341413, 1.01591229, 1.01811612,
       1.0188297 , 1.01514328, 1.00814164, 1.01591229, 1.01289177])

are the same on both architectures

[nsmith-]

How about normal? Note that technically, the exact algorithm of std::normal_distribution in the C++ standard library is implementation-defined so may differ from compiler to compiler. This is why I added a self-coded normal distribution implementaiton.

[ikrommyd]

It doesn't look like this particular correction allows normal from inspecting it (see screenshot at the top). That was my guess as well, that the standard library is being seeded differently in each architecture.

[nsmith-]

Just change the HashPRNG node in the json from "stdnormal" to "normal"

[ikrommyd]

Oh I didn't know that. It took me more time but I did this (created a correction myself based on the tutorial). Indeed normal is the same on both architectures

(higgs-dna) ➜  ~ ipython
Python 3.11.10 | packaged by conda-forge | (main, Oct 16 2024, 01:26:25) [Clang 17.0.6 ]
Type 'copyright', 'credits' or 'license' for more information
IPython 8.30.0 -- An enhanced Interactive Python. Type '?' for help.

In [1]: !uname -m
arm64

In [2]: import correctionlib.schemav2 as cs
   ...: import numpy as np
   ...: rng = np.random.default_rng(1337)
   ...:
   ...: resrng = cs.Correction(
   ...:     name="resrng",
   ...:     description="Deterministic smearing value generator",
   ...:     version=1,
   ...:     inputs=[
   ...:         cs.Variable(name="pt", type="real", description="Unsmeared jet pt"),
   ...:         cs.Variable(name="eta", type="real", description="Jet pseudorapdity"),
   ...:         cs.Variable(name="phi", type="real", description="Jet phi (entropy source)"),
   ...:         cs.Variable(name="evt", type="int", description="Event number (entropy source)"),
   ...:     ],
   ...:     output=cs.Variable(name="rng", type="real"),
   ...:     data=cs.HashPRNG(
   ...:         nodetype="hashprng",
   ...:         inputs=["pt", "eta", "phi", "evt"],
   ...:         distribution="normal",
   ...:     )
   ...: )
   ...:
   ...: ptvals = rng.exponential(scale=10.0, size=10) + rng.exponential(scale=15.0, size=10)
   ...: etavals = rng.normal(scale=1.1, size=10)
   ...: phivals = rng.uniform(-np.pi, np.pi, size=10)
   ...: eventnumber = rng.integers(123456, 123567, size=10)
   ...:
   ...: print(ptvals)
   ...: print(etavals)
   ...: print(phivals)
   ...: print(eventnumber)
   ...:
   ...: print(resrng.to_evaluator().evaluate(ptvals, etavals, phivals, eventnumber))
   ...:
[ 7.70120745 37.88023241 14.80776783 17.79245551 28.9835131   6.08912876
  5.01883723 48.73234253  8.86387115 59.46291897]
[ 0.93256302  0.74482217  0.95995788 -0.37033443  1.11965169  1.21678879
 -0.04246483 -0.89903245 -0.17400397 -1.64892723]
[-2.60203312  2.97591385  1.51105146  3.00818393  1.79806831 -2.30243081
 -2.39144184 -1.74730162  0.79105964  1.26074209]
[123459 123541 123500 123500 123491 123494 123524 123556 123506 123466]
[ 0.50668463  0.84996636  0.8753116  -0.60289565  2.0445962   0.04852718
 -1.16944172  0.64194771  0.12448748  0.53078522]

and

(coffea_0.7.x) ➜  ~ ipython
Python 3.10.13 | packaged by conda-forge | (main, Dec 23 2023, 15:36:39) [GCC 12.3.0]
Type 'copyright', 'credits' or 'license' for more information
IPython 8.21.0 -- An enhanced Interactive Python. Type '?' for help.

In [1]: !uname -m
x86_64

In [2]: import correctionlib.schemav2 as cs
   ...: import numpy as np
   ...: rng = np.random.default_rng(1337)
   ...:
   ...: resrng = cs.Correction(
   ...:     name="resrng",
   ...:     description="Deterministic smearing value generator",
   ...:     version=1,
   ...:     inputs=[
   ...:         cs.Variable(name="pt", type="real", description="Unsmeared jet pt"),
   ...:         cs.Variable(name="eta", type="real", description="Jet pseudorapdity"),
   ...:         cs.Variable(name="phi", type="real", description="Jet phi (entropy source)"),
   ...:         cs.Variable(name="evt", type="int", description="Event number (entropy source)"),
   ...:     ],
   ...:     output=cs.Variable(name="rng", type="real"),
   ...:     data=cs.HashPRNG(
   ...:         nodetype="hashprng",
   ...:         inputs=["pt", "eta", "phi", "evt"],
   ...:         distribution="normal",
   ...:     )
   ...: )
   ...:
   ...: ptvals = rng.exponential(scale=10.0, size=10) + rng.exponential(scale=15.0, size=10)
   ...: etavals = rng.normal(scale=1.1, size=10)
   ...: phivals = rng.uniform(-np.pi, np.pi, size=10)
   ...: eventnumber = rng.integers(123456, 123567, size=10)
   ...:
   ...: print(ptvals)
   ...: print(etavals)
   ...: print(phivals)
   ...: print(eventnumber)
   ...:
   ...: print(resrng.to_evaluator().evaluate(ptvals, etavals, phivals, eventnumber))
   ...:
[ 7.70120745 37.88023241 14.80776783 17.79245551 28.9835131   6.08912876
  5.01883723 48.73234253  8.86387115 59.46291897]
[ 0.93256302  0.74482217  0.95995788 -0.37033443  1.11965169  1.21678879
 -0.04246483 -0.89903245 -0.17400397 -1.64892723]
[-2.60203312  2.97591385  1.51105146  3.00818393  1.79806831 -2.30243081
 -2.39144184 -1.74730162  0.79105964  1.26074209]
[123459 123541 123500 123500 123491 123494 123524 123556 123506 123466]
[ 0.50668463  0.84996636  0.8753116  -0.60289565  2.0445962   0.04852718
 -1.16944172  0.64194771  0.12448748  0.53078522]

[ikrommyd]

For stdnormal the output is this below. stdflat is the same on both architectures as well. So numpy seeds the same way but stdnormal is the only one that differs.

[ 7.70120745 37.88023241 14.80776783 17.79245551 28.9835131   6.08912876
  5.01883723 48.73234253  8.86387115 59.46291897]
[ 0.93256302  0.74482217  0.95995788 -0.37033443  1.11965169  1.21678879
 -0.04246483 -0.89903245 -0.17400397 -1.64892723]
[-2.60203312  2.97591385  1.51105146  3.00818393  1.79806831 -2.30243081
 -2.39144184 -1.74730162  0.79105964  1.26074209]
[123459 123541 123500 123500 123491 123494 123524 123556 123506 123466]
[-0.51854573  0.64308321 -0.46808235 -0.25659287  0.40421679 -0.35262017
 -1.27378241 -0.79503632  0.60463459 -0.29245815]

In [2]: !uname -m
arm64

versus

[ 7.70120745 37.88023241 14.80776783 17.79245551 28.9835131   6.08912876
  5.01883723 48.73234253  8.86387115 59.46291897]
[ 0.93256302  0.74482217  0.95995788 -0.37033443  1.11965169  1.21678879
 -0.04246483 -0.89903245 -0.17400397 -1.64892723]
[-2.60203312  2.97591385  1.51105146  3.00818393  1.79806831 -2.30243081
 -2.39144184 -1.74730162  0.79105964  1.26074209]
[123459 123541 123500 123500 123491 123494 123524 123556 123506 123466]
[ 0.80676139  0.0052043  -1.43167971 -0.53652504  1.88115216  0.6743095
  1.50188694  0.05565888  0.11107543  1.38171626]

In [2]: !uname -m
x86_64

[ikrommyd]

@nsmith- I guess the only thing maybe is to just remove stdnormal all together and use the custom normal implementation only?

[nsmith-]

Yep seems so. We can deprecate it for now and remove later.