Skip to content

P
pytensor
提交 efc098de authored 作者: Frederic's avatar Frederic
浏览文件
操作

pep8

上级 2ef495d9
显示空白字符变更
内嵌 并排
正在显示 包含 131 行增加 和 81 行删除
theano/tensor/raw_random.py
浏览文件 @ efc098de
...@@ -2,14 +2,17 @@ ...@@ -2,14 +2,17 @@
__docformat__ = "restructuredtext en" __docformat__ = "restructuredtext en"
import sys import sys
from copy import copy from copy import copy
import numpy import numpy
#local imports #local imports
import theano
import basic as tensor import basic as tensor
import opt, theano import opt
from theano import gof from theano import gof
from theano.compile import optdb from theano.compile import optdb
class RandomStateType(gof.Type): class RandomStateType(gof.Type):
"""A Type wrapper for numpy.RandomState """A Type wrapper for numpy.RandomState
...@@ -157,8 +160,8 @@ class RandomFunction(gof.Op): ...@@ -157,8 +160,8 @@ class RandomFunction(gof.Op):
print >> sys.stderr, 'WARNING: RandomState instances should be in RandomStateType' print >> sys.stderr, 'WARNING: RandomState instances should be in RandomStateType'
if 0: if 0:
raise TypeError('r must be RandomStateType instance', r) raise TypeError('r must be RandomStateType instance', r)
# the following doesn't work because we want to ignore the broadcastable flags in # the following doesn't work because we want to ignore the
# shape.type # broadcastable flags in shape.type
# assert shape.type == tensor.lvector # assert shape.type == tensor.lvector
# convert args to TensorType instances # convert args to TensorType instances
...@@ -173,7 +176,7 @@ class RandomFunction(gof.Op): ...@@ -173,7 +176,7 @@ class RandomFunction(gof.Op):
r, shp = node.inputs[0:2] r, shp = node.inputs[0:2]
#if shp is a constant array of len 0, then it means 'automatic shape' #if shp is a constant array of len 0, then it means 'automatic shape'
unknown_shape = len(getattr(shp, 'data', [0,1,2])) == 0 unknown_shape = len(getattr(shp, 'data', [0, 1, 2])) == 0
# if ndim_added == 0 and shape != () then shape # if ndim_added == 0 and shape != () then shape
if self.ndim_added == 0 and not unknown_shape: if self.ndim_added == 0 and not unknown_shape:
...@@ -188,8 +191,8 @@ class RandomFunction(gof.Op): ...@@ -188,8 +191,8 @@ class RandomFunction(gof.Op):
def perform(self, node, inputs, out_): def perform(self, node, inputs, out_):
rout, out = out_ rout, out = out_
# Use self.fn to draw shape worth of random numbers. # Use self.fn to draw shape worth of random numbers.
# Numbers are drawn from r if self.inplace is True, and from a copy of r if # Numbers are drawn from r if self.inplace is True, and from a
# self.inplace is False # copy of r if self.inplace is False
r, shape, args = inputs[0], inputs[1], inputs[2:] r, shape, args = inputs[0], inputs[1], inputs[2:]
assert type(r) == numpy.random.RandomState, (type(r), r) assert type(r) == numpy.random.RandomState, (type(r), r)
r_orig = r r_orig = r
...@@ -203,34 +206,44 @@ class RandomFunction(gof.Op): ...@@ -203,34 +206,44 @@ class RandomFunction(gof.Op):
else: else:
shape = tuple(shape) shape = tuple(shape)
if shape is not None and self.outtype.ndim != len(shape) + self.ndim_added: if (shape is not None and
raise ValueError('Shape mismatch: self.outtype.ndim (%i) != len(shape) (%i) + self.ndim_added (%i)'\ self.outtype.ndim != len(shape) + self.ndim_added):
%(self.outtype.ndim, len(shape), self.ndim_added)) raise ValueError('Shape mismatch: self.outtype.ndim (%i) !='
' len(shape) (%i) + self.ndim_added (%i)'
% (self.outtype.ndim, len(shape), self.ndim_added))
if not self.inplace: if not self.inplace:
r = copy(r) r = copy(r)
rout[0] = r rout[0] = r
rval = self.fn(r, *(args + [shape])) rval = self.fn(r, *(args + [shape]))
if not isinstance(rval, numpy.ndarray) \ if not isinstance(rval, numpy.ndarray) \
or str(rval.dtype) != node.outputs[1].type.dtype: or str(rval.dtype) != node.outputs[1].type.dtype:
rval = theano._asarray(rval, dtype = node.outputs[1].type.dtype) rval = theano._asarray(rval, dtype=node.outputs[1].type.dtype)
# When shape is None, numpy has a tendency to unexpectedly # When shape is None, numpy has a tendency to unexpectedly
# return a scalar instead of a higher-dimension array containing # return a scalar instead of a higher-dimension array containing
# only one element. This value should be reshaped # only one element. This value should be reshaped
if shape is None and rval.ndim == 0 and self.outtype.ndim > 0: if shape is None and rval.ndim == 0 and self.outtype.ndim > 0:
rval = rval.reshape([1]*self.outtype.ndim) rval = rval.reshape([1] * self.outtype.ndim)
if len(rval.shape) != self.outtype.ndim: if len(rval.shape) != self.outtype.ndim:
raise ValueError('Shape mismatch: "out" should have dimension %i, but the value produced by "perform" has dimension %i'\ raise ValueError('Shape mismatch: "out" should have dimension %i,'
' but the value produced by "perform" has'
' dimension %i'
% (self.outtype.ndim, len(rval.shape))) % (self.outtype.ndim, len(rval.shape)))
# Check the output has the right shape # Check the output has the right shape
if shape is not None: if shape is not None:
if self.ndim_added == 0 and shape != rval.shape: if self.ndim_added == 0 and shape != rval.shape:
raise ValueError('Shape mismatch: "out" should have shape %s, but the value produced by "perform" has shape %s'\ raise ValueError(
'Shape mismatch: "out" should have shape %s, but the'
' value produced by "perform" has shape %s'
% (shape, rval.shape)) % (shape, rval.shape))
elif self.ndim_added > 0 and shape != rval.shape[:-self.ndim_added]: elif (self.ndim_added > 0 and
raise ValueError('Shape mismatch: "out" should have shape starting with %s (plus %i extra dimensions), but the value produced by "perform" has shape %s'\ shape != rval.shape[:-self.ndim_added]):
raise ValueError(
'Shape mismatch: "out" should have shape starting with'
' %s (plus %i extra dimensions), but the value produced'
' by "perform" has shape %s'
% (shape, self.ndim_added, rval.shape)) % (shape, self.ndim_added, rval.shape))
out[0] = rval out[0] = rval
...@@ -260,9 +273,11 @@ def _infer_ndim_bcast(ndim, shape, *args): ...@@ -260,9 +273,11 @@ def _infer_ndim_bcast(ndim, shape, *args):
# there is a convention that -1 means the corresponding shape of a # there is a convention that -1 means the corresponding shape of a
# potentially-broadcasted symbolic arg # potentially-broadcasted symbolic arg
if (isinstance(shape, (tuple, list)) if (isinstance(shape, (tuple, list))
and numpy.all(numpy.asarray(shape)>=0)): and numpy.all(numpy.asarray(shape) >= 0)):
bcast = [(s==1) for s in shape] bcast = [(s == 1) for s in shape]
v_shape = tensor.TensorConstant(type=tensor.lvector, data=theano._asarray(shape, dtype='int64')) v_shape = tensor.TensorConstant(type=tensor.lvector,
data=theano._asarray(shape,
dtype='int64'))
shape_ndim = len(shape) shape_ndim = len(shape)
if ndim is None: if ndim is None:
ndim = shape_ndim ndim = shape_ndim
...@@ -278,21 +293,21 @@ def _infer_ndim_bcast(ndim, shape, *args): ...@@ -278,21 +293,21 @@ def _infer_ndim_bcast(ndim, shape, *args):
# This case combines together symbolic and non-symbolic shape # This case combines together symbolic and non-symbolic shape
# information # information
if ndim is None: if ndim is None:
ndim=args_ndim ndim = args_ndim
else: else:
ndim = max(args_ndim, ndim) ndim = max(args_ndim, ndim)
ndim = max(args_ndim, len(shape)) ndim = max(args_ndim, len(shape))
shape = [-1]*(ndim - len(shape))+list(shape) shape = [-1] * (ndim - len(shape)) + list(shape)
bcast = [] bcast = []
pre_v_shape = [] pre_v_shape = []
for i,s in enumerate(shape): for i, s in enumerate(shape):
if hasattr(s, 'type'): # s is symbolic if hasattr(s, 'type'): # s is symbolic
bcast.append(False) # todo - introspect further bcast.append(False) # todo - introspect further
pre_v_shape.append(s) pre_v_shape.append(s)
else: else:
if s >= 0: if s >= 0:
pre_v_shape.append(tensor.as_tensor_variable(s)) pre_v_shape.append(tensor.as_tensor_variable(s))
bcast.append((s==1)) bcast.append((s == 1))
elif s == -1: elif s == -1:
n_a_i = 0 n_a_i = 0
for a in args: for a in args:
...@@ -301,7 +316,7 @@ def _infer_ndim_bcast(ndim, shape, *args): ...@@ -301,7 +316,7 @@ def _infer_ndim_bcast(ndim, shape, *args):
# i # i
if i >= ndim - a.ndim: if i >= ndim - a.ndim:
n_a_i += 1 n_a_i += 1
a_i = i + a.ndim -ndim a_i = i + a.ndim - ndim
if not a.broadcastable[a_i]: if not a.broadcastable[a_i]:
pre_v_shape.append(a.shape[a_i]) pre_v_shape.append(a.shape[a_i])
bcast.append(False) bcast.append(False)
...@@ -316,7 +331,8 @@ def _infer_ndim_bcast(ndim, shape, *args): ...@@ -316,7 +331,8 @@ def _infer_ndim_bcast(ndim, shape, *args):
bcast.append(True) bcast.append(True)
else: else:
ValueError('negative shape', s) ValueError('negative shape', s)
# post-condition: shape may still contain both symbolic and non-symbolic things # post-condition: shape may still contain both symbolic and
# non-symbolic things
v_shape = tensor.stack(*pre_v_shape) v_shape = tensor.stack(*pre_v_shape)
elif shape is None: elif shape is None:
...@@ -325,7 +341,7 @@ def _infer_ndim_bcast(ndim, shape, *args): ...@@ -325,7 +341,7 @@ def _infer_ndim_bcast(ndim, shape, *args):
if not args: if not args:
raise TypeError(('_infer_ndim_bcast cannot infer shape without' raise TypeError(('_infer_ndim_bcast cannot infer shape without'
' either shape or args')) ' either shape or args'))
template = reduce(lambda a,b:a+b, args) template = reduce(lambda a, b: a + b, args)
v_shape = template.shape v_shape = template.shape
bcast = template.broadcastable bcast = template.broadcastable
ndim = template.ndim ndim = template.ndim
...@@ -333,18 +349,22 @@ def _infer_ndim_bcast(ndim, shape, *args): ...@@ -333,18 +349,22 @@ def _infer_ndim_bcast(ndim, shape, *args):
v_shape = tensor.as_tensor_variable(shape) v_shape = tensor.as_tensor_variable(shape)
if ndim is None: if ndim is None:
ndim = tensor.get_vector_length(v_shape) ndim = tensor.get_vector_length(v_shape)
bcast = [False]*ndim bcast = [False] * ndim
if not (v_shape.dtype.startswith('int') or v_shape.dtype.startswith('uint')): if (not (v_shape.dtype.startswith('int') or
raise TypeError('shape must be an integer vector or list', v_shape.dtype) v_shape.dtype.startswith('uint'))):
raise TypeError('shape must be an integer vector or list',
v_shape.dtype)
if args_ndim > ndim: if args_ndim > ndim:
raise ValueError('ndim should be at least as big as required by args value', raise ValueError(
'ndim should be at least as big as required by args value',
(ndim, args_ndim), args) (ndim, args_ndim), args)
assert ndim == len(bcast) assert ndim == len(bcast)
return ndim, tensor.cast(v_shape, 'int32'), tuple(bcast) return ndim, tensor.cast(v_shape, 'int32'), tuple(bcast)
def _generate_broadcasting_indices(out_shape, *shapes): def _generate_broadcasting_indices(out_shape, *shapes):
''' '''
Return indices over each shape that broadcast them to match out_shape. Return indices over each shape that broadcast them to match out_shape.
...@@ -359,11 +379,11 @@ def _generate_broadcasting_indices(out_shape, *shapes): ...@@ -359,11 +379,11 @@ def _generate_broadcasting_indices(out_shape, *shapes):
''' '''
all_shapes = (out_shape,) + shapes all_shapes = (out_shape,) + shapes
# Will contain the return value: a list of indices for each argument # Will contain the return value: a list of indices for each argument
ret_indices = [ [()] for shape in all_shapes ] ret_indices = [[()] for shape in all_shapes]
for dim in xrange(len(out_shape)): for dim in xrange(len(out_shape)):
# Temporary list to generate the indices # Temporary list to generate the indices
_ret_indices = [ [] for shape in all_shapes ] _ret_indices = [[] for shape in all_shapes]
out_range = range(out_shape[dim]) out_range = range(out_shape[dim])
...@@ -373,11 +393,14 @@ def _generate_broadcasting_indices(out_shape, *shapes): ...@@ -373,11 +393,14 @@ def _generate_broadcasting_indices(out_shape, *shapes):
for shape in shapes: for shape in shapes:
if shape[dim] == out_shape[dim]: if shape[dim] == out_shape[dim]:
ranges.append(out_range) ranges.append(out_range)
elif shape[dim] == 1: #broadcast elif shape[dim] == 1: # broadcast
ranges.append([0] * out_shape[dim]) ranges.append([0] * out_shape[dim])
else: else:
raise ValueError('shape[%i] (%i) should be equal to out_shape[%i] (%i) or to 1'\ raise ValueError(
% (dim, shape[dim], dim, out_shape[dim]), shape, out_shape, shapes) 'shape[%i] (%i) should be equal to out_shape[%i] (%i) or'
' to 1'
% (dim, shape[dim], dim, out_shape[dim]), shape,
out_shape, shapes)
for prev_index in zip(*ret_indices): for prev_index in zip(*ret_indices):
for dim_index in zip(*ranges): for dim_index in zip(*ranges):
...@@ -435,7 +458,8 @@ def normal(random_state, size=None, avg=0.0, std=1.0, ndim=None, dtype=None): ...@@ -435,7 +458,8 @@ def normal(random_state, size=None, avg=0.0, std=1.0, ndim=None, dtype=None):
return op(random_state, size, avg, std) return op(random_state, size, avg, std)
def binomial(random_state, size=None, n=1, p=0.5, ndim=None, dtype='int64', prob=None): def binomial(random_state, size=None, n=1, p=0.5, ndim=None,
dtype='int64', prob=None):
""" """
Sample n times with probability of success prob for each trial, Sample n times with probability of success prob for each trial,
return the number of successes. return the number of successes.
...@@ -452,7 +476,7 @@ def binomial(random_state, size=None, n=1, p=0.5, ndim=None, dtype='int64', prob ...@@ -452,7 +476,7 @@ def binomial(random_state, size=None, n=1, p=0.5, ndim=None, dtype='int64', prob
n = tensor.as_tensor_variable(n) n = tensor.as_tensor_variable(n)
p = tensor.as_tensor_variable(p) p = tensor.as_tensor_variable(p)
ndim, size, bcast = _infer_ndim_bcast(ndim, size, n, p) ndim, size, bcast = _infer_ndim_bcast(ndim, size, n, p)
if n.dtype=='int64': if n.dtype == 'int64':
### THIS WORKS AROUND A NUMPY BUG on 32bit machine ### THIS WORKS AROUND A NUMPY BUG on 32bit machine
### Erase when the following works on a 32bit machine: ### Erase when the following works on a 32bit machine:
### numpy.random.binomial( ### numpy.random.binomial(
...@@ -460,9 +484,10 @@ def binomial(random_state, size=None, n=1, p=0.5, ndim=None, dtype='int64', prob ...@@ -460,9 +484,10 @@ def binomial(random_state, size=None, n=1, p=0.5, ndim=None, dtype='int64', prob
# p=numpy.asarray([.1, .2, .3], dtype='float64')) # p=numpy.asarray([.1, .2, .3], dtype='float64'))
n = tensor.cast(n, 'int32') n = tensor.cast(n, 'int32')
op = RandomFunction('binomial', op = RandomFunction('binomial',
tensor.TensorType(dtype = dtype, broadcastable = (False,)*ndim) ) tensor.TensorType(dtype=dtype, broadcastable=(False,) * ndim))
return op(random_state, size, n, p) return op(random_state, size, n, p)
def random_integers_helper(random_state, low, high, size): def random_integers_helper(random_state, low, high, size):
''' '''
Helper function to draw random integers. Helper function to draw random integers.
...@@ -477,16 +502,19 @@ def random_integers_helper(random_state, low, high, size): ...@@ -477,16 +502,19 @@ def random_integers_helper(random_state, low, high, size):
out_ndim = max(low.ndim, high.ndim) out_ndim = max(low.ndim, high.ndim)
# broadcast low and high to out_ndim dimensions # broadcast low and high to out_ndim dimensions
if low.ndim > out_ndim: if low.ndim > out_ndim:
raise ValueError('low.ndim (%i) should not be larger than len(size) (%i)' % (low.ndim, out_ndim), raise ValueError(
'low.ndim (%i) should not be larger than len(size) (%i)'
% (low.ndim, out_ndim),
low, size) low, size)
if low.ndim < out_ndim: if low.ndim < out_ndim:
low = low.reshape((1,)*(out_ndim-low.ndim) + low.shape) low = low.reshape((1,) * (out_ndim - low.ndim) + low.shape)
if high.ndim > out_ndim: if high.ndim > out_ndim:
raise ValueError('high.ndim (%i) should not be larger than len(size) (%i)' % (high.ndim, out_ndim), raise ValueError(
high, size) 'high.ndim (%i) should not be larger than len(size) (%i)'
% (high.ndim, out_ndim), high, size)
if high.ndim < out_ndim: if high.ndim < out_ndim:
high = high.reshape((1,)*(out_ndim-high.ndim) + high.shape) high = high.reshape((1,) * (out_ndim - high.ndim) + high.shape)
if size is not None: if size is not None:
out_size = tuple(size) out_size = tuple(size)
...@@ -498,14 +526,17 @@ def random_integers_helper(random_state, low, high, size): ...@@ -498,14 +526,17 @@ def random_integers_helper(random_state, low, high, size):
# Build the indices over which to loop # Build the indices over which to loop
out = numpy.ndarray(out_size) out = numpy.ndarray(out_size)
broadcast_ind = _generate_broadcasting_indices(out_size, low.shape, high.shape) broadcast_ind = _generate_broadcasting_indices(out_size, low.shape,
high.shape)
# Iterate over these indices, drawing one sample at a time from numpy # Iterate over these indices, drawing one sample at a time from numpy
for oi, li, hi in zip(*broadcast_ind): for oi, li, hi in zip(*broadcast_ind):
out[oi] = random_state.random_integers(low = low[li], high = high[hi]) out[oi] = random_state.random_integers(low=low[li], high=high[hi])
return out return out
def random_integers(random_state, size=None, low=0, high=1, ndim=None, dtype='int64'):
def random_integers(random_state, size=None, low=0, high=1, ndim=None,
dtype='int64'):
""" """
Sample a random integer between low and high, both inclusive. Sample a random integer between low and high, both inclusive.
...@@ -522,6 +553,7 @@ def random_integers(random_state, size=None, low=0, high=1, ndim=None, dtype='in ...@@ -522,6 +553,7 @@ def random_integers(random_state, size=None, low=0, high=1, ndim=None, dtype='in
tensor.TensorType(dtype=dtype, broadcastable=bcast)) tensor.TensorType(dtype=dtype, broadcastable=bcast))
return op(random_state, size, low, high) return op(random_state, size, low, high)
def permutation_helper(random_state, n, shape): def permutation_helper(random_state, n, shape):
"""Helper function to generate permutations from integers. """Helper function to generate permutations from integers.
...@@ -552,6 +584,7 @@ def permutation_helper(random_state, n, shape): ...@@ -552,6 +584,7 @@ def permutation_helper(random_state, n, shape):
#print 'RETURNING', out.shape #print 'RETURNING', out.shape
return out return out
def permutation(random_state, size=None, n=1, ndim=None, dtype='int64'): def permutation(random_state, size=None, n=1, ndim=None, dtype='int64'):
""" """
Returns permutations of the integers between 0 and n-1, as many times Returns permutations of the integers between 0 and n-1, as many times
...@@ -569,10 +602,11 @@ def permutation(random_state, size=None, n=1, ndim=None, dtype='int64'): ...@@ -569,10 +602,11 @@ def permutation(random_state, size=None, n=1, ndim=None, dtype='int64'):
ndim, size, bcast = _infer_ndim_bcast(ndim, size) ndim, size, bcast = _infer_ndim_bcast(ndim, size)
#print "NDIM", ndim, size #print "NDIM", ndim, size
op = RandomFunction(permutation_helper, op = RandomFunction(permutation_helper,
tensor.TensorType(dtype=dtype, broadcastable=bcast+(False,)), tensor.TensorType(dtype=dtype, broadcastable=bcast + (False,)),
ndim_added=1) ndim_added=1)
return op(random_state, size, n) return op(random_state, size, n)
def multinomial_helper(random_state, n, pvals, size): def multinomial_helper(random_state, n, pvals, size):
''' '''
Helper function drawing from multinomial distributions. Helper function drawing from multinomial distributions.
...@@ -586,21 +620,25 @@ def multinomial_helper(random_state, n, pvals, size): ...@@ -586,21 +620,25 @@ def multinomial_helper(random_state, n, pvals, size):
if size is not None: if size is not None:
ndim = len(size) ndim = len(size)
else: else:
ndim = max(n.ndim, pvals.ndim-1) ndim = max(n.ndim, pvals.ndim - 1)
out_ndim = ndim+1 out_ndim = ndim + 1
# broadcast n to ndim dimensions and pvals to ndim+1 # broadcast n to ndim dimensions and pvals to ndim+1
if n.ndim > ndim: if n.ndim > ndim:
raise ValueError('n.ndim (%i) should not be larger than len(size) (%i)' % (n.ndim, ndim), raise ValueError(
'n.ndim (%i) should not be larger than len(size) (%i)'
% (n.ndim, ndim),
n, size) n, size)
if n.ndim < ndim: if n.ndim < ndim:
n = n.reshape((1,)*(ndim-n.ndim) + n.shape) n = n.reshape((1,) * (ndim - n.ndim) + n.shape)
if pvals.ndim-1 > ndim: if pvals.ndim - 1 > ndim:
raise ValueError('pvals.ndim-1 (%i) should not be larger than len(size) (%i)' % (pvals.ndim-1, ndim), raise ValueError(
'pvals.ndim-1 (%i) should not be larger than len(size) (%i)'
% (pvals.ndim - 1, ndim),
pvals, size) pvals, size)
if pvals.ndim-1 < ndim: if pvals.ndim - 1 < ndim:
pvals = pvals.reshape((1,)*(ndim-pvals.ndim+1) + pvals.shape) pvals = pvals.reshape((1,) * (ndim - pvals.ndim + 1) + pvals.shape)
if size is not None: if size is not None:
size = tuple(size) size = tuple(size)
...@@ -609,14 +647,16 @@ def multinomial_helper(random_state, n, pvals, size): ...@@ -609,14 +647,16 @@ def multinomial_helper(random_state, n, pvals, size):
for dim in xrange(ndim): for dim in xrange(ndim):
dim_len = max(n.shape[dim], pvals.shape[dim]) dim_len = max(n.shape[dim], pvals.shape[dim])
size = size + (dim_len,) size = size + (dim_len,)
out_size = size+(pvals.shape[-1],) out_size = size + (pvals.shape[-1],)
# Build the indices over which to loop # Build the indices over which to loop
# Note that here, the rows (inner-most 1D subtensors) of pvals and out # Note that here, the rows (inner-most 1D subtensors) of pvals and out
# are indexed, not their individual elements # are indexed, not their individual elements
out = numpy.ndarray(out_size) out = numpy.ndarray(out_size)
broadcast_ind = _generate_broadcasting_indices(size, n.shape, pvals.shape[:-1]) broadcast_ind = _generate_broadcasting_indices(size, n.shape,
# Iterate over these indices, drawing from one multinomial at a time from numpy pvals.shape[:-1])
# Iterate over these indices, drawing from one multinomial at a
# time from numpy
assert pvals.min() >= 0 assert pvals.min() >= 0
for mi, ni, pi in zip(*broadcast_ind): for mi, ni, pi in zip(*broadcast_ind):
pvi = pvals[pi] pvi = pvals[pi]
...@@ -627,24 +667,24 @@ def multinomial_helper(random_state, n, pvals, size): ...@@ -627,24 +667,24 @@ def multinomial_helper(random_state, n, pvals, size):
# In perfect arithmetic this would be correct, but in float32 or # In perfect arithmetic this would be correct, but in float32 or
# float64 it is too strict. # float64 it is too strict.
pisum = numpy.sum(pvi) pisum = numpy.sum(pvi)
if 1.0 < pisum < 1.0+1e-5:#correct if we went a little over if 1.0 < pisum < 1.0 + 1e-5: # correct if we went a little over
# because mtrand.pyx has a ValueError that will trigger if # because mtrand.pyx has a ValueError that will trigger if
# sum(pvals[:-1]) > 1.0 # sum(pvals[:-1]) > 1.0
pvi = pvi * (1.0 - 5e-5) pvi = pvi * (1.0 - 5e-5)
#pvi = pvi * .9 #pvi = pvi * .9
pisum = numpy.sum(pvi) pisum = numpy.sum(pvi)
elif pvi[-1]<5e-5: #will this even work? elif pvi[-1] < 5e-5: # will this even work?
pvi = pvi * (1.0 - 5e-5) pvi = pvi * (1.0 - 5e-5)
pisum = numpy.sum(pvi) pisum = numpy.sum(pvi)
assert pisum<=1.0, pisum assert pisum <= 1.0, pisum
out[mi] = random_state.multinomial(n=n[ni], out[mi] = random_state.multinomial(n=n[ni],
pvals=pvi.astype('float64')) pvals=pvi.astype('float64'))
return out return out
def multinomial(random_state, size=None, n=1, pvals=[0.5, 0.5], def multinomial(random_state, size=None, n=1, pvals=[0.5, 0.5],
ndim=None, dtype='int64'): ndim=None, dtype='int64'):
""" """Sample from one or more multinomial distributions defined by
Sample from one or more multinomial distributions defined by
one-dimensional slices in pvals. one-dimensional slices in pvals.
:param pvals: a tensor of shape "nmulti+(L,)" describing each multinomial :param pvals: a tensor of shape "nmulti+(L,)" describing each multinomial
...@@ -657,15 +697,17 @@ def multinomial(random_state, size=None, n=1, pvals=[0.5, 0.5], ...@@ -657,15 +697,17 @@ def multinomial(random_state, size=None, n=1, pvals=[0.5, 0.5],
right in nmulti. (See examples below.) right in nmulti. (See examples below.)
Default ``None`` means size=nmulti. Default ``None`` means size=nmulti.
:param n: the number of experiments to simulate for each multinomial. This :param n: the number of experiments to simulate for each
can be a scalar, or tensor, it will be broadcasted to have shape "nmulti". multinomial. This can be a scalar, or tensor, it will be
broadcasted to have shape "nmulti".
:param dtype: the dtype of the return value (which will represent counts) :param dtype: the dtype of the return value (which will represent counts)
:returns: tensor of len(size)+1 dimensions, and shape[-1]==L, with the specified ``dtype``, :returns: tensor of len(size)+1 dimensions, and shape[-1]==L, with
with the experiment counts. See examples to understand the shape of the the specified ``dtype``, with the experiment counts. See
return value, which is derived from both size and pvals.shape. examples to understand the shape of the return value, which is
In return value rval, "numpy.allclose(rval.sum(axis=-1), n)" will be true. derived from both size and pvals.shape. In return value rval,
"numpy.allclose(rval.sum(axis=-1), n)" will be true.
For example, to simulate n experiments from each multinomial in a batch of For example, to simulate n experiments from each multinomial in a batch of
size B: size B:
...@@ -685,11 +727,12 @@ def multinomial(random_state, size=None, n=1, pvals=[0.5, 0.5], ...@@ -685,11 +727,12 @@ def multinomial(random_state, size=None, n=1, pvals=[0.5, 0.5],
Using size for broadcasting of pvals: Using size for broadcasting of pvals:
size=(10,1,-1), pvals.shape=(A,B,L) size=(10, 1, -1), pvals.shape=(A, B, L)
--> rval.shape=[10,1,B,L], and requires that A==1. --> rval.shape=[10,1,B,L], and requires that A==1.
rval[l,k,i,j] is the count of possibility j in the distribution specified rval[l,k,i,j] is the count of possibility j in the
by pvals[k,i], in the l'th of 10 draws. distribution specified by pvals[k,i], in the l'th of 10
draws.
""" """
n = tensor.as_tensor_variable(n) n = tensor.as_tensor_variable(n)
...@@ -697,9 +740,9 @@ def multinomial(random_state, size=None, n=1, pvals=[0.5, 0.5], ...@@ -697,9 +740,9 @@ def multinomial(random_state, size=None, n=1, pvals=[0.5, 0.5],
# until ellipsis is implemented (argh) # until ellipsis is implemented (argh)
tmp = pvals.T[0].T tmp = pvals.T[0].T
ndim, size, bcast = _infer_ndim_bcast(ndim, size, n, tmp) ndim, size, bcast = _infer_ndim_bcast(ndim, size, n, tmp)
bcast = bcast+(pvals.type.broadcastable[-1],) bcast = bcast + (pvals.type.broadcastable[-1],)
op = RandomFunction(multinomial_helper, op = RandomFunction(multinomial_helper,
tensor.TensorType(dtype = dtype, broadcastable = bcast), tensor.TensorType(dtype=dtype, broadcastable=bcast),
ndim_added=1) ndim_added=1)
return op(random_state, size, n, pvals) return op(random_state, size, n, pvals)
...@@ -708,17 +751,20 @@ def multinomial(random_state, size=None, n=1, pvals=[0.5, 0.5], ...@@ -708,17 +751,20 @@ def multinomial(random_state, size=None, n=1, pvals=[0.5, 0.5],
def random_make_inplace(node): def random_make_inplace(node):
op = node.op op = node.op
if isinstance(op, RandomFunction) and not op.inplace: if isinstance(op, RandomFunction) and not op.inplace:
new_op = RandomFunction(op.fn, op.outtype, inplace=True, ndim_added=op.ndim_added) new_op = RandomFunction(op.fn, op.outtype, inplace=True,
ndim_added=op.ndim_added)
return new_op.make_node(*node.inputs).outputs return new_op.make_node(*node.inputs).outputs
return False return False
optdb.register('random_make_inplace', opt.in2out(random_make_inplace, ignore_newtrees=True), 99, 'fast_run', 'inplace') optdb.register('random_make_inplace', opt.in2out(random_make_inplace,
ignore_newtrees=True),
99, 'fast_run', 'inplace')
class RandomStreamsBase(object): class RandomStreamsBase(object):
def binomial(self, size=None, n=1, p=0.5, ndim=None, dtype='int64', prob=None): def binomial(self, size=None, n=1, p=0.5, ndim=None, dtype='int64',
prob=None):
""" """
Sample n times with probability of success prob for each trial, Sample n times with probability of success prob for each trial,
return the number of successes. return the number of successes.
...@@ -754,7 +800,8 @@ class RandomStreamsBase(object): ...@@ -754,7 +800,8 @@ class RandomStreamsBase(object):
""" """
return self.gen(normal, size, avg, std, ndim=ndim, dtype=dtype) return self.gen(normal, size, avg, std, ndim=ndim, dtype=dtype)
def random_integers(self, size=None, low=0, high=1, ndim=None, dtype='int64'): def random_integers(self, size=None, low=0, high=1, ndim=None,
dtype='int64'):
""" """
Sample a random integer between low and high, both inclusive. Sample a random integer between low and high, both inclusive.
...@@ -762,7 +809,8 @@ class RandomStreamsBase(object): ...@@ -762,7 +809,8 @@ class RandomStreamsBase(object):
ndim may be a plain integer to supplement the missing ndim may be a plain integer to supplement the missing
information. information.
""" """
return self.gen(random_integers, size, low, high, ndim=ndim, dtype=dtype) return self.gen(random_integers, size, low, high, ndim=ndim,
dtype=dtype)
def permutation(self, size=None, n=1, ndim=None, dtype='int64'): def permutation(self, size=None, n=1, ndim=None, dtype='int64'):
""" """
...@@ -780,7 +828,8 @@ class RandomStreamsBase(object): ...@@ -780,7 +828,8 @@ class RandomStreamsBase(object):
""" """
return self.gen(permutation, size, n, ndim=ndim, dtype=dtype) return self.gen(permutation, size, n, ndim=ndim, dtype=dtype)
def multinomial(self, size=None, n=1, pvals=[0.5, 0.5], ndim=None, dtype='int64'): def multinomial(self, size=None, n=1, pvals=[0.5, 0.5], ndim=None,
dtype='int64'):
""" """
Sample n times from a multinomial distribution defined by Sample n times from a multinomial distribution defined by
probabilities pvals, as many times as required by size. For probabilities pvals, as many times as required by size. For
...@@ -802,7 +851,8 @@ class RandomStreamsBase(object): ...@@ -802,7 +851,8 @@ class RandomStreamsBase(object):
This uses permutation random variable internally, available via This uses permutation random variable internally, available via
the ``.permutation`` attribute of the return value. the ``.permutation`` attribute of the return value.
""" """
perm = self.permutation(size=input.shape[:-1], n=input.shape[-1], ndim=input.ndim-1) perm = self.permutation(size=input.shape[:-1], n=input.shape[-1],
ndim=input.ndim - 1)
shuffled = tensor.permute_row_elements(input, perm) shuffled = tensor.permute_row_elements(input, perm)
shuffled.permutation = perm shuffled.permutation = perm
return shuffled return shuffled
Markdown 格式
0% 或
您添加了 0 人 到此讨论。请谨慎行事。
请先完成此评论的编辑！
请 注册 或者 后发表评论