Skip to content
Permalink
Browse files

Updates the backend to get ready for full refactor (#287)

* extend tensor backend to include new methods among them

gather: index tensor with another tensor
*_logpdf: log pdf methods for both poisson and normal pdfs
  • Loading branch information...
kratsg authored and lukasheinrich committed Oct 1, 2018
1 parent 7486cac commit fdcfb976dc5287590279b0ece7fc700b8208a337
@@ -11,7 +11,7 @@ install:
- pip freeze
script:
- pyflakes pyhf
- pytest -r sx --ignore binder/ --ignore tests/benchmarks/ --ignore tests/test_notebooks.py
- pytest -r sx --ignore tests/benchmarks/ --ignore tests/test_notebooks.py
after_success: coveralls

# always test (on both 'push' and 'pr' builds in Travis)
@@ -21,7 +21,7 @@ def extract_error(h):
Returns:
list: The uncertainty for each bin in the histogram
"""
err = h.variances if h.variances else h.numpy[0]
err = h.variances if h.variances else h.numpy()[0]
return np.sqrt(err).tolist()

def import_root_histogram(rootdir, filename, path, name):
@@ -37,7 +37,7 @@ def import_root_histogram(rootdir, filename, path, name):
h = f[os.path.join(path, name)]
except KeyError:
raise KeyError('Both {0:s} and {1:s} were tried and not found in {2:s}'.format(name, os.path.join(path, name), os.path.join(rootdir, filename)))
return h.numpy[0].tolist(), extract_error(h)
return h.numpy()[0].tolist(), extract_error(h)

def process_sample(sample,rootdir,inputfile, histopath, channelname, track_progress=False):
if 'InputFile' in sample.attrib:
@@ -1,4 +1,5 @@
from mxnet import nd

import logging
import math # Required for normal()
from numbers import Number # Required for normal()
@@ -77,7 +78,10 @@ def outer(self, tensor_in_1, tensor_in_2):
tensor_in_2.reshape((1, rows2, 1, cols2))),
(rows1 * cols1, rows2 * cols2))

def astensor(self, tensor_in):
def gather(self,tensor,indices):
return tensor[indices]

def astensor(self, tensor_in, dtype = 'float'):
"""
Convert to a MXNet NDArray.
@@ -87,10 +91,12 @@ def astensor(self, tensor_in):
Returns:
MXNet NDArray: A multi-dimensional, fixed-size homogenous array.
"""
dtypemap = {'float': 'float32', 'int': 'int32', 'bool': 'uint8'}
dtype = dtypemap[dtype]
try:
tensor = nd.array(tensor_in)
tensor = nd.array(tensor_in, dtype = dtype)
except ValueError:
tensor = nd.array([tensor_in])
tensor = nd.array([tensor_in], dtype = dtype)
return tensor

def sum(self, tensor_in, axis=None):
@@ -324,6 +330,15 @@ def simple_broadcast(self, *args):
for arg in args]
return nd.stack(*broadcast)

def shape(self, tensor):
"""
NB: Returns a tuple of longs
"""
return tensor.shape

def reshape(self, tensor, newshape):
return nd.reshape(tensor,newshape)

def einsum(self, subscripts, *operands):
"""
A generalized contraction between tensors of arbitrary dimension.
@@ -340,6 +355,11 @@ def einsum(self, subscripts, *operands):
raise NotImplementedError("mxnet::einsum is not implemented.")
return self.astensor([])

def poisson_logpdf(self, n, lam):
n = self.astensor(n)
lam = self.astensor(lam)
return n * nd.log(lam) - lam - nd.gammaln(n + 1.)

def poisson(self, n, lam):
r"""
The continous approximation, using :math:`n! = \Gamma\left(n+1\right)`,
@@ -372,6 +392,14 @@ def poisson(self, n, lam):
# https://github.com/pytorch/pytorch/blob/39520ffec15ab7e97691fed048de1832e83785e8/torch/distributions/poisson.py#L59-L63
return nd.exp((nd.log(lam) * n) - lam - nd.gammaln(n + 1.))

def normal_logpdf(self, x, mu, sigma):
# This is currently copied directly from PyTorch's source until a better
# way can be found to do this in MXNet
# https://github.com/pytorch/pytorch/blob/39520ffec15ab7e97691fed048de1832e83785e8/torch/distributions/normal.py#L70-L76
variance = sigma ** 2
log_scale = math.log(sigma) if isinstance(sigma, Number) else sigma.log()
return -((x - mu) ** 2) / (2 * variance) - log_scale - math.log(math.sqrt(2 * math.pi))

def normal(self, x, mu, sigma):
r"""
The probability density function of the Normal distribution evaluated
@@ -402,12 +430,7 @@ def normal(self, x, mu, sigma):
# This is currently copied directly from PyTorch's source until a better
# way can be found to do this in MXNet
# https://github.com/pytorch/pytorch/blob/39520ffec15ab7e97691fed048de1832e83785e8/torch/distributions/normal.py#L70-L76
def log_prob(value, loc, scale):
variance = scale ** 2
log_scale = math.log(scale) if isinstance(
scale, Number) else scale.log()
return -((value - loc) ** 2) / (2 * variance) - log_scale - math.log(math.sqrt(2 * math.pi))
return self.exp(log_prob(x, mu, sigma))
return self.exp(self.normal_logpdf(x, mu, sigma))

def normal_cdf(self, x, mu=0, sigma=1):
"""
@@ -21,7 +21,7 @@ def clip(self, tensor_in, min, max):
>>> pyhf.set_backend(pyhf.tensor.numpy_backend())
>>> a = pyhf.tensorlib.astensor([-2, -1, 0, 1, 2])
>>> pyhf.tensorlib.clip(a, -1, 1)
array([-1, -1, 0, 1, 1])
array([-1., -1., 0., 1., 1.])
Args:
tensor_in (`tensor`): The input tensor object
@@ -31,19 +31,23 @@ def clip(self, tensor_in, min, max):
Returns:
NumPy ndarray: A clipped `tensor`
"""
tensor_in = self.astensor(tensor_in)
return np.clip(tensor_in, min, max)

def tolist(self,tensor_in):
tensor_in = self.astensor(tensor_in)
return tensor_in.tolist()

def outer(self, tensor_in_1, tensor_in_2):
tensor_in_1 = self.astensor(tensor_in_1)
tensor_in_2 = self.astensor(tensor_in_2)
return np.outer(tensor_in_1,tensor_in_2)

def astensor(self, tensor_in):
def gather(self,tensor,indices):
return tensor[indices]

def isfinite(self, tensor):
return np.isfinite(tensor)

def astensor(self, tensor_in, dtype = 'float'):
"""
Convert to a NumPy array.
@@ -53,18 +57,17 @@ def astensor(self, tensor_in):
Returns:
`numpy.ndarray`: A multi-dimensional, fixed-size homogenous array.
"""
return np.asarray(tensor_in)
dtypemap = {'float': np.float64, 'int': np.int64, 'bool': np.bool_}
dtype = dtypemap[dtype]
return np.asarray(tensor_in, dtype = dtype)

def sum(self, tensor_in, axis=None):
tensor_in = self.astensor(tensor_in)
return np.sum(tensor_in, axis=axis)

def product(self, tensor_in, axis=None):
tensor_in = self.astensor(tensor_in)
return np.product(tensor_in, axis = axis)

def abs(self, tensor):
tensor = self.astensor(tensor)
return np.abs(tensor)

def ones(self,shape):
@@ -74,34 +77,24 @@ def zeros(self,shape):
return np.zeros(shape)

def power(self,tensor_in_1, tensor_in_2):
tensor_in_1 = self.astensor(tensor_in_1)
tensor_in_2 = self.astensor(tensor_in_2)
return np.power(tensor_in_1, tensor_in_2)

def sqrt(self,tensor_in):
tensor_in = self.astensor(tensor_in)
return np.sqrt(tensor_in)

def divide(self,tensor_in_1, tensor_in_2):
tensor_in_1 = self.astensor(tensor_in_1)
tensor_in_2 = self.astensor(tensor_in_2)
return np.divide(tensor_in_1, tensor_in_2)

def log(self,tensor_in):
tensor_in = self.astensor(tensor_in)
return np.log(tensor_in)

def exp(self,tensor_in):
tensor_in = self.astensor(tensor_in)
return np.exp(tensor_in)

def stack(self, sequence, axis = 0):
return np.stack(sequence,axis = axis)

def where(self, mask, tensor_in_1, tensor_in_2):
mask = self.astensor(mask)
tensor_in_1 = self.astensor(tensor_in_1)
tensor_in_2 = self.astensor(tensor_in_2)
return np.where(mask, tensor_in_1, tensor_in_2)

def concatenate(self, sequence, axis=0):
@@ -130,7 +123,7 @@ def simple_broadcast(self, *args):
... pyhf.tensorlib.astensor([1]),
... pyhf.tensorlib.astensor([2, 3, 4]),
... pyhf.tensorlib.astensor([5, 6, 7]))
[array([1, 1, 1]), array([2, 3, 4]), array([5, 6, 7])]
[array([1., 1., 1.]), array([2., 3., 4.]), array([5., 6., 7.])]
Args:
args (Array of Tensors): Sequence of arrays
@@ -140,6 +133,12 @@ def simple_broadcast(self, *args):
"""
return np.broadcast_arrays(*args)

def shape(self, tensor):
return tensor.shape

def reshape(self, tensor, newshape):
return np.reshape(tensor,newshape)

def einsum(self, subscripts, *operands):
"""
Evaluates the Einstein summation convention on the operands.
@@ -159,6 +158,11 @@ def einsum(self, subscripts, *operands):
"""
return np.einsum(subscripts, *operands)

def poisson_logpdf(self, n, lam):
n = np.asarray(n)
lam = np.asarray(lam)
return n * np.log(lam) - lam - gammaln(n + 1.)

def poisson(self, n, lam):
r"""
The continous approximation, using :math:`n! = \Gamma\left(n+1\right)`,
@@ -185,6 +189,19 @@ def poisson(self, n, lam):
lam = np.asarray(lam)
return np.exp(n * np.log(lam) - lam - gammaln(n + 1.))

