tpgillam/mt2

Stransverse mass computation, compatible with numpy.

mt2's Issues

Build wheels

Submitting wheels to PyPI is required to make life easier for users, when building from source might be difficult (or impossible, e.g. due to missing python header files).

numpy version mismatch?

Failure here: https://travis-ci.org/github/tpgillam/mt2/jobs/758603656

ImportError while importing test module '/home/travis/build/tpgillam/mt2/tests/test_mt2.py'.
Hint: make sure your test modules/packages have valid Python names.
Traceback:
tests/test_mt2.py:7: in <module>
    from mt2 import mt2
mt2/__init__.py:6: in <module>
    from _mt2 import mt2_ufunc
E   ImportError: numpy.core.multiarray failed to import
------------------------------- Captured stderr --------------------------------
RuntimeError: module compiled against API version 0xe but this version of numpy is 0xd

Use MANIFEST.in to include header files

This seems to be the non-disgusting way to do it, e.g. https://packaging.python.org/guides/using-manifest-in/

Decide how to approach absolute tolerance

Decide on an approach to absolute tolerance in bisection termination.

mt2_bisect_impl has a hard coded absolute tolerance of the type epsilon (ie min x st 1 + x > 1).

This would be unhelpful for exa-electron volt users (ANITA?).

Setting it to 0 might (as far as I know) cause infinite loops when mt2 is denormal.

I suggest the following change:

    auto epsilon = std::numeric_limits<T>::epsilon();
    auto rel_tol = epsilon < precision ? precision : epsilon;
-    auto abs_tol = epsilon;
+    auto abs_tol = std::numeric_limits<T>::min();

Input validation

Can mt2_bisect_impl reach bad states from bad inputs?

Should more input validation be added?

e.g. large negative masses could give us negative bounds here.

Add more complete examples

We should add a few more examples, ideally with a nice plot or two, demonstrating broadcasting.

Add python wrapper for Rupert's bisector

Linux build fail:

• mt2 version: 0.1.4
• Python version: 3.6.8
• Operating System: CentOS Linux 7 (Core)

Description

(tmp4) pclu--bash$ python --version
Python 3.6.8

(tmp4) pclu--bash$ uname -a
Linux pclu 3.10.0-1160.6.1.el7.x86_64 #1 SMP Tue Nov 17 13:59:11 UTC 2020 x86_64 x86_64 x86_64 GNU/Linux

(tmp4) pclu--bash$ cat /etc/os-release 
NAME="CentOS Linux"
VERSION="7 (Core)"
ID="centos"
ID_LIKE="rhel fedora"
VERSION_ID="7"
PRETTY_NAME="CentOS Linux 7 (Core)"
ANSI_COLOR="0;31"
CPE_NAME="cpe:/o:centos:centos:7"
HOME_URL="http://cern.ch/linux/"
BUG_REPORT_URL="http://cern.ch/linux/"
CENTOS_MANTISBT_PROJECT="CentOS-7"
CENTOS_MANTISBT_PROJECT_VERSION="7"
REDHAT_SUPPORT_PRODUCT="centos"
REDHAT_SUPPORT_PRODUCT_VERSION="7"

(tmp4) pclu--bash$ pip list
Package      Version
------------ ----------
click        7.1.2
itsdangerous 1.1.0
Jinja2       2.11.3
MarkupSafe   1.1.1
pip          21.0.1
Python-dev   2.0.0.dev0
setuptools   39.2.0
Werkzeug     1.0.1
wheel        0.36.2

(tmp4) pclu--bash$ pwd
/usera/lester/tmp4

(tmp4) pclu--bash$ pip install mt2
Collecting mt2
  Using cached mt2-0.1.4.tar.gz (17 kB)
  Installing build dependencies ... done
  Getting requirements to build wheel ... done
  Installing backend dependencies ... done
    Preparing wheel metadata ... done
Collecting numpy
  Using cached numpy-1.19.5-cp36-cp36m-manylinux2010_x86_64.whl (14.8 MB)
