scikit-hep / pyhf

pure-Python HistFactory implementation with tensors and autodiff
https://pyhf.readthedocs.io/
Apache License 2.0
276 stars 82 forks source link

Minimization problem with pyhf.infer.mle.fit(data, model, return_uncertainties=True) #1577

Open tancredicarli opened 2 years ago

tancredicarli commented 2 years ago

Description

Dear pyhf developpers,

I have a rather complex workspace and I wanted to make a pull plot based on the tutorial. I get a crash in the minimization step.

$ python test.py
open filename=  ex2.txt
Calculate the observed and expected limit: 
/home/tcarli/pyhftest/pyhf-env/lib/python3.8/site-packages/pyhf/tensor/numpy_backend.py:353: RuntimeWarning: invalid value encountered in log
  return n * np.log(lam) - lam - gammaln(n + 1.0)
      Observed CLs: 0.0009
Expected CLs(-2 σ): 0.0000
Expected CLs(-1 σ): 0.0004
Expected CLs( 0 σ): 0.0043
Expected CLs( 1 σ): 0.0380
Expected CLs( 2 σ): 0.2017
     corr: None
      fun: 827.5196362158953
 hess_inv: None
  message: 'Optimization failed. Call limit was reached.'
   minuit: <FMin algorithm='Migrad' edm=0.00018604103640955825 edm_goal=2e-06 errordef=1.0 fval=827.5196362158953 has_accurate_covar=False has_covariance=True has_made_posdef_covar=False has_parameters_at_limit=True has_posdef_covar=True has_reached_call_limit=True has_valid_parameters=False hesse_failed=False is_above_max_edm=False is_valid=False nfcn=100392 ngrad=0 reduced_chi2=nan>
...
554527922071 x335=7.865288459575842e-05 x336=0.0003142268985029573 x337=0.00010332798675705222 x338=0.009556615091181594>
Traceback (most recent call last):
  File "/home/tcarli/pyhftest/pyhf-env/lib/python3.8/site-packages/pyhf/optimize/mixins.py", line 49, in _internal_minimize
    assert result.success
AssertionError
Traceback (most recent call last):
  File "/home/tcarli/pyhftest/pyhf-env/lib/python3.8/site-packages/pyhf/optimize/mixins.py", line 49, in _internal_minimize
    assert result.success
