diff options
Diffstat (limited to 'numpy/oldnumeric')
24 files changed, 4465 insertions, 0 deletions
diff --git a/numpy/oldnumeric/__init__.py b/numpy/oldnumeric/__init__.py new file mode 100644 index 000000000..05712c02c --- /dev/null +++ b/numpy/oldnumeric/__init__.py @@ -0,0 +1,45 @@ +# Don't add these to the __all__ variable though +from numpy import * + +def _move_axis_to_0(a, axis): + if axis == 0: + return a + n = len(a.shape) + if axis < 0: + axis += n + axes = range(1, axis+1) + [0,] + range(axis+1, n) + return transpose(a, axes) + +# Add these +from compat import * +from functions import * +from precision import * +from ufuncs import * +from misc import * + +import compat +import precision +import functions +import misc +import ufuncs + +import numpy +__version__ = numpy.__version__ +del numpy + +__all__ = ['__version__'] +__all__ += compat.__all__ +__all__ += precision.__all__ +__all__ += functions.__all__ +__all__ += ufuncs.__all__ +__all__ += misc.__all__ + +del compat +del functions +del precision +del ufuncs +del misc + +from numpy.testing import Tester +test = Tester().test +bench = Tester().bench diff --git a/numpy/oldnumeric/alter_code1.py b/numpy/oldnumeric/alter_code1.py new file mode 100644 index 000000000..87538a855 --- /dev/null +++ b/numpy/oldnumeric/alter_code1.py @@ -0,0 +1,240 @@ +""" +This module converts code written for Numeric to run with numpy + +Makes the following changes: + * Changes import statements (warns of use of from Numeric import *) + * Changes import statements (using numerix) ... + * Makes search and replace changes to: + - .typecode() + - .iscontiguous() + - .byteswapped() + - .itemsize() + - .toscalar() + * Converts .flat to .ravel() except for .flat = xxx or .flat[xxx] + * Replace xxx.spacesaver() with True + * Convert xx.savespace(?) to pass + ## xx.savespace(?) + + * Converts uses of 'b' to 'B' in the typecode-position of + functions: + eye, tri (in position 4) + ones, zeros, identity, empty, array, asarray, arange, + fromstring, indices, array_constructor (in position 2) + + and methods: + astype --- only argument + -- converts uses of '1', 's', 'w', and 'u' to + -- 'b', 'h', 'H', and 'I' + + * Converts uses of type(...) is <type> + isinstance(..., <type>) +""" +__all__ = ['convertfile', 'convertall', 'converttree', 'convertsrc'] + +import sys +import os +import re +import glob + + +_func4 = ['eye', 'tri'] +_meth1 = ['astype'] +_func2 = ['ones', 'zeros', 'identity', 'fromstring', 'indices', + 'empty', 'array', 'asarray', 'arange', 'array_constructor'] + +_chars = {'1':'b','s':'h','w':'H','u':'I'} + +func_re = {} +meth_re = {} + +for name in _func2: + _astr = r"""(%s\s*[(][^,]*?[,][^'"]*?['"])b(['"][^)]*?[)])"""%name + func_re[name] = re.compile(_astr, re.DOTALL) + +for name in _func4: + _astr = r"""(%s\s*[(][^,]*?[,][^,]*?[,][^,]*?[,][^'"]*?['"])b(['"][^)]*?[)])"""%name + func_re[name] = re.compile(_astr, re.DOTALL) + +for name in _meth1: + _astr = r"""(.%s\s*[(][^'"]*?['"])b(['"][^)]*?[)])"""%name + func_re[name] = re.compile(_astr, re.DOTALL) + +for char in _chars.keys(): + _astr = r"""(.astype\s*[(][^'"]*?['"])%s(['"][^)]*?[)])"""%char + meth_re[char] = re.compile(_astr, re.DOTALL) + +def fixtypechars(fstr): + for name in _func2 + _func4 + _meth1: + fstr = func_re[name].sub('\\1B\\2',fstr) + for char in _chars.keys(): + fstr = meth_re[char].sub('\\1%s\\2'%_chars[char], fstr) + return fstr + +flatindex_re = re.compile('([.]flat(\s*?[[=]))') + +def changeimports(fstr, name, newname): + importstr = 'import %s' % name + importasstr = 'import %s as ' % name + fromstr = 'from %s import ' % name + fromall=0 + + fstr = re.sub(r'(import\s+[^,\n\r]+,\s*)(%s)' % name, + "\\1%s as %s" % (newname, name), fstr) + fstr = fstr.replace(importasstr, 'import %s as ' % newname) + fstr = fstr.replace(importstr, 'import %s as %s' % (newname,name)) + + ind = 0 + Nlen = len(fromstr) + Nlen2 = len("from %s import " % newname) + while 1: + found = fstr.find(fromstr,ind) + if (found < 0): + break + ind = found + Nlen + if fstr[ind] == '*': + continue + fstr = "%sfrom %s import %s" % (fstr[:found], newname, fstr[ind:]) + ind += Nlen2 - Nlen + return fstr, fromall + +istest_re = {} +_types = ['float', 'int', 'complex', 'ArrayType', 'FloatType', + 'IntType', 'ComplexType'] +for name in _types: + _astr = r'type\s*[(]([^)]*)[)]\s+(?:is|==)\s+(.*?%s)'%name + istest_re[name] = re.compile(_astr) +def fixistesting(astr): + for name in _types: + astr = istest_re[name].sub('isinstance(\\1, \\2)', astr) + return astr + +def replaceattr(astr): + astr = astr.replace(".typecode()",".dtype.char") + astr = astr.replace(".iscontiguous()",".flags.contiguous") + astr = astr.replace(".byteswapped()",".byteswap()") + astr = astr.replace(".toscalar()", ".item()") + astr = astr.replace(".itemsize()",".itemsize") + # preserve uses of flat that should be o.k. + tmpstr = flatindex_re.sub(r"@@@@\2",astr) + # replace other uses of flat + tmpstr = tmpstr.replace(".flat",".ravel()") + # put back .flat where it was valid + astr = tmpstr.replace("@@@@", ".flat") + return astr + +svspc2 = re.compile(r'([^,(\s]+[.]spacesaver[(][)])') +svspc3 = re.compile(r'(\S+[.]savespace[(].*[)])') +#shpe = re.compile(r'(\S+\s*)[.]shape\s*=[^=]\s*(.+)') +def replaceother(astr): + astr = svspc2.sub('True',astr) + astr = svspc3.sub(r'pass ## \1', astr) + #astr = shpe.sub('\\1=\\1.reshape(\\2)', astr) + return astr + +import datetime +def fromstr(filestr): + savestr = filestr[:] + filestr = fixtypechars(filestr) + filestr = fixistesting(filestr) + filestr, fromall1 = changeimports(filestr, 'Numeric', 'numpy.oldnumeric') + filestr, fromall1 = changeimports(filestr, 'multiarray','numpy.oldnumeric') + filestr, fromall1 = changeimports(filestr, 'umath', 'numpy.oldnumeric') + filestr, fromall1 = changeimports(filestr, 'Precision', 'numpy.oldnumeric.precision') + filestr, fromall1 = changeimports(filestr, 'UserArray', 'numpy.oldnumeric.user_array') + filestr, fromall1 = changeimports(filestr, 'ArrayPrinter', 'numpy.oldnumeric.array_printer') + filestr, fromall2 = changeimports(filestr, 'numerix', 'numpy.oldnumeric') + filestr, fromall3 = changeimports(filestr, 'scipy_base', 'numpy.oldnumeric') + filestr, fromall3 = changeimports(filestr, 'Matrix', 'numpy.oldnumeric.matrix') + filestr, fromall3 = changeimports(filestr, 'MLab', 'numpy.oldnumeric.mlab') + filestr, fromall3 = changeimports(filestr, 'LinearAlgebra', 'numpy.oldnumeric.linear_algebra') + filestr, fromall3 = changeimports(filestr, 'RNG', 'numpy.oldnumeric.rng') + filestr, fromall3 = changeimports(filestr, 'RNG.Statistics', 'numpy.oldnumeric.rng_stats') + filestr, fromall3 = changeimports(filestr, 'RandomArray', 'numpy.oldnumeric.random_array') + filestr, fromall3 = changeimports(filestr, 'FFT', 'numpy.oldnumeric.fft') + filestr, fromall3 = changeimports(filestr, 'MA', 'numpy.oldnumeric.ma') + fromall = fromall1 or fromall2 or fromall3 + filestr = replaceattr(filestr) + filestr = replaceother(filestr) + if savestr != filestr: + today = datetime.date.today().strftime('%b %d, %Y') + name = os.path.split(sys.argv[0])[-1] + filestr = '## Automatically adapted for '\ + 'numpy.oldnumeric %s by %s\n\n%s' % (today, name, filestr) + return filestr, 1 + return filestr, 0 + +def makenewfile(name, filestr): + fid = file(name, 'w') + fid.write(filestr) + fid.close() + +def convertfile(filename, orig=1): + """Convert the filename given from using Numeric to using NumPy + + Copies the file to filename.orig and then over-writes the file + with the updated code + """ + fid = open(filename) + filestr = fid.read() + fid.close() + filestr, changed = fromstr(filestr) + if changed: + if orig: + base, ext = os.path.splitext(filename) + os.rename(filename, base+".orig") + else: + os.remove(filename) + makenewfile(filename, filestr) + +def fromargs(args): + filename = args[1] + converttree(filename) + +def convertall(direc=os.path.curdir, orig=1): + """Convert all .py files to use numpy.oldnumeric (from Numeric) in the directory given + + For each changed file, a backup of <usesnumeric>.py is made as + <usesnumeric>.py.orig. A new file named <usesnumeric>.py + is then written with the updated code. + """ + files = glob.glob(os.path.join(direc,'*.py')) + for afile in files: + if afile[-8:] == 'setup.py': continue # skip these + convertfile(afile, orig) + +header_re = re.compile(r'(Numeric/arrayobject.h)') + +def convertsrc(direc=os.path.curdir, ext=None, orig=1): + """Replace Numeric/arrayobject.h with numpy/oldnumeric.h in all files in the + directory with extension give by list ext (if ext is None, then all files are + replaced).""" + if ext is None: + files = glob.glob(os.path.join(direc,'*')) + else: + files = [] + for aext in ext: + files.extend(glob.glob(os.path.join(direc,"*.%s" % aext))) + for afile in files: + fid = open(afile) + fstr = fid.read() + fid.close() + fstr, n = header_re.subn(r'numpy/oldnumeric.h',fstr) + if n > 0: + if orig: + base, ext = os.path.splitext(afile) + os.rename(afile, base+".orig") + else: + os.remove(afile) + makenewfile(afile, fstr) + +def _func(arg, dirname, fnames): + convertall(dirname, orig=0) + convertsrc(dirname, ext=['h','c'], orig=0) + +def converttree(direc=os.path.curdir): + """Convert all .py files and source code files in the tree given + """ + os.path.walk(direc, _func, None) + + +if __name__ == '__main__': + fromargs(sys.argv) diff --git a/numpy/oldnumeric/alter_code2.py b/numpy/oldnumeric/alter_code2.py new file mode 100644 index 000000000..baa6b9d26 --- /dev/null +++ b/numpy/oldnumeric/alter_code2.py @@ -0,0 +1,146 @@ +""" +This module converts code written for numpy.oldnumeric to work +with numpy + +FIXME: Flesh this out. + +Makes the following changes: + * Converts typecharacters '1swu' to 'bhHI' respectively + when used as typecodes + * Changes import statements + * Change typecode= to dtype= + * Eliminates savespace=xxx keyword arguments + * Removes it when keyword is not given as well + * replaces matrixmultiply with dot + * converts functions that don't give axis= keyword that have changed + * converts functions that don't give typecode= keyword that have changed + * converts use of capitalized type-names + * converts old function names in oldnumeric.linear_algebra, + oldnumeric.random_array, and oldnumeric.fft + +""" +#__all__ = ['convertfile', 'convertall', 'converttree'] +__all__ = [] + +import warnings +warnings.warn("numpy.oldnumeric.alter_code2 is not working yet.") + +import sys +import os +import re +import glob + +# To convert typecharacters we need to +# Not very safe. Disabled for now.. +def replacetypechars(astr): + astr = astr.replace("'s'","'h'") + astr = astr.replace("'b'","'B'") + astr = astr.replace("'1'","'b'") + astr = astr.replace("'w'","'H'") + astr = astr.replace("'u'","'I'") + return astr + +def changeimports(fstr, name, newname): + importstr = 'import %s' % name + importasstr = 'import %s as ' % name + fromstr = 'from %s import ' % name + fromall=0 + + fstr = fstr.replace(importasstr, 'import %s as ' % newname) + fstr = fstr.replace(importstr, 'import %s as %s' % (newname,name)) + + ind = 0 + Nlen = len(fromstr) + Nlen2 = len("from %s import " % newname) + while 1: + found = fstr.find(fromstr,ind) + if (found < 0): + break + ind = found + Nlen + if fstr[ind] == '*': + continue + fstr = "%sfrom %s import %s" % (fstr[:found], newname, fstr[ind:]) + ind += Nlen2 - Nlen + return fstr, fromall + +def replaceattr(astr): + astr = astr.replace("matrixmultiply","dot") + return astr + +def replaceother(astr): + astr = re.sub(r'typecode\s*=', 'dtype=', astr) + astr = astr.replace('ArrayType', 'ndarray') + astr = astr.replace('NewAxis', 'newaxis') + return astr + +import datetime +def fromstr(filestr): + #filestr = replacetypechars(filestr) + filestr, fromall1 = changeimports(filestr, 'numpy.oldnumeric', 'numpy') + filestr, fromall1 = changeimports(filestr, 'numpy.core.multiarray', 'numpy') + filestr, fromall1 = changeimports(filestr, 'numpy.core.umath', 'numpy') + filestr, fromall3 = changeimports(filestr, 'LinearAlgebra', + 'numpy.linalg.old') + filestr, fromall3 = changeimports(filestr, 'RNG', 'numpy.random.oldrng') + filestr, fromall3 = changeimports(filestr, 'RNG.Statistics', 'numpy.random.oldrngstats') + filestr, fromall3 = changeimports(filestr, 'RandomArray', 'numpy.random.oldrandomarray') + filestr, fromall3 = changeimports(filestr, 'FFT', 'numpy.fft.old') + filestr, fromall3 = changeimports(filestr, 'MA', 'numpy.core.ma') + fromall = fromall1 or fromall2 or fromall3 + filestr = replaceattr(filestr) + filestr = replaceother(filestr) + today = datetime.date.today().strftime('%b %d, %Y') + name = os.path.split(sys.argv[0])[-1] + filestr = '## Automatically adapted for '\ + 'numpy %s by %s\n\n%s' % (today, name, filestr) + return filestr + +def makenewfile(name, filestr): + fid = file(name, 'w') + fid.write(filestr) + fid.close() + +def getandcopy(name): + fid = file(name) + filestr = fid.read() + fid.close() + base, ext = os.path.splitext(name) + makenewfile(base+'.orig', filestr) + return filestr + +def convertfile(filename): + """Convert the filename given from using Numeric to using NumPy + + Copies the file to filename.orig and then over-writes the file + with the updated code + """ + filestr = getandcopy(filename) + filestr = fromstr(filestr) + makenewfile(filename, filestr) + +def fromargs(args): + filename = args[1] + convertfile(filename) + +def convertall(direc=os.path.curdir): + """Convert all .py files to use NumPy (from Numeric) in the directory given + + For each file, a backup of <usesnumeric>.py is made as + <usesnumeric>.py.orig. A new file named <usesnumeric>.py + is then written with the updated code. + """ + files = glob.glob(os.path.join(direc,'*.py')) + for afile in files: + convertfile(afile) + +def _func(arg, dirname, fnames): + convertall(dirname) + +def converttree(direc=os.path.curdir): + """Convert all .py files in the tree given + + """ + os.path.walk(direc, _func, None) + +if __name__ == '__main__': + fromargs(sys.argv) diff --git a/numpy/oldnumeric/array_printer.py b/numpy/oldnumeric/array_printer.py new file mode 100644 index 000000000..95f3f42c7 --- /dev/null +++ b/numpy/oldnumeric/array_printer.py @@ -0,0 +1,16 @@ + +__all__ = ['array2string'] + +from numpy import array2string as _array2string + +def array2string(a, max_line_width=None, precision=None, + suppress_small=None, separator=' ', + array_output=0): + if array_output: + prefix="array(" + style=repr + else: + prefix = "" + style=str + return _array2string(a, max_line_width, precision, + suppress_small, separator, prefix, style) diff --git a/numpy/oldnumeric/arrayfns.py b/numpy/oldnumeric/arrayfns.py new file mode 100644 index 000000000..dbb910770 --- /dev/null +++ b/numpy/oldnumeric/arrayfns.py @@ -0,0 +1,97 @@ +"""Backward compatible with arrayfns from Numeric +""" + +__all__ = ['array_set', 'construct3', 'digitize', 'error', 'find_mask', + 'histogram', 'index_sort', 'interp', 'nz', 'reverse', 'span', + 'to_corners', 'zmin_zmax'] + +import numpy as np +from numpy import asarray + +class error(Exception): + pass + +def array_set(vals1, indices, vals2): + indices = asarray(indices) + if indices.ndim != 1: + raise ValueError, "index array must be 1-d" + if not isinstance(vals1, np.ndarray): + raise TypeError, "vals1 must be an ndarray" + vals1 = asarray(vals1) + vals2 = asarray(vals2) + if vals1.ndim != vals2.ndim or vals1.ndim < 1: + raise error, "vals1 and vals2 must have same number of dimensions (>=1)" + vals1[indices] = vals2 + +from numpy import digitize +from numpy import bincount as histogram + +def index_sort(arr): + return asarray(arr).argsort(kind='heap') + +def interp(y, x, z, typ=None): + """y(z) interpolated by treating y(x) as piecewise function + """ + res = np.interp(z, x, y) + if typ is None or typ == 'd': + return res + if typ == 'f': + return res.astype('f') + + raise error, "incompatible typecode" + +def nz(x): + x = asarray(x,dtype=np.ubyte) + if x.ndim != 1: + raise TypeError, "intput must have 1 dimension." + indxs = np.flatnonzero(x != 0) + return indxs[-1].item()+1 + +def reverse(x, n): + x = asarray(x,dtype='d') + if x.ndim != 2: + raise ValueError, "input must be 2-d" + y = np.empty_like(x) + if n == 0: + y[...] = x[::-1,:] + elif n == 1: + y[...] = x[:,::-1] + return y + +def span(lo, hi, num, d2=0): + x = np.linspace(lo, hi, num) + if d2 <= 0: + return x + else: + ret = np.empty((d2,num),x.dtype) + ret[...] = x + return ret + +def zmin_zmax(z, ireg): + z = asarray(z, dtype=float) + ireg = asarray(ireg, dtype=int) + if z.shape != ireg.shape or z.ndim != 2: + raise ValueError, "z and ireg must be the same shape and 2-d" + ix, iy = np.nonzero(ireg) + # Now, add more indices + x1m = ix - 1 + y1m = iy-1 + i1 = x1m>=0 + i2 = y1m>=0 + i3 = i1 & i2 + nix = np.r_[ix, x1m[i1], x1m[i1], ix[i2] ] + niy = np.r_[iy, iy[i1], y1m[i3], y1m[i2]] + # remove any negative indices + zres = z[nix,niy] + return zres.min().item(), zres.max().item() + + +def find_mask(fs, node_edges): + raise NotImplementedError + +def to_corners(arr, nv, nvsum): + raise NotImplementedError + + +def construct3(mask, itype): + raise NotImplementedError diff --git a/numpy/oldnumeric/compat.py b/numpy/oldnumeric/compat.py new file mode 100644 index 000000000..3e1d53d0e --- /dev/null +++ b/numpy/oldnumeric/compat.py @@ -0,0 +1,109 @@ +# Compatibility module containing deprecated names + +__all__ = ['NewAxis', + 'UFuncType', 'UfuncType', 'ArrayType', 'arraytype', + 'LittleEndian', 'arrayrange', 'matrixmultiply', + 'array_constructor', 'pickle_array', + 'DumpArray', 'LoadArray', 'multiarray', + # from cPickle + 'dump', 'dumps', 'load', 'loads', + 'Unpickler', 'Pickler' + ] + +import numpy.core.multiarray as multiarray +import numpy.core.umath as um +from numpy.core.numeric import array +import functions +import sys + +from cPickle import dump, dumps + +mu = multiarray + +#Use this to add a new axis to an array +#compatibility only +NewAxis = None + +#deprecated +UFuncType = type(um.sin) +UfuncType = type(um.sin) +ArrayType = mu.ndarray +arraytype = mu.ndarray + +LittleEndian = (sys.byteorder == 'little') + +from numpy import deprecate + +# backward compatibility +arrayrange = deprecate(functions.arange, 'arrayrange', 'arange') + +# deprecated names +matrixmultiply = deprecate(mu.dot, 'matrixmultiply', 'dot') + +def DumpArray(m, fp): + m.dump(fp) + +def LoadArray(fp): + import cPickle + return cPickle.load(fp) + +def array_constructor(shape, typecode, thestr, Endian=LittleEndian): + if typecode == "O": + x = array(thestr, "O") + else: + x = mu.fromstring(thestr, typecode) + x.shape = shape + if LittleEndian != Endian: + return x.byteswap(True) + else: + return x + +def pickle_array(a): + if a.dtype.hasobject: + return (array_constructor, + a.shape, a.dtype.char, a.tolist(), LittleEndian) + else: + return (array_constructor, + (a.shape, a.dtype.char, a.tostring(), LittleEndian)) + +def loads(astr): + import cPickle + arr = cPickle.loads(astr.replace('Numeric', 'numpy.oldnumeric')) + return arr + +def load(fp): + return loads(fp.read()) + +def _LoadArray(fp): + import typeconv + ln = fp.readline().split() + if ln[0][0] == 'A': ln[0] = ln[0][1:] + typecode = ln[0][0] + endian = ln[0][1] + itemsize = int(ln[0][2:]) + shape = [int(x) for x in ln[1:]] + sz = itemsize + for val in shape: + sz *= val + dstr = fp.read(sz) + m = mu.fromstring(dstr, typeconv.convtypecode(typecode)) + m.shape = shape + + if (LittleEndian and endian == 'B') or (not LittleEndian and endian == 'L'): + return m.byteswap(True) + else: + return m + +import pickle, copy +class Unpickler(pickle.Unpickler): + def load_array(self): + self.stack.append(_LoadArray(self)) + + dispatch = copy.copy(pickle.Unpickler.dispatch) + dispatch['A'] = load_array + +class Pickler(pickle.Pickler): + def __init__(self, *args, **kwds): + raise NotImplementedError, "Don't pickle new arrays with this" + def save_array(self, object): + raise NotImplementedError, "Don't pickle new arrays with this" diff --git a/numpy/oldnumeric/fft.py b/numpy/oldnumeric/fft.py new file mode 100644 index 000000000..67f30c750 --- /dev/null +++ b/numpy/oldnumeric/fft.py @@ -0,0 +1,21 @@ + +__all__ = ['fft', 'fft2d', 'fftnd', 'hermite_fft', 'inverse_fft', + 'inverse_fft2d', 'inverse_fftnd', + 'inverse_hermite_fft', 'inverse_real_fft', + 'inverse_real_fft2d', 'inverse_real_fftnd', + 'real_fft', 'real_fft2d', 'real_fftnd'] + +from numpy.fft import fft +from numpy.fft import fft2 as fft2d +from numpy.fft import fftn as fftnd +from numpy.fft import hfft as hermite_fft +from numpy.fft import ifft as inverse_fft +from numpy.fft import ifft2 as inverse_fft2d +from numpy.fft import ifftn as inverse_fftnd +from numpy.fft import ihfft as inverse_hermite_fft +from numpy.fft import irfft as inverse_real_fft +from numpy.fft import irfft2 as inverse_real_fft2d +from numpy.fft import irfftn as inverse_real_fftnd +from numpy.fft import rfft as real_fft +from numpy.fft import rfft2 as real_fft2d +from numpy.fft import rfftn as real_fftnd diff --git a/numpy/oldnumeric/fix_default_axis.py b/numpy/oldnumeric/fix_default_axis.py new file mode 100644 index 000000000..8483de85e --- /dev/null +++ b/numpy/oldnumeric/fix_default_axis.py @@ -0,0 +1,291 @@ +""" +This module adds the default axis argument to code which did not specify it +for the functions where the default was changed in NumPy. + +The functions changed are + +add -1 ( all second argument) +====== +nansum +nanmax +nanmin +nanargmax +nanargmin +argmax +argmin +compress 3 + + +add 0 +====== +take 3 +repeat 3 +sum # might cause problems with builtin. +product +sometrue +alltrue +cumsum +cumproduct +average +ptp +cumprod +prod +std +mean +""" +__all__ = ['convertfile', 'convertall', 'converttree'] + +import sys +import os +import re +import glob + + +_args3 = ['compress', 'take', 'repeat'] +_funcm1 = ['nansum', 'nanmax', 'nanmin', 'nanargmax', 'nanargmin', + 'argmax', 'argmin', 'compress'] +_func0 = ['take', 'repeat', 'sum', 'product', 'sometrue', 'alltrue', + 'cumsum', 'cumproduct', 'average', 'ptp', 'cumprod', 'prod', + 'std', 'mean'] + +_all = _func0 + _funcm1 +func_re = {} + +for name in _all: + _astr = r"""%s\s*[(]"""%name + func_re[name] = re.compile(_astr) + + +import string +disallowed = '_' + string.uppercase + string.lowercase + string.digits + +def _add_axis(fstr, name, repl): + alter = 0 + if name in _args3: + allowed_comma = 1 + else: + allowed_comma = 0 + newcode = "" + last = 0 + for obj in func_re[name].finditer(fstr): + nochange = 0 + start, end = obj.span() + if fstr[start-1] in disallowed: + continue + if fstr[start-1] == '.' \ + and fstr[start-6:start-1] != 'numpy' \ + and fstr[start-2:start-1] != 'N' \ + and fstr[start-9:start-1] != 'numarray' \ + and fstr[start-8:start-1] != 'numerix' \ + and fstr[start-8:start-1] != 'Numeric': + continue + if fstr[start-1] in ['\t',' ']: + k = start-2 + while fstr[k] in ['\t',' ']: + k -= 1 + if fstr[k-2:k+1] == 'def' or \ + fstr[k-4:k+1] == 'class': + continue + k = end + stack = 1 + ncommas = 0 + N = len(fstr) + while stack: + if k>=N: + nochange =1 + break + if fstr[k] == ')': + stack -= 1 + elif fstr[k] == '(': + stack += 1 + elif stack == 1 and fstr[k] == ',': + ncommas += 1 + if ncommas > allowed_comma: + nochange = 1 + break + k += 1 + if nochange: + continue + alter += 1 + newcode = "%s%s,%s)" % (newcode, fstr[last:k-1], repl) + last = k + if not alter: + newcode = fstr + else: + newcode = "%s%s" % (newcode, fstr[last:]) + return newcode, alter + +def _import_change(fstr, names): + # Four possibilities + # 1.) import numpy with subsequent use of numpy.<name> + # change this to import numpy.oldnumeric as numpy + # 2.) import numpy as XXXX with subsequent use of + # XXXX.<name> ==> import numpy.oldnumeric as XXXX + # 3.) from numpy import * + # with subsequent use of one of the names + # 4.) from numpy import ..., <name>, ... (could span multiple + # lines. ==> remove all names from list and + # add from numpy.oldnumeric import <name> + + num = 0 + # case 1 + importstr = "import numpy" + ind = fstr.find(importstr) + if (ind > 0): + found = 0 + for name in names: + ind2 = fstr.find("numpy.%s" % name, ind) + if (ind2 > 0): + found = 1 + break + if found: + fstr = "%s%s%s" % (fstr[:ind], "import numpy.oldnumeric as numpy", + fstr[ind+len(importstr):]) + num += 1 + + # case 2 + importre = re.compile("""import numpy as ([A-Za-z0-9_]+)""") + modules = importre.findall(fstr) + if len(modules) > 0: + for module in modules: + found = 0 + for name in names: + ind2 = fstr.find("%s.%s" % (module, name)) + if (ind2 > 0): + found = 1 + break + if found: + importstr = "import numpy as %s" % module + ind = fstr.find(importstr) + fstr = "%s%s%s" % (fstr[:ind], + "import numpy.oldnumeric as %s" % module, + fstr[ind+len(importstr):]) + num += 1 + + # case 3 + importstr = "from numpy import *" + ind = fstr.find(importstr) + if (ind > 0): + found = 0 + for name in names: + ind2 = fstr.find(name, ind) + if (ind2 > 0) and fstr[ind2-1] not in disallowed: + found = 1 + break + if found: + fstr = "%s%s%s" % (fstr[:ind], + "from numpy.oldnumeric import *", + fstr[ind+len(importstr):]) + num += 1 + + # case 4 + ind = 0 + importstr = "from numpy import" + N = len(importstr) + while 1: + ind = fstr.find(importstr, ind) + if (ind < 0): + break + ind += N + ptr = ind+1 + stack = 1 + while stack: + if fstr[ptr] == '\\': + stack += 1 + elif fstr[ptr] == '\n': + stack -= 1 + ptr += 1 + substr = fstr[ind:ptr] + found = 0 + substr = substr.replace('\n',' ') + substr = substr.replace('\\','') + importnames = [x.strip() for x in substr.split(',')] + # determine if any of names are in importnames + addnames = [] + for name in names: + if name in importnames: + importnames.remove(name) + addnames.append(name) + if len(addnames) > 0: + fstr = "%s%s\n%s\n%s" % \ + (fstr[:ind], + "from numpy import %s" % \ + ", ".join(importnames), + "from numpy.oldnumeric import %s" % \ + ", ".join(addnames), + fstr[ptr:]) + num += 1 + + return fstr, num + +def add_axis(fstr, import_change=False): + total = 0 + if not import_change: + for name in _funcm1: + fstr, num = _add_axis(fstr, name, 'axis=-1') + total += num + for name in _func0: + fstr, num = _add_axis(fstr, name, 'axis=0') + total += num + return fstr, total + else: + fstr, num = _import_change(fstr, _funcm1+_func0) + return fstr, num + + +def makenewfile(name, filestr): + fid = file(name, 'w') + fid.write(filestr) + fid.close() + +def getfile(name): + fid = file(name) + filestr = fid.read() + fid.close() + return filestr + +def copyfile(name, fstr): + base, ext = os.path.splitext(name) + makenewfile(base+'.orig', fstr) + return + +def convertfile(filename, import_change=False): + """Convert the filename given from using Numeric to using NumPy + + Copies the file to filename.orig and then over-writes the file + with the updated code + """ + filestr = getfile(filename) + newstr, total = add_axis(filestr, import_change) + if total > 0: + print "Changing ", filename + copyfile(filename, filestr) + makenewfile(filename, newstr) + sys.stdout.flush() + +def fromargs(args): + filename = args[1] + convertfile(filename) + +def convertall(direc=os.path.curdir, import_change=False): + """Convert all .py files in the directory given + + For each file, a backup of <usesnumeric>.py is made as + <usesnumeric>.py.orig. A new file named <usesnumeric>.py + is then written with the updated code. + """ + files = glob.glob(os.path.join(direc,'*.py')) + for afile in files: + convertfile(afile, import_change) + +def _func(arg, dirname, fnames): + convertall(dirname, import_change=arg) + +def converttree(direc=os.path.curdir, import_change=False): + """Convert all .py files in the tree given + + """ + os.path.walk(direc, _func, import_change) + +if __name__ == '__main__': + fromargs(sys.argv) diff --git a/numpy/oldnumeric/functions.py b/numpy/oldnumeric/functions.py new file mode 100644 index 000000000..5b2b1a8bf --- /dev/null +++ b/numpy/oldnumeric/functions.py @@ -0,0 +1,124 @@ +# Functions that should behave the same as Numeric and need changing + +import numpy as np +import numpy.core.multiarray as mu +import numpy.core.numeric as nn +from typeconv import convtypecode, convtypecode2 + +__all__ = ['take', 'repeat', 'sum', 'product', 'sometrue', 'alltrue', + 'cumsum', 'cumproduct', 'compress', 'fromfunction', + 'ones', 'empty', 'identity', 'zeros', 'array', 'asarray', + 'nonzero', 'reshape', 'arange', 'fromstring', 'ravel', 'trace', + 'indices', 'where','sarray','cross_product', 'argmax', 'argmin', + 'average'] + +def take(a, indicies, axis=0): + return np.take(a, indicies, axis) + +def repeat(a, repeats, axis=0): + return np.repeat(a, repeats, axis) + +def sum(x, axis=0): + return np.sum(x, axis) + +def product(x, axis=0): + return np.product(x, axis) + +def sometrue(x, axis=0): + return np.sometrue(x, axis) + +def alltrue(x, axis=0): + return np.alltrue(x, axis) + +def cumsum(x, axis=0): + return np.cumsum(x, axis) + +def cumproduct(x, axis=0): + return np.cumproduct(x, axis) + +def argmax(x, axis=-1): + return np.argmax(x, axis) + +def argmin(x, axis=-1): + return np.argmin(x, axis) + +def compress(condition, m, axis=-1): + return np.compress(condition, m, axis) + +def fromfunction(args, dimensions): + return np.fromfunction(args, dimensions, dtype=int) + +def ones(shape, typecode='l', savespace=0, dtype=None): + """ones(shape, dtype=int) returns an array of the given + dimensions which is initialized to all ones. + """ + dtype = convtypecode(typecode,dtype) + a = mu.empty(shape, dtype) + a.fill(1) + return a + +def zeros(shape, typecode='l', savespace=0, dtype=None): + """zeros(shape, dtype=int) returns an array of the given + dimensions which is initialized to all zeros + """ + dtype = convtypecode(typecode,dtype) + return mu.zeros(shape, dtype) + +def identity(n,typecode='l', dtype=None): + """identity(n) returns the identity 2-d array of shape n x n. + """ + dtype = convtypecode(typecode, dtype) + return nn.identity(n, dtype) + +def empty(shape, typecode='l', dtype=None): + dtype = convtypecode(typecode, dtype) + return mu.empty(shape, dtype) + +def array(sequence, typecode=None, copy=1, savespace=0, dtype=None): + dtype = convtypecode2(typecode, dtype) + return mu.array(sequence, dtype, copy=copy) + +def sarray(a, typecode=None, copy=False, dtype=None): + dtype = convtypecode2(typecode, dtype) + return mu.array(a, dtype, copy) + +def asarray(a, typecode=None, dtype=None): + dtype = convtypecode2(typecode, dtype) + return mu.array(a, dtype, copy=0) + +def nonzero(a): + res = np.nonzero(a) + if len(res) == 1: + return res[0] + else: + raise ValueError, "Input argument must be 1d" + +def reshape(a, shape): + return np.reshape(a, shape) + +def arange(start, stop=None, step=1, typecode=None, dtype=None): + dtype = convtypecode2(typecode, dtype) + return mu.arange(start, stop, step, dtype) + +def fromstring(string, typecode='l', count=-1, dtype=None): + dtype = convtypecode(typecode, dtype) + return mu.fromstring(string, dtype, count=count) + +def ravel(m): + return np.ravel(m) + +def trace(a, offset=0, axis1=0, axis2=1): + return np.trace(a, offset=0, axis1=0, axis2=1) + +def indices(dimensions, typecode=None, dtype=None): + dtype = convtypecode(typecode, dtype) + return np.indices(dimensions, dtype) + +def where(condition, x, y): + return np.where(condition, x, y) + +def cross_product(a, b, axis1=-1, axis2=-1): + return np.cross(a, b, axis1, axis2) + +def average(a, axis=0, weights=None, returned=False): + return np.average(a, axis, weights, returned) diff --git a/numpy/oldnumeric/linear_algebra.py b/numpy/oldnumeric/linear_algebra.py new file mode 100644 index 000000000..2e7a264fe --- /dev/null +++ b/numpy/oldnumeric/linear_algebra.py @@ -0,0 +1,83 @@ +"""Backward compatible with LinearAlgebra from Numeric +""" +# This module is a lite version of the linalg.py module in SciPy which contains +# high-level Python interface to the LAPACK library. The lite version +# only accesses the following LAPACK functions: dgesv, zgesv, dgeev, +# zgeev, dgesdd, zgesdd, dgelsd, zgelsd, dsyevd, zheevd, dgetrf, dpotrf. + + +__all__ = ['LinAlgError', 'solve_linear_equations', + 'inverse', 'cholesky_decomposition', 'eigenvalues', + 'Heigenvalues', 'generalized_inverse', + 'determinant', 'singular_value_decomposition', + 'eigenvectors', 'Heigenvectors', + 'linear_least_squares' + ] + +from numpy.core import transpose +import numpy.linalg as linalg + +# Linear equations + +LinAlgError = linalg.LinAlgError + +def solve_linear_equations(a, b): + return linalg.solve(a,b) + +# Matrix inversion + +def inverse(a): + return linalg.inv(a) + +# Cholesky decomposition + +def cholesky_decomposition(a): + return linalg.cholesky(a) + +# Eigenvalues + +def eigenvalues(a): + return linalg.eigvals(a) + +def Heigenvalues(a, UPLO='L'): + return linalg.eigvalsh(a,UPLO) + +# Eigenvectors + +def eigenvectors(A): + w, v = linalg.eig(A) + return w, transpose(v) + +def Heigenvectors(A): + w, v = linalg.eigh(A) + return w, transpose(v) + +# Generalized inverse + +def generalized_inverse(a, rcond = 1.e-10): + return linalg.pinv(a, rcond) + +# Determinant + +def determinant(a): + return linalg.det(a) + +# Linear Least Squares + +def linear_least_squares(a, b, rcond=1.e-10): + """returns x,resids,rank,s +where x minimizes 2-norm(|b - Ax|) + resids is the sum square residuals + rank is the rank of A + s is the rank of the singular values of A in descending order + +If b is a matrix then x is also a matrix with corresponding columns. +If the rank of A is less than the number of columns of A or greater than +the number of rows, then residuals will be returned as an empty array +otherwise resids = sum((b-dot(A,x)**2). +Singular values less than s[0]*rcond are treated as zero. +""" + return linalg.lstsq(a,b,rcond) + +def singular_value_decomposition(A, full_matrices=0): + return linalg.svd(A, full_matrices) diff --git a/numpy/oldnumeric/ma.py b/numpy/oldnumeric/ma.py new file mode 100644 index 000000000..532e4f905 --- /dev/null +++ b/numpy/oldnumeric/ma.py @@ -0,0 +1,2267 @@ +"""MA: a facility for dealing with missing observations +MA is generally used as a numpy.array look-alike. +by Paul F. Dubois. + +Copyright 1999, 2000, 2001 Regents of the University of California. +Released for unlimited redistribution. +Adapted for numpy_core 2005 by Travis Oliphant and +(mainly) Paul Dubois. + +""" +import types, sys + +import numpy.core.umath as umath +import numpy.core.fromnumeric as fromnumeric +from numpy.core.numeric import newaxis, ndarray, inf +from numpy.core.fromnumeric import amax, amin +from numpy.core.numerictypes import bool_, typecodes +import numpy.core.numeric as numeric +import warnings + +# Ufunc domain lookup for __array_wrap__ +ufunc_domain = {} +# Ufunc fills lookup for __array__ +ufunc_fills = {} + +MaskType = bool_ +nomask = MaskType(0) +divide_tolerance = 1.e-35 + +class MAError (Exception): + def __init__ (self, args=None): + "Create an exception" + + # The .args attribute must be a tuple. + if not isinstance(args, tuple): + args = (args,) + self.args = args + def __str__(self): + "Calculate the string representation" + return str(self.args[0]) + __repr__ = __str__ + +class _MaskedPrintOption: + "One instance of this class, masked_print_option, is created." + def __init__ (self, display): + "Create the masked print option object." + self.set_display(display) + self._enabled = 1 + + def display (self): + "Show what prints for masked values." + return self._display + + def set_display (self, s): + "set_display(s) sets what prints for masked values." + self._display = s + + def enabled (self): + "Is the use of the display value enabled?" + return self._enabled + + def enable(self, flag=1): + "Set the enabling flag to flag." + self._enabled = flag + + def __str__ (self): + return str(self._display) + + __repr__ = __str__ + +#if you single index into a masked location you get this object. +masked_print_option = _MaskedPrintOption('--') + +# Use single element arrays or scalars. +default_real_fill_value = 1.e20 +default_complex_fill_value = 1.e20 + 0.0j +default_character_fill_value = '-' +default_integer_fill_value = 999999 +default_object_fill_value = '?' + +def default_fill_value (obj): + "Function to calculate default fill value for an object." + if isinstance(obj, types.FloatType): + return default_real_fill_value + elif isinstance(obj, types.IntType) or isinstance(obj, types.LongType): + return default_integer_fill_value + elif isinstance(obj, types.StringType): + return default_character_fill_value + elif isinstance(obj, types.ComplexType): + return default_complex_fill_value + elif isinstance(obj, MaskedArray) or isinstance(obj, ndarray): + x = obj.dtype.char + if x in typecodes['Float']: + return default_real_fill_value + if x in typecodes['Integer']: + return default_integer_fill_value + if x in typecodes['Complex']: + return default_complex_fill_value + if x in typecodes['Character']: + return default_character_fill_value + if x in typecodes['UnsignedInteger']: + return umath.absolute(default_integer_fill_value) + return default_object_fill_value + else: + return default_object_fill_value + +def minimum_fill_value (obj): + "Function to calculate default fill value suitable for taking minima." + if isinstance(obj, types.FloatType): + return numeric.inf + elif isinstance(obj, types.IntType) or isinstance(obj, types.LongType): + return sys.maxint + elif isinstance(obj, MaskedArray) or isinstance(obj, ndarray): + x = obj.dtype.char + if x in typecodes['Float']: + return numeric.inf + if x in typecodes['Integer']: + return sys.maxint + if x in typecodes['UnsignedInteger']: + return sys.maxint + else: + raise TypeError, 'Unsuitable type for calculating minimum.' + +def maximum_fill_value (obj): + "Function to calculate default fill value suitable for taking maxima." + if isinstance(obj, types.FloatType): + return -inf + elif isinstance(obj, types.IntType) or isinstance(obj, types.LongType): + return -sys.maxint + elif isinstance(obj, MaskedArray) or isinstance(obj, ndarray): + x = obj.dtype.char + if x in typecodes['Float']: + return -inf + if x in typecodes['Integer']: + return -sys.maxint + if x in typecodes['UnsignedInteger']: + return 0 + else: + raise TypeError, 'Unsuitable type for calculating maximum.' + +def set_fill_value (a, fill_value): + "Set fill value of a if it is a masked array." + if isMaskedArray(a): + a.set_fill_value (fill_value) + +def getmask (a): + """Mask of values in a; could be nomask. + Returns nomask if a is not a masked array. + To get an array for sure use getmaskarray.""" + if isinstance(a, MaskedArray): + return a.raw_mask() + else: + return nomask + +def getmaskarray (a): + """Mask of values in a; an array of zeros if mask is nomask + or not a masked array, and is a byte-sized integer. + Do not try to add up entries, for example. + """ + m = getmask(a) + if m is nomask: + return make_mask_none(shape(a)) + else: + return m + +def is_mask (m): + """Is m a legal mask? Does not check contents, only type. + """ + try: + return m.dtype.type is MaskType + except AttributeError: + return False + +def make_mask (m, copy=0, flag=0): + """make_mask(m, copy=0, flag=0) + return m as a mask, creating a copy if necessary or requested. + Can accept any sequence of integers or nomask. Does not check + that contents must be 0s and 1s. + if flag, return nomask if m contains no true elements. + """ + if m is nomask: + return nomask + elif isinstance(m, ndarray): + if m.dtype.type is MaskType: + if copy: + result = numeric.array(m, dtype=MaskType, copy=copy) + else: + result = m + else: + result = m.astype(MaskType) + else: + result = filled(m, True).astype(MaskType) + + if flag and not fromnumeric.sometrue(fromnumeric.ravel(result)): + return nomask + else: + return result + +def make_mask_none (s): + "Return a mask of all zeros of shape s." + result = numeric.zeros(s, dtype=MaskType) + result.shape = s + return result + +def mask_or (m1, m2): + """Logical or of the mask candidates m1 and m2, treating nomask as false. + Result may equal m1 or m2 if the other is nomask. + """ + if m1 is nomask: return make_mask(m2) + if m2 is nomask: return make_mask(m1) + if m1 is m2 and is_mask(m1): return m1 + return make_mask(umath.logical_or(m1, m2)) + +def filled (a, value = None): + """a as a contiguous numeric array with any masked areas replaced by value + if value is None or the special element "masked", get_fill_value(a) + is used instead. + + If a is already a contiguous numeric array, a itself is returned. + + filled(a) can be used to be sure that the result is numeric when + passing an object a to other software ignorant of MA, in particular to + numeric itself. + """ + if isinstance(a, MaskedArray): + return a.filled(value) + elif isinstance(a, ndarray) and a.flags['CONTIGUOUS']: + return a + elif isinstance(a, types.DictType): + return numeric.array(a, 'O') + else: + return numeric.array(a) + +def get_fill_value (a): + """ + The fill value of a, if it has one; otherwise, the default fill value + for that type. + """ + if isMaskedArray(a): + result = a.fill_value() + else: + result = default_fill_value(a) + return result + +def common_fill_value (a, b): + "The common fill_value of a and b, if there is one, or None" + t1 = get_fill_value(a) + t2 = get_fill_value(b) + if t1 == t2: return t1 + return None + +# Domain functions return 1 where the argument(s) are not in the domain. +class domain_check_interval: + "domain_check_interval(a,b)(x) = true where x < a or y > b" + def __init__(self, y1, y2): + "domain_check_interval(a,b)(x) = true where x < a or y > b" + self.y1 = y1 + self.y2 = y2 + + def __call__ (self, x): + "Execute the call behavior." + return umath.logical_or(umath.greater (x, self.y2), + umath.less(x, self.y1) + ) + +class domain_tan: + "domain_tan(eps) = true where abs(cos(x)) < eps)" + def __init__(self, eps): + "domain_tan(eps) = true where abs(cos(x)) < eps)" + self.eps = eps + + def __call__ (self, x): + "Execute the call behavior." + return umath.less(umath.absolute(umath.cos(x)), self.eps) + +class domain_greater: + "domain_greater(v)(x) = true where x <= v" + def __init__(self, critical_value): + "domain_greater(v)(x) = true where x <= v" + self.critical_value = critical_value + + def __call__ (self, x): + "Execute the call behavior." + return umath.less_equal (x, self.critical_value) + +class domain_greater_equal: + "domain_greater_equal(v)(x) = true where x < v" + def __init__(self, critical_value): + "domain_greater_equal(v)(x) = true where x < v" + self.critical_value = critical_value + + def __call__ (self, x): + "Execute the call behavior." + return umath.less (x, self.critical_value) + +class masked_unary_operation: + def __init__ (self, aufunc, fill=0, domain=None): + """ masked_unary_operation(aufunc, fill=0, domain=None) + aufunc(fill) must be defined + self(x) returns aufunc(x) + with masked values where domain(x) is true or getmask(x) is true. + """ + self.f = aufunc + self.fill = fill + self.domain = domain + self.__doc__ = getattr(aufunc, "__doc__", str(aufunc)) + self.__name__ = getattr(aufunc, "__name__", str(aufunc)) + ufunc_domain[aufunc] = domain + ufunc_fills[aufunc] = fill, + + def __call__ (self, a, *args, **kwargs): + "Execute the call behavior." +# numeric tries to return scalars rather than arrays when given scalars. + m = getmask(a) + d1 = filled(a, self.fill) + if self.domain is not None: + m = mask_or(m, self.domain(d1)) + result = self.f(d1, *args, **kwargs) + return masked_array(result, m) + + def __str__ (self): + return "Masked version of " + str(self.f) + + +class domain_safe_divide: + def __init__ (self, tolerance=divide_tolerance): + self.tolerance = tolerance + def __call__ (self, a, b): + return umath.absolute(a) * self.tolerance >= umath.absolute(b) + +class domained_binary_operation: + """Binary operations that have a domain, like divide. These are complicated + so they are a separate class. They have no reduce, outer or accumulate. + """ + def __init__ (self, abfunc, domain, fillx=0, filly=0): + """abfunc(fillx, filly) must be defined. + abfunc(x, filly) = x for all x to enable reduce. + """ + self.f = abfunc + self.domain = domain + self.fillx = fillx + self.filly = filly + self.__doc__ = getattr(abfunc, "__doc__", str(abfunc)) + self.__name__ = getattr(abfunc, "__name__", str(abfunc)) + ufunc_domain[abfunc] = domain + ufunc_fills[abfunc] = fillx, filly + + def __call__(self, a, b): + "Execute the call behavior." + ma = getmask(a) + mb = getmask(b) + d1 = filled(a, self.fillx) + d2 = filled(b, self.filly) + t = self.domain(d1, d2) + + if fromnumeric.sometrue(t, None): + d2 = where(t, self.filly, d2) + mb = mask_or(mb, t) + m = mask_or(ma, mb) + result = self.f(d1, d2) + return masked_array(result, m) + + def __str__ (self): + return "Masked version of " + str(self.f) + +class masked_binary_operation: + def __init__ (self, abfunc, fillx=0, filly=0): + """abfunc(fillx, filly) must be defined. + abfunc(x, filly) = x for all x to enable reduce. + """ + self.f = abfunc + self.fillx = fillx + self.filly = filly + self.__doc__ = getattr(abfunc, "__doc__", str(abfunc)) + ufunc_domain[abfunc] = None + ufunc_fills[abfunc] = fillx, filly + + def __call__ (self, a, b, *args, **kwargs): + "Execute the call behavior." + m = mask_or(getmask(a), getmask(b)) + d1 = filled(a, self.fillx) + d2 = filled(b, self.filly) + result = self.f(d1, d2, *args, **kwargs) + if isinstance(result, ndarray) \ + and m.ndim != 0 \ + and m.shape != result.shape: + m = mask_or(getmaskarray(a), getmaskarray(b)) + return masked_array(result, m) + + def reduce (self, target, axis=0, dtype=None): + """Reduce target along the given axis with this function.""" + m = getmask(target) + t = filled(target, self.filly) + if t.shape == (): + t = t.reshape(1) + if m is not nomask: + m = make_mask(m, copy=1) + m.shape = (1,) + if m is nomask: + t = self.f.reduce(t, axis) + else: + t = masked_array (t, m) + # XXX: "or t.dtype" below is a workaround for what appears + # XXX: to be a bug in reduce. + t = self.f.reduce(filled(t, self.filly), axis, + dtype=dtype or t.dtype) + m = umath.logical_and.reduce(m, axis) + if isinstance(t, ndarray): + return masked_array(t, m, get_fill_value(target)) + elif m: + return masked + else: + return t + + def outer (self, a, b): + "Return the function applied to the outer product of a and b." + ma = getmask(a) + mb = getmask(b) + if ma is nomask and mb is nomask: + m = nomask + else: + ma = getmaskarray(a) + mb = getmaskarray(b) + m = logical_or.outer(ma, mb) + d = self.f.outer(filled(a, self.fillx), filled(b, self.filly)) + return masked_array(d, m) + + def accumulate (self, target, axis=0): + """Accumulate target along axis after filling with y fill value.""" + t = filled(target, self.filly) + return masked_array (self.f.accumulate (t, axis)) + def __str__ (self): + return "Masked version of " + str(self.f) + +sqrt = masked_unary_operation(umath.sqrt, 0.0, domain_greater_equal(0.0)) +log = masked_unary_operation(umath.log, 1.0, domain_greater(0.0)) +log10 = masked_unary_operation(umath.log10, 1.0, domain_greater(0.0)) +exp = masked_unary_operation(umath.exp) +conjugate = masked_unary_operation(umath.conjugate) +sin = masked_unary_operation(umath.sin) +cos = masked_unary_operation(umath.cos) +tan = masked_unary_operation(umath.tan, 0.0, domain_tan(1.e-35)) +arcsin = masked_unary_operation(umath.arcsin, 0.0, domain_check_interval(-1.0, 1.0)) +arccos = masked_unary_operation(umath.arccos, 0.0, domain_check_interval(-1.0, 1.0)) +arctan = masked_unary_operation(umath.arctan) +# Missing from numeric +arcsinh = masked_unary_operation(umath.arcsinh) +arccosh = masked_unary_operation(umath.arccosh, 1.0, domain_greater_equal(1.0)) +arctanh = masked_unary_operation(umath.arctanh, 0.0, domain_check_interval(-1.0+1e-15, 1.0-1e-15)) +sinh = masked_unary_operation(umath.sinh) +cosh = masked_unary_operation(umath.cosh) +tanh = masked_unary_operation(umath.tanh) +absolute = masked_unary_operation(umath.absolute) +fabs = masked_unary_operation(umath.fabs) +negative = masked_unary_operation(umath.negative) + +def nonzero(a): + """returns the indices of the elements of a which are not zero + and not masked + """ + return numeric.asarray(filled(a, 0).nonzero()) + +around = masked_unary_operation(fromnumeric.round_) +floor = masked_unary_operation(umath.floor) +ceil = masked_unary_operation(umath.ceil) +logical_not = masked_unary_operation(umath.logical_not) + +add = masked_binary_operation(umath.add) +subtract = masked_binary_operation(umath.subtract) +subtract.reduce = None +multiply = masked_binary_operation(umath.multiply, 1, 1) +divide = domained_binary_operation(umath.divide, domain_safe_divide(), 0, 1) +true_divide = domained_binary_operation(umath.true_divide, domain_safe_divide(), 0, 1) +floor_divide = domained_binary_operation(umath.floor_divide, domain_safe_divide(), 0, 1) +remainder = domained_binary_operation(umath.remainder, domain_safe_divide(), 0, 1) +fmod = domained_binary_operation(umath.fmod, domain_safe_divide(), 0, 1) +hypot = masked_binary_operation(umath.hypot) +arctan2 = masked_binary_operation(umath.arctan2, 0.0, 1.0) +arctan2.reduce = None +equal = masked_binary_operation(umath.equal) +equal.reduce = None +not_equal = masked_binary_operation(umath.not_equal) +not_equal.reduce = None +less_equal = masked_binary_operation(umath.less_equal) +less_equal.reduce = None +greater_equal = masked_binary_operation(umath.greater_equal) +greater_equal.reduce = None +less = masked_binary_operation(umath.less) +less.reduce = None +greater = masked_binary_operation(umath.greater) +greater.reduce = None +logical_and = masked_binary_operation(umath.logical_and) +alltrue = masked_binary_operation(umath.logical_and, 1, 1).reduce +logical_or = masked_binary_operation(umath.logical_or) +sometrue = logical_or.reduce +logical_xor = masked_binary_operation(umath.logical_xor) +bitwise_and = masked_binary_operation(umath.bitwise_and) +bitwise_or = masked_binary_operation(umath.bitwise_or) +bitwise_xor = masked_binary_operation(umath.bitwise_xor) + +def rank (object): + return fromnumeric.rank(filled(object)) + +def shape (object): + return fromnumeric.shape(filled(object)) + +def size (object, axis=None): + return fromnumeric.size(filled(object), axis) + +class MaskedArray (object): + """Arrays with possibly masked values. + Masked values of 1 exclude the corresponding element from + any computation. + + Construction: + x = array(data, dtype=None, copy=True, order=False, + mask = nomask, fill_value=None) + + If copy=False, every effort is made not to copy the data: + If data is a MaskedArray, and argument mask=nomask, + then the candidate data is data.data and the + mask used is data.mask. If data is a numeric array, + it is used as the candidate raw data. + If dtype is not None and + is != data.dtype.char then a data copy is required. + Otherwise, the candidate is used. + + If a data copy is required, raw data stored is the result of: + numeric.array(data, dtype=dtype.char, copy=copy) + + If mask is nomask there are no masked values. Otherwise mask must + be convertible to an array of booleans with the same shape as x. + + fill_value is used to fill in masked values when necessary, + such as when printing and in method/function filled(). + The fill_value is not used for computation within this module. + """ + __array_priority__ = 10.1 + def __init__(self, data, dtype=None, copy=True, order=False, + mask=nomask, fill_value=None): + """array(data, dtype=None, copy=True, order=False, mask=nomask, fill_value=None) + If data already a numeric array, its dtype becomes the default value of dtype. + """ + if dtype is None: + tc = None + else: + tc = numeric.dtype(dtype) + need_data_copied = copy + if isinstance(data, MaskedArray): + c = data.data + if tc is None: + tc = c.dtype + elif tc != c.dtype: + need_data_copied = True + if mask is nomask: + mask = data.mask + elif mask is not nomask: #attempting to change the mask + need_data_copied = True + + elif isinstance(data, ndarray): + c = data + if tc is None: + tc = c.dtype + elif tc != c.dtype: + need_data_copied = True + else: + need_data_copied = False #because I'll do it now + c = numeric.array(data, dtype=tc, copy=True, order=order) + tc = c.dtype + + if need_data_copied: + if tc == c.dtype: + self._data = numeric.array(c, dtype=tc, copy=True, order=order) + else: + self._data = c.astype(tc) + else: + self._data = c + + if mask is nomask: + self._mask = nomask + self._shared_mask = 0 + else: + self._mask = make_mask (mask) + if self._mask is nomask: + self._shared_mask = 0 + else: + self._shared_mask = (self._mask is mask) + nm = size(self._mask) + nd = size(self._data) + if nm != nd: + if nm == 1: + self._mask = fromnumeric.resize(self._mask, self._data.shape) + self._shared_mask = 0 + elif nd == 1: + self._data = fromnumeric.resize(self._data, self._mask.shape) + self._data.shape = self._mask.shape + else: + raise MAError, "Mask and data not compatible." + elif nm == 1 and shape(self._mask) != shape(self._data): + self.unshare_mask() + self._mask.shape = self._data.shape + + self.set_fill_value(fill_value) + + def __array__ (self, t=None, context=None): + "Special hook for numeric. Converts to numeric if possible." + if self._mask is not nomask: + if fromnumeric.ravel(self._mask).any(): + if context is None: + warnings.warn("Cannot automatically convert masked array to "\ + "numeric because data\n is masked in one or "\ + "more locations."); + return self._data + #raise MAError, \ + # """Cannot automatically convert masked array to numeric because data + # is masked in one or more locations. + # """ + else: + func, args, i = context + fills = ufunc_fills.get(func) + if fills is None: + raise MAError, "%s not known to ma" % func + return self.filled(fills[i]) + else: # Mask is all false + # Optimize to avoid future invocations of this section. + self._mask = nomask + self._shared_mask = 0 + if t: + return self._data.astype(t) + else: + return self._data + + def __array_wrap__ (self, array, context=None): + """Special hook for ufuncs. + + Wraps the numpy array and sets the mask according to + context. + """ + if context is None: + return MaskedArray(array, copy=False, mask=nomask) + func, args = context[:2] + domain = ufunc_domain[func] + m = reduce(mask_or, [getmask(a) for a in args]) + if domain is not None: + m = mask_or(m, domain(*[getattr(a, '_data', a) + for a in args])) + if m is not nomask: + try: + shape = array.shape + except AttributeError: + pass + else: + if m.shape != shape: + m = reduce(mask_or, [getmaskarray(a) for a in args]) + + return MaskedArray(array, copy=False, mask=m) + + def _get_shape(self): + "Return the current shape." + return self._data.shape + + def _set_shape (self, newshape): + "Set the array's shape." + self._data.shape = newshape + if self._mask is not nomask: + self._mask = self._mask.copy() + self._mask.shape = newshape + + def _get_flat(self): + """Calculate the flat value. + """ + if self._mask is nomask: + return masked_array(self._data.ravel(), mask=nomask, + fill_value = self.fill_value()) + else: + return masked_array(self._data.ravel(), + mask=self._mask.ravel(), + fill_value = self.fill_value()) + + def _set_flat (self, value): + "x.flat = value" + y = self.ravel() + y[:] = value + + def _get_real(self): + "Get the real part of a complex array." + if self._mask is nomask: + return masked_array(self._data.real, mask=nomask, + fill_value = self.fill_value()) + else: + return masked_array(self._data.real, mask=self._mask, + fill_value = self.fill_value()) + + def _set_real (self, value): + "x.real = value" + y = self.real + y[...] = value + + def _get_imaginary(self): + "Get the imaginary part of a complex array." + if self._mask is nomask: + return masked_array(self._data.imag, mask=nomask, + fill_value = self.fill_value()) + else: + return masked_array(self._data.imag, mask=self._mask, + fill_value = self.fill_value()) + + def _set_imaginary (self, value): + "x.imaginary = value" + y = self.imaginary + y[...] = value + + def __str__(self): + """Calculate the str representation, using masked for fill if + it is enabled. Otherwise fill with fill value. + """ + if masked_print_option.enabled(): + f = masked_print_option + # XXX: Without the following special case masked + # XXX: would print as "[--]", not "--". Can we avoid + # XXX: checks for masked by choosing a different value + # XXX: for the masked singleton? 2005-01-05 -- sasha + if self is masked: + return str(f) + m = self._mask + if m is not nomask and m.shape == () and m: + return str(f) + # convert to object array to make filled work + self = self.astype(object) + else: + f = self.fill_value() + res = self.filled(f) + return str(res) + + def __repr__(self): + """Calculate the repr representation, using masked for fill if + it is enabled. Otherwise fill with fill value. + """ + with_mask = """\ +array(data = + %(data)s, + mask = + %(mask)s, + fill_value=%(fill)s) +""" + with_mask1 = """\ +array(data = %(data)s, + mask = %(mask)s, + fill_value=%(fill)s) +""" + without_mask = """array( + %(data)s)""" + without_mask1 = """array(%(data)s)""" + + n = len(self.shape) + if self._mask is nomask: + if n <= 1: + return without_mask1 % {'data':str(self.filled())} + return without_mask % {'data':str(self.filled())} + else: + if n <= 1: + return with_mask % { + 'data': str(self.filled()), + 'mask': str(self._mask), + 'fill': str(self.fill_value()) + } + return with_mask % { + 'data': str(self.filled()), + 'mask': str(self._mask), + 'fill': str(self.fill_value()) + } + without_mask1 = """array(%(data)s)""" + if self._mask is nomask: + return without_mask % {'data':str(self.filled())} + else: + return with_mask % { + 'data': str(self.filled()), + 'mask': str(self._mask), + 'fill': str(self.fill_value()) + } + + def __float__(self): + "Convert self to float." + self.unmask() + if self._mask is not nomask: + raise MAError, 'Cannot convert masked element to a Python float.' + return float(self.data.item()) + + def __int__(self): + "Convert self to int." + self.unmask() + if self._mask is not nomask: + raise MAError, 'Cannot convert masked element to a Python int.' + return int(self.data.item()) + + def __getitem__(self, i): + "Get item described by i. Not a copy as in previous versions." + self.unshare_mask() + m = self._mask + dout = self._data[i] + if m is nomask: + try: + if dout.size == 1: + return dout + else: + return masked_array(dout, fill_value=self._fill_value) + except AttributeError: + return dout + mi = m[i] + if mi.size == 1: + if mi: + return masked + else: + return dout + else: + return masked_array(dout, mi, fill_value=self._fill_value) + +# -------- +# setitem and setslice notes +# note that if value is masked, it means to mask those locations. +# setting a value changes the mask to match the value in those locations. + + def __setitem__(self, index, value): + "Set item described by index. If value is masked, mask those locations." + d = self._data + if self is masked: + raise MAError, 'Cannot alter masked elements.' + if value is masked: + if self._mask is nomask: + self._mask = make_mask_none(d.shape) + self._shared_mask = False + else: + self.unshare_mask() + self._mask[index] = True + return + m = getmask(value) + value = filled(value).astype(d.dtype) + d[index] = value + if m is nomask: + if self._mask is not nomask: + self.unshare_mask() + self._mask[index] = False + else: + if self._mask is nomask: + self._mask = make_mask_none(d.shape) + self._shared_mask = True + else: + self.unshare_mask() + self._mask[index] = m + + def __nonzero__(self): + """returns true if any element is non-zero or masked + + """ + # XXX: This changes bool conversion logic from MA. + # XXX: In MA bool(a) == len(a) != 0, but in numpy + # XXX: scalars do not have len + m = self._mask + d = self._data + return bool(m is not nomask and m.any() + or d is not nomask and d.any()) + + def __len__ (self): + """Return length of first dimension. This is weird but Python's + slicing behavior depends on it.""" + return len(self._data) + + def __and__(self, other): + "Return bitwise_and" + return bitwise_and(self, other) + + def __or__(self, other): + "Return bitwise_or" + return bitwise_or(self, other) + + def __xor__(self, other): + "Return bitwise_xor" + return bitwise_xor(self, other) + + __rand__ = __and__ + __ror__ = __or__ + __rxor__ = __xor__ + + def __abs__(self): + "Return absolute(self)" + return absolute(self) + + def __neg__(self): + "Return negative(self)" + return negative(self) + + def __pos__(self): + "Return array(self)" + return array(self) + + def __add__(self, other): + "Return add(self, other)" + return add(self, other) + + __radd__ = __add__ + + def __mod__ (self, other): + "Return remainder(self, other)" + return remainder(self, other) + + def __rmod__ (self, other): + "Return remainder(other, self)" + return remainder(other, self) + + def __lshift__ (self, n): + return left_shift(self, n) + + def __rshift__ (self, n): + return right_shift(self, n) + + def __sub__(self, other): + "Return subtract(self, other)" + return subtract(self, other) + + def __rsub__(self, other): + "Return subtract(other, self)" + return subtract(other, self) + + def __mul__(self, other): + "Return multiply(self, other)" + return multiply(self, other) + + __rmul__ = __mul__ + + def __div__(self, other): + "Return divide(self, other)" + return divide(self, other) + + def __rdiv__(self, other): + "Return divide(other, self)" + return divide(other, self) + + def __truediv__(self, other): + "Return divide(self, other)" + return true_divide(self, other) + + def __rtruediv__(self, other): + "Return divide(other, self)" + return true_divide(other, self) + + def __floordiv__(self, other): + "Return divide(self, other)" + return floor_divide(self, other) + + def __rfloordiv__(self, other): + "Return divide(other, self)" + return floor_divide(other, self) + + def __pow__(self, other, third=None): + "Return power(self, other, third)" + return power(self, other, third) + + def __sqrt__(self): + "Return sqrt(self)" + return sqrt(self) + + def __iadd__(self, other): + "Add other to self in place." + t = self._data.dtype.char + f = filled(other, 0) + t1 = f.dtype.char + if t == t1: + pass + elif t in typecodes['Integer']: + if t1 in typecodes['Integer']: + f = f.astype(t) + else: + raise TypeError, 'Incorrect type for in-place operation.' + elif t in typecodes['Float']: + if t1 in typecodes['Integer']: + f = f.astype(t) + elif t1 in typecodes['Float']: + f = f.astype(t) + else: + raise TypeError, 'Incorrect type for in-place operation.' + elif t in typecodes['Complex']: + if t1 in typecodes['Integer']: + f = f.astype(t) + elif t1 in typecodes['Float']: + f = f.astype(t) + elif t1 in typecodes['Complex']: + f = f.astype(t) + else: + raise TypeError, 'Incorrect type for in-place operation.' + else: + raise TypeError, 'Incorrect type for in-place operation.' + + if self._mask is nomask: + self._data += f + m = getmask(other) + self._mask = m + self._shared_mask = m is not nomask + else: + result = add(self, masked_array(f, mask=getmask(other))) + self._data = result.data + self._mask = result.mask + self._shared_mask = 1 + return self + + def __imul__(self, other): + "Add other to self in place." + t = self._data.dtype.char + f = filled(other, 0) + t1 = f.dtype.char + if t == t1: + pass + elif t in typecodes['Integer']: + if t1 in typecodes['Integer']: + f = f.astype(t) + else: + raise TypeError, 'Incorrect type for in-place operation.' + elif t in typecodes['Float']: + if t1 in typecodes['Integer']: + f = f.astype(t) + elif t1 in typecodes['Float']: + f = f.astype(t) + else: + raise TypeError, 'Incorrect type for in-place operation.' + elif t in typecodes['Complex']: + if t1 in typecodes['Integer']: + f = f.astype(t) + elif t1 in typecodes['Float']: + f = f.astype(t) + elif t1 in typecodes['Complex']: + f = f.astype(t) + else: + raise TypeError, 'Incorrect type for in-place operation.' + else: + raise TypeError, 'Incorrect type for in-place operation.' + + if self._mask is nomask: + self._data *= f + m = getmask(other) + self._mask = m + self._shared_mask = m is not nomask + else: + result = multiply(self, masked_array(f, mask=getmask(other))) + self._data = result.data + self._mask = result.mask + self._shared_mask = 1 + return self + + def __isub__(self, other): + "Subtract other from self in place." + t = self._data.dtype.char + f = filled(other, 0) + t1 = f.dtype.char + if t == t1: + pass + elif t in typecodes['Integer']: + if t1 in typecodes['Integer']: + f = f.astype(t) + else: + raise TypeError, 'Incorrect type for in-place operation.' + elif t in typecodes['Float']: + if t1 in typecodes['Integer']: + f = f.astype(t) + elif t1 in typecodes['Float']: + f = f.astype(t) + else: + raise TypeError, 'Incorrect type for in-place operation.' + elif t in typecodes['Complex']: + if t1 in typecodes['Integer']: + f = f.astype(t) + elif t1 in typecodes['Float']: + f = f.astype(t) + elif t1 in typecodes['Complex']: + f = f.astype(t) + else: + raise TypeError, 'Incorrect type for in-place operation.' + else: + raise TypeError, 'Incorrect type for in-place operation.' + + if self._mask is nomask: + self._data -= f + m = getmask(other) + self._mask = m + self._shared_mask = m is not nomask + else: + result = subtract(self, masked_array(f, mask=getmask(other))) + self._data = result.data + self._mask = result.mask + self._shared_mask = 1 + return self + + + + def __idiv__(self, other): + "Divide self by other in place." + t = self._data.dtype.char + f = filled(other, 0) + t1 = f.dtype.char + if t == t1: + pass + elif t in typecodes['Integer']: + if t1 in typecodes['Integer']: + f = f.astype(t) + else: + raise TypeError, 'Incorrect type for in-place operation.' + elif t in typecodes['Float']: + if t1 in typecodes['Integer']: + f = f.astype(t) + elif t1 in typecodes['Float']: + f = f.astype(t) + else: + raise TypeError, 'Incorrect type for in-place operation.' + elif t in typecodes['Complex']: + if t1 in typecodes['Integer']: + f = f.astype(t) + elif t1 in typecodes['Float']: + f = f.astype(t) + elif t1 in typecodes['Complex']: + f = f.astype(t) + else: + raise TypeError, 'Incorrect type for in-place operation.' + else: + raise TypeError, 'Incorrect type for in-place operation.' + mo = getmask(other) + result = divide(self, masked_array(f, mask=mo)) + self._data = result.data + dm = result.raw_mask() + if dm is not self._mask: + self._mask = dm + self._shared_mask = 1 + return self + + def __eq__(self, other): + return equal(self,other) + + def __ne__(self, other): + return not_equal(self,other) + + def __lt__(self, other): + return less(self,other) + + def __le__(self, other): + return less_equal(self,other) + + def __gt__(self, other): + return greater(self,other) + + def __ge__(self, other): + return greater_equal(self,other) + + def astype (self, tc): + "return self as array of given type." + d = self._data.astype(tc) + return array(d, mask=self._mask) + + def byte_swapped(self): + """Returns the raw data field, byte_swapped. Included for consistency + with numeric but doesn't make sense in this context. + """ + return self._data.byte_swapped() + + def compressed (self): + "A 1-D array of all the non-masked data." + d = fromnumeric.ravel(self._data) + if self._mask is nomask: + return array(d) + else: + m = 1 - fromnumeric.ravel(self._mask) + c = fromnumeric.compress(m, d) + return array(c, copy=0) + + def count (self, axis = None): + "Count of the non-masked elements in a, or along a certain axis." + m = self._mask + s = self._data.shape + ls = len(s) + if m is nomask: + if ls == 0: + return 1 + if ls == 1: + return s[0] + if axis is None: + return reduce(lambda x, y:x*y, s) + else: + n = s[axis] + t = list(s) + del t[axis] + return ones(t) * n + if axis is None: + w = fromnumeric.ravel(m).astype(int) + n1 = size(w) + if n1 == 1: + n2 = w[0] + else: + n2 = umath.add.reduce(w) + return n1 - n2 + else: + n1 = size(m, axis) + n2 = sum(m.astype(int), axis) + return n1 - n2 + + def dot (self, other): + "s.dot(other) = innerproduct(s, other)" + return innerproduct(self, other) + + def fill_value(self): + "Get the current fill value." + return self._fill_value + + def filled (self, fill_value=None): + """A numeric array with masked values filled. If fill_value is None, + use self.fill_value(). + + If mask is nomask, copy data only if not contiguous. + Result is always a contiguous, numeric array. +# Is contiguous really necessary now? + """ + d = self._data + m = self._mask + if m is nomask: + if d.flags['CONTIGUOUS']: + return d + else: + return d.copy() + else: + if fill_value is None: + value = self._fill_value + else: + value = fill_value + + if self is masked: + result = numeric.array(value) + else: + try: + result = numeric.array(d, dtype=d.dtype, copy=1) + result[m] = value + except (TypeError, AttributeError): + #ok, can't put that value in here + value = numeric.array(value, dtype=object) + d = d.astype(object) + result = fromnumeric.choose(m, (d, value)) + return result + + def ids (self): + """Return the ids of the data and mask areas""" + return (id(self._data), id(self._mask)) + + def iscontiguous (self): + "Is the data contiguous?" + return self._data.flags['CONTIGUOUS'] + + def itemsize(self): + "Item size of each data item." + return self._data.itemsize + + + def outer(self, other): + "s.outer(other) = outerproduct(s, other)" + return outerproduct(self, other) + + def put (self, values): + """Set the non-masked entries of self to filled(values). + No change to mask + """ + iota = numeric.arange(self.size) + d = self._data + if self._mask is nomask: + ind = iota + else: + ind = fromnumeric.compress(1 - self._mask, iota) + d[ind] = filled(values).astype(d.dtype) + + def putmask (self, values): + """Set the masked entries of self to filled(values). + Mask changed to nomask. + """ + d = self._data + if self._mask is not nomask: + d[self._mask] = filled(values).astype(d.dtype) + self._shared_mask = 0 + self._mask = nomask + + def ravel (self): + """Return a 1-D view of self.""" + if self._mask is nomask: + return masked_array(self._data.ravel()) + else: + return masked_array(self._data.ravel(), self._mask.ravel()) + + def raw_data (self): + """ Obsolete; use data property instead. + The raw data; portions may be meaningless. + May be noncontiguous. Expert use only.""" + return self._data + data = property(fget=raw_data, + doc="The data, but values at masked locations are meaningless.") + + def raw_mask (self): + """ Obsolete; use mask property instead. + May be noncontiguous. Expert use only. + """ + return self._mask + mask = property(fget=raw_mask, + doc="The mask, may be nomask. Values where mask true are meaningless.") + + def reshape (self, *s): + """This array reshaped to shape s""" + d = self._data.reshape(*s) + if self._mask is nomask: + return masked_array(d) + else: + m = self._mask.reshape(*s) + return masked_array(d, m) + + def set_fill_value (self, v=None): + "Set the fill value to v. Omit v to restore default." + if v is None: + v = default_fill_value (self.raw_data()) + self._fill_value = v + + def _get_ndim(self): + return self._data.ndim + ndim = property(_get_ndim, doc=numeric.ndarray.ndim.__doc__) + + def _get_size (self): + return self._data.size + size = property(fget=_get_size, doc="Number of elements in the array.") +## CHECK THIS: signature of numeric.array.size? + + def _get_dtype(self): + return self._data.dtype + dtype = property(fget=_get_dtype, doc="type of the array elements.") + + def item(self, *args): + "Return Python scalar if possible" + if self._mask is not nomask: + m = self._mask.item(*args) + try: + if m[0]: + return masked + except IndexError: + return masked + return self._data.item(*args) + + def itemset(self, *args): + "Set Python scalar into array" + item = args[-1] + args = args[:-1] + self[args] = item + + def tolist(self, fill_value=None): + "Convert to list" + return self.filled(fill_value).tolist() + + def tostring(self, fill_value=None): + "Convert to string" + return self.filled(fill_value).tostring() + + def unmask (self): + "Replace the mask by nomask if possible." + if self._mask is nomask: return + m = make_mask(self._mask, flag=1) + if m is nomask: + self._mask = nomask + self._shared_mask = 0 + + def unshare_mask (self): + "If currently sharing mask, make a copy." + if self._shared_mask: + self._mask = make_mask (self._mask, copy=1, flag=0) + self._shared_mask = 0 + + def _get_ctypes(self): + return self._data.ctypes + + def _get_T(self): + if (self.ndim < 2): + return self + return self.transpose() + + shape = property(_get_shape, _set_shape, + doc = 'tuple giving the shape of the array') + + flat = property(_get_flat, _set_flat, + doc = 'Access array in flat form.') + + real = property(_get_real, _set_real, + doc = 'Access the real part of the array') + + imaginary = property(_get_imaginary, _set_imaginary, + doc = 'Access the imaginary part of the array') + + imag = imaginary + + ctypes = property(_get_ctypes, None, doc="ctypes") + + T = property(_get_T, None, doc="get transpose") + +#end class MaskedArray + +array = MaskedArray + +def isMaskedArray (x): + "Is x a masked array, that is, an instance of MaskedArray?" + return isinstance(x, MaskedArray) + +isarray = isMaskedArray +isMA = isMaskedArray #backward compatibility + +def allclose (a, b, fill_value=1, rtol=1.e-5, atol=1.e-8): + """ Returns true if all components of a and b are equal + subject to given tolerances. + If fill_value is 1, masked values considered equal. + If fill_value is 0, masked values considered unequal. + The relative error rtol should be positive and << 1.0 + The absolute error atol comes into play for those elements + of b that are very small or zero; it says how small a must be also. + """ + m = mask_or(getmask(a), getmask(b)) + d1 = filled(a) + d2 = filled(b) + x = filled(array(d1, copy=0, mask=m), fill_value).astype(float) + y = filled(array(d2, copy=0, mask=m), 1).astype(float) + d = umath.less_equal(umath.absolute(x-y), atol + rtol * umath.absolute(y)) + return fromnumeric.alltrue(fromnumeric.ravel(d)) + +def allequal (a, b, fill_value=1): + """ + True if all entries of a and b are equal, using + fill_value as a truth value where either or both are masked. + """ + m = mask_or(getmask(a), getmask(b)) + if m is nomask: + x = filled(a) + y = filled(b) + d = umath.equal(x, y) + return fromnumeric.alltrue(fromnumeric.ravel(d)) + elif fill_value: + x = filled(a) + y = filled(b) + d = umath.equal(x, y) + dm = array(d, mask=m, copy=0) + return fromnumeric.alltrue(fromnumeric.ravel(filled(dm, 1))) + else: + return 0 + +def masked_values (data, value, rtol=1.e-5, atol=1.e-8, copy=1): + """ + masked_values(data, value, rtol=1.e-5, atol=1.e-8) + Create a masked array; mask is nomask if possible. + If copy==0, and otherwise possible, result + may share data values with original array. + Let d = filled(data, value). Returns d + masked where abs(data-value)<= atol + rtol * abs(value) + if d is of a floating point type. Otherwise returns + masked_object(d, value, copy) + """ + abs = umath.absolute + d = filled(data, value) + if issubclass(d.dtype.type, numeric.floating): + m = umath.less_equal(abs(d-value), atol+rtol*abs(value)) + m = make_mask(m, flag=1) + return array(d, mask = m, copy=copy, + fill_value=value) + else: + return masked_object(d, value, copy=copy) + +def masked_object (data, value, copy=1): + "Create array masked where exactly data equal to value" + d = filled(data, value) + dm = make_mask(umath.equal(d, value), flag=1) + return array(d, mask=dm, copy=copy, fill_value=value) + +def arange(start, stop=None, step=1, dtype=None): + """Just like range() except it returns a array whose type can be specified + by the keyword argument dtype. + """ + return array(numeric.arange(start, stop, step, dtype)) + +arrayrange = arange + +def fromstring (s, t): + "Construct a masked array from a string. Result will have no mask." + return masked_array(numeric.fromstring(s, t)) + +def left_shift (a, n): + "Left shift n bits" + m = getmask(a) + if m is nomask: + d = umath.left_shift(filled(a), n) + return masked_array(d) + else: + d = umath.left_shift(filled(a, 0), n) + return masked_array(d, m) + +def right_shift (a, n): + "Right shift n bits" + m = getmask(a) + if m is nomask: + d = umath.right_shift(filled(a), n) + return masked_array(d) + else: + d = umath.right_shift(filled(a, 0), n) + return masked_array(d, m) + +def resize (a, new_shape): + """resize(a, new_shape) returns a new array with the specified shape. + The original array's total size can be any size.""" + m = getmask(a) + if m is not nomask: + m = fromnumeric.resize(m, new_shape) + result = array(fromnumeric.resize(filled(a), new_shape), mask=m) + result.set_fill_value(get_fill_value(a)) + return result + +def new_repeat(a, repeats, axis=None): + """repeat elements of a repeats times along axis + repeats is a sequence of length a.shape[axis] + telling how many times to repeat each element. + """ + af = filled(a) + if isinstance(repeats, types.IntType): + if axis is None: + num = af.size + else: + num = af.shape[axis] + repeats = tuple([repeats]*num) + + m = getmask(a) + if m is not nomask: + m = fromnumeric.repeat(m, repeats, axis) + d = fromnumeric.repeat(af, repeats, axis) + result = masked_array(d, m) + result.set_fill_value(get_fill_value(a)) + return result + + + +def identity(n): + """identity(n) returns the identity matrix of shape n x n. + """ + return array(numeric.identity(n)) + +def indices (dimensions, dtype=None): + """indices(dimensions,dtype=None) returns an array representing a grid + of indices with row-only, and column-only variation. + """ + return array(numeric.indices(dimensions, dtype)) + +def zeros (shape, dtype=float): + """zeros(n, dtype=float) = + an array of all zeros of the given length or shape.""" + return array(numeric.zeros(shape, dtype)) + +def ones (shape, dtype=float): + """ones(n, dtype=float) = + an array of all ones of the given length or shape.""" + return array(numeric.ones(shape, dtype)) + +def count (a, axis = None): + "Count of the non-masked elements in a, or along a certain axis." + a = masked_array(a) + return a.count(axis) + +def power (a, b, third=None): + "a**b" + if third is not None: + raise MAError, "3-argument power not supported." + ma = getmask(a) + mb = getmask(b) + m = mask_or(ma, mb) + fa = filled(a, 1) + fb = filled(b, 1) + if fb.dtype.char in typecodes["Integer"]: + return masked_array(umath.power(fa, fb), m) + md = make_mask(umath.less(fa, 0), flag=1) + m = mask_or(m, md) + if m is nomask: + return masked_array(umath.power(fa, fb)) + else: + fa = numeric.where(m, 1, fa) + return masked_array(umath.power(fa, fb), m) + +def masked_array (a, mask=nomask, fill_value=None): + """masked_array(a, mask=nomask) = + array(a, mask=mask, copy=0, fill_value=fill_value) + """ + return array(a, mask=mask, copy=0, fill_value=fill_value) + +def sum (target, axis=None, dtype=None): + if axis is None: + target = ravel(target) + axis = 0 + return add.reduce(target, axis, dtype) + +def product (target, axis=None, dtype=None): + if axis is None: + target = ravel(target) + axis = 0 + return multiply.reduce(target, axis, dtype) + +def new_average (a, axis=None, weights=None, returned = 0): + """average(a, axis=None, weights=None) + Computes average along indicated axis. + If axis is None, average over the entire array + Inputs can be integer or floating types; result is of type float. + + If weights are given, result is sum(a*weights,axis=0)/(sum(weights,axis=0)*1.0) + weights must have a's shape or be the 1-d with length the size + of a in the given axis. + + If returned, return a tuple: the result and the sum of the weights + or count of values. Results will have the same shape. + + masked values in the weights will be set to 0.0 + """ + a = masked_array(a) + mask = a.mask + ash = a.shape + if ash == (): + ash = (1,) + if axis is None: + if mask is nomask: + if weights is None: + n = add.reduce(a.raw_data().ravel()) + d = reduce(lambda x, y: x * y, ash, 1.0) + else: + w = filled(weights, 0.0).ravel() + n = umath.add.reduce(a.raw_data().ravel() * w) + d = umath.add.reduce(w) + del w + else: + if weights is None: + n = add.reduce(a.ravel()) + w = fromnumeric.choose(mask, (1.0, 0.0)).ravel() + d = umath.add.reduce(w) + del w + else: + w = array(filled(weights, 0.0), float, mask=mask).ravel() + n = add.reduce(a.ravel() * w) + d = add.reduce(w) + del w + else: + if mask is nomask: + if weights is None: + d = ash[axis] * 1.0 + n = umath.add.reduce(a.raw_data(), axis) + else: + w = filled(weights, 0.0) + wsh = w.shape + if wsh == (): + wsh = (1,) + if wsh == ash: + w = numeric.array(w, float, copy=0) + n = add.reduce(a*w, axis) + d = add.reduce(w, axis) + del w + elif wsh == (ash[axis],): + r = [newaxis]*len(ash) + r[axis] = slice(None, None, 1) + w = eval ("w["+ repr(tuple(r)) + "] * ones(ash, float)") + n = add.reduce(a*w, axis) + d = add.reduce(w, axis) + del w, r + else: + raise ValueError, 'average: weights wrong shape.' + else: + if weights is None: + n = add.reduce(a, axis) + w = numeric.choose(mask, (1.0, 0.0)) + d = umath.add.reduce(w, axis) + del w + else: + w = filled(weights, 0.0) + wsh = w.shape + if wsh == (): + wsh = (1,) + if wsh == ash: + w = array(w, float, mask=mask, copy=0) + n = add.reduce(a*w, axis) + d = add.reduce(w, axis) + elif wsh == (ash[axis],): + r = [newaxis]*len(ash) + r[axis] = slice(None, None, 1) + w = eval ("w["+ repr(tuple(r)) + "] * masked_array(ones(ash, float), mask)") + n = add.reduce(a*w, axis) + d = add.reduce(w, axis) + else: + raise ValueError, 'average: weights wrong shape.' + del w + #print n, d, repr(mask), repr(weights) + if n is masked or d is masked: return masked + result = divide (n, d) + del n + + if isinstance(result, MaskedArray): + result.unmask() + if returned: + if not isinstance(d, MaskedArray): + d = masked_array(d) + if not d.shape == result.shape: + d = ones(result.shape, float) * d + d.unmask() + if returned: + return result, d + else: + return result + +def where (condition, x, y): + """where(condition, x, y) is x where condition is nonzero, y otherwise. + condition must be convertible to an integer array. + Answer is always the shape of condition. + The type depends on x and y. It is integer if both x and y are + the value masked. + """ + fc = filled(not_equal(condition, 0), 0) + xv = filled(x) + xm = getmask(x) + yv = filled(y) + ym = getmask(y) + d = numeric.choose(fc, (yv, xv)) + md = numeric.choose(fc, (ym, xm)) + m = getmask(condition) + m = make_mask(mask_or(m, md), copy=0, flag=1) + return masked_array(d, m) + +def choose (indices, t, out=None, mode='raise'): + "Returns array shaped like indices with elements chosen from t" + def fmask (x): + if x is masked: return 1 + return filled(x) + def nmask (x): + if x is masked: return 1 + m = getmask(x) + if m is nomask: return 0 + return m + c = filled(indices, 0) + masks = [nmask(x) for x in t] + a = [fmask(x) for x in t] + d = numeric.choose(c, a) + m = numeric.choose(c, masks) + m = make_mask(mask_or(m, getmask(indices)), copy=0, flag=1) + return masked_array(d, m) + +def masked_where(condition, x, copy=1): + """Return x as an array masked where condition is true. + Also masked where x or condition masked. + """ + cm = filled(condition,1) + m = mask_or(getmask(x), cm) + return array(filled(x), copy=copy, mask=m) + +def masked_greater(x, value, copy=1): + "masked_greater(x, value) = x masked where x > value" + return masked_where(greater(x, value), x, copy) + +def masked_greater_equal(x, value, copy=1): + "masked_greater_equal(x, value) = x masked where x >= value" + return masked_where(greater_equal(x, value), x, copy) + +def masked_less(x, value, copy=1): + "masked_less(x, value) = x masked where x < value" + return masked_where(less(x, value), x, copy) + +def masked_less_equal(x, value, copy=1): + "masked_less_equal(x, value) = x masked where x <= value" + return masked_where(less_equal(x, value), x, copy) + +def masked_not_equal(x, value, copy=1): + "masked_not_equal(x, value) = x masked where x != value" + d = filled(x, 0) + c = umath.not_equal(d, value) + m = mask_or(c, getmask(x)) + return array(d, mask=m, copy=copy) + +def masked_equal(x, value, copy=1): + """masked_equal(x, value) = x masked where x == value + For floating point consider masked_values(x, value) instead. + """ + d = filled(x, 0) + c = umath.equal(d, value) + m = mask_or(c, getmask(x)) + return array(d, mask=m, copy=copy) + +def masked_inside(x, v1, v2, copy=1): + """x with mask of all values of x that are inside [v1,v2] + v1 and v2 can be given in either order. + """ + if v2 < v1: + t = v2 + v2 = v1 + v1 = t + d = filled(x, 0) + c = umath.logical_and(umath.less_equal(d, v2), umath.greater_equal(d, v1)) + m = mask_or(c, getmask(x)) + return array(d, mask = m, copy=copy) + +def masked_outside(x, v1, v2, copy=1): + """x with mask of all values of x that are outside [v1,v2] + v1 and v2 can be given in either order. + """ + if v2 < v1: + t = v2 + v2 = v1 + v1 = t + d = filled(x, 0) + c = umath.logical_or(umath.less(d, v1), umath.greater(d, v2)) + m = mask_or(c, getmask(x)) + return array(d, mask = m, copy=copy) + +def reshape (a, *newshape): + "Copy of a with a new shape." + m = getmask(a) + d = filled(a).reshape(*newshape) + if m is nomask: + return masked_array(d) + else: + return masked_array(d, mask=numeric.reshape(m, *newshape)) + +def ravel (a): + "a as one-dimensional, may share data and mask" + m = getmask(a) + d = fromnumeric.ravel(filled(a)) + if m is nomask: + return masked_array(d) + else: + return masked_array(d, mask=numeric.ravel(m)) + +def concatenate (arrays, axis=0): + "Concatenate the arrays along the given axis" + d = [] + for x in arrays: + d.append(filled(x)) + d = numeric.concatenate(d, axis) + for x in arrays: + if getmask(x) is not nomask: break + else: + return masked_array(d) + dm = [] + for x in arrays: + dm.append(getmaskarray(x)) + dm = numeric.concatenate(dm, axis) + return masked_array(d, mask=dm) + +def swapaxes (a, axis1, axis2): + m = getmask(a) + d = masked_array(a).data + if m is nomask: + return masked_array(data=numeric.swapaxes(d, axis1, axis2)) + else: + return masked_array(data=numeric.swapaxes(d, axis1, axis2), + mask=numeric.swapaxes(m, axis1, axis2),) + + +def new_take (a, indices, axis=None, out=None, mode='raise'): + "returns selection of items from a." + m = getmask(a) + # d = masked_array(a).raw_data() + d = masked_array(a).data + if m is nomask: + return masked_array(numeric.take(d, indices, axis)) + else: + return masked_array(numeric.take(d, indices, axis), + mask = numeric.take(m, indices, axis)) + +def transpose(a, axes=None): + "reorder dimensions per tuple axes" + m = getmask(a) + d = filled(a) + if m is nomask: + return masked_array(numeric.transpose(d, axes)) + else: + return masked_array(numeric.transpose(d, axes), + mask = numeric.transpose(m, axes)) + + +def put(a, indices, values, mode='raise'): + """sets storage-indexed locations to corresponding values. + + Values and indices are filled if necessary. + + """ + d = a.raw_data() + ind = filled(indices) + v = filled(values) + numeric.put (d, ind, v) + m = getmask(a) + if m is not nomask: + a.unshare_mask() + numeric.put(a.raw_mask(), ind, 0) + +def putmask(a, mask, values): + "putmask(a, mask, values) sets a where mask is true." + if mask is nomask: + return + numeric.putmask(a.raw_data(), mask, values) + m = getmask(a) + if m is nomask: return + a.unshare_mask() + numeric.putmask(a.raw_mask(), mask, 0) + +def inner(a, b): + """inner(a,b) returns the dot product of two arrays, which has + shape a.shape[:-1] + b.shape[:-1] with elements computed by summing the + product of the elements from the last dimensions of a and b. + Masked elements are replace by zeros. + """ + fa = filled(a, 0) + fb = filled(b, 0) + if len(fa.shape) == 0: fa.shape = (1,) + if len(fb.shape) == 0: fb.shape = (1,) + return masked_array(numeric.inner(fa, fb)) + +innerproduct = inner + +def outer(a, b): + """outer(a,b) = {a[i]*b[j]}, has shape (len(a),len(b))""" + fa = filled(a, 0).ravel() + fb = filled(b, 0).ravel() + d = numeric.outer(fa, fb) + ma = getmask(a) + mb = getmask(b) + if ma is nomask and mb is nomask: + return masked_array(d) + ma = getmaskarray(a) + mb = getmaskarray(b) + m = make_mask(1-numeric.outer(1-ma, 1-mb), copy=0) + return masked_array(d, m) + +outerproduct = outer + +def dot(a, b): + """dot(a,b) returns matrix-multiplication between a and b. The product-sum + is over the last dimension of a and the second-to-last dimension of b. + Masked values are replaced by zeros. See also innerproduct. + """ + return innerproduct(filled(a, 0), numeric.swapaxes(filled(b, 0), -1, -2)) + +def compress(condition, x, dimension=-1, out=None): + """Select those parts of x for which condition is true. + Masked values in condition are considered false. + """ + c = filled(condition, 0) + m = getmask(x) + if m is not nomask: + m = numeric.compress(c, m, dimension) + d = numeric.compress(c, filled(x), dimension) + return masked_array(d, m) + +class _minimum_operation: + "Object to calculate minima" + def __init__ (self): + """minimum(a, b) or minimum(a) + In one argument case returns the scalar minimum. + """ + pass + + def __call__ (self, a, b=None): + "Execute the call behavior." + if b is None: + m = getmask(a) + if m is nomask: + d = amin(filled(a).ravel()) + return d + ac = a.compressed() + if len(ac) == 0: + return masked + else: + return amin(ac.raw_data()) + else: + return where(less(a, b), a, b) + + def reduce (self, target, axis=0): + """Reduce target along the given axis.""" + m = getmask(target) + if m is nomask: + t = filled(target) + return masked_array (umath.minimum.reduce (t, axis)) + else: + t = umath.minimum.reduce(filled(target, minimum_fill_value(target)), axis) + m = umath.logical_and.reduce(m, axis) + return masked_array(t, m, get_fill_value(target)) + + def outer (self, a, b): + "Return the function applied to the outer product of a and b." + ma = getmask(a) + mb = getmask(b) + if ma is nomask and mb is nomask: + m = nomask + else: + ma = getmaskarray(a) + mb = getmaskarray(b) + m = logical_or.outer(ma, mb) + d = umath.minimum.outer(filled(a), filled(b)) + return masked_array(d, m) + +minimum = _minimum_operation () + +class _maximum_operation: + "Object to calculate maxima" + def __init__ (self): + """maximum(a, b) or maximum(a) + In one argument case returns the scalar maximum. + """ + pass + + def __call__ (self, a, b=None): + "Execute the call behavior." + if b is None: + m = getmask(a) + if m is nomask: + d = amax(filled(a).ravel()) + return d + ac = a.compressed() + if len(ac) == 0: + return masked + else: + return amax(ac.raw_data()) + else: + return where(greater(a, b), a, b) + + def reduce (self, target, axis=0): + """Reduce target along the given axis.""" + m = getmask(target) + if m is nomask: + t = filled(target) + return masked_array (umath.maximum.reduce (t, axis)) + else: + t = umath.maximum.reduce(filled(target, maximum_fill_value(target)), axis) + m = umath.logical_and.reduce(m, axis) + return masked_array(t, m, get_fill_value(target)) + + def outer (self, a, b): + "Return the function applied to the outer product of a and b." + ma = getmask(a) + mb = getmask(b) + if ma is nomask and mb is nomask: + m = nomask + else: + ma = getmaskarray(a) + mb = getmaskarray(b) + m = logical_or.outer(ma, mb) + d = umath.maximum.outer(filled(a), filled(b)) + return masked_array(d, m) + +maximum = _maximum_operation () + +def sort (x, axis = -1, fill_value=None): + """If x does not have a mask, return a masked array formed from the + result of numeric.sort(x, axis). + Otherwise, fill x with fill_value. Sort it. + Set a mask where the result is equal to fill_value. + Note that this may have unintended consequences if the data contains the + fill value at a non-masked site. + + If fill_value is not given the default fill value for x's type will be + used. + """ + if fill_value is None: + fill_value = default_fill_value (x) + d = filled(x, fill_value) + s = fromnumeric.sort(d, axis) + if getmask(x) is nomask: + return masked_array(s) + return masked_values(s, fill_value, copy=0) + +def diagonal(a, k = 0, axis1=0, axis2=1): + """diagonal(a,k=0,axis1=0, axis2=1) = the k'th diagonal of a""" + d = fromnumeric.diagonal(filled(a), k, axis1, axis2) + m = getmask(a) + if m is nomask: + return masked_array(d, m) + else: + return masked_array(d, fromnumeric.diagonal(m, k, axis1, axis2)) + +def trace (a, offset=0, axis1=0, axis2=1, dtype=None, out=None): + """trace(a,offset=0, axis1=0, axis2=1) returns the sum along diagonals + (defined by the last two dimenions) of the array. + """ + return diagonal(a, offset, axis1, axis2).sum(dtype=dtype) + +def argsort (x, axis = -1, out=None, fill_value=None): + """Treating masked values as if they have the value fill_value, + return sort indices for sorting along given axis. + if fill_value is None, use get_fill_value(x) + Returns a numpy array. + """ + d = filled(x, fill_value) + return fromnumeric.argsort(d, axis) + +def argmin (x, axis = -1, out=None, fill_value=None): + """Treating masked values as if they have the value fill_value, + return indices for minimum values along given axis. + if fill_value is None, use get_fill_value(x). + Returns a numpy array if x has more than one dimension. + Otherwise, returns a scalar index. + """ + d = filled(x, fill_value) + return fromnumeric.argmin(d, axis) + +def argmax (x, axis = -1, out=None, fill_value=None): + """Treating masked values as if they have the value fill_value, + return sort indices for maximum along given axis. + if fill_value is None, use -get_fill_value(x) if it exists. + Returns a numpy array if x has more than one dimension. + Otherwise, returns a scalar index. + """ + if fill_value is None: + fill_value = default_fill_value (x) + try: + fill_value = - fill_value + except: + pass + d = filled(x, fill_value) + return fromnumeric.argmax(d, axis) + +def fromfunction (f, s): + """apply f to s to create array as in umath.""" + return masked_array(numeric.fromfunction(f, s)) + +def asarray(data, dtype=None): + """asarray(data, dtype) = array(data, dtype, copy=0) + """ + if isinstance(data, MaskedArray) and \ + (dtype is None or dtype == data.dtype): + return data + return array(data, dtype=dtype, copy=0) + +# Add methods to support ndarray interface +# XXX: I is better to to change the masked_*_operation adaptors +# XXX: to wrap ndarray methods directly to create ma.array methods. +from types import MethodType +def _m(f): + return MethodType(f, None, array) +def not_implemented(*args, **kwds): + raise NotImplementedError, "not yet implemented for numpy.ma arrays" +array.all = _m(alltrue) +array.any = _m(sometrue) +array.argmax = _m(argmax) +array.argmin = _m(argmin) +array.argsort = _m(argsort) +array.base = property(_m(not_implemented)) +array.byteswap = _m(not_implemented) + +def _choose(self, *args, **kwds): + return choose(self, args) +array.choose = _m(_choose) +del _choose + +def _clip(self,a_min,a_max,out=None): + return MaskedArray(data = self.data.clip(asarray(a_min).data, + asarray(a_max).data), + mask = mask_or(self.mask, + mask_or(getmask(a_min),getmask(a_max)))) +array.clip = _m(_clip) + +def _compress(self, cond, axis=None, out=None): + return compress(cond, self, axis) +array.compress = _m(_compress) +del _compress + +array.conj = array.conjugate = _m(conjugate) +array.copy = _m(not_implemented) + +def _cumprod(self, axis=None, dtype=None, out=None): + m = self.mask + if m is not nomask: + m = umath.logical_or.accumulate(self.mask, axis) + return MaskedArray(data = self.filled(1).cumprod(axis, dtype), mask=m) +array.cumprod = _m(_cumprod) + +def _cumsum(self, axis=None, dtype=None, out=None): + m = self.mask + if m is not nomask: + m = umath.logical_or.accumulate(self.mask, axis) + return MaskedArray(data=self.filled(0).cumsum(axis, dtype), mask=m) +array.cumsum = _m(_cumsum) + +array.diagonal = _m(diagonal) +array.dump = _m(not_implemented) +array.dumps = _m(not_implemented) +array.fill = _m(not_implemented) +array.flags = property(_m(not_implemented)) +array.flatten = _m(ravel) +array.getfield = _m(not_implemented) + +def _max(a, axis=None, out=None): + if out is not None: + raise TypeError("Output arrays Unsupported for masked arrays") + if axis is None: + return maximum(a) + else: + return maximum.reduce(a, axis) +array.max = _m(_max) +del _max +def _min(a, axis=None, out=None): + if out is not None: + raise TypeError("Output arrays Unsupported for masked arrays") + if axis is None: + return minimum(a) + else: + return minimum.reduce(a, axis) +array.min = _m(_min) +del _min +array.mean = _m(new_average) +array.nbytes = property(_m(not_implemented)) +array.newbyteorder = _m(not_implemented) +array.nonzero = _m(nonzero) +array.prod = _m(product) + +def _ptp(a,axis=None,out=None): + return a.max(axis,out)-a.min(axis) +array.ptp = _m(_ptp) +array.repeat = _m(new_repeat) +array.resize = _m(resize) +array.searchsorted = _m(not_implemented) +array.setfield = _m(not_implemented) +array.setflags = _m(not_implemented) +array.sort = _m(not_implemented) # NB: ndarray.sort is inplace + +def _squeeze(self): + try: + result = MaskedArray(data = self.data.squeeze(), + mask = self.mask.squeeze()) + except AttributeError: + result = _wrapit(self, 'squeeze') + return result +array.squeeze = _m(_squeeze) + +array.strides = property(_m(not_implemented)) +array.sum = _m(sum) +def _swapaxes(self,axis1,axis2): + return MaskedArray(data = self.data.swapaxes(axis1, axis2), + mask = self.mask.swapaxes(axis1, axis2)) +array.swapaxes = _m(_swapaxes) +array.take = _m(new_take) +array.tofile = _m(not_implemented) +array.trace = _m(trace) +array.transpose = _m(transpose) + +def _var(self,axis=None,dtype=None, out=None): + if axis is None: + return numeric.asarray(self.compressed()).var() + a = self.swapaxes(axis,0) + a = a - a.mean(axis=0) + a *= a + a /= a.count(axis=0) + return a.swapaxes(0,axis).sum(axis) +def _std(self,axis=None, dtype=None, out=None): + return (self.var(axis,dtype))**0.5 +array.var = _m(_var) +array.std = _m(_std) + +array.view = _m(not_implemented) +array.round = _m(around) +del _m, MethodType, not_implemented + + +masked = MaskedArray(0, int, mask=1) + +def repeat(a, repeats, axis=0): + return new_repeat(a, repeats, axis) + +def average(a, axis=0, weights=None, returned=0): + return new_average(a, axis, weights, returned) + +def take(a, indices, axis=0): + return new_take(a, indices, axis) diff --git a/numpy/oldnumeric/matrix.py b/numpy/oldnumeric/matrix.py new file mode 100644 index 000000000..5f8c1ca5e --- /dev/null +++ b/numpy/oldnumeric/matrix.py @@ -0,0 +1,67 @@ +# This module is for compatibility only. + +__all__ = ['UserArray', 'squeeze', 'Matrix', 'asarray', 'dot', 'k', 'Numeric', 'LinearAlgebra', 'identity', 'multiply', 'types', 'string'] + +import types +from user_array import UserArray, asarray +import numpy.oldnumeric as Numeric +from numpy.oldnumeric import dot, identity, multiply +import numpy.oldnumeric.linear_algebra as LinearAlgebra +from numpy import matrix as Matrix, squeeze + +# Hidden names that will be the same. + +_table = [None]*256 +for k in range(256): + _table[k] = chr(k) +_table = ''.join(_table) + +_numchars = '0123456789.-+jeEL' +_todelete = [] +for k in _table: + if k not in _numchars: + _todelete.append(k) +_todelete = ''.join(_todelete) + + +def _eval(astr): + return eval(astr.translate(_table,_todelete)) + +def _convert_from_string(data): + data.find + rows = data.split(';') + newdata = [] + count = 0 + for row in rows: + trow = row.split(',') + newrow = [] + for col in trow: + temp = col.split() + newrow.extend(map(_eval,temp)) + if count == 0: + Ncols = len(newrow) + elif len(newrow) != Ncols: + raise ValueError, "Rows not the same size." + count += 1 + newdata.append(newrow) + return newdata + + +_lkup = {'0':'000', + '1':'001', + '2':'010', + '3':'011', + '4':'100', + '5':'101', + '6':'110', + '7':'111'} + +def _binary(num): + ostr = oct(num) + bin = '' + for ch in ostr[1:]: + bin += _lkup[ch] + ind = 0 + while bin[ind] == '0': + ind += 1 + return bin[ind:] diff --git a/numpy/oldnumeric/misc.py b/numpy/oldnumeric/misc.py new file mode 100644 index 000000000..d6e0f2430 --- /dev/null +++ b/numpy/oldnumeric/misc.py @@ -0,0 +1,29 @@ +# Functions that already have the correct syntax or miscellaneous functions + + +__all__ = ['sort', 'copy_reg', 'clip', 'rank', + 'sign', 'shape', 'types', 'allclose', 'size', + 'choose', 'swapaxes', 'array_str', + 'pi', 'math', 'concatenate', 'putmask', 'put', + 'around', 'vdot', 'transpose', 'array2string', 'diagonal', + 'searchsorted', 'copy', 'resize', + 'array_repr', 'e', 'StringIO', 'pickle', + 'argsort', 'convolve', 'cross_correlate', + 'dot', 'outerproduct', 'innerproduct', 'insert'] + +import types +import StringIO +import pickle +import math +import copy +import copy_reg + +from numpy import sort, clip, rank, sign, shape, putmask, allclose, size,\ + choose, swapaxes, array_str, array_repr, e, pi, put, \ + resize, around, concatenate, vdot, transpose, \ + diagonal, searchsorted, argsort, convolve, dot, \ + outer as outerproduct, inner as innerproduct, \ + correlate as cross_correlate, \ + place as insert + +from array_printer import array2string diff --git a/numpy/oldnumeric/mlab.py b/numpy/oldnumeric/mlab.py new file mode 100644 index 000000000..c11e34c1f --- /dev/null +++ b/numpy/oldnumeric/mlab.py @@ -0,0 +1,126 @@ +# This module is for compatibility only. All functions are defined elsewhere. + +__all__ = ['rand', 'tril', 'trapz', 'hanning', 'rot90', 'triu', 'diff', 'angle', + 'roots', 'ptp', 'kaiser', 'randn', 'cumprod', 'diag', 'msort', + 'LinearAlgebra', 'RandomArray', 'prod', 'std', 'hamming', 'flipud', + 'max', 'blackman', 'corrcoef', 'bartlett', 'eye', 'squeeze', 'sinc', + 'tri', 'cov', 'svd', 'min', 'median', 'fliplr', 'eig', 'mean'] + +import numpy.oldnumeric.linear_algebra as LinearAlgebra +import numpy.oldnumeric.random_array as RandomArray +from numpy import tril, trapz as _Ntrapz, hanning, rot90, triu, diff, \ + angle, roots, ptp as _Nptp, kaiser, cumprod as _Ncumprod, \ + diag, msort, prod as _Nprod, std as _Nstd, hamming, flipud, \ + amax as _Nmax, amin as _Nmin, blackman, bartlett, \ + squeeze, sinc, median, fliplr, mean as _Nmean, transpose + +from numpy.linalg import eig, svd +from numpy.random import rand, randn +import numpy as np + +from typeconv import convtypecode + +def eye(N, M=None, k=0, typecode=None, dtype=None): + """ eye returns a N-by-M 2-d array where the k-th diagonal is all ones, + and everything else is zeros. + """ + dtype = convtypecode(typecode, dtype) + if M is None: M = N + m = np.equal(np.subtract.outer(np.arange(N), np.arange(M)),-k) + if m.dtype != dtype: + return m.astype(dtype) + +def tri(N, M=None, k=0, typecode=None, dtype=None): + """ returns a N-by-M array where all the diagonals starting from + lower left corner up to the k-th are all ones. + """ + dtype = convtypecode(typecode, dtype) + if M is None: M = N + m = np.greater_equal(np.subtract.outer(np.arange(N), np.arange(M)),-k) + if m.dtype != dtype: + return m.astype(dtype) + +def trapz(y, x=None, axis=-1): + return _Ntrapz(y, x, axis=axis) + +def ptp(x, axis=0): + return _Nptp(x, axis) + +def cumprod(x, axis=0): + return _Ncumprod(x, axis) + +def max(x, axis=0): + return _Nmax(x, axis) + +def min(x, axis=0): + return _Nmin(x, axis) + +def prod(x, axis=0): + return _Nprod(x, axis) + +def std(x, axis=0): + N = asarray(x).shape[axis] + return _Nstd(x, axis)*sqrt(N/(N-1.)) + +def mean(x, axis=0): + return _Nmean(x, axis) + +# This is exactly the same cov function as in MLab +def cov(m, y=None, rowvar=0, bias=0): + if y is None: + y = m + else: + y = y + if rowvar: + m = transpose(m) + y = transpose(y) + if (m.shape[0] == 1): + m = transpose(m) + if (y.shape[0] == 1): + y = transpose(y) + N = m.shape[0] + if (y.shape[0] != N): + raise ValueError, "x and y must have the same number "\ + "of observations" + m = m - _Nmean(m,axis=0) + y = y - _Nmean(y,axis=0) + if bias: + fact = N*1.0 + else: + fact = N-1.0 + return squeeze(dot(transpose(m), conjugate(y)) / fact) + +from numpy import sqrt, multiply +def corrcoef(x, y=None): + c = cov(x,y) + d = diag(c) + return c/sqrt(multiply.outer(d,d)) + +from compat import * +from functions import * +from precision import * +from ufuncs import * +from misc import * + +import compat +import precision +import functions +import misc +import ufuncs + +import numpy +__version__ = numpy.__version__ +del numpy + +__all__ += ['__version__'] +__all__ += compat.__all__ +__all__ += precision.__all__ +__all__ += functions.__all__ +__all__ += ufuncs.__all__ +__all__ += misc.__all__ + +del compat +del functions +del precision +del ufuncs +del misc diff --git a/numpy/oldnumeric/precision.py b/numpy/oldnumeric/precision.py new file mode 100644 index 000000000..c095ceb19 --- /dev/null +++ b/numpy/oldnumeric/precision.py @@ -0,0 +1,168 @@ +# Lifted from Precision.py. This is for compatibility only. +# +# The character strings are still for "new" NumPy +# which is the only Incompatibility with Numeric + +__all__ = ['Character', 'Complex', 'Float', + 'PrecisionError', 'PyObject', 'Int', 'UInt', + 'UnsignedInt', 'UnsignedInteger', 'string', 'typecodes', 'zeros'] + +from functions import zeros +import string # for backwards compatibility + +typecodes = {'Character':'c', 'Integer':'bhil', 'UnsignedInteger':'BHIL', 'Float':'fd', 'Complex':'FD'} + +def _get_precisions(typecodes): + lst = [] + for t in typecodes: + lst.append( (zeros( (1,), t ).itemsize*8, t) ) + return lst + +def _fill_table(typecodes, table={}): + for key, value in typecodes.items(): + table[key] = _get_precisions(value) + return table + +_code_table = _fill_table(typecodes) + +class PrecisionError(Exception): + pass + +def _lookup(table, key, required_bits): + lst = table[key] + for bits, typecode in lst: + if bits >= required_bits: + return typecode + raise PrecisionError, key+" of "+str(required_bits)+" bits not available on this system" + +Character = 'c' + +try: + UnsignedInt8 = _lookup(_code_table, "UnsignedInteger", 8) + UInt8 = UnsignedInt8 + __all__.extend(['UnsignedInt8', 'UInt8']) +except(PrecisionError): + pass +try: + UnsignedInt16 = _lookup(_code_table, "UnsignedInteger", 16) + UInt16 = UnsignedInt16 + __all__.extend(['UnsignedInt16', 'UInt16']) +except(PrecisionError): + pass +try: + UnsignedInt32 = _lookup(_code_table, "UnsignedInteger", 32) + UInt32 = UnsignedInt32 + __all__.extend(['UnsignedInt32', 'UInt32']) +except(PrecisionError): + pass +try: + UnsignedInt64 = _lookup(_code_table, "UnsignedInteger", 64) + UInt64 = UnsignedInt64 + __all__.extend(['UnsignedInt64', 'UInt64']) +except(PrecisionError): + pass +try: + UnsignedInt128 = _lookup(_code_table, "UnsignedInteger", 128) + UInt128 = UnsignedInt128 + __all__.extend(['UnsignedInt128', 'UInt128']) +except(PrecisionError): + pass +UInt = UnsignedInt = UnsignedInteger = 'u' + +try: + Int0 = _lookup(_code_table, 'Integer', 0) + __all__.append('Int0') +except(PrecisionError): + pass +try: + Int8 = _lookup(_code_table, 'Integer', 8) + __all__.append('Int8') +except(PrecisionError): + pass +try: + Int16 = _lookup(_code_table, 'Integer', 16) + __all__.append('Int16') +except(PrecisionError): + pass +try: + Int32 = _lookup(_code_table, 'Integer', 32) + __all__.append('Int32') +except(PrecisionError): + pass +try: + Int64 = _lookup(_code_table, 'Integer', 64) + __all__.append('Int64') +except(PrecisionError): + pass +try: + Int128 = _lookup(_code_table, 'Integer', 128) + __all__.append('Int128') +except(PrecisionError): + pass +Int = 'l' + +try: + Float0 = _lookup(_code_table, 'Float', 0) + __all__.append('Float0') +except(PrecisionError): + pass +try: + Float8 = _lookup(_code_table, 'Float', 8) + __all__.append('Float8') +except(PrecisionError): + pass +try: + Float16 = _lookup(_code_table, 'Float', 16) + __all__.append('Float16') +except(PrecisionError): + pass +try: + Float32 = _lookup(_code_table, 'Float', 32) + __all__.append('Float32') +except(PrecisionError): + pass +try: + Float64 = _lookup(_code_table, 'Float', 64) + __all__.append('Float64') +except(PrecisionError): + pass +try: + Float128 = _lookup(_code_table, 'Float', 128) + __all__.append('Float128') +except(PrecisionError): + pass +Float = 'd' + +try: + Complex0 = _lookup(_code_table, 'Complex', 0) + __all__.append('Complex0') +except(PrecisionError): + pass +try: + Complex8 = _lookup(_code_table, 'Complex', 16) + __all__.append('Complex8') +except(PrecisionError): + pass +try: + Complex16 = _lookup(_code_table, 'Complex', 32) + __all__.append('Complex16') +except(PrecisionError): + pass +try: + Complex32 = _lookup(_code_table, 'Complex', 64) + __all__.append('Complex32') +except(PrecisionError): + pass +try: + Complex64 = _lookup(_code_table, 'Complex', 128) + __all__.append('Complex64') +except(PrecisionError): + pass +try: + Complex128 = _lookup(_code_table, 'Complex', 256) + __all__.append('Complex128') +except(PrecisionError): + pass +Complex = 'D' + +PyObject = 'O' diff --git a/numpy/oldnumeric/random_array.py b/numpy/oldnumeric/random_array.py new file mode 100644 index 000000000..e84aedf1e --- /dev/null +++ b/numpy/oldnumeric/random_array.py @@ -0,0 +1,266 @@ +# Backward compatible module for RandomArray + +__all__ = ['ArgumentError','F','beta','binomial','chi_square', 'exponential', + 'gamma', 'get_seed', 'mean_var_test', 'multinomial', + 'multivariate_normal', 'negative_binomial', 'noncentral_F', + 'noncentral_chi_square', 'normal', 'permutation', 'poisson', + 'randint', 'random', 'random_integers', 'seed', 'standard_normal', + 'uniform'] + +ArgumentError = ValueError + +import numpy.random.mtrand as mt +import numpy as np + +def seed(x=0, y=0): + if (x == 0 or y == 0): + mt.seed() + else: + mt.seed((x,y)) + +def get_seed(): + raise NotImplementedError, \ + "If you want to save the state of the random number generator.\n"\ + "Then you should use obj = numpy.random.get_state() followed by.\n"\ + "numpy.random.set_state(obj)." + +def random(shape=[]): + "random(n) or random([n, m, ...]) returns array of random numbers" + if shape == []: + shape = None + return mt.random_sample(shape) + +def uniform(minimum, maximum, shape=[]): + """uniform(minimum, maximum, shape=[]) returns array of given shape of random reals + in given range""" + if shape == []: + shape = None + return mt.uniform(minimum, maximum, shape) + +def randint(minimum, maximum=None, shape=[]): + """randint(min, max, shape=[]) = random integers >=min, < max + If max not given, random integers >= 0, <min""" + if not isinstance(minimum, int): + raise ArgumentError, "randint requires first argument integer" + if maximum is None: + maximum = minimum + minimum = 0 + if not isinstance(maximum, int): + raise ArgumentError, "randint requires second argument integer" + a = ((maximum-minimum)* random(shape)) + if isinstance(a, np.ndarray): + return minimum + a.astype(np.int) + else: + return minimum + int(a) + +def random_integers(maximum, minimum=1, shape=[]): + """random_integers(max, min=1, shape=[]) = random integers in range min-max inclusive""" + return randint(minimum, maximum+1, shape) + +def permutation(n): + "permutation(n) = a permutation of indices range(n)" + return mt.permutation(n) + +def standard_normal(shape=[]): + """standard_normal(n) or standard_normal([n, m, ...]) returns array of + random numbers normally distributed with mean 0 and standard + deviation 1""" + if shape == []: + shape = None + return mt.standard_normal(shape) + +def normal(mean, std, shape=[]): + """normal(mean, std, n) or normal(mean, std, [n, m, ...]) returns + array of random numbers randomly distributed with specified mean and + standard deviation""" + if shape == []: + shape = None + return mt.normal(mean, std, shape) + +def multivariate_normal(mean, cov, shape=[]): + """multivariate_normal(mean, cov) or multivariate_normal(mean, cov, [m, n, ...]) + returns an array containing multivariate normally distributed random numbers + with specified mean and covariance. + + mean must be a 1 dimensional array. cov must be a square two dimensional + array with the same number of rows and columns as mean has elements. + + The first form returns a single 1-D array containing a multivariate + normal. + + The second form returns an array of shape (m, n, ..., cov.shape[0]). + In this case, output[i,j,...,:] is a 1-D array containing a multivariate + normal.""" + if shape == []: + shape = None + return mt.multivariate_normal(mean, cov, shape) + +def exponential(mean, shape=[]): + """exponential(mean, n) or exponential(mean, [n, m, ...]) returns array + of random numbers exponentially distributed with specified mean""" + if shape == []: + shape = None + return mt.exponential(mean, shape) + +def beta(a, b, shape=[]): + """beta(a, b) or beta(a, b, [n, m, ...]) returns array of beta distributed random numbers.""" + if shape == []: + shape = None + return mt.beta(a, b, shape) + +def gamma(a, r, shape=[]): + """gamma(a, r) or gamma(a, r, [n, m, ...]) returns array of gamma distributed random numbers.""" + if shape == []: + shape = None + return mt.gamma(a, r, shape) + +def F(dfn, dfd, shape=[]): + """F(dfn, dfd) or F(dfn, dfd, [n, m, ...]) returns array of F distributed random numbers with dfn degrees of freedom in the numerator and dfd degrees of freedom in the denominator.""" + if shape == []: + shape = None + return mt.f(dfn, dfd, shape) + +def noncentral_F(dfn, dfd, nconc, shape=[]): + """noncentral_F(dfn, dfd, nonc) or noncentral_F(dfn, dfd, nonc, [n, m, ...]) returns array of noncentral F distributed random numbers with dfn degrees of freedom in the numerator and dfd degrees of freedom in the denominator, and noncentrality parameter nconc.""" + if shape == []: + shape = None + return mt.noncentral_f(dfn, dfd, nconc, shape) + +def chi_square(df, shape=[]): + """chi_square(df) or chi_square(df, [n, m, ...]) returns array of chi squared distributed random numbers with df degrees of freedom.""" + if shape == []: + shape = None + return mt.chisquare(df, shape) + +def noncentral_chi_square(df, nconc, shape=[]): + """noncentral_chi_square(df, nconc) or chi_square(df, nconc, [n, m, ...]) returns array of noncentral chi squared distributed random numbers with df degrees of freedom and noncentrality parameter.""" + if shape == []: + shape = None + return mt.noncentral_chisquare(df, nconc, shape) + +def binomial(trials, p, shape=[]): + """binomial(trials, p) or binomial(trials, p, [n, m, ...]) returns array of binomially distributed random integers. + + trials is the number of trials in the binomial distribution. + p is the probability of an event in each trial of the binomial distribution.""" + if shape == []: + shape = None + return mt.binomial(trials, p, shape) + +def negative_binomial(trials, p, shape=[]): + """negative_binomial(trials, p) or negative_binomial(trials, p, [n, m, ...]) returns + array of negative binomially distributed random integers. + + trials is the number of trials in the negative binomial distribution. + p is the probability of an event in each trial of the negative binomial distribution.""" + if shape == []: + shape = None + return mt.negative_binomial(trials, p, shape) + +def multinomial(trials, probs, shape=[]): + """multinomial(trials, probs) or multinomial(trials, probs, [n, m, ...]) returns + array of multinomial distributed integer vectors. + + trials is the number of trials in each multinomial distribution. + probs is a one dimensional array. There are len(prob)+1 events. + prob[i] is the probability of the i-th event, 0<=i<len(prob). + The probability of event len(prob) is 1.-np.sum(prob). + + The first form returns a single 1-D array containing one multinomially + distributed vector. + + The second form returns an array of shape (m, n, ..., len(probs)). + In this case, output[i,j,...,:] is a 1-D array containing a multinomially + distributed integer 1-D array.""" + if shape == []: + shape = None + return mt.multinomial(trials, probs, shape) + +def poisson(mean, shape=[]): + """poisson(mean) or poisson(mean, [n, m, ...]) returns array of poisson + distributed random integers with specified mean.""" + if shape == []: + shape = None + return mt.poisson(mean, shape) + + +def mean_var_test(x, type, mean, var, skew=[]): + n = len(x) * 1.0 + x_mean = np.sum(x,axis=0)/n + x_minus_mean = x - x_mean + x_var = np.sum(x_minus_mean*x_minus_mean,axis=0)/(n-1.0) + print "\nAverage of ", len(x), type + print "(should be about ", mean, "):", x_mean + print "Variance of those random numbers (should be about ", var, "):", x_var + if skew != []: + x_skew = (np.sum(x_minus_mean*x_minus_mean*x_minus_mean,axis=0)/9998.)/x_var**(3./2.) + print "Skewness of those random numbers (should be about ", skew, "):", x_skew + +def test(): + obj = mt.get_state() + mt.set_state(obj) + obj2 = mt.get_state() + if (obj2[1] - obj[1]).any(): + raise SystemExit, "Failed seed test." + print "First random number is", random() + print "Average of 10000 random numbers is", np.sum(random(10000),axis=0)/10000. + x = random([10,1000]) + if len(x.shape) != 2 or x.shape[0] != 10 or x.shape[1] != 1000: + raise SystemExit, "random returned wrong shape" + x.shape = (10000,) + print "Average of 100 by 100 random numbers is", np.sum(x,axis=0)/10000. + y = uniform(0.5,0.6, (1000,10)) + if len(y.shape) !=2 or y.shape[0] != 1000 or y.shape[1] != 10: + raise SystemExit, "uniform returned wrong shape" + y.shape = (10000,) + if np.minimum.reduce(y) <= 0.5 or np.maximum.reduce(y) >= 0.6: + raise SystemExit, "uniform returned out of desired range" + print "randint(1, 10, shape=[50])" + print randint(1, 10, shape=[50]) + print "permutation(10)", permutation(10) + print "randint(3,9)", randint(3,9) + print "random_integers(10, shape=[20])" + print random_integers(10, shape=[20]) + s = 3.0 + x = normal(2.0, s, [10, 1000]) + if len(x.shape) != 2 or x.shape[0] != 10 or x.shape[1] != 1000: + raise SystemExit, "standard_normal returned wrong shape" + x.shape = (10000,) + mean_var_test(x, "normally distributed numbers with mean 2 and variance %f"%(s**2,), 2, s**2, 0) + x = exponential(3, 10000) + mean_var_test(x, "random numbers exponentially distributed with mean %f"%(s,), s, s**2, 2) + x = multivariate_normal(np.array([10,20]), np.array(([1,2],[2,4]))) + print "\nA multivariate normal", x + if x.shape != (2,): raise SystemExit, "multivariate_normal returned wrong shape" + x = multivariate_normal(np.array([10,20]), np.array([[1,2],[2,4]]), [4,3]) + print "A 4x3x2 array containing multivariate normals" + print x + if x.shape != (4,3,2): raise SystemExit, "multivariate_normal returned wrong shape" + x = multivariate_normal(np.array([-100,0,100]), np.array([[3,2,1],[2,2,1],[1,1,1]]), 10000) + x_mean = np.sum(x,axis=0)/10000. + print "Average of 10000 multivariate normals with mean [-100,0,100]" + print x_mean + x_minus_mean = x - x_mean + print "Estimated covariance of 10000 multivariate normals with covariance [[3,2,1],[2,2,1],[1,1,1]]" + print np.dot(np.transpose(x_minus_mean),x_minus_mean)/9999. + x = beta(5.0, 10.0, 10000) + mean_var_test(x, "beta(5.,10.) random numbers", 0.333, 0.014) + x = gamma(.01, 2., 10000) + mean_var_test(x, "gamma(.01,2.) random numbers", 2*100, 2*100*100) + x = chi_square(11., 10000) + mean_var_test(x, "chi squared random numbers with 11 degrees of freedom", 11, 22, 2*np.sqrt(2./11.)) + x = F(5., 10., 10000) + mean_var_test(x, "F random numbers with 5 and 10 degrees of freedom", 1.25, 1.35) + x = poisson(50., 10000) + mean_var_test(x, "poisson random numbers with mean 50", 50, 50, 0.14) + print "\nEach element is the result of 16 binomial trials with probability 0.5:" + print binomial(16, 0.5, 16) + print "\nEach element is the result of 16 negative binomial trials with probability 0.5:" + print negative_binomial(16, 0.5, [16,]) + print "\nEach row is the result of 16 multinomial trials with probabilities [0.1, 0.5, 0.1 0.3]:" + x = multinomial(16, [0.1, 0.5, 0.1], 8) + print x + print "Mean = ", np.sum(x,axis=0)/8. + +if __name__ == '__main__': + test() diff --git a/numpy/oldnumeric/rng.py b/numpy/oldnumeric/rng.py new file mode 100644 index 000000000..b4c72a68c --- /dev/null +++ b/numpy/oldnumeric/rng.py @@ -0,0 +1,135 @@ +# This module re-creates the RNG interface from Numeric +# Replace import RNG with import numpy.oldnumeric.rng as RNG +# +# It is for backwards compatibility only. + + +__all__ = ['CreateGenerator','ExponentialDistribution','LogNormalDistribution', + 'NormalDistribution', 'UniformDistribution', 'error', 'ranf', + 'default_distribution', 'random_sample', 'standard_generator'] + +import numpy.random.mtrand as mt +import math + +class error(Exception): + pass + +class Distribution(object): + def __init__(self, meth, *args): + self._meth = meth + self._args = args + + def density(self,x): + raise NotImplementedError + + def __call__(self, x): + return self.density(x) + + def _onesample(self, rng): + return getattr(rng, self._meth)(*self._args) + + def _sample(self, rng, n): + kwds = {'size' : n} + return getattr(rng, self._meth)(*self._args, **kwds) + + +class ExponentialDistribution(Distribution): + def __init__(self, lambda_): + if (lambda_ <= 0): + raise error, "parameter must be positive" + Distribution.__init__(self, 'exponential', lambda_) + + def density(x): + if x < 0: + return 0.0 + else: + lambda_ = self._args[0] + return lambda_*math.exp(-lambda_*x) + +class LogNormalDistribution(Distribution): + def __init__(self, m, s): + m = float(m) + s = float(s) + if (s <= 0): + raise error, "standard deviation must be positive" + Distribution.__init__(self, 'lognormal', m, s) + sn = math.log(1.0+s*s/(m*m)); + self._mn = math.log(m)-0.5*sn + self._sn = math.sqrt(sn) + self._fac = 1.0/math.sqrt(2*math.pi)/self._sn + + def density(x): + m,s = self._args + y = (math.log(x)-self._mn)/self._sn + return self._fac*math.exp(-0.5*y*y)/x + + +class NormalDistribution(Distribution): + def __init__(self, m, s): + m = float(m) + s = float(s) + if (s <= 0): + raise error, "standard deviation must be positive" + Distribution.__init__(self, 'normal', m, s) + self._fac = 1.0/math.sqrt(2*math.pi)/s + + def density(x): + m,s = self._args + y = (x-m)/s + return self._fac*math.exp(-0.5*y*y) + +class UniformDistribution(Distribution): + def __init__(self, a, b): + a = float(a) + b = float(b) + width = b-a + if (width <=0): + raise error, "width of uniform distribution must be > 0" + Distribution.__init__(self, 'uniform', a, b) + self._fac = 1.0/width + + def density(x): + a, b = self._args + if (x < a) or (x >= b): + return 0.0 + else: + return self._fac + +default_distribution = UniformDistribution(0.0,1.0) + +class CreateGenerator(object): + def __init__(self, seed, dist=None): + if seed <= 0: + self._rng = mt.RandomState() + elif seed > 0: + self._rng = mt.RandomState(seed) + if dist is None: + dist = default_distribution + if not isinstance(dist, Distribution): + raise error, "Not a distribution object" + self._dist = dist + + def ranf(self): + return self._dist._onesample(self._rng) + + def sample(self, n): + return self._dist._sample(self._rng, n) + + +standard_generator = CreateGenerator(-1) + +def ranf(): + "ranf() = a random number from the standard generator." + return standard_generator.ranf() + +def random_sample(*n): + """random_sample(n) = array of n random numbers; + + random_sample(n1, n2, ...)= random array of shape (n1, n2, ..)""" + + if not n: + return standard_generator.ranf() + m = 1 + for i in n: + m = m * i + return standard_generator.sample(m).reshape(*n) diff --git a/numpy/oldnumeric/rng_stats.py b/numpy/oldnumeric/rng_stats.py new file mode 100644 index 000000000..8c7fec433 --- /dev/null +++ b/numpy/oldnumeric/rng_stats.py @@ -0,0 +1,35 @@ + +__all__ = ['average', 'histogram', 'standardDeviation', 'variance'] + +import numpy.oldnumeric as Numeric + +def average(data): + data = Numeric.array(data) + return Numeric.add.reduce(data)/len(data) + +def variance(data): + data = Numeric.array(data) + return Numeric.add.reduce((data-average(data,axis=0))**2)/(len(data)-1) + +def standardDeviation(data): + data = Numeric.array(data) + return Numeric.sqrt(variance(data)) + +def histogram(data, nbins, range = None): + data = Numeric.array(data, Numeric.Float) + if range is None: + min = Numeric.minimum.reduce(data) + max = Numeric.maximum.reduce(data) + else: + min, max = range + data = Numeric.repeat(data, + Numeric.logical_and(Numeric.less_equal(data, max), + Numeric.greater_equal(data, + min)),axis=0) + bin_width = (max-min)/nbins + data = Numeric.floor((data - min)/bin_width).astype(Numeric.Int) + histo = Numeric.add.reduce(Numeric.equal( + Numeric.arange(nbins)[:,Numeric.NewAxis], data), -1) + histo[-1] = histo[-1] + Numeric.add.reduce(Numeric.equal(nbins, data)) + bins = min + bin_width*(Numeric.arange(nbins)+0.5) + return Numeric.transpose(Numeric.array([bins, histo])) diff --git a/numpy/oldnumeric/setup.py b/numpy/oldnumeric/setup.py new file mode 100644 index 000000000..31b5ff3cc --- /dev/null +++ b/numpy/oldnumeric/setup.py @@ -0,0 +1,10 @@ + +def configuration(parent_package='',top_path=None): + from numpy.distutils.misc_util import Configuration + config = Configuration('oldnumeric',parent_package,top_path) + config.add_data_dir('tests') + return config + +if __name__ == '__main__': + from numpy.distutils.core import setup + setup(configuration=configuration) diff --git a/numpy/oldnumeric/setupscons.py b/numpy/oldnumeric/setupscons.py new file mode 100644 index 000000000..82e8a6201 --- /dev/null +++ b/numpy/oldnumeric/setupscons.py @@ -0,0 +1,8 @@ + +def configuration(parent_package='',top_path=None): + from numpy.distutils.misc_util import Configuration + return Configuration('oldnumeric',parent_package,top_path) + +if __name__ == '__main__': + from numpy.distutils.core import setup + setup(configuration=configuration) diff --git a/numpy/oldnumeric/tests/test_oldnumeric.py b/numpy/oldnumeric/tests/test_oldnumeric.py new file mode 100644 index 000000000..24d709d2c --- /dev/null +++ b/numpy/oldnumeric/tests/test_oldnumeric.py @@ -0,0 +1,94 @@ +import unittest + +from numpy.testing import * + +from numpy import array +from numpy.oldnumeric import * +from numpy.core.numeric import float32, float64, complex64, complex128, int8, \ + int16, int32, int64, uint, uint8, uint16, uint32, uint64 + +class test_oldtypes(unittest.TestCase): + def test_oldtypes(self, level=1): + a1 = array([0,1,0], Float) + a2 = array([0,1,0], float) + assert_array_equal(a1, a2) + a1 = array([0,1,0], Float8) + a2 = array([0,1,0], float) + assert_array_equal(a1, a2) + a1 = array([0,1,0], Float16) + a2 = array([0,1,0], float) + assert_array_equal(a1, a2) + a1 = array([0,1,0], Float32) + a2 = array([0,1,0], float32) + assert_array_equal(a1, a2) + a1 = array([0,1,0], Float64) + a2 = array([0,1,0], float64) + assert_array_equal(a1, a2) + a1 = array([0,1,0], Complex) + a2 = array([0,1,0], complex) + assert_array_equal(a1, a2) + a1 = array([0,1,0], Complex8) + a2 = array([0,1,0], complex) + assert_array_equal(a1, a2) + a1 = array([0,1,0], Complex16) + a2 = array([0,1,0], complex) + assert_array_equal(a1, a2) + a1 = array([0,1,0], Complex32) + a2 = array([0,1,0], complex64) + assert_array_equal(a1, a2) + a1 = array([0,1,0], Complex64) + a2 = array([0,1,0], complex128) + assert_array_equal(a1, a2) + a1 = array([0,1,0], Int) + a2 = array([0,1,0], int) + assert_array_equal(a1, a2) + a1 = array([0,1,0], Int8) + a2 = array([0,1,0], int8) + assert_array_equal(a1, a2) + a1 = array([0,1,0], Int16) + a2 = array([0,1,0], int16) + assert_array_equal(a1, a2) + a1 = array([0,1,0], Int32) + a2 = array([0,1,0], int32) + assert_array_equal(a1, a2) + try: + a1 = array([0,1,0], Int64) + a2 = array([0,1,0], int64) + assert_array_equal(a1, a2) + except NameError: + # Not all systems have 64-bit integers. + pass + a1 = array([0,1,0], UnsignedInt) + a2 = array([0,1,0], UnsignedInteger) + a3 = array([0,1,0], uint) + assert_array_equal(a1, a3) + assert_array_equal(a2, a3) + a1 = array([0,1,0], UInt8) + a2 = array([0,1,0], UnsignedInt8) + a3 = array([0,1,0], uint8) + assert_array_equal(a1, a3) + assert_array_equal(a2, a3) + a1 = array([0,1,0], UInt16) + a2 = array([0,1,0], UnsignedInt16) + a3 = array([0,1,0], uint16) + assert_array_equal(a1, a3) + assert_array_equal(a2, a3) + a1 = array([0,1,0], UInt32) + a2 = array([0,1,0], UnsignedInt32) + a3 = array([0,1,0], uint32) + assert_array_equal(a1, a3) + assert_array_equal(a2, a3) + try: + a1 = array([0,1,0], UInt64) + a2 = array([0,1,0], UnsignedInt64) + a3 = array([0,1,0], uint64) + assert_array_equal(a1, a3) + assert_array_equal(a2, a3) + except NameError: + # Not all systems have 64-bit integers. + pass + + +if __name__ == "__main__": + import nose + nose.main() diff --git a/numpy/oldnumeric/typeconv.py b/numpy/oldnumeric/typeconv.py new file mode 100644 index 000000000..4e203d4ae --- /dev/null +++ b/numpy/oldnumeric/typeconv.py @@ -0,0 +1,60 @@ +__all__ = ['oldtype2dtype', 'convtypecode', 'convtypecode2', 'oldtypecodes'] + +import numpy as np + +oldtype2dtype = {'1': np.dtype(np.byte), + 's': np.dtype(np.short), +# 'i': np.dtype(np.intc), +# 'l': np.dtype(int), +# 'b': np.dtype(np.ubyte), + 'w': np.dtype(np.ushort), + 'u': np.dtype(np.uintc), +# 'f': np.dtype(np.single), +# 'd': np.dtype(float), +# 'F': np.dtype(np.csingle), +# 'D': np.dtype(complex), +# 'O': np.dtype(object), +# 'c': np.dtype('c'), + None: np.dtype(int) + } + +# converts typecode=None to int +def convtypecode(typecode, dtype=None): + if dtype is None: + try: + return oldtype2dtype[typecode] + except: + return np.dtype(typecode) + else: + return dtype + +#if both typecode and dtype are None +# return None +def convtypecode2(typecode, dtype=None): + if dtype is None: + if typecode is None: + return None + else: + try: + return oldtype2dtype[typecode] + except: + return np.dtype(typecode) + else: + return dtype + +_changedtypes = {'B': 'b', + 'b': '1', + 'h': 's', + 'H': 'w', + 'I': 'u'} + +class _oldtypecodes(dict): + def __getitem__(self, obj): + char = np.dtype(obj).char + try: + return _changedtypes[char] + except KeyError: + return char + + +oldtypecodes = _oldtypecodes() diff --git a/numpy/oldnumeric/ufuncs.py b/numpy/oldnumeric/ufuncs.py new file mode 100644 index 000000000..c26050f55 --- /dev/null +++ b/numpy/oldnumeric/ufuncs.py @@ -0,0 +1,19 @@ +__all__ = ['less', 'cosh', 'arcsinh', 'add', 'ceil', 'arctan2', 'floor_divide', + 'fmod', 'hypot', 'logical_and', 'power', 'sinh', 'remainder', 'cos', + 'equal', 'arccos', 'less_equal', 'divide', 'bitwise_or', + 'bitwise_and', 'logical_xor', 'log', 'subtract', 'invert', + 'negative', 'log10', 'arcsin', 'arctanh', 'logical_not', + 'not_equal', 'tanh', 'true_divide', 'maximum', 'arccosh', + 'logical_or', 'minimum', 'conjugate', 'tan', 'greater', + 'bitwise_xor', 'fabs', 'floor', 'sqrt', 'arctan', 'right_shift', + 'absolute', 'sin', 'multiply', 'greater_equal', 'left_shift', + 'exp', 'divide_safe'] + +from numpy import less, cosh, arcsinh, add, ceil, arctan2, floor_divide, \ + fmod, hypot, logical_and, power, sinh, remainder, cos, \ + equal, arccos, less_equal, divide, bitwise_or, bitwise_and, \ + logical_xor, log, subtract, invert, negative, log10, arcsin, \ + arctanh, logical_not, not_equal, tanh, true_divide, maximum, \ + arccosh, logical_or, minimum, conjugate, tan, greater, bitwise_xor, \ + fabs, floor, sqrt, arctan, right_shift, absolute, sin, \ + multiply, greater_equal, left_shift, exp, divide as divide_safe diff --git a/numpy/oldnumeric/user_array.py b/numpy/oldnumeric/user_array.py new file mode 100644 index 000000000..375c4013b --- /dev/null +++ b/numpy/oldnumeric/user_array.py @@ -0,0 +1,9 @@ + + +from numpy.oldnumeric import * +from numpy.lib.user_array import container as UserArray + +import numpy.oldnumeric as nold +__all__ = nold.__all__[:] +__all__ += ['UserArray'] +del nold |