Saving result from mle fit

[ssaha1234]

Dear pyhf members,
I am trying to review some post mle fit results, in a proper format. I am saving the result in a .pkl file as below -

result, result_obj = pyhf.infer.mle.fixed_poi_fit(0.0, data, model, return_uncertainties=True, return_result_obj=True)
retval = [result[:,0],result[:,1],result_obj.corr,float(sys.argv[2])]
pickle.dump(retval, open("Bonly.pkl","wb"))

I do get a list of post-fit values when I print the contents of the .pkl file, but I would like to have access to the exact labels I have used in the input file, to understand the numbers. Is there an elegant way to do the same ?

Thanks!
Best,
Shreya

[alexander-held]

I have not tried this, but I guess you could add model.config.par_names to your retval list to get those labels also saved?

[ssaha1234]

Hi @alexander-held,
Thanks for the good point! I will try it now :))
Best,
Shreya

[ssaha1234]

Hi @alexander-held,
It does print out the par_names, as you suggested, but I am not completely sure if the printed values are in the same order as the order of the par_names. I will try to make it a bit more readable, in the meanwhile.
Thanks!
Best,
Shreya

[alexander-held]

The order in model.config.par_names matches exactly the order of the parameters reported from the fit, it is what pyhf uses internally for ordering them (#1747 is related). Perhaps it would be useful to add the piece of information "this matches the order in which your fit results are returned" explicitly to the documentation as well though.

[ssaha1234]

Hi @alexander-held,
Okay that is perfect! I will add a note in my script, thank you very much for the clarification.
Best,
Shreya

[matthewfeickert]

@ssaha1234 as Alex has already answered your question very nicely, I will just make a suggestion that you avoid using pickle in your workflows if possible. You may already be very well aware of all the reasons not to use pickle and have intentionally chosen it because it is the best tool for your needs, and in that case please carry on. Regardless, my apologies if this comes across as telling you what to do with your own code — not my intent.

The Python developers have a nice big warning right at the top of the pickle module docs on why not to use it

Warning: The pickle module is not secure. Only unpickle data you trust.

It is possible to construct malicious pickle data which will execute arbitrary code during unpickling. Never unpickle data that could have come from an untrusted source, or that could have been tampered with.

Consider signing data with hmac if you need to ensure that it has not been tampered with.

Safer serialization formats such as json may be more appropriate if you are processing untrusted data. See Comparison with json.

and while I'm sure you trust yourself and your colleagues, I believe [citation needed] that there are instances where pickle is not stable across different Python runtimes (edit: That was incorrect of me. I was apparently thinking back to the old days of Python 2 and Python 3 as pickle is now stable and you need to now work to use the unstable C optimizer version: import _pickle as cPickle. Source: "What’s New In Python 3.0" library changes section). So in general I would suggest only using pickle if you are going to be using the exact same version of Python to unpickle that file on the exact same computer later.

To try to be more constructive, I appreciate that there are a huge amount of data serialization options out there depending on what you want to do. But if you're already using pickle and not trying to store arrays (as you'd have gone with .npz) then JSON might be an easy alternative to use until you find a need to go to something more extensive or robust.

Example:

import json

import pyhf

pyhf.set_backend("numpy", "minuit")

model = pyhf.simplemodels.uncorrelated_background(
    signal=[12.0, 11.0], bkg=[50.0, 52.0], bkg_uncertainty=[3.0, 7.0]
)
data = [51, 48] + model.config.auxdata
test_poi = 0.0

result, result_obj = pyhf.infer.mle.fixed_poi_fit(
    test_poi, data, model, return_uncertainties=True, return_result_obj=True
)

# I think your sys.argv[2] is maybe similar to https://github.com/scikit-hep/pyhf/discussions/2004 ?
# so I'll guess here and make it at least runnable
region_number = 2

# Just using .tolist here as we know these are numpy
results = {
    "parameter_names": model.config.par_names,
    "parameters": result[:, 0].tolist(),
    "uncertainties": result[:, 1].tolist(),
    "correlations": result_obj.corr.tolist(),
    "region_number": region_number,
}

with open("results.json", "w", encoding="utf-8") as write_file:
    write_file.write(json.dumps(results, indent=4))

$ cat results.json
{
    "parameter_names": [
        "mu",
        "uncorr_bkguncrt[0]",
        "uncorr_bkguncrt[1]"
    ],
    "parameters": [
        0.0,
        1.0030477468945027,
        0.9622145339358428
    ],
    "uncertainties": [
        0.0,
        0.05531543447327775,
        0.09471935230388406
    ],
    "correlations": [
        [
            0.0,
            0.0,
            0.0
        ],
        [
            0.0,
            0.9999999999999999,
            0.0
        ],
        [
            0.0,
            0.0,
            1.0
        ]
    ],
    "region_number": 2
}% 

Anyway, this is off topic already, but thank you Alex for responding so quickly and thank you @ssaha1234 for asking good questions here and for using pyhf!