Building wheels for collected packages: mt2
  Building wheel for mt2 (PEP 517) ... error
  ERROR: Command errored out with exit status 1:
   command: /usera/lester/tmp4/bin/python3 /usera/lester/tmp4/lib64/python3.6/site-packages/pip/_vendor/pep517/_in_process.py build_wheel /tmp/lester/tmpxxxygywb
       cwd: /tmp/lester/pip-install-9zxp06u5/mt2_54ceb8da895b4a8485d5ee3ba8343f01
  Complete output (23 lines):
  running bdist_wheel
  running build
  running build_ext
  gcc -pthread -Wno-unused-result -Wsign-compare -DNDEBUG -O2 -g -pipe -Wall -Wp,-D_FORTIFY_SOURCE=2 -fexceptions -fstack-protector-strong --param=ssp-buffer-size=4 -grecord-gcc-switches -m64 -mtune=generic -D_GNU_SOURCE -fPIC -fwrapv -fPIC -I/usera/lester/tmp4/include -I/usr/include/python3.6m -c flagcheck.cpp -o flagcheck.o -std=c++17
  gcc: error: unrecognized command line option '-std=c++17'
  gcc -pthread -Wno-unused-result -Wsign-compare -DNDEBUG -O2 -g -pipe -Wall -Wp,-D_FORTIFY_SOURCE=2 -fexceptions -fstack-protector-strong --param=ssp-buffer-size=4 -grecord-gcc-switches -m64 -mtune=generic -D_GNU_SOURCE -fPIC -fwrapv -fPIC -I/usera/lester/tmp4/include -I/usr/include/python3.6m -c flagcheck.cpp -o flagcheck.o -std=c++14
  gcc: error: unrecognized command line option '-std=c++14'
  gcc -pthread -Wno-unused-result -Wsign-compare -DNDEBUG -O2 -g -pipe -Wall -Wp,-D_FORTIFY_SOURCE=2 -fexceptions -fstack-protector-strong --param=ssp-buffer-size=4 -grecord-gcc-switches -m64 -mtune=generic -D_GNU_SOURCE -fPIC -fwrapv -fPIC -I/usera/lester/tmp4/include -I/usr/include/python3.6m -c flagcheck.cpp -o flagcheck.o -std=c++11
  building 'mt2' extension
  creating build
  creating build/temp.linux-x86_64-3.6
  creating build/temp.linux-x86_64-3.6/src
  gcc -pthread -Wno-unused-result -Wsign-compare -DNDEBUG -O2 -g -pipe -Wall -Wp,-D_FORTIFY_SOURCE=2 -fexceptions -fstack-protector-strong --param=ssp-buffer-size=4 -grecord-gcc-switches -m64 -mtune=generic -D_GNU_SOURCE -fPIC -fwrapv -fPIC -DVERSION_INFO=0.1.4 -I/tmp/lester/pip-build-env-zmcyf3yy/overlay/lib/python3.6/site-packages/pybind11/include -I/usera/lester/tmp4/include -I/usr/include/python3.6m -c src/main.cpp -o build/temp.linux-x86_64-3.6/src/main.o -std=c++11 -fvisibility=hidden -g0
  In file included from /tmp/lester/pip-build-env-zmcyf3yy/overlay/lib/python3.6/site-packages/pybind11/include/pybind11/pytypes.h:12:0,
                   from /tmp/lester/pip-build-env-zmcyf3yy/overlay/lib/python3.6/site-packages/pybind11/include/pybind11/cast.h:13,
                   from /tmp/lester/pip-build-env-zmcyf3yy/overlay/lib/python3.6/site-packages/pybind11/include/pybind11/attr.h:13,
                   from /tmp/lester/pip-build-env-zmcyf3yy/overlay/lib/python3.6/site-packages/pybind11/include/pybind11/pybind11.h:45,
                   from src/main.cpp:1:
  /tmp/lester/pip-build-env-zmcyf3yy/overlay/lib/python3.6/site-packages/pybind11/include/pybind11/detail/common.h:124:20: fatal error: Python.h: No such file or directory
   #include <Python.h>
                      ^
  compilation terminated.
  error: command 'gcc' failed with exit status 1
  ----------------------------------------
  ERROR: Failed building wheel for mt2
Failed to build mt2
ERROR: Could not build wheels for mt2 which use PEP 517 and cannot be installed directly

CI site not aware of numpy dependency

Show copyright message only on installation

• mt2 version: 0.1.4
• Python version: N/A
• Operating System: N/A

Description

Currently we show a copyright message for every invocation of the mt2 root finding. We would prefer it if the message was shown once on installation, and not thereafter.

This will ideally require forking the header file from the versions available on arxiv, if you don't object @kesterlester ? (Alternatively I can make do with the existing API, but there'll be a small amount of overhead for doing so that we could otherwise avoid)

On that topic, as my basis I'll be using "v4" from arxiv, since this doesn't have the various inlining optimisations removed. No ROOT shall be used here!

Move CI to github actions from Travis

Public repositories (like this one) get unlimited github actions usage. My previous experience with github actions is that stuff runs very quickly, with limited time in the queue. Travis seems to have long wait times (>1 hour), so moving would benefit its utility for checking PRs.

Tom requested this new issue

Aha! That's interesting. I think this failure mode is caused by your version of numpy being too old -- at some point the second and third arguments went from having to be long* to const long*. The latter is definitely more correct, but the documentation still used the former.

Can you create an issue for this please, ideally with the version of numpy that you're having problems with? I reckon there'll be some kind of macro switch I can insert to support both. Thanks!

