scikit-hep / pyhf

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

simple_broadcast is broken in numpy? #158

Closed kratsg closed 6 years ago

kratsg commented 6 years ago

Description

simple_broadcast does not work the same in different backends

Expected Behavior

Arrays got broadcasted correctly in numpy

Actual Behavior

Arrays did not get broadcast correctly in numpy

Steps to Reproduce

>>> import pyhf
>>> pyhf.set_backend(pyhf.tensor.mxnet_backend())
>>> pyhf.tensorlib.simple_broadcast(pyhf.tensorlib.astensor([1]), pyhf.tensorlib.astensor([2, 2]), pyhf.tensorlib.astensor([3, 3, 3]))

[[1. 1. 1.]
 [2. 2. 2.]
 [3. 3. 3.]]
<NDArray 3x3 @cpu(0)>
>>> import pyhf
>>> pyhf.set_backend(pyhf.tensor.numpy_backend())
>>> pyhf.tensorlib.simple_broadcast(pyhf.tensorlib.astensor([1]), pyhf.tensorlib.astensor([2, 2]), pyhf.tensorlib.astensor([3, 3, 3]))
>>> pyhf.tensorlib.simple_broadcast(pyhf.tensorlib.astensor([1]), pyhf.tensorlib.astensor([2, 2]), pyhf.tensorlib.astensor([3, 3, 3]))
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/Users/kratsg/pyhf/pyhf/tensor/numpy_backend.py", line 90, in simple_broadcast
    return np.broadcast_arrays(*args)
  File "/Users/kratsg/.virtualenvs/pyhf/lib/python2.7/site-packages/numpy/lib/stride_tricks.py", line 249, in broadcast_arrays
    shape = _broadcast_shape(*args)
  File "/Users/kratsg/.virtualenvs/pyhf/lib/python2.7/site-packages/numpy/lib/stride_tricks.py", line 184, in _broadcast_shape
    b = np.broadcast(*args[:32])
ValueError: shape mismatch: objects cannot be broadcast to a single shape

Checklist

matthewfeickert commented 6 years ago

I think that the MXNet backend is the only one that does this correctly at the moment as the other backends only broadcast properly if a single entry array is given (when I was writing the MXNet backend I remember possibly noticing this but I guess I never went and checked the other backends until now).

>>> import pyhf
>>> import tensorflow as tf
>>> pyhf.set_backend(pyhf.tensor.tensorflow_backend(session=tf.Session()))
>>> a = pyhf.tensorlib.simple_broadcast(pyhf.tensorlib.astensor([1]), pyhf.tensorlib.astensor([2, 2]), pyhf.tensorlib.astensor([3, 3, 3]))
>>> tf.Session().run(a)
[array([1., 1., 1.], dtype=float32), array([2., 2.], dtype=float32), array([3., 3., 3.], dtype=float32)]
>>> b = pyhf.tensorlib.simple_broadcast(pyhf.tensorlib.astensor([1]), pyhf.tensorlib.astensor([2]), pyhf.tensorlib.astensor([3, 3, 3]))
>>> tf.Session().run(b)
[array([1., 1., 1.], dtype=float32), array([2., 2., 2.], dtype=float32), array([3., 3., 3.], dtype=float32)]
>>>
>>> pyhf.set_backend(pyhf.tensor.pytorch_backend())
>>> pyhf.tensorlib.simple_broadcast(pyhf.tensorlib.astensor([1]), pyhf.tensorlib.astensor([2, 2]), pyhf.tensorlib.astensor([3, 3, 3]))
[tensor([ 1.,  1.,  1.]), tensor([ 2.,  2.]), tensor([ 3.,  3.,  3.])]
>>> pyhf.tensorlib.simple_broadcast(pyhf.tensorlib.astensor([1]), pyhf.tensorlib.astensor([2]), pyhf.tensorlib.astensor([3, 3, 3]))
[tensor([ 1.,  1.,  1.]), tensor([ 2.,  2.,  2.]), tensor([ 3.,  3.,  3.])]

I can look into fixing this unless @kratsg has already started.

kratsg commented 6 years ago

I can look into fixing this unless @kratsg has already started.

I haven't. Just reporting it for now.

lukasheinrich commented 6 years ago

Hi, thanks for the report. I think the canonical one is the one by numpy (since there, the tensorlib call is a simple shim around np.broadcast_arrays, which unfortunately is not present in other libraries). The full broadcasting rules are specified in https://docs.scipy.org/doc/numpy-1.13.0/user/basics.broadcasting.html but we don't need those (yet). THe only type of broadcasting we need is the scenario of taking a set of (1,)- or (N,)-shape arrays and normalizing them to (N,)-shape arrays.