AssertionError

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "test.py", line 55, in <module>
    result = pyhf.infer.mle.fit(data, model, return_uncertainties=True)
  File "/home/tcarli/pyhftest/pyhf-env/lib/python3.8/site-packages/pyhf/infer/mle.py", line 131, in fit
    return opt.minimize(
  File "/home/tcarli/pyhftest/pyhf-env/lib/python3.8/site-packages/pyhf/optimize/mixins.py", line 160, in minimize
    result = self._internal_minimize(**minimizer_kwargs, options=kwargs)
  File "/home/tcarli/pyhftest/pyhf-env/lib/python3.8/site-packages/pyhf/optimize/mixins.py", line 52, in _internal_minimize
    raise exceptions.FailedMinimization(result)

This might not be a bug, but rather a problem in the set-up. May be there is some better setting for the minimize ?

Steps to Reproduce

You just need to do $ python test.py and download my workspace ex2.txt

I use version: pyhf, version 0.6.2

Regards, Tancredi

matthewfeickert commented 2 years ago

@tancredicarli Can you please provide us with the full contents of your test.py as well as your virtual environments contents (so the output of python -m pip freeze)? It seems that you're using the iminuit optimizer but we'll need to see what your code is to give any advice.

alexander-held commented 2 years ago

Hi @tancredicarli, the default number of allowed iterations in the minimization is 100k. If that is reached, the above error occurs. You can pass a maxiter argument to pyhf.infer.mle.fit to increase the allowed number of iterations:

pyhf.infer.mle.fit(data, model, maxiter=1000000)

When setting the minimizer, you can also set the verbosity to see more MINUIT output, including updates for the number of calls used:

pyhf.set_backend("numpy", pyhf.optimize.minuit_optimizer(verbose=2))
tancredicarli commented 2 years ago

Hello, thank you for the fast answer. I did not notice that I did not included the test.py. Here it is:

if __name__ == "__main__":  

 filename='ex2.txt'
 print ('open filename= ',filename)
 with open(filename) as serialized:
  spec = json.load(serialized)

 workspace = pyhf.Workspace(spec)
 model = workspace.model()
 observations=workspace.data(model)
 init_pars = model.config.suggested_init()

 print ('Calculate the observed and expected limit: ')
 mu_test = 1.0
 CLs_obs, CLs_exp = pyhf.infer.hypotest(
    mu_test,  # null hypothesis
    #observations+ model.config.auxdata,
    observations,    
    model,
    init_pars,
    test_stat="qtilde",
    #qtilde=True,
    return_expected_set=True   
 )
 print(f"      Observed CLs: {CLs_obs:.4f}")
 for expected_value, n_sigma in zip(CLs_exp, np.arange(-2, 3)):
  print(f"Expected CLs({n_sigma:2d} σ): {expected_value:.4f}")

 model = workspace.model(poi_name=None)  # background-only!
 data = workspace.data(model)

 #pyhf.set_backend("numpy", "minuit")
 pyhf.set_backend(pyhf.tensorlib, pyhf.optimize.minuit_optimizer(tolerance=1e-3)) 
 pyhf.optimizer.tolerance

 result = pyhf.infer.mle.fit(data, model, return_uncertainties=True,maxiter=1000000)
 print('result= ',result)

exit()

My code is based on: https://pyhf.readthedocs.io/en/v0.3.3/_generated/pyhf.pdf.Model.html?highlight=model#pyhf.pdf.Model

I now included the suggestions of @alexander-held (see new output below) @matthewfeickert Can you give me more details how I do the python -m pip freeze ? I am not familiar with settings python environments. I have nothing special. I just followed the tutorials. By the way: Is there an easy way to attache *.py files ?

Regards, Tancredi

P.S.: The new output with maxiter=1000000

python test.py 
open filename=  ex2.txt
Calculate the observed and expected limit: 
/home/tcarli/pyhftest/pyhf-env/lib/python3.8/site-packages/pyhf/tensor/numpy_backend.py:353: RuntimeWarning: invalid value encountered in log
  return n * np.log(lam) - lam - gammaln(n + 1.0)
      Observed CLs: 0.0009
Expected CLs(-2 σ): 0.0000
Expected CLs(-1 σ): 0.0004
Expected CLs( 0 σ): 0.0043
Expected CLs( 1 σ): 0.0380
Expected CLs( 2 σ): 0.2017
/home/tcarli/pyhftest/pyhf-env/lib/python3.8/site-packages/pyhf/interpolators/code4p.py:60: RuntimeWarning: invalid value encountered in greater
  alphasets > 1, self.mask_on, self.mask_off
/home/tcarli/pyhftest/pyhf-env/lib/python3.8/site-packages/pyhf/interpolators/code4p.py:64: RuntimeWarning: invalid value encountered in less
  alphasets < -1, self.mask_on, self.mask_off
     corr: [[ 1.00000000e+00 -6.67016779e-05  1.03798186e-05 ...  5.06273880e-05
   4.36985866e-05  8.71480409e-05]
 [-6.67016779e-05  1.00000000e+00  1.18333310e-02 ...  4.28055962e-03
   3.93097455e-03  1.24747483e-02]
 [ 1.03798186e-05  1.18333310e-02  1.00000000e+00 ... -5.11578528e-03
  -4.72067943e-03 -1.50296118e-02]
 ...[ 5.12354745e-09  3.93725562e-05 -1.05791456e-05 ...  2.54216202e-02
  -4.25475123e-05 -6.10509665e-05]
 [ 4.31473915e-09  3.52772894e-05 -9.52455539e-06 ... -4.25475123e-05
   2.41995369e-02 -5.50561014e-05]
 [ 3.92689321e-09  5.10894322e-05 -1.38386037e-05 ... -6.10509665e-05
  -5.50561014e-05  5.03982748e-03]]
  message: 'Optimization terminated successfully.'
   minuit: <FMin algorithm='Migrad' edm=0.00011951129548291648 edm_goal=2e-06 errordef=1.0 fval=827.515715667355 has_accurate_covar=False has_covariance=True has_made_posdef_covar=True has_parameters_at_limit=True has_posdef_covar=False has_reached_call_limit=False has_valid_parameters=False hesse_failed=False is_above_max_edm=True is_valid=False nfcn=506546 ngrad=0 reduced_chi2=nan>
(Param(number=0, name='x0', value=5.110156913547053e-07, error=0.19580235187648984, merror=None, is_const=False, is_fixed=False, lower_limit=0.0, upper_limit=10.0), Param(number=1, name='x1', value=-0.01872309165651593, error=0.057687512416157674, merror=None, is_const=False, is_fixed=False, lower_limit=-5.0, upper_limit=5.0), Param(number=2, name='x2', value=0.044267184801822256, error=0.012969886841663922, merror=None, is_const=False, is_fixed=False, lower_limit=-5.0, upper_limit=5.0), Param(number=3, name='x3', value=-0.01976563421454359, error=0.029824702532300615, 
...
78822224 x332=0.00012842188552101988 x333=-0.000983523601681538 x334=-0.0004962854111166893 x335=5.8972914359304435e-05 x336=0.00044492138759491496 x337=0.00022334046093944312 x338=0.009691767651235163>
Traceback (most recent call last):
  File "/home/tcarli/pyhftest/pyhf-env/lib/python3.8/site-packages/pyhf/optimize/mixins.py", line 49, in _internal_minimize
    assert result.success
AssertionError
Traceback (most recent call last):
  File "/home/tcarli/pyhftest/pyhf-env/lib/python3.8/site-packages/pyhf/optimize/mixins.py", line 49, in _internal_minimize
    assert result.success
AssertionError

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "test.py", line 55, in <module>
    result = pyhf.infer.mle.fit(data, model, return_uncertainties=True,maxiter=1000000)
  File "/home/tcarli/pyhftest/pyhf-env/lib/python3.8/site-packages/pyhf/infer/mle.py", line 131, in fit
    return opt.minimize(
  File "/home/tcarli/pyhftest/pyhf-env/lib/python3.8/site-packages/pyhf/optimize/mixins.py", line 160, in minimize
    result = self._internal_minimize(**minimizer_kwargs, options=kwargs)
  File "/home/tcarli/pyhftest/pyhf-env/lib/python3.8/site-packages/pyhf/optimize/mixins.py", line 52, in _internal_minimize
    raise exceptions.FailedMinimization(result)
pyhf.exceptions.FailedMinimization: Optimization terminated successfully.
tancredicarli commented 2 years ago

Here is the output of

python -m pip freeze
absl-py==0.13.0
altair==4.1.0
anyio==3.3.0
argon2-cffi==20.1.0
astunparse==1.6.3
attrs==21.2.0
awkward0==0.15.5
Babel==2.9.1
backcall==0.2.0
bleach==4.0.0
cachetools==4.2.2
certifi==2021.5.30
cffi==1.14.6
charset-normalizer==2.0.4
click==7.1.2
cloudpickle==1.3.0
cycler==0.10.0
debugpy==1.4.1
decorator==5.0.9
defusedxml==0.7.1
entrypoints==0.3
gast==0.3.3
google-auth==1.35.0
google-auth-oauthlib==0.4.5
google-pasta==0.2.0
grpcio==1.39.0
h5py==2.10.0
idna==3.2
iminuit==2.8.2
ipykernel==6.2.0
ipython==7.26.0
ipython-genutils==0.2.0
ipywidgets==7.6.3
jedi==0.18.0
Jinja2==3.0.1
json5==0.9.6
jsonpatch==1.32
jsonpointer==2.1
jsonschema==3.2.0
jupyter==1.0.0
jupyter-client==6.1.12
jupyter-console==6.4.0
jupyter-core==4.7.1
jupyter-server==1.10.2
jupyterlab==3.1.7
jupyterlab-pygments==0.1.2
jupyterlab-server==2.7.1
jupyterlab-widgets==1.0.0
Keras-Preprocessing==1.1.2
kiwisolver==1.3.1
Markdown==3.3.4
MarkupSafe==2.0.1
matplotlib==3.4.3
matplotlib-inline==0.1.2
mistune==0.8.4
nbclassic==0.3.1
nbclient==0.5.4
nbconvert==6.1.0
nbformat==5.1.3
nest-asyncio==1.5.1
notebook==6.4.3
numpy==1.18.5
oauthlib==3.1.1
opt-einsum==3.3.0
packaging==21.0
pandas==1.3.2
pandocfilters==1.4.3
parso==0.8.2
pexpect==4.8.0
pickleshare==0.7.5
Pillow==8.3.1
pkg_resources==0.0.0
prometheus-client==0.11.0
prompt-toolkit==3.0.19
protobuf==3.17.3
ptyprocess==0.7.0
pyasn1==0.4.8
pyasn1-modules==0.2.8
pycparser==2.20
Pygments==2.10.0
pyhf==0.6.2
pyparsing==2.4.7
PyQt5==5.15.4
PyQt5-Qt5==5.15.2
PyQt5-sip==12.9.0
pyrsistent==0.18.0
python-dateutil==2.8.2
pytz==2021.1
PyYAML==5.4.1
pyzmq==22.2.1
qtconsole==5.1.1
QtPy==1.10.0
requests==2.26.0
requests-oauthlib==1.3.0
requests-unixsocket==0.2.0
rsa==4.7.2
scipy==1.7.1
Send2Trash==1.8.0
six==1.16.0
sniffio==1.2.0
tensorboard==2.2.2
tensorboard-plugin-wit==1.8.0
tensorflow==2.2.3
tensorflow-estimator==2.2.0
tensorflow-probability==0.10.1
termcolor==1.1.0
terminado==0.11.0
testpath==0.5.0
toolz==0.11.1
torch==1.9.0
tornado==6.1
tqdm==4.62.1
traitlets==5.0.5
typing-extensions==3.10.0.0
uproot==4.0.11
uproot3==3.14.4
uproot3-methods==0.10.1
urllib3==1.26.6
wcwidth==0.2.5
webencodings==0.5.1
websocket-client==1.2.1
Werkzeug==2.0.1
widgetsnbextension==3.5.1
wrapt==1.12.1
alexander-held commented 2 years ago

Hi @tancredicarli, I made a bit of progress debugging this setup and have an idea for what may be going wrong.

The original code uses

model = workspace.model(poi_name=None)  # background-only!

which however does not do what is intended according to this comment. The poi_name=None means that no normalization factor is treated as POI (relevant e.g. for limit calculations), but it does not keep the parameter defined as POI constant. For a background-only fit, the following can be used in the measurement section of the workspace specification:

"measurements": [
    {
        "config": {
            "parameters": [{"bounds": [[-5, 5 ]], "fixed": true, "name": "mu", "inits": [0.0]}],
            "poi": "mu"
        },
        "name": "Measurement"
    }
]

This explicitly sets the parameter mu to 0.0 and keeps it constant. That addresses one of the things flagged in the errors of this fit (a non-constant parameter that has a best-fit value close to its allowed boundary), but the fit will still fail regardless, one reason is a non-positive-definite covariance matrix.

Something also seems to go wrong with the final error message:

pyhf.exceptions.FailedMinimization: Optimization terminated successfully.

suggests that the process is successful, but the boolean for success gets set to false. This may not be a pyhf issue but possibly come from iminuit, but I have not yet tracked it down.

I believe there may be a misunderstanding in the definition of the histosys template histograms. The numbers used as values in hi_data and lo_data for these modifiers are the bin yields when applying the up/down systematic variation. That means that a systematic with no effect at all should have the same numbers there as it has for the nominal yield (under the data field for the corresponding sample). The presence of many modifiers with bins where the yield is 0.0 looks to me as if the histograms were built to express the difference to nominal. Other histograms in histosys variations that have negative yields in the bins also point to that. I would suggest updating the workspace to use the absolute yield under systematic variations in the hi_data and lo_data fields (effectively adding the nominal sample yields to all histosys modifier histograms) and believe that may resolve the issues observed here.

By the way: Is there an easy way to attache *.py files ?

I am not aware of one. It is possible to drag files into messages and they get uploaded, but this seems not supported for *.py. One possibility would be changing the file ending to .txt and then dragging it into the comment field.

kratsg commented 2 years ago

You can use https://gist.github.com/ for those files and link us to your gist of files :)

tancredicarli commented 2 years ago

Hello @alexander-held

thank you very for your help. There was indeed some misunderstanding from my side about the meaning of the modifiers. I corrected this exupdatederrorconvention.txt However, I still have a problem (output testoutputexupdatederrorconvention.txt). May be it related to my extreme set-up.I will investigate a bit more where the Chi=nan might come from.

For your second suggestion: Could you help me to set-up the fit correctly fixing mu=0 ? I followed your suggestion in exupdatedpoi.txt

With this work-space I get:

Calculate the observed and expected limit: 
Traceback (most recent call last):
  `File "test.py", line 29, in <module>
    CLs_obs, CLs_exp = pyhf.infer.hypotest(
  File "/home/tcarli/pyhftest/pyhf-env/lib/python3.8/site-packages/pyhf/infer/__init__.py", line 147, in hypotest
    _check_hypotest_prerequisites(pdf, data, init_pars, par_bounds, fixed_params)
  File "/home/tcarli/pyhftest/pyhf-env/lib/python3.8/site-packages/pyhf/infer/__init__.py", line 15, in _check_hypotest_prerequisites
    raise exceptions.InvalidModel(
pyhf.exceptions.InvalidModel: POI at index [0] is set as fixed, which makes inference impossible. Please unfix the POI to continue.

May be have now some inconsistency with the function call ?

The output of python test.py with updated workspace exupdatederrorconvention.txt: testoutputexupdatederrorconvention.txt The updated work-spaces: 1) Updated errors: Modifiers are now varied cross-sections not absolute errors: exupdatederrorconvention.txt

2) Update the errors and treatment of the parameter of interest according to Alex suggestions: exupdatedpoi.txt

Thank you for all the help. Tancredi