Originally posted by @tpgillam in #16 (comment)

Drop python 3.5 compatibility

Python 3.5 reached end-of-life last year, and support will soon be dropped from pip. Its presence is also causing failure on travis.

Resolve Lally vs LesterNachman differences in test-suite

Reminder to self that it would be a good idea for @kesterlester to check out the events and MT2 discrepancies talked about in Appendix B of Lally's https://arxiv.org/pdf/1509.01831.pdf ....

Decide whether to fix small bug in Lally code ..

Copied here, for the record, are the code-related bits of my discussion with Lally about a potential bug in his imp:

Lally's reply:

OK, bottom line is it looks like an error. I would have to have been thinking of something overly complicated to have put an assignment like this in an If statement. So I'm pretty sure it should have been a simple comparison.

The good news (and the reason I guess why I never picked up on it) is that, as far as I can see, it has no material effect. Just to quickly explain, even though I knew I had to watch out for the Newton-Raphson finding either higher-than-lowest-value roots, or getting caught in a local minimum, I still wanted to at least be able to output what would normally be a new, lower, upperbound from this function, which could then be used in the Regula Falsi function instead of the kinematic upperbound i.e. to gain some value / speed from going through the Newton-Raphson method first, even if it didn't actually find the correct root. What this bit of the code is doing is selecting a suitable new upperbound from the last few iterations. What the error does is pretty much ALWAYS (because x1 and xNR are unlikely to be both zero valued and so the If statement will always be true) increase the second-to-last iterated value a little (rather than only increase it if it is EQUAL to the final iterated value). This has the potential (in practice very rarely) to output a slightly higher new upperbound than would have been the case without the error, but the effect on the speed of the Regula Falsi function will be de minimis as the higher upperbound is not going to be that much higher than the next lower iterated value as we are in a local minimum type situation.

Anyway, it's a good spot! I presume it was spotted more from a 'forensic' reading of my code i.e. it just 'looked wrong' given the If statement seemed to only require a simple comparison, and not because the overall algorithm was doing something wrong? If this isn't the case I would need to look in more detail, as maybe I was actually trying to do something fancy (unlikely!).

The above followed my message to him:

I have a question about your MT2 implementation from your most recent paper. My question relates to line 421. You can see line 421 by going to the link below and scrolling down far enough.

a1148d6#diff-be879a9f457b24439e1e1e360f477a90c04c7734a61de906c28589dd25f86da1

The line in question says:

and I am querying why it says:

