"""NumArray Universal Function Module
$Id: ufunc.py,v 1.157 2006/05/12 17:41:01 jaytmiller Exp $
"""
import types
import memory
import _sort, _bytes
from operator import and_
import _ufuncall
import generic as _gen
import _ufunc
import _converter
import numerictypes as _nt
import numarraycore as _nc
import safethread
import warnings as _warnings
# note more imports are at end of file
#import pywin.debugger
"""DEBUGGING NOTE:
A reasonable amount of care went into maintaining the Python
prototype for _callOverDimensions, _UnaryUFunc, _BinaryUFunc,
_Converter, and _Operator. Thus, problems in the C
implementation can sometimes be isolated by substituting the
original Python code back in, either in whole or part.
Substituting in the whole is easy: set _PROTOTYPE to 1 below.
"""
_PROTOTYPE = 0
_FPE_DIVIDE_BY_ZERO = 1
_FPE_OVERFLOW = 2
_FPE_UNDERFLOW = 4
_FPE_INVALID = 8
class _NumErrorMode:
_errorset = ("ignore", "warn", "raise")
def _num_mode_test(self, val, all, default):
if val not in self._errorset and val is not None:
raise ValueError("Error mode value must be in " +
repr(self._errorset) + " or None.")
if all is not None:
if val is not default:
return 1
else:
return 0
else:
return 1
def __init__(self, all=None, overflow="warn", underflow="ignore",
dividebyzero="warn", invalid="warn"):
if all in self._errorset:
self.overflow = all
self.underflow = all
self.dividebyzero = all
self.invalid = all
if self._num_mode_test(overflow, all, "warn"):
self.overflow = overflow
if self._num_mode_test(underflow, all, "ignore"):
self.underflow = underflow
if self._num_mode_test(dividebyzero, all, "warn"):
self.dividebyzero = dividebyzero
if self._num_mode_test(invalid, all, "warn"):
self.invalid = invalid
def __repr__(self):
return ("_NumErrorMode(overflow='%s', underflow='%s', dividebyzero='%s', invalid='%s')" % (self.overflow, self.underflow, self.dividebyzero, self.invalid))
class NumError:
"""Defines how numeric errors should be handled.
Error configuration is managed as a mapping (by thread) of stacks.
Thus it is possible for each thread to have an independent error
configuration, and also for multiple routines within the same
thread to save and restore error configurations as needed.
"""
def __init__(self, all=None,
overflow="warn", underflow="ignore", dividebyzero="warn",
invalid="warn"):
self._modestack = {} # map of stacks indexed by thread id
self._defaultmode = _NumErrorMode()
self.setMode(all=all, underflow=underflow, overflow=overflow,
dividebyzero=dividebyzero, invalid=invalid)
def setDefaultMode(self, all=None, overflow="warn", underflow="ignore",
dividebyzero="warn", invalid="warn"):
self._defaultmode = _NumErrorMode(all, overflow, underflow,
dividebyzero, invalid)
def setMode(self, all=None, overflow="warn", underflow="ignore",
dividebyzero="warn", invalid="warn"):
self._setmodes( _NumErrorMode(all, overflow, underflow,
dividebyzero, invalid) )
def pushMode(self, all=None, overflow="warn", underflow="ignore",
dividebyzero="warn", invalid="warn"):
self._pushmodes( _NumErrorMode(all, overflow, underflow,
dividebyzero, invalid) )
def _getmodestack(self):
id = safethread.get_ident()
try:
l = self._modestack[id]
except KeyError:
l = [ self._defaultmode ]
self._modestack[id] = l
if l == []:
l = [ self._defaultmode ]
return l
def _setmodes(self, modes):
l = self._getmodestack()
l[-1] = modes
def getMode(self):
l = self._getmodestack()
return l[-1]
def _pushmodes(self, modes):
l = self._getmodestack()
l.append(modes)
def popMode(self):
l = self._getmodestack()
return l.pop()
Error = NumError()
class MathDomainError(ArithmeticError): pass
class UnderflowError(ArithmeticError): pass
class NumOverflowError(OverflowError, ArithmeticError): pass
def handleError(errorStatus, sourcemsg):
"""Take error status and use error mode to handle it."""
modes = Error.getMode()
if errorStatus & _FPE_INVALID:
if modes.invalid == "warn":
print "Warning: Encountered invalid numeric result(s)", sourcemsg
if modes.invalid == "raise":
raise MathDomainError(sourcemsg)
if errorStatus & _FPE_DIVIDE_BY_ZERO:
if modes.dividebyzero == "warn":
print "Warning: Encountered divide by zero(s)", sourcemsg
if modes.dividebyzero == "raise":
raise ZeroDivisionError(sourcemsg)
if errorStatus & _FPE_OVERFLOW:
if modes.overflow == "warn":
print "Warning: Encountered overflow(s)", sourcemsg
if modes.overflow == "raise":
raise NumOverflowError(sourcemsg)
if errorStatus & _FPE_UNDERFLOW:
if modes.underflow == "warn":
print "Warning: Encountered underflow(s)", sourcemsg
if modes.underflow == "raise":
raise UnderflowError(sourcemsg)
def _nIOArgs(cfuncs):
"""Determine number of arguments from signature
Returns tuple of number of input and output arguments.
Also check for consistency in numbers of arguments.
"""
# The number of items in the optype key implies the number of inputs,
# hence the number of inputs = # of underscores+1
noptypeinputs = []
ninputs = []
noutputs = []
for key in cfuncs.keys():
if 'R' in key or 'A' in key: # reduce, accumulate
continue
noptypeinputs.append(key.count('_'))
for informkey in cfuncs[key].keys():
ninputs.append(len(informkey))
noutputs.append(len(cfuncs[key][informkey][0]))
return max(ninputs), max(noutputs)
if _PROTOTYPE:
def _getBlockingParameters(shape, niter, overlap=0):
try:
return _blockingParametersCache[ (shape, niter, overlap) ]
except KeyError:
return _getBlockingMiss(shape, niter, overlap)
def _getBlockingMiss(shape, niter, overlap=0):
"""Figure out how to split the input array into subnumarray
close (but less than) the size of the block buffers.
overlap is used in cases where it is necessary to overlap the
last dimension in computations (for example for reduce and
accumulate.)
returns:
level: index at which the c functions should be called
blockingparameters: tuple consisting of
number of full subnumarray at the given level
subarray shape
a flag indicating if there is a leftover part of the array
which is smaller than the regular subarray (leftover)
the shape of the leftover subarray
"""
if shape == ():
retval = 0, (1, (), 0, None)
_blockingParametersCache[(shape,niter,overlap)] = retval
return retval
ndim = len(shape)
csize = long(shape[-1]) # cumulative shape
if csize > niter:
# last dimension is too big for buffer
# special overlap handling in this case
nblocks, leftover = divmod(csize, niter-overlap)
nblocks, leftover = int(nblocks), int(leftover)
if overlap:
if leftover:
retval = ndim-1, (nblocks, (niter,), 1, (leftover,))
else:
retval = ndim-1, (nblocks-1, (niter,), 1, (niter-1,))
else:
if leftover:
retval = ndim-1, (nblocks, (niter,), 1, (leftover,))
else:
retval = ndim-1, (nblocks, (niter,), 0, None)
else:
for i in xrange(2,ndim+1):
csize *= shape[-i]
if csize >= niter:
# need to break up this dimension
tdim = long(shape[-i])
blockdim = niter*tdim/csize
nblocks, leftover = divmod(tdim, blockdim)
nblocks, leftover = int(nblocks), int(leftover)
stail = shape[-i+1:]
if leftover:
retval = ndim-i, (nblocks, (blockdim,)+stail,
1, (leftover,) + stail)
else:
retval = ndim-i, (nblocks, (blockdim,)+stail, 0, None)
break
else:
# don't need to break up
retval = 0, (1, shape, 0, None)
_blockingParametersCache[(shape,niter,overlap)] = retval
return retval
# _doOverDimensions has been re-implemented in C.
def _doOverDimensions(objects, outshape,
dims, indexlevel, blockingparameters, overlap=0, level=0):
"""This version handles alignment, byteswapping and n-dimensions
This is a recursive method dims starts out as an empty list, and
accumulates the indices for each dimension. Eg. if the array has
shape (7, 3, 5) on the initial call, the next level call will
specify the the 'first' index (corresponding to the 7), the next
will have the first two indices (and that would be the last level
since the lowest dimension is taken care of by the cfuncs.)
"""
if level == indexlevel:
nregShapeIters, shape, leftover, leftoverShape, = blockingparameters
if shape != ():
dimval = shape[0]-overlap
for i in xrange(nregShapeIters):
for o in objects:
o.compute(dims+[i*dimval], shape)
if leftover:
for o in objects:
o.compute(dims+[(i+1)*dimval], leftoverShape)
else:
for o in objects:
o.compute([0], shape)
else:
# recurse
for i in xrange(outshape[level]):
_doOverDimensions(objects, outshape,
dims+[i], indexlevel, blockingparameters,
overlap, level+1)
def _callOverDimensions(objects, outshape, index, blocking, overlap=0, level=0):
return _doOverDimensions(objects, outshape, [], index, blocking, overlap, level)
# Replace the Python version of _callOverDimensions with one implemented in C.
# Although _doOverDimensions is not replaced here, it has been re-implemented
# in C as well
if not _PROTOTYPE:
from _ufunc import _callOverDimensions
def _copyFromAndConvertMiss(inarr, outarr):
key = (_digest(inarr), _digest(outarr), safethread.get_ident())
# Type conversion done using inputconverter
input = _InputConverter(inarr, outarr._type.name)
# Ensure output step uses output buffer of input step
output = _OutputConverter(outarr, outarr._type.name,
inbuffer=input.result_buff,
forcestride=1)
maxitemsize = max(inarr._type.bytes, outarr._type.bytes)
niter = _ufunc.getBufferSize()/maxitemsize
if (isinstance(inarr, _nc.NumArray) and
isinstance(outarr, _nc.NumArray)):
_copyCache[ key ] = input, output, niter
return input, output, niter
if _PROTOTYPE:
def _copyFromAndConvert(inarr, outarr):
"""Copy from one array to another handling strides, form & conv"""
if inarr._shape != outarr._shape:
raise ValueError("Arrays must have the same shape")
key = (_digest(inarr), _digest(outarr), safethread.get_ident())
try:
input, output, niter = _copyCache[ key ]
except KeyError:
input, output, niter = _copyFromAndConvertMiss(inarr, outarr)
input.rebuffer(inarr)
output.rebuffer(outarr, input.result_buff)
indexlevel, blockingparameters = \
_ufunc._getBlockingParameters(outarr._shape, niter)
objects = (input, output)
_ufunc.CheckFPErrors()
_callOverDimensions(objects, outarr._shape, indexlevel,
blockingparameters)
errorstatus = _ufunc.CheckFPErrors()
input.clean(inarr)
output.clean(inarr)
output.clean(outarr)
if errorstatus:
handleError(errorstatus, " during type conversion")
else:
from _ufunc import _copyFromAndConvert
def _noConversionNeeded(x, xt):
"""returns 1 iff x is a c_array of type xt."""
if isinstance(x, _nt.scalarTypes):
return True
return x.is_c_array() and (x._type.name == xt or x._type == xt)
def _inputcheck(*inargs):
"""Check input args for type, convert sequences to numnumarray
Return list of 1) input arg tuple with scalars and numnumarray only
2) signature string indicating which args are scalars
and which are numnumarray
3) flag indicating if only scalar arguments are present
"""
retargs = []
retsig = []
scalar = 1
for inarg in inargs:
if isinstance(inarg, _nc.NumArray):
scalar = 0
retsig.append('v')
elif _nc.PyNUMERIC_TYPES.has_key(type(inarg)):
retsig.append('s')
else:
scalar = 0
retsig.append('v')
# see if it can be made into an array
try:
inarg = _nc.array(inarg)
except TypeError:
raise TypeError(
"UFunc arguments must be numarray, scalars or numeric sequences")
retargs.append(inarg)
return retargs, ''.join(retsig), scalar
def _pseudoArrayFromScalars(scalarvalues, type):
"""Wrap a scalar in a buffer so it can be used as an array"""
arr = _bufferPool.getBuffer()
arr._check_overflow = 1
newtype = type # _numtypedict[type]
arr._strides = (newtype.bytes,)
arr._type = newtype
arr._itemsize = newtype.bytes
arr._strides = None
if isinstance(scalarvalues, (list, tuple)):
arr._shape = (len(scalarvalues),)
for i in xrange(len(scalarvalues)):
arr[i] = scalarvalues[i]
else:
arr._shape = ()
arr[()] = scalarvalues
# Modify block buffer attributes to look like vector/vector setup.
return arr
def _normalize_results(inputs, outputs, results, return_rank1=False):
"""normalize_results() does standard return value handling for
a given set of ufunc actual parameters and results. Ufunc inputs
are assumed to be arrays or scalars, ufunc outputs can be arrays
or unspecified, and ufunc "preliminary results" must be arrays.
The standard processing assumptions and rules are summarized as
follows:
1. All ufunc inputs must be specified with *some* value.
2. If any outputs were specified, all outputs must be specified.
If no outputs were specified, outputs is None or (). Specified
outputs are always arrays.
3. If any outputs were specified, the return value is supressed
(is None).
4. If no outputs were specified, at least one array value will be
returned. All results are assumed to have the same shape.
a. If len(results)==1, the tuple is discarded and the return value
is a single array or scalar.
b. If len(results) > 1, the return value is a tuple.
5. If the result array(s) are rank-0, either a scalar, rank-0, or
rank-1 array is returned.
a. If at least 1 input is an array, the result will be an array.
b. If all inputs are scalars, the result(s) is(are) a scalar.
c. If return_rank1 is True, rank-0 results are converted to rank-1-len-1.
6. Arrays of dimension > 0 have no special shape handling.
"""
# catch "void" ufuncs here just in case, also specified outputs => None
if len(results) == 0 or (outputs is not None and outputs is not ()):
return None
if results[0].rank == 0: # scalar, rank-0, or rank-1 results
for i in inputs: # any rank0 input implies rank0 output
if isinstance(i, _nc.NumArray) and i.rank==0:
if return_rank1: # convert rank0 arrays to rank1
for result in results:
result.shape = (1,)
break
else: # all scalar inputs implies scalar outputs
results = list(results)
for i in range(len(results)):
results[i] = results[i][()]
results = tuple(results)
else: # rank >= 1 arrays are returned unaltered
pass
if len(results) == 1:
# len-1 results converts to the result itself
results = results[0]
return results
if not _PROTOTYPE:
from numarray._ufunc import _normalize_results
# Portions of the UFunc classes have been re-implemented in C.
class _UFunc(_ufunc._ufunc):
"""Class to handle all element-by-element mathematical functions
This is the base class -- subclasses must implement __call__ and
(if allowed) the outer and _cumulative methods.
"""
def __init__(self, operator, ufuncs, inputs, outputs, identity):
self._cfuncs = self._organize_cfuncs(ufuncs)
self._coercions = _nt.genericCoercions
self._promotionExclusions = _nt.genericPromotionExclusions
self._typerank = _nt.genericTypeRank
self._cache = {}
def _evaltypes(self, intypes):
return tuple([ eval("_nt." + t) for t in intypes ])
def _sort_typesigs(self, cfuncs):
l = []
for intypes, rest in cfuncs.items(): # input signatures
it = self._evaltypes( intypes )
l.append((it,)+rest)
l.sort() # intype tuples should sort in promotion order!?
return l
def _organize_cfuncs(self, ufuncs):
"""Replace unordered cfunc dictionary with list ordered by
input signature type comparison."""
ufs = {}
for form in ufuncs:
ufs[form] = self._sort_typesigs(ufuncs[form])
return ufs
def _getidentity(self):
return _nc.array([self._identity])
if _PROTOTYPE:
def _cache_flush(self):
# print "flushing cache for", self
self._cache = {}
def __call__(self, *args):
raise NotImplementedError("__call__ is not implemented by base UFunc class")
def __repr__(self):
return "<UFunc: '%s'>" % self.operator
def _scan_sigs(self, types, binTypeList):
"""scan_sigs matches a list of input types against a sorted list of
signature,cfunc tuples.
"""
inplen = len(types)
for bt in binTypeList:
btsig = bt[0]
if inplen != len(btsig): # same parameter count
continue
for i in range(inplen): # find first sig >= types
t = types[i]
if (isinstance(t, _nt.NumericType) and
(t > btsig[i] or
btsig[i].name in self._promotionExclusions[t])):
break
else: # entire sig "fits", first match.
return bt
raise TypeError("Couldn't find a function which matches the inputs")
def _typematch_N(self, types, form):
"""Determine which C function is most appropriate for the given
types. Returns a new form (possibly the all-vector default form),
the input types to be converted to, the output types, and the
c function.
"""
for t in types:
assert isinstance(t, _nt.NumericType) or t in _nt.scalarTypes
assert isinstance(form, str)
try: # Look for an exact match of the current form
form2 = form
binTypeList = self._cfuncs[ form2 ]
except KeyError: # Look for an all "vector" standby
form2 = "v".join(form.split("s"))
binTypeList = self._cfuncs[ form2 ]
try: # Do "binary" ufunc type coercions
intype = _nt.getType(self._coercions[ tuple(types) ])
types1 = (intype,)*len(types)
except KeyError:
types1 = types
try: # Attempt with coerced types...
return (form2,) + self._scan_sigs(types1, binTypeList)
except TypeError: # fail back to non-coerced match
if form == "s"*len(types): # all scalars --> scalars must match
types = [ eval("_nt."+_nt.scalarTypeMap[t]) for t in types ]
return (form2,) + self._scan_sigs(types, binTypeList)
def outer(self, inarr1, inarr2, outarr=None):
raise ValueError("Outer only available for binary functions")
_blockingParametersCache = {}
_copyCache = {}
def flush_caches():
"""flush_cache elminates all ufunc cache entries."""
import sys
global _blockingParametersCache, _copyCache
module = sys.modules[__name__]
for n in module.__dict__.keys():
f = module.__dict__[n]
if isinstance(f, _UFunc):
f._cache_flush()
_blockingParametersCache = {}
_copyCache = {}
class _BufferPool:
"""A Pool of available buffers available for block compuations
Buffers can be gotten. When they are deleted they are automatically
returned to the buffer pool so they can be reused.
"""
def __init__(self, buffersize=100000):
# buffers are initially allocated as needed
self.buffers = []
self.setBufferSize(buffersize)
def getBuffer(self):
try:
return _nc._UBuffer(self.buffers.pop())
except IndexError:
# out of buffers, create a new one
return _nc._UBuffer(memory.new_memory(_ufunc.getBufferSize()))
def getBufferSize(self):
return _ufunc.getBufferSize()
def setBufferSize(self, buffersize):
"""Change buffer size (useful mainly for testing)"""
oldsize = _ufunc.getBufferSize()
if buffersize != oldsize:
_ufunc.setBufferSize(buffersize)
# delete existing (wrong-sized) buffers
flush_caches()
self.buffers = []
return oldsize
_bufferPool = _BufferPool()
def _sequence(x):
return isinstance(x, list) or isinstance(x, tuple)
if _PROTOTYPE:
def _digest(x):
if isinstance(x, _nc.NumArray):
return ("array",
x.iscontiguous(),
x.isaligned(),
x.isbyteswapped(),
x.type())
elif _isScalar(x):
return type(x), x # for demo only
elif x is None:
return None
else:
raise KeyError, "_digest force cache miss"
else:
_digest = _ufunc.digest
def _restuff_pseudo(p, x):
if _isScalar(x):
if p is None:
raise RuntimeError, "scalar value with no pseudo array"
else: # restuff p with x
p[0] = x
return p
else:
if p is not None:
return p
else:
return x
_restuff_pseudo = _ufunc.restuff_pseudo
def _cache(original_input, pseudo_array):
if _isScalar(original_input):
return pseudo_array
else:
return None
def _rank(seq, r=0):
if isinstance(seq, _nc.NumArray):
return r + len(seq._shape)
try:
x = len(seq)
except:
return r
else:
if x == 0:
return r
else:
return _rank(seq[0], r+1)
class _UnaryUFunc(_UFunc):
"""Class for ufuncs with 1 input and 1 output argument"""
if _PROTOTYPE:
def __call__(self, inarr1, outarr=None):
"""The standard calling interface for UFuncs"""
return self._cached_dispatch1(inarr1, outarr)
def _fast_exec1(self, in1, out, cached):
mode, outtype, cfunc, ufargs_junk, inform, cin1 = cached
ufargs = (out.nelements(), 1, 1,
((in1._data, in1._byteoffset),
(out._data, out._byteoffset)))
apply(cfunc, ufargs)
def _slow_exec1(self, in1, out, cached):
mode, outtype, cfunc, ufargs, inform, cin1 = cached
inputs, outputs, maxitemsize = ufargs
niter = _ufunc.getBufferSize()/maxitemsize
indexlevel, blockingparameters = \
_ufunc._getBlockingParameters(out._shape, niter)
operator = _Operator(cfunc,
[inputs[0].rebuffer(in1)],
[outputs[0].rebuffer(out)])
objects = (inputs[0],operator,outputs[0])
_callOverDimensions(objects, out._shape, indexlevel, blockingparameters)
inputs[0].clean(in1)
outputs[0].clean(out)
def _cache_exec1(self, in1, out, cached):
if out.nelements(): # skip 0-element arrays
_ufunc.CheckFPErrors()
if cached[0] == "fast":
self._fast_exec1(in1, out, cached)
else:
self._slow_exec1(in1, out, cached)
error = _ufunc.CheckFPErrors()
if error:
handleError(error, " in "+self.operator)
return out
def _cache_lookup1(self, in1, out):
"""_cache_lookup1 checks the ufunc cache for an entry corresponding to
(in1, in2, out) and returns it if it is found. If the cache misses,
_cache_lookup1 performs a full ufunc setup and returns it. In either
case, vector_vector inputs are dual broadcast to match shapes.
"""
win1 = in1
try:
key = (_digest(win1), _digest(out), safethread.get_ident())
cached = self._cache[ key ]
except KeyError:
return self._cache_miss1(win1, out)
else:
mode, otype, cfunc, ufargs, inform, cin1 = cached
cin1 = _restuff_pseudo(cin1, in1)
if out is None:
if inform == "v":
wout = win1.new(otype)
else: # scalar
wout = _nc.zeros((), type=otype)
else:
if inform == "v":
shape = win1._shape
else: # scalar
shape = ()
if shape != out._shape:
raise ValueError("Supplied output array does not have appropriate shape")
wout = out
return cin1, wout, cached
def _cached_dispatch1(self, inarr, outarr):
params = self._cache_lookup1(inarr, outarr)
result = self._cache_exec1(*params)
return _normalize_results((inarr,),outarr,(result,),False)
def _cache_miss1(self, n1, out):
(in1,), inform, scalar = _inputcheck(n1)
mode, win1, wout, cfunc, ufargs = \
self._setup(in1, inform, out)
cached = mode, wout._type, cfunc, ufargs, inform, _cache(n1, win1)
try:
key = (_digest(n1), _digest(out), safethread.get_ident())
except KeyError:
pass
else:
if _PROTOTYPE:
self._cache[ key ] = cached
else:
self._cache_insert( cached, n1, None, out)
return win1, wout, cached
def _setup(self, inarr, inform, outarr):
"""Setup for unary ufunc"""
if inform == "s": # scalar
t = _nc.array(inarr) # create rank-0
t._check_overflow = 1 # enable overflow checking for t
t[()] = inarr # force any overflow
inarr = t
shape = inarr._shape
fform, convtypes, outtypes, cfunc = self._typematch_N((inarr._type,), inform)
intype = convtypes[0]
if outarr is None:
outarr = inarr.new(_numtypedict[outtypes[0]])
if shape != outarr._shape:
raise ValueError("Supplied output array does not have appropriate shape")
# If fast case, return the arguments for an immediate call to the
# C function
if (_noConversionNeeded(inarr, intype) and
_noConversionNeeded(outarr, outtypes[0])):
return ("fast", inarr, outarr, cfunc, None)
# slower, general case
inputs = (_InputConverter(inarr, intype), )
outputs = (_OutputConverter(arr=outarr, type=outtypes[0]),)
# find largest itemsize involved, it will determine how many
# iterations can be done on the temporary buffers (i.e., niter)
maxitemsize = max(_numtypedict[inarr._type.name].bytes,
_numtypedict[intype].bytes,
_numtypedict[outtypes[0]].bytes,
_numtypedict[outarr._type.name].bytes)
preprocessing_outbuffers = (inputs[0].result_buff,)
postprocessing_inbuffers = (outputs[0].result_buff,)
return ("slow", inarr, outarr, cfunc,
(inputs, outputs, maxitemsize))
def _firstcol(arr):
rval = arr.view()
rval._shape = rval._shape[:-1]
rval._strides = rval._strides[:-1]
return rval
def _moveToLast(dim, iseq):
oseq = []
for i in range(len(iseq)):
if i != dim:
oseq.append(iseq[i])
else:
s = iseq[i]
oseq.append(s)
return oseq
def _fixdim(axis, dim):
"""support the deprecation of the 'dim' keyword in favor of 'axis'"""
if dim is not None:
if axis != 0:
raise RuntimeError("Specify 'axis' or 'dim', but not both. 'dim' is deprecated.")
_warnings.warn("The 'dim' keyword is deprecated. Specify 'axis' instead.",
DeprecationWarning,
stacklevel=3)
return dim
return axis
class _BinaryUFunc(_UFunc):
"""Class for ufuncs with 2 input and 1 output arguments"""
if _PROTOTYPE:
def _fast_exec2(self, in1, in2, out, cached):
mode, outtype, cfunc, ufargs_junk, inform, cin1, cin2 = cached
ufargs = (out.nelements(), 2, 1,
((in1._data, in1._byteoffset),
(in2._data, in2._byteoffset),
(out._data, out._byteoffset)))
apply(cfunc, ufargs)
def _slow_exec2(self, in1, in2, out, cached):
mode, outtype, cfunc, ufargs, inform, cin1, cin2 = cached
inputs, outputs, maxitemsize = ufargs
niter = _ufunc.getBufferSize()/maxitemsize
indexlevel, blockingparameters = \
_ufunc._getBlockingParameters(out._shape, niter)
operator = _Operator(cfunc,[inputs[0].rebuffer(in1),
inputs[1].rebuffer(in2)],
[outputs[0].rebuffer(out)])
objects = inputs + (operator,) + outputs
_callOverDimensions(objects, out._shape, indexlevel,
blockingparameters)
inputs[0].clean(in1)
inputs[1].clean(in2)
outputs[0].clean(out)
def _cache_exec2(self, in1, in2, out, cached):
if out.nelements(): # skip 0-element arrays
_ufunc.CheckFPErrors()
if cached[0] == "fast":
self._fast_exec2(in1, in2, out, cached)
else:
self._slow_exec2(in1, in2, out, cached)
error = _ufunc.CheckFPErrors()
if error:
handleError(error, " in "+self.operator)
return out
def _cache_lookup2(self, in1, in2, out):
"""_cache_lookup2 checks the ufunc cache for an entry corresponding to
(in1, in2, out) and returns it if it is found. If the cache misses,
_cache_lookup2 performs a full ufunc setup and returns it. In either
case, vector_vector inputs are dual broadcast to match shapes.
"""
if (isinstance(in1, _nc.NumArray) and
isinstance(in2, _nc.NumArray)):
win1, win2 = in1._dualbroadcast(in2)
else:
win1, win2 = in1, in2
try:
key = (_digest(win1), _digest(win2), _digest(out), safethread.get_ident())
cached = self._cache[ key ]
except KeyError:
return self._cache_miss2(win1, win2, out)
else:
mode, otype, cfunc, ufargs, inform, cin1, cin2 = cached
cin1 = _restuff_pseudo(cin1, win1)
cin2 = _restuff_pseudo(cin2, win2)
if out is None:
if inform in ["vs", "vv"]: # vector_scalar, vector_vector
wout = win1.new(cached[1])
elif inform == "sv": # scalar_vector
wout = win2.new(cached[1])
else: # scalar_scalar
wout = _nc.zeros((), type=otype)
else:
if inform in ["vs", "vv"]:
shape = win1._shape
elif inform == "sv":
shape = win2._shape
else: # scalar_scalar
shape = ()
if shape != out._shape:
raise ValueError("Supplied output array does not have appropriate shape")
wout = out
return cin1, cin2, wout, cached
def _cached_dispatch2(self, inarr1, inarr2, outarr):
params = self._cache_lookup2(inarr1, inarr2, outarr)
result = self._cache_exec2(*params)
return _normalize_results((inarr1, inarr2), outarr, (result,), False)
def __call__(self, inarr1, inarr2, outarr=None):
"""The standard calling interface for UFuncs"""
return self._cached_dispatch2(inarr1, inarr2, outarr)
def _cache_miss2(self, n1, n2, out):
(in1, in2), inform, scalar = _inputcheck(n1, n2)
mode, win1, win2, wout, cfunc, ufargs = \
self._setup(in1, in2, inform, out)
cached = mode, wout._type, cfunc, ufargs, inform, \
_cache(n1, win1), _cache(n2, win2)
try:
key = (_digest(n1), _digest(n2), _digest(out),
safethread.get_ident())
except KeyError:
pass
else:
if _PROTOTYPE:
self._cache[ key ] = cached
else:
self._cache_insert(cached, n1, n2, out)
return win1, win2, wout, cached
def _setup(self, in1, in2, inform, out):
"""Setup for binary ufunc"""
wout = out
if inform == "ss":
# do it as two 0-d numarray.
in1, in2 = _nc.array(in1), _nc.array(in2)
inform = "vv"
if out is None: wout = in1.copy()
elif inform == "vs":
intypes = (in1._type, type(in2))
fform, convtypes, outtypes, cfunc = self._typematch_N(intypes, inform)
if fform == "vs":
inarr1, inarr2 = in1, _pseudoArrayFromScalars(in2, convtypes[1])
else:
in1, in2 = in1, _nc.array(in2, type=convtypes[1])
inform = "vv"
if out is None: wout = in1.new(outtypes[0])
elif inform == "sv":
intypes = (type(in1), in2._type)
fform, convtypes, outtypes, cfunc = self._typematch_N(intypes, inform)
if fform == "sv":
inarr1, inarr2 = _pseudoArrayFromScalars(in1, convtypes[0]), in2
else:
in1, in2 = _nc.array(in1, type=convtypes[0]), in2
inform = "vv"
if out is None: wout = in2.new(outtypes[0])
if inform == "vv":
intypes = (in1._type, in2._type)
inarr1, inarr2 = in1._dualbroadcast(in2)
fform, convtypes, outtypes, cfunc = self._typematch_N(intypes, inform)
if out is None: wout = inarr1.new(outtypes[0])
#If fast case, return the arguments for an immediate call to the
# C function
if (_noConversionNeeded(inarr1, convtypes[0]) and
_noConversionNeeded(inarr2, convtypes[1]) and
_noConversionNeeded(wout, _numtypedict[outtypes[0]])):
return ("fast", inarr1, inarr2, wout, cfunc, None)
# slower, general case
inputs = (_InputConverter(inarr1, convtypes[0]),
_InputConverter(inarr2, convtypes[1]))
outputs = (_OutputConverter(arr=wout, type=outtypes[0]),)
# find largest itemsize involved, it will determine how many
# iterations can be done on the temporary buffers (i.e., niter)
maxitemsize = max(_numtypedict[inarr1._type.name].bytes,
_numtypedict[inarr2._type.name].bytes,
_numtypedict[convtypes[0]].bytes,
_numtypedict[convtypes[1]].bytes,
_numtypedict[outtypes[0]].bytes,
_numtypedict[wout._type.name].bytes)
return ("slow", inarr1, inarr2, wout, cfunc,
(inputs, outputs, maxitemsize))
def outer(self, inarr1, inarr2, outarr=None):
"""Return outer product of 2 numarray"""
(inarr1, inarr2), inform, scalar = _inputcheck(inarr1, inarr2)
ndim2 = len(inarr2._shape)
indexarg = (slice(None, None, None),)*len(inarr1._shape) + (None,)*ndim2
if outarr is not None:
self(inarr1.__getitem__(indexarg), inarr2, outarr)
else:
return self(inarr1.__getitem__(indexarg), inarr2)
def zreduce(self, array, axis=0, out=None, type=None, dim=None):
"""zreduce returns a rank-0 array as the result of the
reduction of array <= rank-1 array, similar to Numeric.
"""
axis = _fixdim(axis, dim)
in1 = _nc.asarray(array)
r = self._cum_swapped(in1, axis, out, "R", type)
if in1.rank <= 1:
r.shape = ()
r.strides = ()
return r
if _PROTOTYPE:
def accumulate(self, array, axis=0, out=None, type=None, dim=None):
"""accumulate applies the binary operator 'self' at successive
pairs of elements of 'array' along dimension 'dim', storing the
result in 'out'. If no 'out' is provided, the result of the
accumulation will be returned. The result of accumulating
'array' has the same shape as 'array'.
"""
axis = _fixdim(axis, dim)
in1 = _nc.asarray(array)
return self._cum_swapped(in1, axis, out, "A", type)
def areduce(self, array, axis=0, out=None, type=None,
dim=None):
"""areduce applies the operator of ufunc 'self' to reduce
'array' along its 'dim' axis. If no 'out' is provided,
the result of the reduction is returned. The result of
reducing an N-dimensional array is an N-1 dimensional
array. The result of reducing a rank-1 array is *still* a
rank-1 array.
see also 'reduce'.
"""
axis = _fixdim(axis, dim)
in1 = _nc.asarray(array)
return self._cum_swapped(in1, axis, out, "R", type)
def reduce(self, array, axis=0, out=None, type=None, dim=None):
"""areduce applies the operator of ufunc 'self' to reduce
'array' along its 'dim' axis. If no 'out' is provided,
the result of the reduction is returned. The result of
reducing an N-dimensional array is an N-1 dimensional array.
The result of reducing a rank-1 array is a scalar.
see also 'areduce'.
"""
axis = _fixdim(axis, dim)
in1 = _nc.asarray(array)
r = self._cum_swapped(in1, axis, out, "R", type)
if r is not None and in1.rank <= 1:
if r._shape == (1,):
r = r[0]
elif r._shape == ():
r = r[()]
return r
def _cum_cached(self, cumop, in1, out, type):
if (out is not None):
if not isinstance(out, _nc.NumArray):
raise TypeError("output array must be a NumArray")
if (not out.isaligned() or out.isbyteswapped()):
raise ValueError("Reduce/Accumulate: no support for misaligned / byteswapped output numarray")
if cumop not in ["R", "A"]:
raise ValueError("Unknown cumulative option")
if in1.rank == 0:
if out is None:
out = in1.astype(type)
out[()] = in1[()]
return out
params = self._cum_lookup(cumop, in1, out, type)
wout = self._cum_exec(*params)
if cumop == "R":
wout._shape = in1._shape[:-1]
wout._strides = wout._stridesFromShape()
if wout._shape == ():
wout._shape = (1,)
wout._strides = (wout._itemsize,)
if out is None:
return wout
i, o, cached = params
mode, otype, cfunc, ufargs = cached
if out._type != otype:
_copyFromAndConvert(wout, out)
return out
def _cum_lookup(self, cumop, in1, out, type):
try:
key = (_digest(in1), _digest(out), cumop, safethread.get_ident(), type)
cached = self._cache[ key ]
except KeyError: # lists always miss
return self._cum_cache_miss(cumop, in1, out, type)
else:
mode, otype, cfunc, ufargs = cached
if cumop == "R":
wout = self._reduce_out(in1, out, otype)
else:
wout = self._accumulate_out(in1, out, otype)
return in1, wout, cached
def _cum_exec(self, in1, out, cached):
mode, otype, cfunc, ufargs = cached
if otype == _nt.Bool:
if in1._type != _nt.Bool:
in1 = in1.astype(_nt.Bool)
if in1.nelements(): # skip 0-element arrays
_ufunc.CheckFPErrors()
if mode == "fast":
self._cum_fast_exec(in1, out, cached)
else:
self._cum_slow_exec(in1, out, cached)
errorstatus = _ufunc.CheckFPErrors()
if errorstatus:
handleError(errorstatus,
"in cumulative application of " + self.operator)
return out
def _cum_fast_exec(self, in1, out, cached):
mode, otype, cfunc, ufargs = cached
cfunc(in1._shape,
in1._data, in1._byteoffset, in1._strides,
out._data, out._byteoffset, out._strides)
def _cum_slow_exec(self, in1, out, cached):
mode, otype, cfunc, ufargs = cached
input, output, maxitemsize = ufargs
niter = _ufunc.getBufferSize()/maxitemsize
if in1._shape[-1] > niter:
overlap = 1
else:
overlap = 0
indexlevel, blockingparameters = \
_ufunc._getBlockingParameters(in1._shape, niter, overlap)
operator = _Operator(cfunc,
[input.rebuffer(in1)],
[output.rebuffer(out)],
otype.bytes)
objects = (input, operator, output)
_callOverDimensions(objects, out._shape,
indexlevel, blockingparameters, overlap)
input.clean(in1)
output.clean(out)
def _accumulate_out(self, inarr, outarr, outtype):
"""_accumulate_out creates output array for accumulate"""
# Create output array if not supplied.
if outarr is None:
toutarr = inarr.new( outtype )
else:
if outarr._shape != inarr._shape:
raise ValueError(
"Supplied output array does not have the appropriate shape")
toutarr = outarr
# It is necessary to initialize the first subarray of the output
# to the first subarray of the input (C functions depend on it).
# toutarr[...,0] = inarr[...,0]
if inarr.nelements():
_firstcol(toutarr)._copyFrom(_firstcol(inarr))
return toutarr
def _reduce_out(self, inarr, outarr, outtype):
"""_reduce_out creates output array for reduce"""
firstcol = _firstcol(inarr) # inarr[..., 0]
# Create output array if not supplied or wrong type
if outarr is None or outarr._type != outtype:
toutarr = firstcol.new(outtype)
else:
toutarr = outarr
if len(inarr._shape) == 1 and toutarr._shape == (1,):
toutarr._shape, toutarr._strides = (), ()
if firstcol._shape != toutarr._shape:
raise ValueError("output array shape mismatch")
if inarr.nelements(): # skip 0-element arrays
toutarr._copyFrom(firstcol) # Assign first subarray of output
else:
toutarr._copyFrom(self._identity)
toutarr._strides += (0,)
toutarr._shape = inarr._shape
diff = len(toutarr._strides) - len(toutarr._shape)
# if diff > 0:
# toutarr._strides = toutarr._strides[diff:]
return toutarr
def _cum_swapped(self, in1, dim, out, cumop, type=None):
if in1.rank == 0:
return in1.copy()
if dim == in1.rank-1:
dim = -1
if dim != -1:
_in1 = _gen.swapaxes(in1, -1, dim)
_out = _gen.swapaxes(out, -1, dim)
else:
_in1, _out = in1, out
_out1 = self._cum_cached(cumop, _in1, _out, type)
if cumop == "A":
if dim != -1:
_out1.swapaxes(-1, dim) # swap axes to "fix" the result
else:
if dim != -1:
if in1.rank:
_out1._shape = _moveToLast(dim, _out1._shape)
_out1._strides = _moveToLast(dim, _out1._strides)
if out is None:
return _out1
def _cum_cache_miss(self, cumop, in1, out, type):
"""computes the setup values and caches them."""
mode, win1, wout, cfunc, ufargs = \
self._cum_setup(cumop, in1, out, type)
cached = mode, wout._type, cfunc, ufargs
try:
key = (_digest(in1), _digest(out), cumop,
safethread.get_ident(), type)
except KeyError:
pass
else:
if _PROTOTYPE:
self._cache[ key ] = cached
else:
self._cache_insert(cached, in1, None, out, cumop, type);
return win1, wout, cached
'''
try:
intype = _numtypedict[self._typePromoter(win1._type.name, typekeys)]
except TypeError:
assert 0
otype, cfunc = ufdict[(optype,)]
'''
def _cum_setup(self, cumop, in1, out, type=None):
"""Used by both reduce and accumulate to compute accumulations.
Assumes that the dimension of accumulation/reduction is 0.
in1 -- The input array to accumulate/reduce
out -- Optional output array
The output array must have a consistent shape with the input
Array (the same for accumulate and minus the accumulated dimension
for reduce).
The output type must match that produced by default
and cannot be nonaligned or byteswapped. (It is a real mess to
remove this restriction. (The original numeric doesn't even support
output numarray.)
"""
(win1,), inform, scalar = _inputcheck(in1) # win1 substitute for lists
if scalar:
raise ValueError("First argument must be an array or sequence")
# The original numeric returns the operations identity for empty numarray.
sourcemsg = "in %s.%s" % (self.operator, cumop) # For error message
# Get appropriate function, figure out any needed input type conversions
# and what the output type will be
intype = win1._type
if out is None:
if type is None:
optype = intype
else:
optype = _nt.getType(type)
else:
if type is not None:
raise ValueError("Can't specify both 'type' and 'out'")
else:
optype = out.type()
try:
fform, convtypes, otypes, cfunc = self._typematch_N(
(optype,), cumop)
except TypeError: # Hack for logical operators
ufdict = self._cfuncs[cumop]
typekeys = [tupl[0] for tupl in ufdict]
if len(typekeys) == 1 and typekeys[0] == (_nt.Bool,):
return self._cum_setup( cumop, win1.astype('Bool'), out)
else:
raise
otype = _numtypedict[otypes[0]]
if out is None and type is None:
optype = otype
# A few words on shape/stride manipulations:
# The C function assumes the last dimension is being accumulated.
# Create the accumulate/reduce specific output array, tout
# type of reduce tout is always otype to avoid output conversion
# with mismatched shapes.
if cumop == "R":
tout = self._reduce_out(win1, out, otype)
else:
tout = self._accumulate_out(win1, out, otype)
# check for fast case (Unlike Ufuncs, contiguous numarray not required)
if (not win1.isbyteswapped() and win1.isaligned()
and otype == intype
and not tout.isbyteswapped()
and tout.isaligned()):
return "fast", win1, tout, cfunc, None
else:
# slow case
# And now a few words about overlap. It is used to allow the
# accumulate C function to work across blocks. Overlap in effect
# makes the first subarray of the next block the last subarray
# of the previous block (or last array value for 1-d numarray) so
# the accumulate can carry on. The computation of the block size
# must take it into account as does computeOverDimensions in its
# computation of the offsets for the blocks. This is only needed
# where the accumulated dimension is larger than a block. It is
# not needed at all for reduce since the output array contains
# the net accumulation.
# check special case where last dimension is bigger than blocksize
maxitemsize = max(intype.bytes, optype.bytes, otype.bytes)
niter = _ufunc.getBufferSize()/maxitemsize
# This is a bit tricky. If the accumulated dimension is bigger
# than the block size, it is necessary to overlap the block
# computations by one element so that the last value of the
# previous block can continue the count into the next block.
if win1._shape[-1] > niter:
overlap = 1
else:
overlap = 0
indexlevel, blockingparameters = \
_ufunc._getBlockingParameters(win1._shape, niter, overlap)
# The usual compute stuff
input = _InputConverter(win1, type=convtypes[0].name)
# Since cfunc is selected by output type, no type conversion
# is required. Since the cfuncs are striding, no striding
# is required. Byteswapping and aligning are disallowed.
# Thus, the output converter is a NULL converter.
output = _OutputConverter(tout, type=optype.name, nonstriding=1)
return "slow", win1, tout, cfunc, (input, output, maxitemsize)
class _CacheEntry:
def __init__(self, **keys):
self.__dict__.update(keys)
class _NaryUFunc(_UFunc):
"""Class for ufuncs with M input and N output arguments"""
def __init__(self, operator, ufuncs, inputs, outputs, identity=None):
_UFunc.__init__(self, operator, ufuncs, inputs, outputs, identity)
self._cfuncs = self._organize_cfuncs(ufuncs)
self._coercions = _nt.genericCoercions
self._promotionExclusions = _nt.genericPromotionExclusions
self._typerank = _nt.genericTypeRank
self._cache = {}
def _cache_flush(self):
self._cache = {}
def __repr__(self):
return "_NaryUfunc(%s, inputs=%d, outputs=%d)" % \
(self.operator, self.n_inputs, self.n_outputs)
def __call__(self, *args):
"""The standard calling interface for UFuncs"""
if (len(args) < self.n_inputs or
self.n_inputs+self.n_outputs < len(args)):
raise ValueError(
"Ufunc %s takes %d inputs and %d outputs "
"but %d parameters given." %
(self.operator, self.n_inputs, self.n_outputs, len(args)))
if self.n_inputs < len(args) < self.n_inputs + self.n_outputs:
raise ValueError("Ufunc %s takes %d inputs and %d outputs; "
"if any outputs are specified, all outputs "
"must be specified." %
(self.operator, self.n_inputs, self.n_outputs))
inputs = tuple(args[:self.n_inputs])
outputs = tuple(args[self.n_inputs:])
params = self._cache_lookup(inputs, outputs)
results = self._cache_exec(*params)
return _normalize_results(inputs, outputs, results, False)
def _fast_path(self, parameters):
"""_fast_path decides if a set of parameters is suitable for
small array performance optimiation.
The fast path is as follows:
1. Well-behaved NumArray and scalar parameters are supported
2. Broadcast arrays are not supported
3. Misbehaved arrays are not supported
4. Type converted arrays are not supported
"""
shape = None
type = None
for p in parameters:
if isinstance(p, _nc.NumArray):
if shape is None:
shape = p._shape
else:
if p._shape != shape:
return False, None
if type is None:
type = p._type
else:
if type != p._type:
return False, None
if not p.is_c_array():
return False, None
elif not _isScalar(p):
return False, None
else:
if shape is None:
shape = ()
return True, shape
def _cache_lookup(self, inputs, outputs):
"""_cache_lookup checks the ufunc cache for an entry
corresponding to (inputs... outputs...) and returns it if it
is found. If the cache misses, _cache_lookup performs a full
ufunc setup and returns it. In either case, inputs are nWay
broadcast to maintain shape.
"""
fast, shape = self._fast_path(inputs + outputs)
if fast:
try:
key = self._key(inputs, outputs) # some inputs fail to key
cached = self._cache[ key ]
except KeyError:
winputs, woutputs, cached = self._cache_miss(inputs, outputs)
key = self._key(inputs, outputs)
self._cache[ key ] = cached
else:
winputs = self._manage_cached_inputs(cached.cinputs, inputs)
woutputs = self._manage_outputs(
shape, cached.outtypes, outputs)
else:
winputs, woutputs, cached = self._cache_miss(inputs, outputs)
return winputs, woutputs, cached
def _key(self, inputs, outputs):
if outputs is None:
outputs = [ None ]
digests = [ _digest(x) for x in inputs + outputs ]
return tuple([safethread.get_ident()] + digests)
def _manage_cached_inputs(self, cinputs, winputs):
# stuff pseudo buffers with current scalar values
cinputs = list(cinputs)
for i in range(len(winputs)):
cinputs[i] = _restuff_pseudo(cinputs[i], winputs[i])
return tuple(cinputs)
def _manage_outputs(self, shape, outtypes, outputs):
if outputs == (): # create output arrays
woutputs = []
for ot in outtypes:
woutputs.append(_nc.NumArray(shape=shape, type=ot))
else: # check specified output arrays
if len(outputs) != len(outtypes):
raise RuntimeError("Wrong number of output arrays. Either supply no output arrays, or supply all output arrays.")
for out in outputs:
if shape != out._shape:
raise ValueError("Supplied output array does not"
" have appropriate shape")
woutputs = outputs
return tuple(woutputs)
def _cache_miss(self, inputs, outputs):
mode, outtypes, winputs, woutputs, cfunc, ufargs = \
self._setup(inputs, outputs)
cached_inputs = [_cache(inputs[i], winputs[i]) \
for i in range(self.n_inputs)]
cached = _CacheEntry(mode=mode, outtypes=outtypes, cfunc=cfunc,
ufargs=ufargs, cinputs=cached_inputs)
return winputs, woutputs, cached
def _setup(self, inputs, outputs):
"""Setup for nary ufunc"""
inputs, inform, scalar = _inputcheck(*tuple(inputs))
intypes = []
for i in inputs:
if isinstance(i, _nc.NumArray):
intypes.append(i._type)
else:
intypes.append(type(i))
inform1, convtypes, outtypes, cfunc = \
self._typematch_N(intypes, inform)
for i in range(len(inform1)):
if isinstance(inputs[i], _nt.scalarTypes):
if inform1[i] == "s":
inputs[i] = _pseudoArrayFromScalars(
[inputs[i]], type=convtypes[i])
elif inform1[i] == "v":
inputs[i] = _nc.array(inputs[i], type=convtypes[i])
else:
raise RuntimeError("Unexpected array form.")
common_shape = _gen._common_shapes(inputs)
winputs = _gen._broadcast_all(inputs, common_shape)
woutputs = self._manage_outputs(common_shape, outtypes, outputs)
#If fast case, return the arguments for an immediate call to the
# C function
types = convtypes + outtypes
parameters = list(winputs) + list(outputs)
for i in range(len(parameters)):
if not _noConversionNeeded(parameters[i], types[i]):
break
else:
return ("fast", outtypes, winputs, woutputs, cfunc, None)
# slower, general case
input_convs = tuple([ _InputConverter(a, ctype) for a,ctype in
zip(winputs, convtypes) ])
output_convs = tuple([_OutputConverter(a, ctype) for a,ctype in
zip(woutputs, outtypes) ])
# find largest itemsize involved, it will determine how many
# iterations can be done on the temporary buffers (i.e., niter)
types = tuple([ a._type.name for a in winputs + woutputs ]) + \
convtypes + outtypes
maxitemsize = 0
for t in types:
tsize = _numtypedict[t].bytes
if tsize > maxitemsize:
maxitemsize = tsize
return ("slow", outtypes, winputs, woutputs, cfunc,
(input_convs, output_convs, maxitemsize))
def _cache_exec(self, inputs, outputs, cached):
if outputs[0].nelements(): # skip 0-element arrays
_ufunc.CheckFPErrors()
if cached.mode == "fast":
self._fast_exec(inputs, outputs, cached)
else:
self._slow_exec(inputs, outputs, cached)
error = _ufunc.CheckFPErrors()
if error:
handleError(error, " in "+self.operator)
return outputs
def _fast_exec(self, inputs, outputs, cached):
arrays = inputs + outputs
buffers = []
ne = 1
for i in range(self.n_inputs+self.n_outputs):
a = arrays[i]
buffers.append((a._data, a._byteoffset))
ne = max(ne, a.nelements())
ufargs = (ne, self.n_inputs, self.n_outputs,
tuple(buffers))
apply(cached.cfunc, ufargs)
def _slow_exec(self, inputs, outputs, cached):
input_convs, output_convs, maxitemsize = cached.ufargs
niter = _ufunc.getBufferSize()/maxitemsize
indexlevel, blockingparameters = _ufunc._getBlockingParameters(
outputs[0]._shape, niter)
rebuffered_inputs = []
for i in range(self.n_inputs):
rebuffered_inputs.append(input_convs[i].rebuffer(inputs[i]))
rebuffered_outputs = []
for i in range(self.n_outputs):
rebuffered_outputs.append(output_convs[i].rebuffer(outputs[i]))
operator = _Operator(cached.cfunc, rebuffered_inputs,
rebuffered_outputs)
objects = input_convs + (operator,) + output_convs
_callOverDimensions(objects, outputs[0]._shape, indexlevel,
blockingparameters)
for i in range(len(input_convs)):
input_convs[i].clean(inputs[i])
for i in range(len(output_convs)):
output_convs[i].clean(outputs[i])
# Portions of the following class are implemented in C. See _convertermodule.c.
# In particular, compute, stride, convert, and rebuffer.
class _Converter(_converter._converter):
def __init__(self, arr, type, inbuffer=None, forcestride=0, nonstriding=0):
self.stridefunction = None
self.convfunction = None
self.type = type
if arr.isbyteswapped():
if isinstance(arr._type, _nt.ComplexType):
fname = "byteswap" + arr._type.name
else:
fname = "byteswap"+`arr._itemsize`+"bytes"
self.stridefunction = _bytes.functionDict[fname]
elif not arr.isaligned():
fname = "align"+`arr._itemsize`+"bytes"
self.stridefunction = (_bytes.functionDict.get(fname) or
_bytes.functionDict["copyNbytes"])
elif not nonstriding and not arr.iscontiguous() or forcestride:
fname = "copy"+`arr._itemsize`+"bytes"
self.stridefunction = (_bytes.functionDict.get(fname) or
_bytes.functionDict["copyNbytes"])
if type != arr._type.name:
self.convfunction = \
_numtypedict[type]._conv.astype[arr._type.name]
if self.stridefunction:
# Since the output array is iterated over dimensions, even
# the input array, always contiguous, needs strides defined.
# But first see if input buffer is a real array.
if (inbuffer is not None) and inbuffer._strides:
if not inbuffer.is_c_array():
raise ValueError("input buffer is not a c_array")
self.bytestrides = [ inbuffer._strides,
list(arr._strides) ]
self.generated = 0
else:
self.bytestrides = [ _stridesFromShape(arr._shape,
arr._itemsize),
list(arr._strides) ]
self.generated = 1
else:
self.bytestrides = [ None, None ]
# If inbuffer not specified, allocate one
if (inbuffer is None) and (self.stridefunction or self.convfunction):
inbuffer = _bufferPool.getBuffer()
self.conversion_required = 1
if self.convfunction:
if self.stridefunction:
buff = _bufferPool.getBuffer()
buffers = [inbuffer, buff, buff, arr]
self.arr_position = 3
self.inb_position = 0
else:
buffers = [inbuffer, arr, None, None]
self.arr_position = 1
self.inb_position = 0
else:
if self.stridefunction:
buffers = [None, None,inbuffer, arr]
self.arr_position = 3
self.inb_position = 2
else:
if inbuffer is not None:
raise ValueError(
"Specified both input and output but no conversion implied")
buffers = [None, None, None, None]
inbuffer = arr
self.conversion_required = 0
self.buffers = buffers
self.result_buff = inbuffer
def clean(self, arr):
l = self.buffers
for i in range(len(l)):
if l[i] is arr:
l[i] = None
self.buffers = l
if self.result_buff is arr:
self.result_buff = None
if _PROTOTYPE:
# The following method has been re-implemented in C.
def rebuffer(self, arr, inbuffer=None):
if not self.conversion_required:
self.result_buff = arr
return arr
l = self.buffers
l[self.arr_position] = arr
if not self.generated and inbuffer is not None:
self.result_buff = inbuffer
l[self.inb_position] = inbuffer
self.buffers = l
if arr is not None and self.bytestrides != [None, None]:
if self.generated:
ustrides = _stridesFromShape(arr._shape, arr._itemsize)
else:
if inbuffer is not None:
ustrides = inbuffer._strides
else:
ustrides = self.bytestrides[ not self.direction ]
if self.direction:
self.bytestrides = [ ustrides, arr._strides ]
else:
self.bytestrides = [ arr._strides, ustrides ]
return self.result_buff
def convert(self, buf, indices, shape):
"""Perform type conversions."""
if self.convfunction:
niter = _gen.product(shape)
inoffset = self.buffers[buf*2+0]._getByteOffset(indices)
offset = self.buffers[buf*2+1]._getByteOffset(indices)
self.convfunction(niter, 1, 1,
((self.buffers[buf*2+0]._data, inoffset),
(self.buffers[buf*2+1]._data, offset)))
def stride(self, buf, indices, shape):
"""Perform copies and re-alignments."""
if self.stridefunction:
inoffset = self.buffers[buf*2+0]._getByteOffset(indices)
offset = self.buffers[buf*2+1]._getByteOffset(indices)
self.stridefunction(shape,
self.buffers[buf*2+0]._data, inoffset,
self.bytestrides[0][-len(shape):],
self.buffers[buf*2+1]._data, offset,
self.bytestrides[1][-len(shape):])
# The following code has been implemented in C.
# Renaming "compute" below affects compute, stride, and convert.
def compute(self, indices, shape):
""" Perform the necessary transformations on the array.
"""
if self.direction:
self.convert(0, indices, shape)
self.stride(1, indices, shape)
else:
self.stride(0, indices, shape)
self.convert(1, indices, shape)
class _InputConverter(_Converter):
"""This class handles alignment, byteswaping, copying and type conversions
Stride here is a bit of a misnomer, it refers to one of
functions (three currently) that deal with strides. These are:
byteswap, align, copy operations. The conv functions are type conversion
operations. They never deal with strides. At most one of each is
necessary; neither is required. Strides and offsets are *byte*!
"""
def __init__(self, arr, type, inbuffer=None, forcestride=0):
_Converter.__init__(self, arr, type, inbuffer, forcestride)
self.direction = 0 # input
if type != arr._type.name:
self.convfunction = arr._type._conv.astype[type]
l = self.bytestrides
l.reverse()
self.bytestrides = l
l = self.buffers
l.reverse()
self.buffers = l
if self.conversion_required:
self.arr_position = {1:2, 3:0}[self.arr_position]
self.inb_position = {0:3, 2:1}[self.inb_position]
class _OutputConverter(_Converter):
"""This class handles alignment, byteswaping, copying and type
conversions processing the output of a ufunc.
"""
def __init__(self, arr, type, inbuffer=None, forcestride=0, nonstriding=0):
_Converter.__init__(self, arr, type, inbuffer, forcestride,
nonstriding)
self.direction = 1 # output
# The following class has been re-implemented in C. See _operatormodule.c.
class _Operator:
"""performs the operation"""
def __init__(self, cfunction, inputs, outputs, striding=0):
self.cfunction = cfunction
self.inputs = inputs
self.outputs = outputs
self.striding = striding
def _buffer_offset(self, buffer, indices):
if buffer._strides:
return buffer._getByteOffset(indices)
else:
return 0
def _buffer_strides(self, buffer, shape, itemsize):
if buffer._strides:
return buffer._strides
else:
return _stridesFromShape(shape, itemsize)
def compute(self, indices, shape):
if self.striding:
assert len(self.inputs) == 1 and len(self.outputs) == 1
input, output = self.inputs[0], self.outputs[0]
inbuffer, outbuffer = input._data, output._data
inoffset = self._buffer_offset(input, indices)
outoffset = self._buffer_offset(output, indices)
instrides = self._buffer_strides(input, shape, self.striding)
outstrides = self._buffer_strides(output, shape, self.striding)
if len(instrides) < len(outstrides): # hack for reductions
outstrides = outstrides[len(outstrides)-len(instrides):]
self.cfunction(shape, inbuffer, inoffset, instrides,
outbuffer, outoffset, outstrides)
else:
niter = _gen.product(shape)
args = []
for item in self.inputs+self.outputs:
args.append((item._data, self._buffer_offset(item, indices)))
self.cfunction(
niter, len(self.inputs), len(self.outputs), tuple(args))
# Replace the Python version of "_Operator" with a C version.
if not _PROTOTYPE:
from _operator import _operator as _Operator
def _stridesFromShape(shape, bytestride):
"""Compute the strides from shape for a contiguous array, sort of"""
if shape != ():
ndim = len(shape)
strides = [0]*ndim
strides[-1] = bytestride
for i in xrange(ndim-2, -1, -1):
strides[i] = strides[i+1] * shape[i+1]
else:
strides = ()
return strides
def _makeCUFuncDict(functionDict):
"""Organize C functions by function name and signature"""
dict = {}
for keystr in functionDict.keys():
operator, optype, signature = keystr
if not dict.has_key(operator):
dict[operator] = {}
if not dict[operator].has_key(optype):
dict[operator][optype] = {}
dict[operator][optype][signature[0]] = (signature[1],functionDict[keystr])
return dict
def _isScalar(x):
return isinstance(x, _nt.scalarTypes)
def _maxPopType(xs):
"""_maxPopType determines the maximum type of a sequence of
numarray and scalars. see _nt.genericCoercions.
"""
if isinstance(xs[0], types.InstanceType) and \
not isinstance(xs[0], _nc.NumArray):
return None # If they're not all NumArrays, quit now.
maxT = None
for x in xs:
if isinstance(x, _nc.NumArray):
t = x._type
elif _nc.PyNUMERIC_TYPES.has_key(type(x)):
t = type(x)
else:
t = _nc._maxtype(x)
if maxT is None:
maxT = t
else:
maxT = _nt.genericCoercions[(maxT, t)]
return maxT
CLIP = 0 # Peg indices > N-1 at N-1, indices < 0 at 0
WRAP = 1 # Index = Index % N
RAISE = 2 # Raise an exception for out of range indices
class _ChooseUFunc(_UFunc):
"""Class for building the "choose" ufunc.
Class for ufuncs with 2 input (S, [~S...]) and 1 output (S) arguments"""
def _doit(self, computation_mode, outarr, cfunc, ufargs):
_ufunc.CheckFPErrors()
if computation_mode == "fast":
apply(cfunc, ufargs)
else:
inputs, outputs, preprocessing_outbuffers, \
postprocessing_inbuffers, maxitemsize = ufargs
niter = _ufunc.getBufferSize()/maxitemsize
outshape = outarr._shape
indexlevel, blockingparameters = \
_ufunc._getBlockingParameters(outshape, niter)
operator = _Operator(cfunc, preprocessing_outbuffers,
postprocessing_inbuffers)
objects = inputs + (operator,)+ outputs
_callOverDimensions(objects, outshape, indexlevel,
blockingparameters)
errorstatus = _ufunc.CheckFPErrors()
if errorstatus:
handleError(errorstatus, " in "+self.operator)
return outarr
def __call__(self, inarr1, inarr2, outarr=None, clipmode=RAISE):
"""The standard calling interface for UFuncs"""
computation_mode, woutarr, cfunc, ufargs = \
self._setup(inarr1, inarr2, outarr, clipmode)
result = self._doit(computation_mode, woutarr, cfunc, ufargs)
if outarr is not None:
outarr = (outarr,)
return _normalize_results((inarr1,)+tuple(inarr2), outarr,
(result,), False)
def _setup(self, in1, in2, outarr=None, clipmode=RAISE):
"""Setup for choose()"""
in1 = _nc.asarray(in1, type=MaybeLong)
in2 = list(in2)
if outarr is None:
convType = _maxPopType(in2)
else:
convType = outarr._type
for i in range(len(in2)):
if isinstance(in2[i], _gen.NDArray):
pass
else:
in2[i] = _nc.array(in2[i], type=convType)
result = _gen._nWayBroadcast( [in1] + in2 )
in1, in2 = result[0], result[1:]
if outarr is None:
outarr = in1.new(convType)
if in1._shape != outarr._shape:
raise ValueError("Supplied output array does not have appropriate shape")
N = outarr._itemsize
if N in [1,2,4,8,16]:
ucfname = "choose" + `N` + "bytes"
cfunc = _bytes.functionDict[ucfname]
else:
cfunc = _bytes.functionDict["chooseNbytes"]
# pars = clipmode, population count, itemsize
pars = _pseudoArrayFromScalars([clipmode, len(in2),N],
_nt.MaybeLong)
args = [pars, in1] + in2 + [outarr]
fastargs = reduce(and_,
[ _noConversionNeeded(x, convType) for x in in2 + [outarr]])
fastargs = fastargs and _noConversionNeeded(in1, _nt.MaybeLong)
if fastargs:
# If fast case, return the arguments for an immediate call to the
# C function
assert outarr._byteoffset == 0
fastparms = tuple([ (x._data, x._byteoffset) for x in args])
return ("fast", outarr, cfunc,
(outarr.nelements(), len(args)-1, 1, fastparms))
# slower, general case
cpars = _InputConverter(pars, _nt.MaybeLong)
selector = _InputConverter(in1, _nt.MaybeLong)
population = [_InputConverter(x, convType) for x in in2 ]
inputs = (cpars, selector,) + tuple(population)
outputs = (_OutputConverter(arr=outarr, type=convType),)
alltypes = [_numtypedict[convType], _nt.MaybeLong] + \
[a.type() for a in in2]
maxitemsize = max([t.bytes for t in alltypes])
preprocessing_outbuffers = tuple([x.result_buff for x in inputs])
postprocessing_inbuffers = tuple([x.result_buff for x in outputs])
return ("slow", outarr, cfunc,
(inputs, outputs, preprocessing_outbuffers,
postprocessing_inbuffers, maxitemsize))
_choose = _ChooseUFunc("choose", [], 0, 0, None)
def choose(selector, population, outarr=None, clipmode=RAISE):
"""
choose() returns a new array shaped like 'selector' with elements
chosen from members of sequence 'population' by the values of
selector. The shape of each member of 'population' must be
broadcastable to the shape of 'selector'. The value of each
member of 'selector' must satisfy: 0 <= value < len(population).
clipmode=RAISE if clipmode is CLIP, out of range selector values
are clipped at [0, shape[i]).
if clipmode is WRAP, out of range selector values
are wrapped around from 0 to shape[i] (<0) or
from shape[i] to 0 (>= shape[i]).
if clipmode is RAISE, selector values out of range
[0, shape[i]) result in an exception.
"""
return _choose(selector, population, outarr, clipmode)
def _scatteredPseudos( scattered ):
scatteredStrides = _pseudoArrayFromScalars(scattered._strides,
type=MaybeLong)
scatteredShape = _pseudoArrayFromScalars(scattered._shape,
type=MaybeLong)
return scatteredStrides, scatteredShape
class _TakeUFunc(_ChooseUFunc):
"""take(scattered, indexArrays, gathered=None)
scattered: input, array from which elements are "gathered"
indexArrays: input, tuple of index numarray or scalars
gathered: input/output result array
take plucks the elements of 'scattered' specified by the tuple of
index numarray, 'indexArrays', and stores the result in the output array
'gathered'.
'indexArrays' may be partial, i.e. incompletely specified,
hence the shape of 'gathered' is derived from both the shape of
'indexArrays' (primarily) and the shape of the trailing unspecified
dimensions of 'scattered'.
"""
def _setup(self, scattered, indexArrays, gathered=None, clipmode=RAISE):
"""Setup for Nary ufunc
"""
scattered = _nc.asarray(scattered)
indexArrays = list(indexArrays)
gatheredWasSpecified = (gathered is not None)
if gatheredWasSpecified:
ctype = gathered._type
else:
ctype = scattered._type
# *All* of the scattered array must be simultaneously available.
# So blocked conversions won't work.
if not _noConversionNeeded(scattered, ctype):
scattered = scattered.astype(ctype)
impliedShape, N = _gen._takeShape(scattered, indexArrays)
scatteredStrides, scatteredShape = _scatteredPseudos(scattered)
if gatheredWasSpecified:
if gathered._shape != impliedShape:
raise ValueError("Inconsistent array shapes...")
else:
gathered = scattered.__class__(shape=impliedShape,
type=scattered._type)
# scattered and gathered should now be identically typed.
# Since N can be arbitrarily large depending on the trailing
# dimensions of scattered, may have to make gathered contiguous
# to ensure a buffer large enough to store at least N.
# if (N > _ufunc.getBufferSize()
if not _noConversionNeeded(gathered, scattered._type):
raise ValueError("take() destination array must be contiguous, aligned, not byteswapped, and of the same type as the source array")
niter = (gathered.nelements()*gathered._itemsize)/N
cfunc = _bytes.functionDict[self.operator + "Nbytes"]
pars = _pseudoArrayFromScalars([clipmode, N], type=MaybeLong)
# Index arrays must be aligned, not byteswapped, and correctly typed.
# Throw in contiguous for now.
for i in range(len(indexArrays)):
x = indexArrays[i]
if not _noConversionNeeded(x, MaybeLong):
indexArrays[i] = x.astype(MaybeLong)
args = [pars, scattered, scatteredStrides, scatteredShape] + \
indexArrays + [gathered]
fastparms = tuple([ (x._data, x._byteoffset) for x in args])
return ("fast", gathered, cfunc, (niter, len(args)-1, 1, fastparms))
_take = _TakeUFunc("take", [], 0, 0, None)
def take(array, indices, axis=0, outarr=None, clipmode=RAISE):
"""take() picks elements of 'array' specified by sequence of
numerical sequences 'indices'.
parameters which must be specified by keyword:
array data to be indexed & collected (taken from).
indices An integer sequence, or tuple of integer sequences
specifying the array coordinates from which data
is to be taken. Partial indices result in entire
inner blocks being taken.
axis=0 selects the axis (or axes) along which the take
should be performed.
clipmode=RAISE if clipmode is CLIP, out of range indices are clipped
at [0, shape[i]).
if clipmode is WRAP, out of range indices are wrapped
around from 0 to shape[i] (<0) or from shape[i] to
0 (>= shape[i]).
if clipmode is RAISE, indices in the range -N..N-1 are
treated like Python sequence indices. Out of range
indices result in an exception.
"""
if axis == 0:
array = _nc.asarray(array)
return array._take((indices,), outarr=outarr, clipmode=clipmode)
elif isinstance(axis, (int, long)):
if isinstance(indices, (int,long,float)):
raise ValueError("indices must be a sequence")
work = _gen.swapaxes(array, 0, axis)
work = work._take((indices,), outarr=outarr, clipmode=clipmode)
return _gen.swapaxes(work, 0, axis)
else:
def_axes = range(array.rank)
for x in axis:
def_axes.remove(x)
axis = list(axis) + def_axes
work = _gen.transpose(array, axis)
return work._take(indices, outarr=outarr, clipmode=clipmode)
def _nonIteratedArray(a):
ni = a.view()
ni._shape=(0,)
ni._strides = (0,)
ni._contiguous = 1
return ni
class _PutUFunc(_TakeUFunc):
"""put(scattered, indexArrays, gathered) is the inverse function of "take",
and scatters the elements of array 'gathered' across array 'scattered' as
specified by tuple of index numarray, 'indexArrays'.
scattered: input/output Array where "stuff" gets put
indexArrays: input Arrays of indexArrays of scattered where
corresponding pieces of gathered go.
gathered: input Source array
gathered can be an array which is broadcastable to the array shape
implied by scattered and indexArrays -or- gathered can be a scalar.
"""
def _setup(self, scattered, indexArrays, gathered, clipmode=RAISE):
"""Setup for Nary ufunc
"""
indexArrays = list(indexArrays)
if ((len(indexArrays) < len(scattered._shape)) and
(scattered.isbyteswapped() or not scattered.iscontiguous())):
raise ValueError("Invalid destination array: partial indices require contiguous non-byteswapped destination")
impliedShape, N = _gen._takeShape(scattered, indexArrays)
scatteredStrides, scatteredShape = _scatteredPseudos(scattered)
scattered = _nonIteratedArray(scattered)
pars = _pseudoArrayFromScalars([clipmode, N], type=MaybeLong)
gathered = _nc.asarray(gathered, type=scattered._type)
gathered = _gen._broadcast(gathered, impliedShape)
if not gathered.iscontiguous():
gathered = gathered.copy()
# gathered and scattered should now be identically typed.
# gathered need not be contiguous since it is "bufferable".
# Since N can be arbitrarily large depending on the trailing
# dimensions of scattered, may have to pre-process gathered.
# if N > _ufunc.getBufferSize() and \
if not _noConversionNeeded(gathered, scattered._type):
gathered = gathered.astype( scattered._type )
niter = (gathered.nelements()*gathered._itemsize)/N
cfunc = _bytes.functionDict[self.operator + "Nbytes"]
for i in range(len(indexArrays)):
x = indexArrays[i]
if not _noConversionNeeded(x, MaybeLong):
indexArrays[i] = x.astype(MaybeLong)
args = [pars, gathered, scatteredStrides, scatteredShape] + \
indexArrays + [scattered]
fastparms = tuple([ (x._data, x._byteoffset) for x in args])
return ("fast", scattered, cfunc, (niter, len(args)-1, 1, fastparms))
_put = _PutUFunc("put", [], 0, 0, None)
def put(array, indices, values, axis=0, clipmode=RAISE):
"""put(array, indices, values, clipmode=RAISE, axis=0)
stores 'values' into 'array' at 'indices...'.
parameters which must be specified by keyword:
array data to be indexed & stuffed (put to).
indices An integer sequence, or tuple of integer sequences
specifying the array coordinates to hich data
is to be put. Partial indices result in entire
inner blocks being overwritten.
values A sequence of values to be written to the specified
indices of array.
axis=0 selects the axis along which the put should be performed.
clipmode=RAISE if clipmode is CLIP, out of range indices are clipped
at [0, shape[i]).
if clipmode is WRAP, out of range indices are wrapped
around from 0 to shape[i] (<0) or from shape[i] to
0 (>= shape[i])
if clipmode is RAISE, indices in the range -N..N-1 are
treated like Python sequence indices. Out of range
indices result in an exception.
"""
if not isinstance(array, _gen.NDArray):
raise TypeError("put only works on NDArray and its subclasses")
work = _nc.asarray(array)
if not work.is_c_array():
work = work.copy()
if axis == 0:
work._put((indices,), values, clipmode=clipmode)
elif isinstance(axis, (int, long)):
if isinstance(indices, (int,long,float)):
raise ValueError("indices must be a sequence")
work = _gen.swapaxes(work, 0, axis)
work._put((indices,), values, clipmode=clipmode)
work = _gen.swapaxes(work, 0, axis)
else:
def_axes = range(work.rank)
for x in axis:
def_axes.remove(x)
axis = list(axis) + def_axes
work = _gen.transpose(work, axis)
work._put(indices, values, clipmode=clipmode)
work = _gen.transpose(work, axis)
if work is not array:
if isinstance(array, _gen.NDArray):
array._copyFrom(work)
class _NonzeroUFunc:
"""nonzero(array)
array: input, array scanned for non-zero elements
Nonzero scans "array" for nonzero elements, and returns a tuple
of coordinate numarray corresponding to the nonzero elements.
"""
def __call__(self, inarr1):
"""The standard calling interface for UFuncs"""
nz = _nc.asarray(not_equal(inarr1, 0))
if nz.rank < 1:
nz.shape = (1,)
nonzeroCount = add.reduce(_gen.ravel(nz).astype(MaybeLong))
outarr = [ _nc.NumArray(shape=(nonzeroCount,), type=Long)
for i in range(nz.rank) ]
sstrides = _nc.array(nz._strides, type=MaybeLong)
sstrides /= nz._bytestride
niter = nz.nelements()
cfunc = _sort.functionDict[repr(("Bool", "nonzeroCoords"))]
args = [ nz, sstrides] + outarr
fastparms = tuple([ (x._data, x._byteoffset) for x in args])
cfunc(niter, 2, len(args)-2, fastparms)
return tuple(outarr)
def searchsorted(bins, values):
"""searchsort(bins, values) returns the array A[j] of
greatest indices 'i' such that each values[j] <= bins[i].
"""
bins = _nc.asarray(bins)
values = _nc.asarray(values)
if len(bins._shape) > 1 or len(values._shape) > 1:
raise ValueError("Multi-dimensional searchsort not supported.")
outarr = _nc.NumArray(shape=values.shape, type=Long)
maxtype = bins._type
if maxtype < values._type:
maxtype = values._type
if not _noConversionNeeded(bins, maxtype):
bins = bins.astype(maxtype)
if not _noConversionNeeded(values, maxtype):
values = values.astype(maxtype)
cfunc = _sort.functionDict[repr((bins._type.name, 'searchsorted'))]
nbins = _pseudoArrayFromScalars([len(bins)], MaybeLong)
args = [(x._data, x._byteoffset) for x in [nbins, bins, values, outarr]]
cfunc(values.nelements(), 3, 1, tuple(args))
if outarr.rank == 0:
return outarr[()]
else:
return outarr
nonzero = _NonzeroUFunc()
def _sort1(inarr1):
"""1D in-place Sort"""
conversionNeeded = not _noConversionNeeded(inarr1, inarr1._type)
contigarr1 = _nc.array(inarr1, copy=conversionNeeded)
cfunc = _sort.functionDict[repr((contigarr1._type.name, 'sort'))]
cfunc(contigarr1.nelements(), 0, 1,
((contigarr1._data, contigarr1._byteoffset),))
if conversionNeeded:
inarr1._copyFrom(contigarr1)
def _argsort1(inarr1, witness=None):
"""1D ArgSort"""
niter = inarr1.nelements()
inarr1 = _nc.array(inarr1) # make a discardable copy
if witness is None:
witness =_nc.arange(niter,type=Long)
elif not _noConversionNeeded(witness, Long):
raise ValueError("No support for converting the witness array")
cfunc = _sort.functionDict[repr((inarr1._type.name,'asort'))]
cfunc(niter, 0, 2,
tuple([ (x._data, x._byteoffset) for x in [inarr1, witness]]))
return witness
## # =====================================================================
def _fbroadcast(f, N, shape, args, params=()):
"""_fbroadcast(f, N, args, shape, params=()) calls 'f' for each of the
'N'-dimensional inner subnumarray of 'args'. Each subarray has
.shape == 'shape'[-N:]. There are a total of product(shape[:-N])
calls to 'f'.
"""
if len(shape) == N:
apply(f, tuple(args)+params)
else:
for i in range(shape[0]):
_fbroadcast(f, N, shape[1:], [x[i] for x in args], params)
def _sortN(a):
"""_sortN implements N-D in-place sort in terms of 1D sort"""
_fbroadcast(_sort1, 1, a._shape, (a,))
def _argsortN(a, w):
"""_argsortN implements N-D in-place argsort in terms of 1D argsort"""
_fbroadcast(_argsort1, 1, a._shape, (a,w))
#
# Direct sort types
#
def _broadcast_direct_sort(a, name):
"""Helper function to broadcast direct sorts.
Need to finish error handling.
"""
if len(a.shape) == 1 :
cfunc = _sort.functionDict[repr((a._type.name, name))]
if a.is_c_array() :
err = cfunc(a.nelements(), 0, 1, ((a._data, a._byteoffset),))
else :
b = _nc.array(a, copy=1)
err = cfunc(b.nelements(), 0, 1, ((b._data, b._byteoffset),))
a._copyFrom(b)
else :
for i in range(a.getshape()[0]) :
_broadcast_direct_sort(a[i], name)
def _direct_sort(a, name, axis=-1):
"""Sorts an array in-place along the specified axis.
It takes care of the sort axis.
"""
if axis==-1:
_broadcast_direct_sort(a, name)
else:
a.swapaxes(axis,-1)
_broadcast_direct_sort(a, name)
a.swapaxes(axis,-1)
#
# Indirect sort types.
#
def _broadcast_indirect_sort(a, w, name) :
"""Helper function to broadcast indirect sorts.
The witness array 'w' is assumed set in _indirect_sort
It must be of the same shape as a, be of c_type, and
contain longs. It is not checked. Need to finish error
handling.
"""
if len(a.shape) == 1 :
cfunc = _sort.functionDict[repr((a._type.name, name))]
if a.is_c_array() :
err = cfunc(a.nelements(), 1, 1, ((a._data, a._byteoffset), (w._data, w._byteoffset)))
else :
b = _nc.array(a, copy=1)
err = cfunc(a.nelements(), 1, 1, ((b._data, b._byteoffset), (w._data, w._byteoffset)))
a._copyFrom(b)
else :
for i in range(a.getshape()[0]) :
_broadcast_indirect_sort(a[i], w[i], name)
def _indirect_sort(a, name, axis=-1) :
"""Helper function to setup indirect sorts.
It takes care of the sort axis and sets
up the witness array.
"""
if axis == -1 :
ashape = a.getshape()
w = _nc.array(shape=ashape, type=_nt.Long)
w[...,:] = _nc.arange(ashape[-1], type=_nt.Long)
_broadcast_indirect_sort(a, w, 'a' + name)
return w
else :
a.swapaxes(axis,-1)
w = _indirect_sort(a, name, axis=-1)
a.swapaxes(axis,-1)
w.swapaxes(axis,-1)
return w
def divide_remainder(a,b):
"""returns (a/b, a%b)."""
a, b = _nc.asarray(a), _nc.asarray(b)
return (a/b,a%b)
Long = _nt.Long
MaybeLong = _nt.MaybeLong
# module-level variables
_numtypedict = _nt.typeDict
# Short term solution -- identity needs to be defined in _ufuncall module XXX
_identities = {"add":0, "subtract":0, "multiply":1, "divide":1,
"bitwise_and":1, "bitwise_or":0, "bitwise_xor":0,
"logical_and":1, "logical_or":0, "logical_xor":0,
}
def ufuncFactory(operator, cfuncs, identity=None):
"""Create UFunc instance based on signature of functions"""
ninputs, noutputs = _nIOArgs(cfuncs)
if ninputs==1 and noutputs==1:
return _UnaryUFunc(operator, cfuncs, 1, 1, identity)
elif ninputs==2 and noutputs==1:
return _BinaryUFunc(operator, cfuncs, 2, 1, identity)
else:
return _NaryUFunc(operator, cfuncs, ninputs, noutputs, identity)
def make_ufuncs(m):
"""Creates a dictionary of UFunc objects from a C module."""
_cufuncs = _makeCUFuncDict(m.functionDict)
_UFuncs = {}
for operator in _cufuncs.keys():
_UFuncs[operator] = ufuncFactory(operator, _cufuncs[operator],
identity=_identities.get(operator))
return _UFuncs
_UFuncs = make_ufuncs(_ufuncall)
globals().update(_UFuncs)
syntax highlighted by Code2HTML, v. 0.9.1