"""numarray: The big enchilada numeric module """ import sys as _sys import types, math, os.path import operator as _operator import copy as _copy import warnings as _warnings from math import pi, e import memory import generic as _gen import _bytes import _numarray import _ufunc import _sort import numerictypes as _nt _PROTOTYPE = 0 # Set to 1 to switch to Python prototype code. # Set to 0 to inherit C code from C basetype. # rename built-in function type so not to conflict with keyword _type = type MAX_LINE_WIDTH = None PRECISION = None SUPPRESS_SMALL = None PyINT_TYPES = { bool: 1, int: 1, long: 1, } PyREAL_TYPES = { bool: 1, int: 1, long: 1, float: 1, } # Python numeric types with values indicating level in type hierarchy PyNUMERIC_TYPES = { bool: 0, int: 1, long: 2, float: 3, complex: 4 } # Mapping back from level to type PyLevel2Type = {} for key, value in PyNUMERIC_TYPES.items(): PyLevel2Type[value] = key del key, value # Mapping from Python to Numeric types Py2NumType = { bool: _nt.Bool, int: _nt.Long, long: _nt.Int64, float: _nt.Float, complex: _nt.Complex } def array2list(arr): return arr.tolist() # array factory functions def _all_arrays(args): for x in args: if not isinstance(x, NumArray): return 0 return len(args) > 0 def _maxtype(args): """Find the maximum scalar numeric type in the arguments. An exception is raised if the types are not python numeric types. """ if not len(args): return None elif isinstance(args, NumArray): return args.type() elif _all_arrays(args): temp = args[0].type() for x in args[1:]: if temp < x.type(): temp = x.type() if isinstance(temp, _nt.BooleanType): return bool elif isinstance(temp, _nt.IntegralType): return int elif isinstance(temp, _nt.FloatingType): return float elif isinstance(temp, _nt.ComplexType): return complex else: return PyLevel2Type[_numarray._maxtype(args)] def _storePyValueInBuffer(buffer, Ctype, index, value): """Store a python value in a buffer, index is in element units, not bytes""" # Do not use for complex scalars! Ctype._conv.fromPyValue(value, buffer._data, index*Ctype.bytes, Ctype.bytes, 0) def _storePyValueListInBuffer(buffer, Ctype, valuelist): # Do not use for complex values! for i in xrange(len(valuelist)): _storePyValueInBuffer(buffer, Ctype, i, valuelist[i]) def _fillarray(size, start, delta, type=None): ptype = _maxtype((start, delta)) if ptype == long: ptype = _nt.Int64 elif PyINT_TYPES.has_key(ptype): ptype = _nt.Long elif PyREAL_TYPES.has_key(ptype): ptype = _nt.Float else: ptype = _nt.Complex if type: outtype = _nt.getType(type) if (isinstance(ptype, _nt.ComplexType) and not isinstance( outtype, _nt.ComplexType)): raise TypeError("outtype must be a complex type") else: outtype = ptype if outtype > ptype: # Hack for Int64/UInt64 on 32-bit platforms. ptype = outtype if isinstance(outtype, _nt.ComplexType): # Not memory efficient at the moment real = _fillarray(size, complex(start).real, complex(delta).real, type = _realtype(ptype)) image = _fillarray(size, complex(start).imag, complex(delta).imag, type = _realtype(ptype)) outarr = NumArray((size,), outtype, real=real, imag=image) else: # save parameters in a buffer parbuffer = ufunc._bufferPool.getBuffer() _storePyValueListInBuffer(parbuffer, ptype, [start, delta]) cfunction = _sort.functionDict[repr((ptype.name, 'fillarray'))] outarr = NumArray((size,), outtype) if ptype == outtype: # no conversion necessary, simple case _ufunc.CheckFPErrors() cfunction(size, 1, 1, ((outarr._data, 0), (parbuffer._data, 0))) errorstatus = _ufunc.CheckFPErrors() if errorstatus: ufunc.handleError(errorstatus, " in fillarray") else: # use buffer loop convbuffer = ufunc._bufferPool.getBuffer() convfunction = ptype._conv.astype[outtype.name] bsize = len(convbuffer._data)/ptype.bytes iters, lastbsize = divmod(size, bsize) _ufunc.CheckFPErrors() outoff = 0 for i in xrange(iters + (lastbsize>0)): if i == iters: bsize = lastbsize cfunction(bsize, 1, 1, ((convbuffer._data, 0), (parbuffer._data, 0))) convfunction(bsize, 1, 1, ((convbuffer._data, 0), (outarr._data, outoff))) outoff += bsize*outtype.bytes start += delta * bsize _storePyValueListInBuffer(parbuffer, ptype, [start, delta]) errorstatus = _ufunc.CheckFPErrors() if errorstatus: ufunc.handleError(errorstatus, " in fillarray") return outarr def _frontseqshape(seq): """Find the length of all the first elements, return as a list""" if not len(seq): return (0,) if isinstance(seq, types.StringType): return (len(seq),) try: shape = [] while 1: shape.append(len(seq)) try: seq = seq[0] if isinstance(seq, types.StringType): return shape except IndexError: return shape except TypeError: return shape except ValueError: if isinstance(seq, NumArray) and seq.rank == 0: return shape def fromlist(seq, type=None, shape=None, check_overflow=0, typecode=None): """fromlist creates a NumArray from the sequence 'seq' which must be a list or tuple of python numeric types. If type is specified, it is as the type of the resulting NumArray. If shape is specified, it becomes the shape of the result and must have the same number of elements as seq. """ type = _typeFromTypeAndTypecode(type, typecode) if not len(seq) and type is None: type = _nt.Long if _all_arrays(seq): a = _gen.concatenate(seq) if shape is not None: a.shape = shape else: a.shape = (len(seq),a._shape[0]/len(seq)) + a._shape[1:] if type is not None and a.type() != type: a = a.astype(type) return a if type is None: highest_type = _maxtype(seq) tshape = _frontseqshape(seq) if shape is not None and _gen.product(shape) != _gen.product(tshape): raise ValueError("shape incompatible with sequence") ndim = len(tshape) if ndim <= 0: raise TypeError("Argument must be a sequence") if type is None: type = Py2NumType.get(highest_type) if type is None: raise TypeError("Cannot create array of type %s" % highest_type.__name__) tshape = tuple(tshape) arr = NumArray(shape=tshape, type=type) arr._check_overflow = check_overflow arr.fromlist(seq) # _numarray._setarray(arr, seq) if shape is not None: arr.setshape(shape) arr._check_overflow = 0 return arr def _typeFromTypeAndTypecode(type, typecode): """returns a type object from a type or typecode specifier (keyword) or returns the type() of any sequence which is an NDArray. """ if type is not None and typecode is not None and type != typecode: raise ValueError("Can't define both 'type' and 'typecode' for an array.") elif type is not None: # Still might be a string or typecode return _nt.getType(type) elif typecode is not None: return _nt.getType(typecode) else: return None def getTypeObject(sequence, type, typecode): """getTypeObject computes the typeObject for 'sequence' if both 'type' and 'typecode' are unspecified. Otherwise, it returns the typeObject specified by 'type' or 'typecode'. """ rtype = _typeFromTypeAndTypecode(type, typecode) if rtype is not None: return rtype elif isinstance(sequence, _gen.NDArray): # handle array([]) return sequence.type() elif hasattr(sequence, "typecode"): # for Numeric/MA return sequence.typecode() elif (isinstance(sequence, (types.ListType, types.TupleType)) and len(sequence) == 0): return _nt.Long else: if isinstance(sequence, (types.IntType, types.LongType, types.FloatType, types.ComplexType)): sequence = [sequence] try: return Py2NumType[ _maxtype(sequence) ] except KeyError: raise TypeError("Can't determine a reasonable type from sequence") def array(sequence=None, typecode=None, copy=1, savespace=0, type=None, shape=None): """array() constructs a NumArray by calling NumArray, one of its factory functions (fromstring, fromfile, fromlist), or by making a copy of an existing array. If copy=0, array() will create a new array only if sequence specifies the contents or storage for the array type and typecode are interchangeable and define the array type using either a type object or string. copy when sequence is an array, copy determines whether a copy or the original is returned. savespace is ignored by numarray. shape defines the shape of the array object and is necessary when not implied by the sequence parameter. """ if sequence is None and shape is None: return None if isinstance(shape, types.IntType): shape = (shape,) if sequence is None and (shape is None or type is None): raise ValueError("Must define shape and type if sequence is None") if isinstance(sequence, _gen.NDArray): if type is None and typecode is None: if copy: a = sequence.copy() else: a = sequence else: type = getTypeObject(sequence, type, typecode) if not copy and type is sequence.type(): a = sequence else: a = sequence.astype(type) # make a re-typed copy if shape is not None: a.setshape(shape) return a type = getTypeObject(sequence, type, typecode) if sequence is None or _gen.SuitableBuffer(sequence): return NumArray(buffer=sequence, shape=shape, type=type) elif isinstance(sequence, types.StringType): return fromstring(sequence, shape=shape, type=type) elif isinstance(sequence, (types.ListType, types.TupleType)): return fromlist(sequence, type, shape) elif PyNUMERIC_TYPES.has_key(_type(sequence)): if shape and shape != (): raise ValueError("Only shape () is valid for a rank 0 array.") return fromlist([sequence], shape=(), type=type or Py2NumType[_type(sequence)]) elif isinstance(sequence, types.FileType): return fromfile(sequence, type=type, shape=shape) else: try: return sequence.__array__(type) except AttributeError: try: sequence[0] except: raise ValueError("Unknown input type") else: return fromlist(sequence, type, shape) def asarray(seq, type=None, typecode=None): """converts scalars, lists and tuples to numarray if possible. passes NumArrays thru making copies only to convert types. """ if isinstance(seq, _gen.NDArray) and type is None and typecode is None: return seq return array(seq, type=type, typecode=typecode, copy=0) inputarray = asarray # Obsolete synonym def fromstring(datastring, type=None, shape=None, typecode=None): """Create an array from binary data contained in a string (by copying)""" type = _typeFromTypeAndTypecode(type, typecode) if not shape: size, rem = divmod(len(datastring), type.bytes) if rem: raise ValueError("Type size inconsistent with string length") else: shape = (size,) # default to a 1-d array elif _type(shape) is types.IntType: shape = (shape,) if len(datastring) != (_gen.product(shape)*type.bytes): raise ValueError("Specified shape and type inconsistent with string length") arr = NumArray(shape=shape, type=type) strbuff = buffer(datastring) nelements = arr.nelements() # Currently uses only the byte-by-byte copy, should be optimized to use # larger element copies when possible. cfunc = _bytes.functionDict["copyNbytes"] cfunc((nelements,), strbuff, 0, (type.bytes,), arr._data, 0, (type.bytes,), type.bytes) if arr._type == _nt.Bool: arr = ufunc.not_equal(arr, 0) return arr def fromfile(file, type, shape=None): """Create an array from binary file data If file is a string then that file is opened, else it is assumed to be a file object. No options at the moment, all file positioning must be done prior to this function call with a file object """ type = _nt.getType(type) name = 0 if _type(file) == _type(""): name = 1 file = open(file, 'rb') size = os.path.getsize(file.name) - file.tell() if not shape: nelements, rem = divmod(size, type.bytes) if rem: raise ValueError( "Type size inconsistent with shape or remaining bytes in file") shape = (nelements,) elif _type(shape) is types.IntType: shape=(shape,) nbytes = _gen.product(shape)*type.bytes if nbytes > size: raise ValueError( "Not enough bytes left in file for specified shape and type") # create the array arr = NumArray(shape=shape, type=type) # Use of undocumented file method! XXX nbytesread = file.readinto(arr._data) if nbytesread != nbytes: raise IOError("Didn't read as many bytes as expected") if name: file.close() if arr._type == _nt.Bool: arr = ufunc.not_equal(arr, 0) return arr class UsesOpPriority(object): """Classes can subclass from UsesOpPriority to signal to numarray that perhaps the class' r-operator hook (e.g. __radd__) should be given preference over NumArray's l-operator hook (e.g. __add__). This would be done so that when different object types are used in an operation (e.g. NumArray + MaskedArray) the type of the result is well defined and independent of the order of the operands (e.g. MaskedArray). Before altering the "normal" behavior of an operator, this scheme (implemented in the operator hook functions of NumArray) first checks to see if the other operand subclasses UsesOpPriority. If it does, the op priorities of both operands are compared, and an appropriate hook function from the one with the highest priority is called. Thus, a subclass of NumArray which wants to ensure that its type dominates in mixed type operations should define a class level op_priority > 0. If several subclasses wind up doing this, op_priority will determine how they relate to one another as well. """ op_priority = 0.0 class NumArray(_numarray._numarray, _gen.NDArray, UsesOpPriority): """Fundamental Numeric Array type The type of each data element, e.g. Int32 byteorder The actual ordering of bytes in buffer: "big" or "little". """ if _PROTOTYPE: def __init__(self, shape=None, type=None, buffer=None, byteoffset=0, bytestride=None, byteorder=_sys.byteorder, aligned=1, real=None, imag=None): type = _nt.getType(type) itemsize = type.bytes _gen.NDArray.__init__(self, shape, itemsize, buffer, byteoffset, bytestride) self._type = type if byteorder in ["little", "big"]: self._byteorder = byteorder else: raise ValueError("byteorder must be 'little' or 'big'") if real is not None: self.real = real if imag is not None: self.imag = imag def _copyFrom(self, arr): """Copy elements from another array. This overrides the _generic NDArray version """ # Test for simple case first if isinstance(arr, NumArray): if (arr.nelements() == 0 and self.nelements() == 0): return if (arr._type == self._type and self._shape == arr._shape and arr._byteorder == self._byteorder and _gen.product(arr._strides) != 0 and arr.isaligned() and self.isaligned()): name = 'copy'+`self._itemsize`+'bytes' cfunc = ( _bytes.functionDict.get(name) or _bytes.functionDict["copyNbytes"]) cfunc(self._shape, arr._data, arr._byteoffset, arr._strides, self._data, self._byteoffset, self._strides, self._itemsize) return elif PyNUMERIC_TYPES.has_key(_type(arr)): # Convert scalar to a one element array for broadcasting arr = array([arr]) else: raise TypeError('argument is not array or number') barr = self._broadcast(arr) ufunc._copyFromAndConvert(barr, self) def view(self): """Returns a new array object which refers to the same data as the original array. The new array object can be manipulated (reshaped, restrided, new attributes, etc.) without affecting the original array. Modifications of the array data *do* affect the original array. """ v = _gen.NDArray.view(self) v._type = self._type v._byteorder = self._byteorder return v def __del__(self): if self._shadows != None: self._shadows._copyFrom(self) self._shadows = None def __getstate__(self): """returns state of NumArray for pickling.""" # assert not hasattr(self, "_shadows") # Not a good idea for pickling. state = _gen.NDArray.__getstate__(self) state["_byteorder"] = self._byteorder state["_type"] = self._type.name return state def __setstate__(self, state): """sets state of NumArray after unpickling.""" _gen.NDArray.__setstate__(self, state) self._byteorder = state["_byteorder"] self._type = _nt.getType(state["_type"]) def _put(self, indices, values, **keywds): ufunc._put(self, indices, values, **keywds) def _take(self, indices, **keywds): return ufunc._take(self, indices, **keywds) def getreal(self): if isinstance(self._type, _nt.ComplexType): t = _realtype(self._type) arr = NumArray(self._shape, t, buffer=self._data, byteoffset=self._byteoffset, bytestride=self._bytestride, byteorder=self._byteorder) arr._strides = self._strides[:] return arr elif isinstance(self._type, _nt.FloatingType): return self else: return self.astype(_nt.Float64) def setreal(self, value): if isinstance(self._type, (_nt.ComplexType, _nt.FloatingType)): self.getreal()[:] = value else: raise TypeError("Can't setreal() on a non-floating-point array") real = property(getreal, setreal, doc="real component of a non-complex numarray") def getimag(self): if isinstance(self._type, _nt.ComplexType): t = _realtype(self._type) arr = NumArray(self._shape, t, buffer=self._data, byteoffset=self._byteoffset+t.bytes, bytestride=self._bytestride, byteorder=self._byteorder) arr._strides = self._strides[:] return arr else: zeros = self.new(_nt.Float64) zeros[:] = 0.0 return zeros def setimag(self, value): if isinstance(self._type, _nt.ComplexType): self.getimag()[:] = value else: raise TypeError("Can't set imaginary component of a non-comlex type") imag = property(getimag, setimag, doc="imaginary component of complex array") real = property(getreal, setreal, doc="real component of complex array") setimaginary = setimag getimaginary = getimag imaginary = imag def conjugate(self): """Returns the conjugate a - bj of complex array a + bj.""" a = self.copy() ufunc.minus(a.getimag(), a.getimag()) return a def togglebyteorder(self): """reverses the state of the _byteorder attribute: little <-> big.""" self._byteorder = {"little":"big","big":"little"}[self._byteorder] def byteswap(self): """Byteswap data in place, leaving the _byteorder attribute untouched. """ if self._itemsize == 1: return if isinstance(self._type, _nt.ComplexType): fname = "byteswap" + self._type.name else: fname = "byteswap"+str(self._itemsize)+"bytes" cfunc = _bytes.functionDict[fname] cfunc(self._shape, self._data, self._byteoffset, self._strides, self._data, self._byteoffset, self._strides) _byteswap = byteswap # alias to keep old code working. def byteswapped(self): """returns a byteswapped copy of self, with adjusted _byteorder.""" b = self.copy() b._byteswap() return b def info(self): """info() prints out the key attributes of a numarray.""" _gen.NDArray.info(self) print "byteorder:", self._byteorder print "byteswap:", self.isbyteswapped() print "type:", repr(self._type) def astype(self, type=None): """Return a copy of the array converted to the given type""" if type is None: type = self._type type = _nt.getType(type) if type == self._type: # always return a copy even if type is same return self.copy() if type._conv: retarr = self.__class__(buffer=None, shape=self._shape, type=type) if retarr.nelements() == 0: return retarr if self.is_c_array(): _ufunc.CheckFPErrors() cfunc = self._type._conv.astype[type.name] cfunc(self.nelements(), 1, 1, ((self._data, self._byteoffset), (retarr._data, 0))) errorstatus = _ufunc.CheckFPErrors() if errorstatus: ufunc.handleError(errorstatus, " during type conversion") else: ufunc._copyFromAndConvert(self, retarr) elif type.fromtype: retarr = type.fromtype(self) else: raise TypeError("Don't know how to convert from %s to %s" % (self._type.name, type.name)) return retarr def __tonumtype__(self, type): """__tonumtype__ supports C-API NA_setFromPythonScalar permitting a rank-0 array to be used in lieu of a numerical scalar value.""" if self.rank != 0: raise ValueError("Can't use non-rank-0 array as a scalar.") if type is not self.type(): s = self.astype(type) else: s = self return s[()] def is_c_array(self): """returns 1 iff an array is c-contiguous, aligned, and not byteswapped.""" return (self.iscontiguous() and self.isaligned() and not self.isbyteswapped()) != 0 def is_f_array(self): """returns 1 iff an array is fortan-contiguous, aligned, and not byteswapped.""" return (self.is_fortran_contiguous() and self.isaligned() and not self.isbyteswapped()) def new(self, type=None): """Return a new array of given type with same shape as this array Note this only creates the array; it does not copy the data. """ if type is None: type = self._type return self.__class__(shape=self._shape, type=type) def type(self): """Return the type object for the array""" return self._type def copy(self): """Returns a native byte order copy of the array.""" c = _gen.NDArray.copy(self) c._byteorder = self._byteorder c._type = self._type if self.isbyteswapped(): c.byteswap() c.togglebyteorder() return c def _stdtype(self, t): return t in [_nt.Long, _nt.Float, _nt.Complex] def __str__(self): return array_str(self) def __repr__(self): return array_repr(self) def __int__(self): if len(self._shape) == 0: return int(self[()]) else: raise TypeError, "Only rank-0 numarray can be cast to integers." def __float__(self): if len(self._shape) == 0: return float(self[()]) else: raise TypeError, "Only rank-0 numarray can be cast to floats." def __complex__(self): if len(self._shape) == 0: return complex(self[()]) else: raise TypeError, "Only rank-0 numarray can be cast to complex." def __abs__(self): return ufunc.abs(self) def __neg__(self): return ufunc.minus(self) def __invert__(self): return ufunc.bitwise_not(self) def __pos__(self): return self def __add__(self, operand): if (isinstance(operand, UsesOpPriority) and self.op_priority < operand.op_priority): return operand.__radd__(self) else: return ufunc.add(self, operand) def __radd__(self, operand): if (isinstance(operand, UsesOpPriority) and self.op_priority < operand.op_priority): return operand.__add__(self) else: r = ufunc.add(operand, self) if isinstance(r, NumArray): r.__class__ = self.__class__ return r def __sub__(self, operand): if (isinstance(operand, UsesOpPriority) and self.op_priority < operand.op_priority): return operand.__rsub__(self) else: return ufunc.subtract(self, operand) def __rsub__(self, operand): if (isinstance(operand, UsesOpPriority) and self.op_priority < operand.op_priority): return operand.__sub__(self) else: r = ufunc.subtract(operand, self) if isinstance(r, NumArray): r.__class__ = self.__class__ return r def __mul__(self, operand): if (isinstance(operand, UsesOpPriority) and self.op_priority < operand.op_priority): return operand.__rmul__(self) else: return ufunc.multiply(self, operand) def __rmul__(self, operand): if (isinstance(operand, UsesOpPriority) and self.op_priority < operand.op_priority): return operand.__mul__(self) else: r = ufunc.multiply(operand, self) if isinstance(r, NumArray): r.__class__ = self.__class__ return r def __div__(self, operand): if (isinstance(operand, UsesOpPriority) and self.op_priority < operand.op_priority): return operand.__rdiv__(self) else: return ufunc.divide(self, operand) def __truediv__(self, operand): if (isinstance(operand, UsesOpPriority) and self.op_priority < operand.op_priority): return operand.__rtruediv__(self) else: return ufunc.true_divide(self, operand) def __floordiv__(self, operand): if (isinstance(operand, UsesOpPriority) and self.op_priority < operand.op_priority): return operand.__rfloordiv__(self) else: return ufunc.floor_divide(self, operand) def __rdiv__(self, operand): if (isinstance(operand, UsesOpPriority) and self.op_priority < operand.op_priority): return operand.__div__(self) else: r = ufunc.divide(operand, self) if isinstance(r, NumArray): r.__class__ = self.__class__ return r def __rtruediv__(self, operand): if (isinstance(operand, UsesOpPriority) and self.op_priority < operand.op_priority): return operand.__truediv__(self) else: r = ufunc.true_divide(operand, self) if isinstance(r, NumArray): r.__class__ = self.__class__ return r def __rfloordiv__(self, operand): if (isinstance(operand, UsesOpPriority) and self.op_priority < operand.op_priority): return operand.__floordiv__(self) else: r = ufunc.floor_divide(operand, self) if isinstance(r, NumArray): r.__class__ = self.__class__ return r def __mod__(self, operand): if (isinstance(operand, UsesOpPriority) and self.op_priority < operand.op_priority): return operand.__rmod__(self) else: return ufunc.remainder(self, operand) def __rmod__(self, operand): if (isinstance(operand, UsesOpPriority) and self.op_priority < operand.op_priority): return operand.__mod__(self) else: r = ufunc.remainder(operand, self) if isinstance(r, NumArray): r.__class__ = self.__class__ return r def __divmod__(self,operand): """returns the tuple (self/operand, self%operand).""" if (isinstance(operand, UsesOpPriority) and self.op_priority < operand.op_priority): return operand.__rdivmod__(self) else: return ufunc.divide_remainder(self, operand) def __rdivmod__(self,operand): """returns the tuple (operand/self, operand%self).""" if (isinstance(operand, UsesOpPriority) and self.op_priority < operand.op_priority): return operand.__divmod__(self) else: r = ufunc.divide_remainder(operand, self) if isinstance(r, NumArray): r.__class__ = self.__class__ return r def __pow__(self, operand): if (isinstance(operand, UsesOpPriority) and self.op_priority < operand.op_priority): return operand.__rpow__(self) else: return ufunc.power(self, operand) def __rpow__(self, operand): if (isinstance(operand, UsesOpPriority) and self.op_priority < operand.op_priority): return operand.__pow__(self) else: r = ufunc.power(operand, self) if isinstance(r, NumArray): r.__class__ = self.__class__ return r def __and__(self, operand): if (isinstance(operand, UsesOpPriority) and self.op_priority < operand.op_priority): return operand.__rand__(self) else: return ufunc.bitwise_and(self, operand) def __rand__(self, operand): if (isinstance(operand, UsesOpPriority) and self.op_priority < operand.op_priority): return operand.__and__(self) else: r = ufunc.bitwise_and(operand, self) if isinstance(r, NumArray): r.__class__ = self.__class__ return r def __or__(self, operand): if (isinstance(operand, UsesOpPriority) and self.op_priority < operand.op_priority): return operand.__ror__(self) else: return ufunc.bitwise_or(self, operand) def __ror__(self, operand): if (isinstance(operand, UsesOpPriority) and self.op_priority < operand.op_priority): return operand.__or__(self) else: r = ufunc.bitwise_or(operand, self) if isinstance(r, NumArray): r.__class__ = self.__class__ return r def __xor__(self, operand): if (isinstance(operand, UsesOpPriority) and self.op_priority < operand.op_priority): return operand.__rxor__(self) else: return ufunc.bitwise_xor(self, operand) def __rxor__(self, operand): if (isinstance(operand, UsesOpPriority) and self.op_priority < operand.op_priority): return operand.__xor__(self) else: r = ufunc.bitwise_xor(operand, self) if isinstance(r, NumArray): r.__class__ = self.__class__ return r def __rshift__(self, operand): if (isinstance(operand, UsesOpPriority) and self.op_priority < operand.op_priority): return operand.__rrshift__(self) else: return ufunc.rshift(self, operand) def __rrshift__(self, operand): if (isinstance(operand, UsesOpPriority) and self.op_priority < operand.op_priority): return operand.__rshift__(self) else: r = ufunc.rshift(operand, self) if isinstance(r, NumArray): r.__class__ = self.__class__ return r def __lshift__(self, operand): if (isinstance(operand, UsesOpPriority) and self.op_priority < operand.op_priority): return operand.__rlshift__(self) else: return ufunc.lshift(self, operand) def __rlshift__(self, operand): if (isinstance(operand, UsesOpPriority) and self.op_priority < operand.op_priority): return operand.__lshift__(self) else: r = ufunc.lshift(operand, self) if isinstance(r, NumArray): r.__class__ = self.__class__ return r # augmented assignment operators def __iadd__(self, operand): ufunc.add(self, operand, self) return self def __isub__(self, operand): ufunc.subtract(self, operand, self) return self def __imul__(self, operand): ufunc.multiply(self, operand, self) return self def __idiv__(self, operand): ufunc.divide(self, operand, self) return self def __ifloordiv__(self, operand): ufunc.floor_divide(self, operand, self) return self def __itruediv__(self, operand): ufunc.true_divide(self, operand, self) return self def __imod__(self, operand): ufunc.remainder(self, operand, self) return self def __ipow__(self, operand): ufunc.power(self, operand, self) return self def __iand__(self, operand): ufunc.bitwise_and(self, operand, self) return self def __ior__(self, operand): ufunc.bitwise_or(self, operand, self) return self def __ixor__(self, operand): ufunc.bitwise_xor(self, operand, self) return self def __irshift__(self, operand): ufunc.rshift(self, operand, self) return self def __ilshift__(self, operand): ufunc.lshift(self, operand, self) return self # rich comparisons (only works in Python 2.1 and later) def __lt__(self, operand): if isinstance(self._type, _nt.ComplexType): raise TypeError("Complex numarray don't support < comparison") else: return ufunc.less(self, operand) def __gt__(self, operand): if isinstance(self._type, _nt.ComplexType): raise TypeError("Complex numarray don't support > comparison") else: return ufunc.greater(self, operand) def __le__(self, operand): if isinstance(self._type, _nt.ComplexType): raise TypeError("Complex numarray don't support <= comparison") else: return ufunc.less_equal(self, operand) def __ge__(self, operand): if isinstance(self._type, _nt.ComplexType): raise TypeError("Complex numarray don't support >= comparison") else: return ufunc.greater_equal(self, operand) def __eq__(self, operand): if operand is None: return 0 else: return ufunc.equal(self, operand) def __ne__(self, operand): if operand is None: return 1 else: return ufunc.not_equal(self, operand) def sort(self, axis=-1): """sorts an array in-place along the specified axis, returning None.""" if axis==-1: ufunc._sortN(self) else: self.swapaxes(axis,-1) ufunc._sortN(self) self.swapaxes(axis,-1) def _argsort(self, axis=-1): if axis==-1: ashape = self.getshape() w = array(shape=ashape, type=_nt.Long) w[...,:] = arange(ashape[-1], type=_nt.Long) ufunc._argsortN(self,w) return w else: self.swapaxes(axis,-1) w = self._argsort() w.swapaxes(axis, -1) return w def argsort(self, axis=-1): """returns the indices of 'self' which if taken from self would return a copy of 'self' sorted along the specified 'axis'""" return self.copy()._argsort(axis) def argmax(self, axis=-1): """returns the index(es) of the greatest element(s) of 'self' along the specified 'axis'.""" import numeric return numeric.argmax(self, axis) def argmin(self, axis=-1): """returns the index(es) of the least element(s) of 'self' along the specified 'axis'.""" import numeric return numeric.argmin(self, axis) def diagonal(self, *args, **keywords): """returns the diagonal elements of the array.""" return diagonal(self, *args, **keywords) def trace(self, *args, **keywords): """returns the sum of the diagonal elements of the array.""" return trace(self, *args, **keywords) def typecode(self): """returns the Numeric typecode of the array.""" return _nt.typecode[self._type] def min(self): """Returns the minimum element in the array.""" return ufunc.minimum.reduce(ufunc.minimum.areduce(self).flat) def max(self): """Returns the maximum element in the array.""" return ufunc.maximum.reduce(ufunc.maximum.areduce(self).flat) def sum(self, type=None): """Returns the sum of all elements in the array.""" if type is None: type = _nt.MaximumType(self._type) return ufunc.add.reduce(ufunc.add.areduce(self, type=type).flat, type=type) def mean(self): """Returns the average of all elements in the array.""" return self.sum()/(self.nelements()*1.0) def stddev(self): """Returns the standard deviation of the array.""" m = self.mean() N = self.nelements() return math.sqrt(((self - m)**2).sum()/(N-1)) def spacesaver(self): """Always False. Supported for Numeric compatibility.""" return 0 class ComplexArray(NumArray): # Deprecated pass def Complex32_fromtype(arr): """Used for converting other types to Complex32. This is used to set an fromtype attribute of the ComplexType objects""" rarr = arr.astype(Float32) retarr = ComplexArray(arr._shape, _nt.Complex32) retarr.getreal()[:] = rarr retarr.getimag()[:] = 0. return retarr def Complex64_fromtype(arr): """Used for converting other types to Complex64. This is used to set an fromtype attribute of the ComplexType objects""" rarr = arr.astype(Float64) retarr = ComplexArray(arr._shape, _nt.Complex64) retarr.getreal()[:] = rarr retarr.getimag()[:] = 0. return retarr # Check whether byte order is big endian or little endian. from sys import byteorder isBigEndian = (byteorder == "big") del byteorder # Add fromtype function to Complex types _nt.Complex32.fromtype = Complex32_fromtype _nt.Complex64.fromtype = Complex64_fromtype # Return type of complex type components def _isComplexType(type): # Deprecated return type in [_nt.Complex32, _nt.Complex64] def _realtype(complextype): if complextype == _nt.Complex32: return _nt.Float32 else: return _nt.Float64 def conjugate(a): """conjugate(a) returns the complex conjugate of 'a'""" a = asarray(a) if not isinstance(a._type, _nt.ComplexType): a = a.astype(_nt.Complex64) return a.conjugate() def zeros(shape, type=None, typecode=None): """Constructs a zero filled array of the specified shape and type.""" type = _typeFromTypeAndTypecode(type, typecode) if type is None: type = _nt.Long retarr = NumArray(shape=shape, type=type) retarr._data.clear() return retarr def ones(shape, type=None, typecode=None): """Constructs an array of the specified shape and type filled with ones.""" type = _typeFromTypeAndTypecode(type, typecode) shape = _gen.getShape(shape) retarr = _fillarray(_gen.product(shape), 1, 0, type) retarr.setshape(shape) return retarr def arange(a1, a2=None, stride=1, type=None, shape=None, typecode=None): """Returns an array of numbers in sequence over the specified range.""" # Return empty range of correct type for single negative paramter. type = _typeFromTypeAndTypecode(type, typecode) if not isinstance(a1, types.ComplexType) and a1 < 0 and a2 is None: t = __builtins__["type"](a1) return zeros(shape=(0,), type=Py2NumType[t]) if a2 == None: start = 0 + (0*a1) # to make it same type as stop stop = a1 else: start = a1 +(0*a2) stop = a2 delta = (stop-start)/stride ## xxx int divide issue if _type(delta) == types.ComplexType: # xxx What make sense here? size = int(math.ceil(delta.real)) else: size = int(math.ceil((stop-start)/float(stride))) if size < 0: size = 0 r = _fillarray(size, start, stride, type) if shape is not None: r.setshape(shape) return r arrayrange = arange # alias arange as arrayrange. def identity(n, type=None, typecode=None): """Returns an array resembling and identity matrix.""" type = _typeFromTypeAndTypecode(type, typecode) a = zeros(shape=(n,n), type=type) i = arange(n) a.put((i, i), 1, axis=(0,)) return a if _PROTOTYPE: def dot(array1, array2): """dot matrix-multiplies array1 by array2. """ return ufunc.innerproduct(array1, _gen.swapaxes(array2, -1, -2)) else: from _numarray import dot matrixmultiply = dot # Deprecated in Numeric def outerproduct(array1, array2): """outerproduct(array1, array2) computes the NxM outerproduct of N vector 'array1' and M vector 'array2', where result[i,j] = array1[i]*array2[j]. """ array1=_gen.reshape(asarray(array1), (-1,1)) # ravel array1 into an Nx1 array2=_gen.reshape(asarray(array2), (1,-1)) # ravel array2 into a 1xM return matrixmultiply(array1,array2) # return NxM result def allclose (array1, array2, rtol=1.e-5, atol=1.e-8): # From Numeric 20.3 """allclose returns true if all components of array1 and array2 are equal subject to given tolerances. The relative error rtol must be positive and << 1.0. The absolute error atol comes into play for those elements of y that are very small or zero; it says how small x must be also. """ x, y = asarray(array1), asarray(array2) d = ufunc.less(ufunc.abs(x-y), atol + rtol * ufunc.abs(y)) return bool(alltrue(_gen.ravel(d))) # JC's revised diagonal for supporting dimensions > 2 correctly. def diagonal(a, offset= 0, axis1=0, axis2=1): """diagonal(a, offset=0, axis1=0, axis2=1) returns the given diagonals defined by the last two dimensions of the array. """ a = inputarray(a) ### do not swap axis1 and axis2, so the offset sign is meaningful ###if axis2 < axis1: axis1, axis2 = axis2, axis1 ###if axis2 > 1: if 1: ### new_axes = range(len(a._shape)) ###del new_axes[axis2]; del new_axes[axis1] try: ### new_axes.remove(axis1) ### new_axes.remove(axis2) ### except: ### raise ValueError, "axis1(=%d) and axis2(=%d) must be different and within range." % (axis1, axis2) ### ###new_axes [0:0] = [axis1, axis2] new_axes = new_axes + [axis1, axis2] ### insert at the end, not the beginning a = _gen.transpose(a, new_axes) s = a._shape rank = len(s) ### if rank == 2: ### n1 = s[0] n2 = s[1] n = n1 * n2 s = (n,) a = _gen.reshape(a, s) if offset < 0: return _gen.take(a, range(- n2 * offset, min(n2, n1+offset) * (n2+1) - n2 * offset, n2+1), axis=0) else: return _gen.take(a, range(offset, min(n1, n2-offset) * (n2+1) + offset, n2+1), axis=0) else : my_diagonal = [] for i in range(s[0]): my_diagonal.append(diagonal(a[i], offset, rank-3, rank-2)) ### return array(my_diagonal) # From Numeric-21.0 def trace(array, offset=0, axis1=0, axis2=1): """trace returns the sum along diagonals (defined by the last two dimenions) of the array. """ return ufunc.add.reduce(diagonal(array, offset, axis1, axis2)) def rank(array): """rank returns the number of dimensions in an array.""" return len(shape(array)) def shape(array): """shape returns the shape tuple of an array.""" try: return array.shape except AttributeError: return asarray(array).getshape() def size(array, axis=None): """size returns the number of elements in an array or along the specified axis.""" array = asarray(array) if axis is None: return array.nelements() else: s = array.shape if axis < 0: axis += len(s) return s[axis] def array_str(array, max_line_width=MAX_LINE_WIDTH, precision=PRECISION, suppress_small=SUPPRESS_SMALL): """Formats and array as a string.""" array = asarray(array) return arrayprint.array2string( array, max_line_width, precision, suppress_small, ' ', "") def array_repr(array, max_line_width=MAX_LINE_WIDTH, precision=PRECISION, suppress_small=SUPPRESS_SMALL): """Returns the repr string of an array.""" array = asarray(array) lst = arrayprint.array2string( array, max_line_width, precision, suppress_small, ', ', "array(") typeless = array._stdtype(array._type) if array.__class__ is not NumArray: cName= array.__class__.__name__ else: cName = "array" if typeless and not hasattr(array, "_explicit_type"): return cName + "(%s)" % lst else: return cName + "(%s, type=%s)" % (lst, array._type.name) def around(array, digits=0, output=None): """rounds 'array' to 'digits' of precision, storing the result in 'output', or returning the result as new array if output=None""" array = asarray(array) scale = 10.0**digits if output is None: wout = array.copy() else: wout = output if digits != 0: wout *= scale # bug in 2.2.1 and earlier causes fail as bad sequence op ufunc._round(wout, wout) if digits != 0: wout /= scale if output is None: return wout def round(*args, **keys): _warnings.warn("round() is deprecated. Switch to around().", DeprecationWarning) return ufunc._round(*args, **keys) def explicit_type(x): """explicit_type(x) returns a view of x which will always show it's type in it's repr. This is useful when the same test is run in two places, one where the type used *is* the default and hence not normally repr'ed, and one where the type used *is not* the default and so is displayed. """ y = x.view() y._explicit_type = 1 return y ArrayType = NumArray # Alias for backwards compatability with Numeric def array_equiv(array1, array2): """array_equiv returns True if 'a' and 'b' are shape consistent (mutually broadcastable) and have all elements equal and False otherwise.""" try: array1, array2 = asarray(array1), asarray(array2) except TypeError: return 0 if not isinstance(array1, NumArray) or not isinstance(array2, NumArray): return 0 try: array1, array2 = array1._dualbroadcast(array2) except ValueError: return 0 return ufunc.logical_and.reduce(_gen.ravel(array1 == array2)) def array_equal(array1, array2): """array_equal returns True if array1 and array2 have identical shapes and all elements equal and False otherwise.""" try: array1, array2 = asarray(array1), asarray(array2) except TypeError: return 0 if not isinstance(array1, NumArray) or not isinstance(array2, NumArray): return 0 if array1._shape != array2._shape: return 0 return ufunc.logical_and.reduce(_gen.ravel(array1 == array2)) class _UBuffer(NumArray): """_UBuffer is used to hold a single "block" of ufunc data during the block-wise processing of all elements in an array. Subclassing the buffer object from numnumarray simplifies (and speeds!) their usage at the C level. They are not intended to be used as public array objects, hence they are private! """ def __init__(self, pybuffer): NumArray.__init__(self, (len(pybuffer),), _nt.Int8, pybuffer) self._strides = None # how it is distinguished from a real array def isbyteswapped(self): return 0 def isaligned(self): return 1 def iscontiguous(self): return 1 def is_c_array(self): return 1 def _getByteOffset(self, shape): return 0 def __del__(self): """On deletion return the data to the buffer pool""" if self._data is not None: if ufunc is not None and ufunc._bufferPool is not None: ufunc._bufferPool.buffers.append(self._data) def __repr__(self): return "<_UBuffer of size:%d>" % self.shape[0] import ufunc import arrayprint sum = ufunc.add.reduce cumsum = ufunc.add.accumulate product = ufunc.multiply.reduce cumproduct = ufunc.multiply.accumulate absolute = ufunc.abs negative = ufunc.minus fmod = ufunc.remainder def alltrue(array, axis=0): """For 1D arrays, returns True IFF all the elements of the array are nonzero. For higher dimensional arrays, returns a reduction. """ array = asarray(array) if array.rank == 0: return array[()] == 1 else: return ufunc.logical_and.reduce(array, axis=axis) def sometrue(array, axis=0): """For 1D arrays, returns True IFF any one of the elements of the array is nonzero. For higher dimensional arrays, returns a reduction. """ array = asarray(array) if array.rank == 0: return array[()] == 1 else: return ufunc.logical_or.reduce(array, axis=axis) NewArray = NumArray # unnecessary following merger of real/complex Arrays def tensormultiply(array1, array2): """tensormultiply returns the product for any rank >=1 arrays, defined as: r_{xxx, yyy} = \sum_k array1_{xxx, k} array2_{k, yyyy} where xxx, yyy denote the rest of the a and b dimensions. """ array1, array2 = asarray(array1), asarray(array2) if array1.shape[-1] != array2.shape[0]: raise ValueError, "Unmatched dimensions" shape = array1.shape[:-1] + array2.shape[1:] return _gen.reshape(dot(_gen.reshape(array1, (-1, array1.shape[-1])), _gen.reshape(array2, (array2.shape[0], -1))), shape) def kroneckerproduct(a,b): '''Computes a otimes b where otimes is the Kronecker product operator. Note: the Kronecker product is also known as the matrix direct product or tensor product. It is defined as follows for 2D arrays a and b where shape(a)=(m,n) and shape(b)=(p,q): c = a otimes b => cij = a[i,j]*b where cij is the ij-th submatrix of c. So shape(c)=(m*p,n*q). >>> print kroneckerproduct([[1,2]],[[3],[4]]) [[3 6] [4 8]] >>> print kroneckerproduct([[1,2]],[[3,4]]) [ [3 4 6 8]] >>> print kroneckerproduct([[1],[2]],[[3],[4]]) [[3] [4] [6] [8]] ''' a, b = asarray(a), asarray(b) if not (len(shape(a))==2 and len(shape(b))==2): raise ValueError, 'Input must be 2D arrays.' if not a.iscontiguous(): a = _gen.reshape(a, a.shape) if not b.iscontiguous(): b = _gen.reshape(b, b.shape) o = outerproduct(a,b) o.shape = a.shape + b.shape return _gen.concatenate(_gen.concatenate(o, axis=1), axis=1)