if (xNR = x1) {

rather than

if (xNR == x1) {

Both are valid C and either could (potentially) be what you want … I’ve not yet figured out … however they have different effects. The first version assigns x1 to xNR and then passes the test if xNR is non zero, while the second version passes the test if xNR is equal to x1.

As I’ve said, it might be that you want to do the former, as written, in which case there is no bug. But if you wanted to do the latter then there is a bug.

Implementation switch?

Will Fawcett gave a short presentation in the ATLAS Susy Background forum yesterday advertising this mt2 library. There some polite interest, which is nice, with people who comments pleased to hear of (a) the vectorisation, and (b) Rupert's speedup. But there was one issue where you could hear a slight fear in the voices which concerned validation of Rupert's alg. Bascially since atlas has used the arXiv.h code for years and years and years without change, it was clear that they were nervous about using a re-implementation even if the alg was basically the same. I got the impression, though, that if there was a switch that would let them choose "fast_rupert" vs "as_on_arXiv" then people would then have no hurdle to adoption.

Minor refactoring of nice 2d example

TODO:

• use numpy
• use matplotlib
• include figure in README

Check it runs on the pyboard under micropython

Manually create a ufunc instead of using pybind11:vectorize

Currently pybind11::vectorize is a very simple way to wrap the function, however what it generates isn't actually a numpy ufunc.

Creating a ufunc by hand would be potentially nice, as it would remove the pybind11 dependency (which sadly is causing confusion and issues during building due to the newness of PEP 517/518, which I'm currently depending upon). It could also be a bit more efficient, as it would allow use of "where", "out" etc.

It may not support default arguments (?) — if not then we can declare a simple wrapper that does.

Surprising test comparison outcome

There is a line in a new test which I created which sets a tolerance in the approx-equal command. Docs say that the default tolerance is 1e-6, but the MT2 calculator should be able to do better in this test, so I pushed it up to 1e-12 and the test still passes. However if I push this tolerance up to 1e-13 then the test fails (on my laptop -- I don't know if it fails on github or elsewhere at 1e-13). This I fund surprising as the printouts suggest that python thinks the numbers it is comparing agree to at least 1e-15 or maybe even higher, EVEN WHEN IT SAYS THE COMPARISON FAILS. This means that something fishy is happening which should be understood.

Originally posted by @kesterlester in #32 (comment)

Python 3.9 support

Currently only testing up to 3.8, but at some point should verify that it works on 3.9 too.

Support both `const long*` and `long*` form of ufunc inner loop

At some point numpy has changed the API for ufunc declaration - ideally we want to support both, and also reflect the version requirements in our dependency structure.

Better bisectors?

@tpgillam and @kesterlester to think about ways of using the ideas behind Lally ( https://arxiv.org/abs/1509.01831 ) to make for a supralinear/quadratically better bisector

Make more pythonic

get_mT2 -> mt2 ?

copyright messages only in install phase?

Decide if/how/whether to implement/add Rupert's implementation

Rupert's notes:

No significant changes to the method itself but, by hoisting heavy lifting out of the hot loop, this version is significantly faster.

Running `make' in the attached directory (perhaps after source setup.sh from atlas space) should print out some timings
based on tests stolen from Tom's python module.

Time is no big issue but the ~1e-10 precision we get is a little disappointing (though a 0 better in the exponent than what Lally achieves)

I also include a version using long double for a few extra bits of accuracy.

I think the main precision loss is in that highest order "thing" in the disjointness condition. If I look at this again, I will try some tricks there!

mt2_email.zip

Switch to TLN as default

Tombs-Lester-Nachman implementation is faster than the original LN, and is scale invariant (see #44 ). Therefore it would be beneficial to switch to this method in the public-facing API.

We should retain the LN code if required for unit tests.

The Lally method should (for) probably be moved into a separate branch, so as not to include it in the next release, but to retain the ability to run experiments on it when desired.

Pip install fail

• mt2 version: 0.1.4
• Python version: 3.7.6
• Operating System: macOS Catalina, version 10.15.7

Description

Pip install fails.

What I Did

Will Fawcett  10:56 AM
okay, I removed the old version of mt2 (from the git package) and tried again, unfortunately I the pip install doesn’t seem to work:
$ pip install mt2
Collecting mt2
  Using cached mt2-0.1.4.tar.gz (17 kB)
  Installing build dependencies ... error
  ERROR: Command errored out with exit status 1:
   command: /Users/wfawcett/opt/anaconda3/bin/python /Users/wfawcett/opt/anaconda3/lib/python3.7/site-packages/pip install --ignore-installed --no-user --prefix /private/var/folders/tg/d7mypqb95lz0r1wrdbxqbs840000gn/T/pip-build-env-5ercqndn/overlay --no-warn-script-location --no-binary :none: --only-binary :none: -i https://pypi.org/simple -- 'setuptools>=42' wheel 'pybind11>=2.6.0'
       cwd: None
  Complete output (14 lines):
  Traceback (most recent call last):
    File "/Users/wfawcett/opt/anaconda3/lib/python3.7/runpy.py", line 193, in _run_module_as_main
      "__main__", mod_spec)
    File "/Users/wfawcett/opt/anaconda3/lib/python3.7/runpy.py", line 85, in _run_code
      exec(code, run_globals)
    File "/Users/wfawcett/opt/anaconda3/lib/python3.7/site-packages/pip/__main__.py", line 23, in <module>
      from pip._internal.cli.main import main as _main  # isort:skip # noqa
    File "/Users/wfawcett/opt/anaconda3/lib/python3.7/site-packages/pip/_internal/cli/main.py", line 5, in <module>
      import locale
    File "/Users/wfawcett/opt/anaconda3/lib/python3.7/locale.py", line 16, in <module>
      import re
    File "/Users/wfawcett/opt/anaconda3/lib/python3.7/re.py", line 143, in <module>
      class RegexFlag(enum.IntFlag):
  AttributeError: module 'enum' has no attribute 'IntFlag'
  ----------------------------------------
ERROR: Command errored out with exit status 1: /Users/wfawcett/opt/anaconda3/bin/python /Users/wfawcett/opt/anaconda3/lib/python3.7/site-packages/pip install --ignore-installed --no-user --prefix /private/var/folders/tg/d7mypqb95lz0r1wrdbxqbs840000gn/T/pip-build-env-5ercqndn/overlay --no-warn-script-location --no-binary :none: --only-binary :none: -i https://pypi.org/simple -- 'setuptools>=42' wheel 'pybind11>=2.6.0' Check the logs for full command output.

Lally's alg should be in the test suite.

Ask Tom nicely if he can FOR USE IN THE TEST-SUITE ONLY -- NOT IN THE PUBLIC LIBRARY -- give https://arxiv.org/src/1509.01831v2/anc/mt2_Lallyver2.h some python bindings so that we can create tests that compare our MT2 values to those it produces. This implementation need not be a fancy ufunc thing -- just anything that "works".