def normal_logpdf(self, x, mu, sigma):
# this is much faster than
# norm.logpdf(x, loc=mu, scale=sigma)
# https://codereview.stackexchange.com/questions/69718/fastest-computation-of-n-likelihoods-on-normal-distributions
root2 = np.sqrt(2)
root2pi = np.sqrt(2*np.pi)
prefactor = -np.log(sigma * root2pi)
summand = -np.square(np.divide((x - mu),(root2 * sigma)))
return prefactor + summand

# def normal_logpdf(self, x, mu, sigma):
# return norm.logpdf(x, loc=mu, scale=sigma)

def normal(self, x, mu, sigma):
r"""
The probability density function of the Normal distribution evaluated
@@ -42,7 +42,7 @@ def outer(self, tensor_in_1, tensor_in_2):
tensor_in_2 = self.astensor(tensor_in_2)
return torch.ger(tensor_in_1,tensor_in_2)

def astensor(self, tensor_in):
def astensor(self, tensor_in, dtype = 'float'):
"""
Convert to a PyTorch Tensor.
@@ -52,13 +52,24 @@ def astensor(self, tensor_in):
Returns:
torch.Tensor: A multi-dimensional matrix containing elements of a single data type.
"""
dtypemap = {'float': torch.float, 'int': torch.int, 'bool': torch.uint8}
dtype = dtypemap[dtype]
if isinstance(tensor_in, torch.Tensor):
v = tensor_in
else:
if not isinstance(tensor_in, list):
tensor_in = [tensor_in]
v = torch.Tensor(tensor_in)
return v.type(torch.FloatTensor)
v = torch.tensor(tensor_in, dtype = dtype)
return v

def gather(self,tensor,indices):
return torch.take(tensor,indices.type(torch.LongTensor))

def reshape(self, tensor, newshape):
return torch.reshape(tensor,newshape)

def shape(self, tensor):
return tuple(map(int,tensor.shape))

def sum(self, tensor_in, axis=None):
tensor_in = self.astensor(tensor_in)
@@ -104,7 +115,7 @@ def stack(self, sequence, axis = 0):
return torch.stack(sequence,dim = axis)

def where(self, mask, tensor_in_1, tensor_in_2):
mask = self.astensor(mask)
mask = self.astensor(mask).type(torch.FloatTensor)
tensor_in_1 = self.astensor(tensor_in_1)
tensor_in_2 = self.astensor(tensor_in_2)
return mask * tensor_in_1 + (1-mask) * tensor_in_2
@@ -123,6 +134,9 @@ def concatenate(self, sequence, axis=0):
"""
return torch.cat(sequence, dim=axis)

def isfinite(self, tensor):
return torch.isfinite(tensor)

def simple_broadcast(self, *args):
"""
Broadcast a sequence of 1 dimensional arrays.
@@ -170,6 +184,11 @@ def einsum(self, subscripts, *operands):
ops = tuple(self.astensor(op) for op in operands)
return torch.einsum(subscripts, ops)

def poisson_logpdf(self, n, lam):
n = self.astensor(n)
lam = self.astensor(lam)
return torch.distributions.Poisson(lam).log_prob(n)

def poisson(self, n, lam):
r"""
The continous approximation, using :math:`n! = \Gamma\left(n+1\right)`,
@@ -196,6 +215,13 @@ def poisson(self, n, lam):
lam = self.astensor(lam)
return torch.exp(torch.distributions.Poisson(lam).log_prob(n))

def normal_logpdf(self, x, mu, sigma):
x = self.astensor(x)
mu = self.astensor(mu)
sigma = self.astensor(sigma)
normal = torch.distributions.Normal(mu, sigma)
return normal.log_prob(x)

def normal(self, x, mu, sigma):
r"""
The probability density function of the Normal distribution evaluated

0 comments on commit fdcfb97

Please sign in to comment.
You can’t perform that action at this time.