summaryrefslogtreecommitdiff
path: root/numpy/lib
diff options
context:
space:
mode:
Diffstat (limited to 'numpy/lib')
-rw-r--r--numpy/lib/convertcode.py147
-rw-r--r--numpy/lib/function_base.py815
-rw-r--r--numpy/lib/getlimits.py118
-rw-r--r--numpy/lib/index_tricks.py291
-rw-r--r--numpy/lib/machar.py268
-rw-r--r--numpy/lib/mlab.py14
-rw-r--r--numpy/lib/polynomial.py554
-rw-r--r--numpy/lib/scimath.py77
-rw-r--r--numpy/lib/shape_base.py539
-rw-r--r--numpy/lib/src/_compiled_base.c453
-rw-r--r--numpy/lib/test_shape_base.py364
-rw-r--r--numpy/lib/tests/test_function_base.py338
-rw-r--r--numpy/lib/tests/test_getlimits.py38
-rw-r--r--numpy/lib/tests/test_index_tricks.py53
-rw-r--r--numpy/lib/tests/test_polynomial.py83
-rw-r--r--numpy/lib/tests/test_twodim_base.py134
-rw-r--r--numpy/lib/twodim_base.py123
-rw-r--r--numpy/lib/type_check.py180
-rw-r--r--numpy/lib/ufunclike.py77
-rw-r--r--numpy/lib/utils.py28
20 files changed, 4694 insertions, 0 deletions
diff --git a/numpy/lib/convertcode.py b/numpy/lib/convertcode.py
new file mode 100644
index 000000000..5c532b394
--- /dev/null
+++ b/numpy/lib/convertcode.py
@@ -0,0 +1,147 @@
+
+# This module converts code written for Numeric to run with scipy.base
+
+# Makes the following changes:
+# * Converts typecharacters
+# * 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()
+# * Converts .flat to .ravel() except for .flat = xxx or .flat[xxx]
+# * Change typecode= to dtype=
+# * Eliminates savespace=xxx
+# * Replace xxx.spacesaver() with True
+# * Convert xx.savespace(?) to pass + ## xx.savespace(?)
+# #### -- not * Convert a.shape = ? to a.reshape(?)
+# * Prints warning for use of bool, int, float, copmlex, object, and unicode
+#
+
+__all__ = ['fromfile', 'fromstr']
+
+import sys
+import os
+import re
+import glob
+
+flatindex_re = re.compile('([.]flat(\s*?[[=]))')
+
+def replacetypechars(astr):
+# astr = astr.replace("'s'","'h'")
+# astr = astr.replace("'c'","'S1'")
+ 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(".typecode()",".dtypechar")
+ 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("@@@@\\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
+
+svspc = re.compile(r'(\S+\s*[(].+),\s*savespace\s*=.+\s*[)]')
+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 = astr.replace("typecode=","dtype=")
+ astr = astr.replace("UserArray","ndarray")
+ astr = svspc.sub('\\1)',astr)
+ astr = svspc2.sub('True',astr)
+ astr = svspc3.sub('pass ## \\1', astr)
+ #astr = shpe.sub('\\1=\\1.reshape(\\2)', astr)
+ return astr
+
+import datetime
+def fromstr(filestr):
+ filestr = replacetypechars(filestr)
+ filestr, fromall1 = changeimports(filestr, 'Numeric', 'scipy')
+ filestr, fromall1 = changeimports(filestr, 'multiarray',
+ 'scipy.base.multiarray')
+ filestr, fromall1 = changeimports(filestr, 'umath',
+ 'scipy.base.umath')
+ filestr, fromall1 = changeimports(filestr, 'Precision', 'scipy.base')
+ filestr, fromall2 = changeimports(filestr, 'numerix', 'scipy.base')
+ filestr, fromall3 = changeimports(filestr, 'scipy_base', 'scipy.base')
+ filestr, fromall3 = changeimports(filestr, 'MLab', 'scipy.base.mlab')
+ filestr, fromall3 = changeimports(filestr, 'LinearAlgebra', 'scipy.corelinalg')
+ filestr, fromall3 = changeimports(filestr, 'RNG', 'scipy.random')
+ filestr, fromall3 = changeimports(filestr, 'RandomArray', 'scipy.random')
+ filestr, fromall3 = changeimports(filestr, 'FFT', 'scipy.corefft')
+ filestr, fromall3 = changeimports(filestr, 'MA', 'scipy.base.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 '\
+ 'scipy %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 fromfile(filename):
+ filestr = getandcopy(filename)
+ filestr = fromstr(filestr)
+ makenewfile(filename, filestr)
+
+def fromargs(args):
+ filename = args[1]
+ fromfile(filename)
+
+def convertall(direc=''):
+ files = glob.glob(os.path.join(direc,'*.py'))
+ for afile in files:
+ fromfile(afile)
+
+if __name__ == '__main__':
+ fromargs(sys.argv)
+
+
+
diff --git a/numpy/lib/function_base.py b/numpy/lib/function_base.py
new file mode 100644
index 000000000..60e4b4be0
--- /dev/null
+++ b/numpy/lib/function_base.py
@@ -0,0 +1,815 @@
+
+l__all__ = ['logspace', 'linspace', 'round_',
+ 'select', 'piecewise', 'trim_zeros',
+ 'copy', 'iterable', 'base_repr', 'binary_repr',
+ 'diff', 'gradient', 'angle', 'unwrap', 'sort_complex', 'disp',
+ 'unique', 'extract', 'insert', 'nansum', 'nanmax', 'nanargmax',
+ 'nanargmin', 'nanmin', 'vectorize', 'asarray_chkfinite', 'average',
+ 'histogram', 'bincount', 'digitize', 'cov', 'corrcoef', 'msort',
+ 'median', 'sinc', 'hamming', 'hanning', 'bartlett', 'blackman',
+ 'kaiser', 'trapz'
+ ]
+
+import types
+import math
+import numeric as _nx
+from numeric import ones, zeros, arange, concatenate, array, asarray, empty
+from numeric import ScalarType, dot, where, newaxis
+from umath import pi, multiply, add, arctan2, maximum, minimum, frompyfunc, \
+ isnan, absolute, cos, less_equal, sqrt, sin, mod
+from oldnumeric import ravel, nonzero, choose, \
+ sometrue, alltrue, reshape, any, all, typecodes, ArrayType, squeeze,\
+ sort
+from type_check import ScalarType, isscalar
+from shape_base import atleast_1d
+from twodim_base import diag
+from _compiled_base import digitize, bincount, _insert
+from ufunclike import sign
+
+_lkup = {'0':'000',
+ '1':'001',
+ '2':'010',
+ '3':'011',
+ '4':'100',
+ '5':'101',
+ '6':'110',
+ '7':'111',
+ 'L':''}
+
+def binary_repr(num):
+ """Return the binary representation of the input number as a string.
+
+ This is equivalent to using base_repr with base 2, but about 25x
+ faster.
+ """
+ ostr = oct(num)
+ bin = ''
+ for ch in ostr[1:]:
+ bin += _lkup[ch]
+ ind = 0
+ while bin[ind] == '0':
+ ind += 1
+ return bin[ind:]
+
+def base_repr (number, base=2, padding=0):
+ """Return the representation of a number in any given base.
+ """
+ chars = '0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ'
+
+ lnb = math.log(base)
+ res = padding*chars[0]
+ if number == 0:
+ return res + chars[0]
+ exponent = int (math.log (number)/lnb)
+ while(exponent >= 0):
+ term = long(base)**exponent
+ lead_digit = int(number / term)
+ res += chars[lead_digit]
+ number -= term*lead_digit
+ exponent -= 1
+ return res
+#end Fernando's utilities
+
+
+def linspace(start, stop, num=50, endpoint=True, retstep=False):
+ """Return evenly spaced numbers.
+
+ Return 'num' evenly spaced samples from 'start' to 'stop'. If
+ 'endpoint' is True, the last sample is 'stop'. If 'retstep' is
+ True then return the step value used.
+ """
+ num = int(num)
+ if num <= 0:
+ return array([])
+ if endpoint:
+ if num == 1:
+ return array([start])
+ step = (stop-start)/float((num-1))
+ else:
+ step = (stop-start)/float(num)
+ y = _nx.arange(0, num) * step + start
+ if retstep:
+ return y, step
+ else:
+ return y
+
+def logspace(start,stop,num=50,endpoint=True,base=10.0):
+ """Evenly spaced numbers on a logarithmic scale.
+
+ Computes int(num) evenly spaced exponents from start to stop.
+ If endpoint=True, then last exponent is stop.
+ Returns base**exponents.
+ """
+ y = linspace(start,stop,num=num,endpoint=endpoint)
+ return _nx.power(base,y)
+
+def iterable(y):
+ try: iter(y)
+ except: return 0
+ return 1
+
+def histogram(a, bins=10, range=None, normed=False):
+ a = asarray(a).ravel()
+ if not iterable(bins):
+ if range is None:
+ range = (a.min(), a.max())
+ mn, mx = [mi+0.0 for mi in range]
+ if mn == mx:
+ mn -= 0.5
+ mx += 0.5
+ bins = linspace(mn, mx, bins, endpoint=False)
+
+ n = sort(a).searchsorted(bins)
+ n = concatenate([n, [len(a)]])
+ n = n[1:]-n[:-1]
+
+ if normed:
+ db = bins[1] - bins[0]
+ return 1.0/(a.size*db) * n, bins
+ else:
+ return n, bins
+
+def average(a, axis=0, weights=None, returned=False):
+ """average(a, axis=0, weights=None, returned=False)
+
+ Average the array over the given axis. If the axis is None, average
+ over all dimensions of the array. Equivalent to a.mean(axis), but
+ with a default axis of 0 instead of None.
+
+ If an integer axis is given, this equals:
+ a.sum(axis) * 1.0 / len(a)
+
+ If axis is None, this equals:
+ a.sum(axis) * 1.0 / product(a.shape)
+
+ If weights are given, result is:
+ sum(a * weights) / sum(weights),
+ where the weights must have a's shape or be 1D with length the
+ size of a in the given axis. Integer weights are converted to
+ Float. Not specifying weights is equivalent to specifying
+ weights that are all 1.
+
+ If 'returned' is True, return a tuple: the result and the sum of
+ the weights or count of values. The shape of these two results
+ will be the same.
+
+ Raises ZeroDivisionError if appropriate. (The version in MA does
+ not -- it returns masked values).
+ """
+ if axis is None:
+ a = array(a).ravel()
+ if weights is None:
+ n = add.reduce(a)
+ d = len(a) * 1.0
+ else:
+ w = array(weights).ravel() * 1.0
+ n = add.reduce(multiply(a, w))
+ d = add.reduce(w)
+ else:
+ a = array(a)
+ ash = a.shape
+ if ash == ():
+ a.shape = (1,)
+ if weights is None:
+ n = add.reduce(a, axis)
+ d = ash[axis] * 1.0
+ if returned:
+ d = ones(n.shape) * d
+ else:
+ w = array(weights, copy=False) * 1.0
+ wsh = w.shape
+ if wsh == ():
+ wsh = (1,)
+ if wsh == ash:
+ n = add.reduce(a*w, axis)
+ d = add.reduce(w, axis)
+ elif wsh == (ash[axis],):
+ ni = ash[axis]
+ r = [newaxis]*ni
+ r[axis] = slice(None, None, 1)
+ w1 = eval("w["+repr(tuple(r))+"]*ones(ash, Float)")
+ n = add.reduce(a*w1, axis)
+ d = add.reduce(w1, axis)
+ else:
+ raise ValueError, 'averaging weights have wrong shape'
+
+ if not isinstance(d, ArrayType):
+ if d == 0.0:
+ raise ZeroDivisionError, 'zero denominator in average()'
+ if returned:
+ return n/d, d
+ else:
+ return n/d
+
+def asarray_chkfinite(a):
+ """Like asarray, but check that no NaNs or Infs are present.
+ """
+ a = asarray(a)
+ if (a.dtypechar in _nx.typecodes['AllFloat']) \
+ and (_nx.isnan(a).any() or _nx.isinf(a).any()):
+ raise ValueError, "array must not contain infs or NaNs"
+ return a
+
+
+
+
+def piecewise(x, condlist, funclist, *args, **kw):
+ """Return a piecewise-defined function.
+
+ x is the domain
+
+ condlist is a list of boolean arrays or a single boolean array
+ The length of the condition list must be n2 or n2-1 where n2
+ is the length of the function list. If len(condlist)==n2-1, then
+ an 'otherwise' condition is formed by |'ing all the conditions
+ and inverting.
+
+ funclist is a list of functions to call of length (n2).
+ Each function should return an array output for an array input
+ Each function can take (the same set) of extra arguments and
+ keyword arguments which are passed in after the function list.
+
+ The output is the same shape and type as x and is found by
+ calling the functions on the appropriate portions of x.
+
+ Note: This is similar to choose or select, except
+ the the functions are only evaluated on elements of x
+ that satisfy the corresponding condition.
+
+ The result is
+ |--
+ | f1(x) for condition1
+ y = --| f2(x) for condition2
+ | ...
+ | fn(x) for conditionn
+ |--
+
+ """
+ n2 = len(funclist)
+ if not isinstance(condlist, type([])):
+ condlist = [condlist]
+ n = len(condlist)
+ if n == n2-1: # compute the "otherwise" condition.
+ totlist = condlist[0]
+ for k in range(1, n):
+ totlist |= condlist
+ condlist.append(~totlist)
+ n += 1
+ if (n != n2):
+ raise ValueError, "function list and condition list must be the same"
+ y = empty(x.shape, x.dtype)
+ for k in range(n):
+ item = funclist[k]
+ if not callable(item):
+ y[condlist[k]] = item
+ else:
+ y[condlist[k]] = item(x[condlist[k]], *args, **kw)
+ return y
+
+def select(condlist, choicelist, default=0):
+ """ Return an array composed of different elements of choicelist
+ depending on the list of conditions.
+
+ condlist is a list of condition arrays containing ones or zeros
+
+ choicelist is a list of choice arrays (of the "same" size as the
+ arrays in condlist). The result array has the "same" size as the
+ arrays in choicelist. If condlist is [c0, ..., cN-1] then choicelist
+ must be of length N. The elements of the choicelist can then be
+ represented as [v0, ..., vN-1]. The default choice if none of the
+ conditions are met is given as the default argument.
+
+ The conditions are tested in order and the first one statisfied is
+ used to select the choice. In other words, the elements of the
+ output array are found from the following tree (notice the order of
+ the conditions matters):
+
+ if c0: v0
+ elif c1: v1
+ elif c2: v2
+ ...
+ elif cN-1: vN-1
+ else: default
+
+ Note that one of the condition arrays must be large enough to handle
+ the largest array in the choice list.
+ """
+ n = len(condlist)
+ n2 = len(choicelist)
+ if n2 != n:
+ raise ValueError, "list of cases must be same length as list of conditions"
+ choicelist.insert(0, default)
+ S = 0
+ pfac = 1
+ for k in range(1, n+1):
+ S += k * pfac * asarray(condlist[k-1])
+ if k < n:
+ pfac *= (1-asarray(condlist[k-1]))
+ # handle special case of a 1-element condition but
+ # a multi-element choice
+ if type(S) in ScalarType or max(asarray(S).shape)==1:
+ pfac = asarray(1)
+ for k in range(n2+1):
+ pfac = pfac + asarray(choicelist[k])
+ S = S*ones(asarray(pfac).shape)
+ return choose(S, tuple(choicelist))
+
+def _asarray1d(arr, copy=False):
+ """Ensure 1D array for one array.
+ """
+ if copy:
+ return asarray(arr).flatten()
+ else:
+ return asarray(arr).ravel()
+
+def copy(a):
+ """Return an array copy of the given object.
+ """
+ return array(a, copy=True)
+
+# Basic operations
+
+def gradient(f, *varargs):
+ """Calculate the gradient of an N-dimensional scalar function.
+
+ Uses central differences on the interior and first differences on boundaries
+ to give the same shape.
+
+ Inputs:
+
+ f -- An N-dimensional array giving samples of a scalar function
+
+ varargs -- 0, 1, or N scalars giving the sample distances in each direction
+
+ Outputs:
+
+ N arrays of the same shape as f giving the derivative of f with respect
+ to each dimension.
+ """
+ N = len(f.shape) # number of dimensions
+ n = len(varargs)
+ if n==0:
+ dx = [1.0]*N
+ elif n==1:
+ dx = [varargs[0]]*N
+ elif n==N:
+ dx = list(varargs)
+ else:
+ raise SyntaxError, "invalid number of arguments"
+
+ # use central differences on interior and first differences on endpoints
+
+ print dx
+ outvals = []
+
+ # create slice objects --- initially all are [:, :, ..., :]
+ slice1 = [slice(None)]*N
+ slice2 = [slice(None)]*N
+ slice3 = [slice(None)]*N
+
+ otype = f.dtypechar
+ if otype not in ['f', 'd', 'F', 'D']:
+ otype = 'd'
+
+ for axis in range(N):
+ # select out appropriate parts for this dimension
+ out = zeros(f.shape, f.dtypechar)
+ slice1[axis] = slice(1, -1)
+ slice2[axis] = slice(2, None)
+ slice3[axis] = slice(None, -2)
+ # 1D equivalent -- out[1:-1] = (f[2:] - f[:-2])/2.0
+ out[slice1] = (f[slice2] - f[slice3])/2.0
+ slice1[axis] = 0
+ slice2[axis] = 1
+ slice3[axis] = 0
+ # 1D equivalent -- out[0] = (f[1] - f[0])
+ out[slice1] = (f[slice2] - f[slice3])
+ slice1[axis] = -1
+ slice2[axis] = -1
+ slice3[axis] = -2
+ # 1D equivalent -- out[-1] = (f[-1] - f[-2])
+ out[slice1] = (f[slice2] - f[slice3])
+
+ # divide by step size
+ outvals.append(out / dx[axis])
+
+ # reset the slice object in this dimension to ":"
+ slice1[axis] = slice(None)
+ slice2[axis] = slice(None)
+ slice3[axis] = slice(None)
+
+ if N == 1:
+ return outvals[0]
+ else:
+ return outvals
+
+
+def diff(a, n=1, axis=-1):
+ """Calculate the nth order discrete difference along given axis.
+ """
+ if n==0:
+ return a
+ if n<0:
+ raise ValueError, 'order must be non-negative but got ' + `n`
+ a = asarray(a)
+ nd = len(a.shape)
+ slice1 = [slice(None)]*nd
+ slice2 = [slice(None)]*nd
+ slice1[axis] = slice(1, None)
+ slice2[axis] = slice(None, -1)
+ slice1 = tuple(slice1)
+ slice2 = tuple(slice2)
+ if n > 1:
+ return diff(a[slice1]-a[slice2], n-1, axis=axis)
+ else:
+ return a[slice1]-a[slice2]
+
+def angle(z, deg=0):
+ """Return the angle of the complex argument z.
+ """
+ if deg:
+ fact = 180/pi
+ else:
+ fact = 1.0
+ z = asarray(z)
+ if (issubclass(z.dtype, _nx.complexfloating)):
+ zimag = z.imag
+ zreal = z.real
+ else:
+ zimag = 0
+ zreal = z
+ return arctan2(zimag, zreal) * fact
+
+def unwrap(p, discont=pi, axis=-1):
+ """Unwrap radian phase p by changing absolute jumps greater than
+ 'discont' to their 2*pi complement along the given axis.
+ """
+ p = asarray(p)
+ nd = len(p.shape)
+ dd = diff(p, axis=axis)
+ slice1 = [slice(None, None)]*nd # full slices
+ slice1[axis] = slice(1, None)
+ ddmod = mod(dd+pi, 2*pi)-pi
+ _nx.putmask(ddmod, (ddmod==-pi) & (dd > 0), pi)
+ ph_correct = ddmod - dd;
+ _nx.putmask(ph_correct, abs(dd)<discont, 0)
+ up = array(p, copy=True, dtype='d')
+ up[slice1] = p[slice1] + ph_correct.cumsum(axis)
+ return up
+
+def sort_complex(a):
+ """ Sort 'a' as a complex array using the real part first and then
+ the imaginary part if the real part is equal (the default sort order
+ for complex arrays). This function is a wrapper ensuring a complex
+ return type.
+ """
+ b = array(a,copy=True)
+ b.sort()
+ if not issubclass(b.dtype, _nx.complexfloating):
+ if b.dtypechar in 'bhBH':
+ return b.astype('F')
+ elif b.dtypechar == 'g':
+ return b.astype('G')
+ else:
+ return b.astype('D')
+ else:
+ return b
+
+def trim_zeros(filt, trim='fb'):
+ """ Trim the leading and trailing zeros from a 1D array.
+
+ Example:
+ >>> import scipy
+ >>> a = array((0, 0, 0, 1, 2, 3, 2, 1, 0))
+ >>> scipy.trim_zeros(a)
+ array([1, 2, 3, 2, 1])
+ """
+ first = 0
+ trim = trim.upper()
+ if 'F' in trim:
+ for i in filt:
+ if i != 0.: break
+ else: first = first + 1
+ last = len(filt)
+ if 'B' in trim:
+ for i in filt[::-1]:
+ if i != 0.: break
+ else: last = last - 1
+ return filt[first:last]
+
+def unique(inseq):
+ """Return unique items from a 1-dimensional sequence.
+ """
+ # Dictionary setting is quite fast.
+ set = {}
+ for item in inseq:
+ set[item] = None
+ return asarray(set.keys())
+
+def extract(condition, arr):
+ """Return the elements of ravel(arr) where ravel(condition) is True
+ (in 1D).
+
+ Equivalent to compress(ravel(condition), ravel(arr)).
+ """
+ return _nx.take(ravel(arr), nonzero(ravel(condition)))
+
+def insert(arr, mask, vals):
+ """Similar to putmask arr[mask] = vals but the 1D array vals has the
+ same number of elements as the non-zero values of mask. Inverse of
+ extract.
+ """
+ return _insert(arr, mask, vals)
+
+def nansum(a, axis=-1):
+ """Sum the array over the given axis, treating NaNs as 0.
+ """
+ y = array(a)
+ if not issubclass(y.dtype, _nx.integer):
+ y[isnan(a)] = 0
+ return y.sum(axis)
+
+def nanmin(a, axis=-1):
+ """Find the minimium over the given axis, ignoring NaNs.
+ """
+ y = array(a)
+ if not issubclass(y.dtype, _nx.integer):
+ y[isnan(a)] = _nx.inf
+ return y.min(axis)
+
+def nanargmin(a, axis=-1):
+ """Find the indices of the minimium over the given axis ignoring NaNs.
+ """
+ y = array(a)
+ if not issubclass(y.dtype, _nx.integer):
+ y[isnan(a)] = _nx.inf
+ return y.argmin(axis)
+
+def nanmax(a, axis=-1):
+ """Find the maximum over the given axis ignoring NaNs.
+ """
+ y = array(a)
+ if not issubclass(y.dtype, _nx.integer):
+ y[isnan(a)] = -_nx.inf
+ return y.max(axis)
+
+def nanargmax(a, axis=-1):
+ """Find the maximum over the given axis ignoring NaNs.
+ """
+ y = array(a)
+ if not issubclass(y.dtype, _nx.integer):
+ y[isnan(a)] = -_nx.inf
+ return y.argmax(axis)
+
+def disp(mesg, device=None, linefeed=True):
+ """Display a message to the given device (default is sys.stdout)
+ with or without a linefeed.
+ """
+ if device is None:
+ import sys
+ device = sys.stdout
+ if linefeed:
+ device.write('%s\n' % mesg)
+ else:
+ device.write('%s' % mesg)
+ device.flush()
+ return
+
+class vectorize(object):
+ """
+ vectorize(somefunction, otypes=None, doc=None)
+ Generalized Function class.
+
+ Description:
+
+ Define a vectorized function which takes nested sequence
+ objects or scipy arrays as inputs and returns a
+ scipy array as output, evaluating the function over successive
+ tuples of the input arrays like the python map function except it uses
+ the broadcasting rules of scipy.
+
+ Input:
+
+ somefunction -- a Python function or method
+
+ Example:
+
+ def myfunc(a, b):
+ if a > b:
+ return a-b
+ else
+ return a+b
+
+ vfunc = vectorize(myfunc)
+
+ >>> vfunc([1, 2, 3, 4], 2)
+ array([3, 4, 1, 2])
+
+ """
+ def __init__(self, pyfunc, otypes='', doc=None):
+ try:
+ fcode = pyfunc.func_code
+ except AttributeError:
+ raise TypeError, "object is not a callable Python object"
+
+ self.thefunc = pyfunc
+ self.ufunc = None
+ self.nin = fcode.co_argcount
+ if pyfunc.func_defaults:
+ self.nin_wo_defaults = self.nin - len(pyfunc.func_defaults)
+ else:
+ self.nin_wo_defaults = self.nin
+ self.nout = None
+ if doc is None:
+ self.__doc__ = pyfunc.__doc__
+ else:
+ self.__doc__ = doc
+ if isinstance(otypes, types.StringType):
+ self.otypes=otypes
+ else:
+ raise ValueError, "output types must be a string"
+ for char in self.otypes:
+ if char not in typecodes['All']:
+ raise ValueError, "invalid typecode specified"
+ self.lastcallargs = 0
+
+ def __call__(self, *args):
+ # get number of outputs and output types by calling
+ # the function on the first entries of args
+ nargs = len(args)
+ if (nargs > self.nin) or (nargs < self.nin_wo_defaults):
+ raise ValueError, "mismatch between python function inputs"\
+ " and received arguments"
+ if self.nout is None or self.otypes == '':
+ newargs = []
+ for arg in args:
+ newargs.append(asarray(arg).flat[0])
+ theout = self.thefunc(*newargs)
+ if isinstance(theout, types.TupleType):
+ self.nout = len(theout)
+ else:
+ self.nout = 1
+ theout = (theout,)
+ if self.otypes == '':
+ otypes = []
+ for k in range(self.nout):
+ otypes.append(asarray(theout[k]).dtypechar)
+ self.otypes = ''.join(otypes)
+
+ if (self.ufunc is None) or (self.lastcallargs != nargs):
+ self.ufunc = frompyfunc(self.thefunc, nargs, self.nout)
+ self.lastcallargs = nargs
+
+ if self.nout == 1:
+ return self.ufunc(*args).astype(self.otypes[0])
+ else:
+ return tuple([x.astype(c) for x, c in zip(self.ufunc(*args), self.otypes)])
+
+
+def round_(a, decimals=0):
+ """Round 'a' to the given number of decimal places. Rounding
+ behaviour is equivalent to Python.
+
+ Return 'a' if the array is not floating point. Round both the real
+ and imaginary parts separately if the array is complex.
+ """
+ a = asarray(a)
+ if not issubclass(a.dtype, _nx.inexact):
+ return a
+ if issubclass(a.dtype, _nx.complexfloating):
+ return round_(a.real, decimals) + 1j*round_(a.imag, decimals)
+ if decimals is not 0:
+ decimals = asarray(decimals)
+ s = sign(a)
+ if decimals is not 0:
+ a = absolute(multiply(a, 10.**decimals))
+ else:
+ a = absolute(a)
+ rem = a-asarray(a).astype(_nx.intp)
+ a = _nx.where(_nx.less(rem, 0.5), _nx.floor(a), _nx.ceil(a))
+ # convert back
+ if decimals is not 0:
+ return multiply(a, s/(10.**decimals))
+ else:
+ return multiply(a, s)
+
+
+def cov(m,y=None, rowvar=0, bias=0):
+ """Estimate the covariance matrix.
+
+ If m is a vector, return the variance. For matrices where each row
+ is an observation, and each column a variable, return the covariance
+ matrix. Note that in this case diag(cov(m)) is a vector of
+ variances for each column.
+
+ cov(m) is the same as cov(m, m)
+
+ Normalization is by (N-1) where N is the number of observations
+ (unbiased estimate). If bias is 1 then normalization is by N.
+
+ If rowvar is zero, then each row is a variable with
+ observations in the columns.
+ """
+ if y is None:
+ y = asarray(m)
+ else:
+ y = asarray(y)
+ m = asarray(m)
+ if rowvar:
+ m = m.transpose()
+ y = y.transpose()
+ if (m.shape[0] == 1):
+ m = m.transpose()
+ if (y.shape[0] == 1):
+ y = y.transpose()
+ N = m.shape[0]
+ if (y.shape[0] != N):
+ raise ValueError, "x and y must have the same number of observations."
+ m = m - m.mean(axis=0)
+ y = y - y.mean(axis=0)
+ if bias:
+ fact = N*1.0
+ else:
+ fact = N-1.0
+
+ val = squeeze(dot(m.transpose(),y.conj()) / fact)
+ return val
+
+def corrcoef(x, y=None):
+ """The correlation coefficients
+ """
+ c = cov(x, y)
+ d = diag(c)
+ return c/sqrt(multiply.outer(d,d))
+
+def blackman(M):
+ """blackman(M) returns the M-point Blackman window.
+ """
+ n = arange(0,M)
+ return 0.42-0.5*cos(2.0*pi*n/(M-1)) + 0.08*cos(4.0*pi*n/(M-1))
+
+def bartlett(M):
+ """bartlett(M) returns the M-point Bartlett window.
+ """
+ n = arange(0,M)
+ return where(less_equal(n,(M-1)/2.0),2.0*n/(M-1),2.0-2.0*n/(M-1))
+
+def hanning(M):
+ """hanning(M) returns the M-point Hanning window.
+ """
+ n = arange(0,M)
+ return 0.5-0.5*cos(2.0*pi*n/(M-1))
+
+def hamming(M):
+ """hamming(M) returns the M-point Hamming window.
+ """
+ n = arange(0,M)
+ return 0.54-0.46*cos(2.0*pi*n/(M-1))
+
+def kaiser(M,beta):
+ """kaiser(M, beta) returns a Kaiser window of length M with shape parameter
+ beta. It depends on scipy.special (in full scipy) for the modified bessel
+ function i0.
+ """
+ from scipy.special import i0
+ n = arange(0,M)
+ alpha = (M-1)/2.0
+ return i0(beta * sqrt(1-((n-alpha)/alpha)**2.0))/i0(beta)
+
+def sinc(x):
+ """sinc(x) returns sin(pi*x)/(pi*x) at all points of array x.
+ """
+ y = pi* where(x == 0, 1.0e-20, x)
+ return sin(y)/y
+
+def msort(a):
+ b = array(a,copy=True)
+ b.sort(0)
+ return b
+
+def median(m):
+ """median(m) returns a median of m along the first dimension of m.
+ """
+ sorted = msort(m)
+ if sorted.shape[0] % 2 == 1:
+ return sorted[int(sorted.shape[0]/2)]
+ else:
+ sorted = msort(m)
+ index=sorted.shape[0]/2
+ return (sorted[index-1]+sorted[index])/2.0
+
+def trapz(y, x=None, dx=1.0, axis=-1):
+ """Integrate y(x) using samples along the given axis and the composite
+ trapezoidal rule. If x is None, spacing given by dx is assumed.
+ """
+ y = asarray(y)
+ if x is None:
+ d = dx
+ else:
+ d = diff(x,axis=axis)
+ nd = len(y.shape)
+ slice1 = [slice(None)]*nd
+ slice2 = [slice(None)]*nd
+ slice1[axis] = slice(1,None)
+ slice2[axis] = slice(None,-1)
+ return add.reduce(d * (y[slice1]+y[slice2])/2.0,axis)
diff --git a/numpy/lib/getlimits.py b/numpy/lib/getlimits.py
new file mode 100644
index 000000000..41030af2d
--- /dev/null
+++ b/numpy/lib/getlimits.py
@@ -0,0 +1,118 @@
+""" Machine limits for Float32 and Float64 and (long double) if available...
+"""
+
+__all__ = ['finfo']
+
+from machar import MachAr
+import numeric
+from numeric import array
+
+def _frz(a):
+ """fix rank-0 --> rank-1"""
+ if a.ndim == 0: a.shape = (1,)
+ return a
+
+_convert_to_float = {
+ numeric.csingle: numeric.single,
+ numeric.complex_: numeric.float_,
+ numeric.clongfloat: numeric.longfloat
+ }
+
+class finfo(object):
+
+ _finfo_cache = {}
+
+ def __new__(cls, dtype):
+ obj = cls._finfo_cache.get(dtype,None)
+ if obj is not None:
+ return obj
+ dtypes = [dtype]
+ newdtype = numeric.obj2dtype(dtype)
+ if newdtype is not dtype:
+ dtypes.append(newdtype)
+ dtype = newdtype
+ if not issubclass(dtype, numeric.inexact):
+ raise ValueError, "data type %r not inexact" % (dtype)
+ obj = cls._finfo_cache.get(dtype,None)
+ if obj is not None:
+ return obj
+ if not issubclass(dtype, numeric.floating):
+ newdtype = _convert_to_float[dtype]
+ if newdtype is not dtype:
+ dtypes.append(newdtype)
+ dtype = newdtype
+ obj = cls._finfo_cache.get(dtype,None)
+ if obj is not None:
+ return obj
+ obj = object.__new__(cls)._init(dtype)
+ for dt in dtypes:
+ cls._finfo_cache[dt] = obj
+ return obj
+
+ def _init(self, dtype):
+ self.dtype = dtype
+ if dtype is numeric.float_:
+ machar = MachAr(lambda v:array([v],'d'),
+ lambda v:_frz(v.astype('i'))[0],
+ lambda v:array(_frz(v)[0],'d'),
+ lambda v:'%24.16e' % array(_frz(v)[0],'d'),
+ 'scipy float precision floating point '\
+ 'number')
+ elif dtype is numeric.single:
+ machar = MachAr(lambda v:array([v],'f'),
+ lambda v:_frz(v.astype('i'))[0],
+ lambda v:array(_frz(v)[0],'f'), #
+ lambda v:'%15.7e' % array(_frz(v)[0],'f'),
+ "scipy single precision floating "\
+ "point number")
+ elif dtype is numeric.longfloat:
+ machar = MachAr(lambda v:array([v],'g'),
+ lambda v:_frz(v.astype('i'))[0],
+ lambda v:array(_frz(v)[0],'g'), #
+ lambda v:str(array(_frz(v)[0],'g')),
+ "scipy longfloat precision floating "\
+ "point number")
+ else:
+ raise ValueError,`dtype`
+
+ for word in ['tiny', 'precision', 'resolution','iexp',
+ 'maxexp','minexp','epsneg','negep',
+ 'machep']:
+ setattr(self,word,getattr(machar, word))
+ self.max = machar.huge
+ self.min = -self.max
+ self.eps = machar.epsilon
+ self.nexp = machar.iexp
+ self.nmant = machar.it
+ self.machar = machar
+ self._str_tiny = machar._str_xmin
+ self._str_max = machar._str_xmax
+ self._str_epsneg = machar._str_epsneg
+ self._str_eps = machar._str_eps
+ self._str_resolution = machar._str_resolution
+ return self
+
+ def __str__(self):
+ return '''\
+Machine parameters for %(dtype)s
+---------------------------------------------------------------------
+precision=%(precision)3s resolution=%(_str_resolution)s
+machep=%(machep)6s eps= %(_str_eps)s
+negep =%(negep)6s epsneg= %(_str_epsneg)s
+minexp=%(minexp)6s tiny= %(_str_tiny)s
+maxexp=%(maxexp)6s max= %(_str_max)s
+nexp =%(nexp)6s min= -max
+---------------------------------------------------------------------
+''' % self.__dict__
+
+if __name__ == '__main__':
+ f = finfo(numeric.single)
+ print 'single epsilon:',f.eps
+ print 'single tiny:',f.tiny
+ f = finfo(numeric.float)
+ print 'float epsilon:',f.eps
+ print 'float tiny:',f.tiny
+ f = finfo(numeric.longfloat)
+ print 'longfloat epsilon:',f.eps
+ print 'longfloat tiny:',f.tiny
+
diff --git a/numpy/lib/index_tricks.py b/numpy/lib/index_tricks.py
new file mode 100644
index 000000000..71d30a387
--- /dev/null
+++ b/numpy/lib/index_tricks.py
@@ -0,0 +1,291 @@
+## Automatically adapted for scipy Sep 19, 2005 by convertcode.py
+
+__all__ = ['mgrid','ogrid','r_', 'c_', 'index_exp', 'ix_','ndenumerate']
+
+import sys
+import types
+import numeric as _nx
+from numeric import asarray
+
+from type_check import ScalarType
+import function_base
+import twodim_base as matrix_base
+import matrix
+makemat = matrix.matrix
+
+def ix_(*args):
+ """ Construct an open mesh from multiple sequences.
+
+ This function takes n 1-d sequences and returns n outputs with n
+ dimensions each such that the shape is 1 in all but one dimension and
+ the dimension with the non-unit shape value cycles through all n
+ dimensions.
+
+ Using ix_() one can quickly construct index arrays that will index
+ the cross product.
+
+ a[ix_([1,3,7],[2,5,8])] returns the array
+
+ a[1,2] a[1,5] a[1,8]
+ a[3,2] a[3,5] a[3,8]
+ a[7,2] a[7,5] a[7,8]
+ """
+ out = []
+ nd = len(args)
+ baseshape = [1]*nd
+ for k in range(nd):
+ new = _nx.array(args[k])
+ if (new.ndim <> 1):
+ raise ValueError, "Cross index must be 1 dimensional"
+ baseshape[k] = len(new)
+ new.shape = tuple(baseshape)
+ out.append(new)
+ baseshape[k] = 1
+ return tuple(out)
+
+class nd_grid(object):
+ """ Construct a "meshgrid" in N-dimensions.
+
+ grid = nd_grid() creates an instance which will return a mesh-grid
+ when indexed. The dimension and number of the output arrays are equal
+ to the number of indexing dimensions. If the step length is not a
+ complex number, then the stop is not inclusive.
+
+ However, if the step length is a COMPLEX NUMBER (e.g. 5j), then the
+ integer part of it's magnitude is interpreted as specifying the
+ number of points to create between the start and stop values, where
+ the stop value IS INCLUSIVE.
+
+ If instantiated with an argument of 1, the mesh-grid is open or not
+ fleshed out so that only one-dimension of each returned argument is
+ greater than 1
+
+ Example:
+
+ >>> mgrid = nd_grid()
+ >>> mgrid[0:5,0:5]
+ array([[[0, 0, 0, 0, 0],
+ [1, 1, 1, 1, 1],
+ [2, 2, 2, 2, 2],
+ [3, 3, 3, 3, 3],
+ [4, 4, 4, 4, 4]],
+ [[0, 1, 2, 3, 4],
+ [0, 1, 2, 3, 4],
+ [0, 1, 2, 3, 4],
+ [0, 1, 2, 3, 4],
+ [0, 1, 2, 3, 4]]])
+ >>> mgrid[-1:1:5j]
+ array([-1. , -0.5, 0. , 0.5, 1. ])
+
+ >>> ogrid = nd_grid(1)
+ >>> ogrid[0:5,0:5]
+ [array([[0],[1],[2],[3],[4]]), array([[0, 1, 2, 3, 4]])]
+ """
+ def __init__(self, sparse=False):
+ self.sparse = sparse
+ def __getitem__(self,key):
+ try:
+ size = []
+ typecode = _nx.Int
+ for k in range(len(key)):
+ step = key[k].step
+ start = key[k].start
+ if start is None: start=0
+ if step is None: step=1
+ if type(step) is type(1j):
+ size.append(int(abs(step)))
+ typecode = _nx.Float
+ else:
+ size.append(int((key[k].stop - start)/(step*1.0)))
+ if isinstance(step,types.FloatType) or \
+ isinstance(start, types.FloatType) or \
+ isinstance(key[k].stop, types.FloatType):
+ typecode = _nx.Float
+ if self.sparse:
+ nn = map(lambda x,t: _nx.arange(x,dtype=t),size,(typecode,)*len(size))
+ else:
+ nn = _nx.indices(size,typecode)
+ for k in range(len(size)):
+ step = key[k].step
+ start = key[k].start
+ if start is None: start=0
+ if step is None: step=1
+ if type(step) is type(1j):
+ step = int(abs(step))
+ step = (key[k].stop - start)/float(step-1)
+ nn[k] = (nn[k]*step+start)
+ if self.sparse:
+ slobj = [_nx.NewAxis]*len(size)
+ for k in range(len(size)):
+ slobj[k] = slice(None,None)
+ nn[k] = nn[k][slobj]
+ slobj[k] = _nx.NewAxis
+ return nn
+ except (IndexError, TypeError):
+ step = key.step
+ stop = key.stop
+ start = key.start
+ if start is None: start = 0
+ if type(step) is type(1j):
+ step = abs(step)
+ length = int(step)
+ step = (key.stop-start)/float(step-1)
+ stop = key.stop+step
+ return _nx.arange(0,length,1,_nx.Float)*step + start
+ else:
+ return _nx.arange(start, stop, step)
+
+ def __getslice__(self,i,j):
+ return _nx.arange(i,j)
+
+ def __len__(self):
+ return 0
+
+mgrid = nd_grid()
+ogrid = nd_grid(1)
+
+class concatenator(object):
+ """ Translates slice objects to concatenation along an axis.
+ """
+ def _retval(self, res):
+ if self.matrix:
+ oldndim = res.ndim
+ res = makemat(res)
+ if oldndim == 1 and self.col:
+ res = res.T
+ self.axis=self._axis
+ self.matrix=self._matrix
+ self.col=0
+ return res
+
+ def __init__(self, axis=0, matrix=False):
+ self._axis = axis
+ self._matrix = matrix
+ self.axis = axis
+ self.matrix = matrix
+ self.col = 0
+
+ def __getitem__(self,key):
+ if isinstance(key,types.StringType):
+ frame = sys._getframe().f_back
+ mymat = matrix.bmat(key,frame.f_globals,frame.f_locals)
+ return mymat
+ if type(key) is not types.TupleType:
+ key = (key,)
+ objs = []
+ for k in range(len(key)):
+ if type(key[k]) is types.SliceType:
+ step = key[k].step
+ start = key[k].start
+ stop = key[k].stop
+ if start is None: start = 0
+ if step is None:
+ step = 1
+ if type(step) is type(1j):
+ size = int(abs(step))
+ newobj = function_base.linspace(start, stop, num=size)
+ else:
+ newobj = _nx.arange(start, stop, step)
+ elif type(key[k]) is types.StringType:
+ if (key[k] in 'rc'):
+ self.matrix = True
+ self.col = (key[k] == 'c')
+ continue
+ try:
+ self.axis = int(key[k])
+ continue
+ except:
+ raise ValueError, "Unknown special directive."
+ elif type(key[k]) in ScalarType:
+ newobj = asarray([key[k]])
+ else:
+ newobj = key[k]
+ objs.append(newobj)
+ res = _nx.concatenate(tuple(objs),axis=self.axis)
+ return self._retval(res)
+
+ def __getslice__(self,i,j):
+ res = _nx.arange(i,j)
+ return self._retval(res)
+
+ def __len__(self):
+ return 0
+
+r_=concatenator(0)
+c_=concatenator(-1)
+#row = concatenator(0,1)
+#col = concatenator(-1,1)
+
+
+# A simple nd index iterator over an array:
+
+class ndenumerate(object):
+ def __init__(self, arr):
+ arr = asarray(arr)
+ self.iter = enumerate(arr.flat)
+ self.ashape = arr.shape
+ self.nd = arr.ndim
+ self.factors = [None]*(self.nd-1)
+ val = self.ashape[-1]
+ for i in range(self.nd-1,0,-1):
+ self.factors[i-1] = val
+ val *= self.ashape[i-1]
+
+ def next(self):
+ res = self.iter.next()
+ indxs = [None]*self.nd
+ val = res[0]
+ for i in range(self.nd-1):
+ indxs[i] = val / self.factors[i]
+ val = val % self.factors[i]
+ indxs[self.nd-1] = val
+ return tuple(indxs), res[1]
+
+ def __iter__(self):
+ return self
+
+
+
+# A nicer way to build up index tuples for arrays.
+#
+# You can do all this with slice() plus a few special objects,
+# but there's a lot to remember. This version is simpler because
+# it uses the standard array indexing syntax.
+#
+# Written by Konrad Hinsen <hinsen@cnrs-orleans.fr>
+# last revision: 1999-7-23
+#
+# Cosmetic changes by T. Oliphant 2001
+#
+#
+# This module provides a convenient method for constructing
+# array indices algorithmically. It provides one importable object,
+# 'index_expression'.
+#
+# For any index combination, including slicing and axis insertion,
+# 'a[indices]' is the same as 'a[index_expression[indices]]' for any
+# array 'a'. However, 'index_expression[indices]' can be used anywhere
+# in Python code and returns a tuple of slice objects that can be
+# used in the construction of complex index expressions.
+
+class _index_expression_class(object):
+ maxint = sys.maxint
+
+ def __getitem__(self, item):
+ if type(item) != type(()):
+ return (item,)
+ else:
+ return item
+
+ def __len__(self):
+ return self.maxint
+
+ def __getslice__(self, start, stop):
+ if stop == self.maxint:
+ stop = None
+ return self[start:stop:None]
+
+index_exp = _index_expression_class()
+
+# End contribution from Konrad.
+
diff --git a/numpy/lib/machar.py b/numpy/lib/machar.py
new file mode 100644
index 000000000..e00d112ef
--- /dev/null
+++ b/numpy/lib/machar.py
@@ -0,0 +1,268 @@
+#
+# Machine arithmetics - determine the parameters of the
+# floating-point arithmetic system
+#
+# Author: Pearu Peterson, September 2003
+#
+
+__all__ = ['MachAr']
+
+from numeric import array
+from oldnumeric import any
+
+# Need to speed this up...especially for longfloat
+
+class MachAr(object):
+ """Diagnosing machine parameters.
+
+ The following attributes are available:
+
+ ibeta - radix in which numbers are represented
+ it - number of base-ibeta digits in the floating point mantissa M
+ machep - exponent of the smallest (most negative) power of ibeta that,
+ added to 1.0,
+ gives something different from 1.0
+ eps - floating-point number beta**machep (floating point precision)
+ negep - exponent of the smallest power of ibeta that, substracted
+ from 1.0, gives something different from 1.0
+ epsneg - floating-point number beta**negep
+ iexp - number of bits in the exponent (including its sign and bias)
+ minexp - smallest (most negative) power of ibeta consistent with there
+ being no leading zeros in the mantissa
+ xmin - floating point number beta**minexp (the smallest (in
+ magnitude) usable floating value)
+ maxexp - smallest (positive) power of ibeta that causes overflow
+ xmax - (1-epsneg)* beta**maxexp (the largest (in magnitude)
+ usable floating value)
+ irnd - in range(6), information on what kind of rounding is done
+ in addition, and on how underflow is handled
+ ngrd - number of 'guard digits' used when truncating the product
+ of two mantissas to fit the representation
+
+ epsilon - same as eps
+ tiny - same as xmin
+ huge - same as xmax
+ precision - int(-log10(eps))
+ resolution - 10**(-precision)
+
+ Reference:
+ Numerical Recipies.
+ """
+ def __init__(self, float_conv=float,int_conv=int,
+ float_to_float=float,
+ float_to_str = lambda v:'%24.16e' % v,
+ title = 'Python floating point number'):
+ """
+ float_conv - convert integer to float (array)
+ int_conv - convert float (array) to integer
+ float_to_float - convert float array to float
+ float_to_str - convert array float to str
+ title - description of used floating point numbers
+ """
+ one = float_conv(1)
+ two = one + one
+ zero = one - one
+
+ # Do we really need to do this? Aren't they 2 and 2.0?
+ # Determine ibeta and beta
+ a = one
+ while 1:
+ a = a + a
+ temp = a + one
+ temp1 = temp - a
+ if any(temp1 - one != zero):
+ break
+ b = one
+ while 1:
+ b = b + b
+ temp = a + b
+ itemp = int_conv(temp-a)
+ if any(itemp != 0):
+ break
+ ibeta = itemp
+ beta = float_conv(ibeta)
+
+ # Determine it and irnd
+ it = -1
+ b = one
+ while 1:
+ it = it + 1
+ b = b * beta
+ temp = b + one
+ temp1 = temp - b
+ if any(temp1 - one != zero):
+ break
+
+ betah = beta / two
+ a = one
+ while 1:
+ a = a + a
+ temp = a + one
+ temp1 = temp - a
+ if any(temp1 - one != zero):
+ break
+ temp = a + betah
+ irnd = 0
+ if any(temp-a != zero):
+ irnd = 1
+ tempa = a + beta
+ temp = tempa + betah
+ if irnd==0 and any(temp-tempa != zero):
+ irnd = 2
+
+ # Determine negep and epsneg
+ negep = it + 3
+ betain = one / beta
+ a = one
+ for i in range(negep):
+ a = a * betain
+ b = a
+ while 1:
+ temp = one - a
+ if any(temp-one != zero):
+ break
+ a = a * beta
+ negep = negep - 1
+ # Prevent infinite loop on PPC with gcc 4.0:
+ if negep < 0:
+ raise RuntimeError, "could not determine machine tolerance " \
+ "for 'negep'"
+ negep = -negep
+ epsneg = a
+
+ # Determine machep and eps
+ machep = - it - 3
+ a = b
+
+ while 1:
+ temp = one + a
+ if any(temp-one != zero):
+ break
+ a = a * beta
+ machep = machep + 1
+ eps = a
+
+ # Determine ngrd
+ ngrd = 0
+ temp = one + eps
+ if irnd==0 and any(temp*one - one != zero):
+ ngrd = 1
+
+ # Determine iexp
+ i = 0
+ k = 1
+ z = betain
+ t = one + eps
+ nxres = 0
+ while 1:
+ y = z
+ z = y*y
+ a = z*one # Check here for underflow
+ temp = z*t
+ if any(a+a == zero) or any(abs(z)>=y):
+ break
+ temp1 = temp * betain
+ if any(temp1*beta == z):
+ break
+ i = i + 1
+ k = k + k
+ if ibeta != 10:
+ iexp = i + 1
+ mx = k + k
+ else:
+ iexp = 2
+ iz = ibeta
+ while k >= iz:
+ iz = iz * ibeta
+ iexp = iexp + 1
+ mx = iz + iz - 1
+
+ # Determine minexp and xmin
+ while 1:
+ xmin = y
+ y = y * betain
+ a = y * one
+ temp = y * t
+ if any(a+a != zero) and any(abs(y) < xmin):
+ k = k + 1
+ temp1 = temp * betain
+ if any(temp1*beta == y) and any(temp != y):
+ nxres = 3
+ xmin = y
+ break
+ else:
+ break
+ minexp = -k
+
+ # Determine maxexp, xmax
+ if mx <= k + k - 3 and ibeta != 10:
+ mx = mx + mx
+ iexp = iexp + 1
+ maxexp = mx + minexp
+ irnd = irnd + nxres
+ if irnd >= 2:
+ maxexp = maxexp - 2
+ i = maxexp + minexp
+ if ibeta == 2 and not i:
+ maxexp = maxexp - 1
+ if i > 20:
+ maxexp = maxexp - 1
+ if any(a != y):
+ maxexp = maxexp - 2
+ xmax = one - epsneg
+ if any(xmax*one != xmax):
+ xmax = one - beta*epsneg
+ xmax = xmax / (xmin*beta*beta*beta)
+ i = maxexp + minexp + 3
+ for j in range(i):
+ if ibeta==2:
+ xmax = xmax + xmax
+ else:
+ xmax = xmax * beta
+
+ self.ibeta = ibeta
+ self.it = it
+ self.negep = negep
+ self.epsneg = float_to_float(epsneg)
+ self._str_epsneg = float_to_str(epsneg)
+ self.machep = machep
+ self.eps = float_to_float(eps)
+ self._str_eps = float_to_str(eps)
+ self.ngrd = ngrd
+ self.iexp = iexp
+ self.minexp = minexp
+ self.xmin = float_to_float(xmin)
+ self._str_xmin = float_to_str(xmin)
+ self.maxexp = maxexp
+ self.xmax = float_to_float(xmax)
+ self._str_xmax = float_to_str(xmax)
+ self.irnd = irnd
+
+ self.title = title
+ # Commonly used parameters
+ self.epsilon = self.eps
+ self.tiny = self.xmin
+ self.huge = self.xmax
+
+ import math
+ self.precision = int(-math.log10(float_to_float(self.eps)))
+ ten = two + two + two + two + two
+ resolution = ten ** (-self.precision)
+ self.resolution = float_to_float(resolution)
+ self._str_resolution = float_to_str(resolution)
+
+ def __str__(self):
+ return '''\
+Machine parameters for %(title)s
+---------------------------------------------------------------------
+ibeta=%(ibeta)s it=%(it)s iexp=%(iexp)s ngrd=%(ngrd)s irnd=%(irnd)s
+machep=%(machep)s eps=%(_str_eps)s (beta**machep == epsilon)
+negep =%(negep)s epsneg=%(_str_epsneg)s (beta**epsneg)
+minexp=%(minexp)s xmin=%(_str_xmin)s (beta**minexp == tiny)
+maxexp=%(maxexp)s xmax=%(_str_xmax)s ((1-epsneg)*beta**maxexp == huge)
+---------------------------------------------------------------------
+''' % self.__dict__
+
+
+if __name__ == '__main__':
+ print MachAr()
diff --git a/numpy/lib/mlab.py b/numpy/lib/mlab.py
new file mode 100644
index 000000000..749600d9b
--- /dev/null
+++ b/numpy/lib/mlab.py
@@ -0,0 +1,14 @@
+# This module is for compatibility only. All functions are defined elsewhere.
+
+from numeric import *
+
+from twodim_base import eye, tri, diag, fliplr, flipud, rot90, tril, triu
+from oldnumeric import amax as max
+from oldnumeric import amin as min
+from function_base import msort, median, trapz, diff, cov, corrcoef, kaiser, blackman, \
+ bartlett, hanning, hamming, sinc, angle
+from oldnumeric import cumsum, ptp, mean, std, prod, cumprod, squeeze
+from polynomial import roots
+
+from scipy.random import rand, randn
+from scipy.corelinalg import eig, svd
diff --git a/numpy/lib/polynomial.py b/numpy/lib/polynomial.py
new file mode 100644
index 000000000..df7013bab
--- /dev/null
+++ b/numpy/lib/polynomial.py
@@ -0,0 +1,554 @@
+"""
+Functions to operate on polynomials.
+"""
+
+__all__ = ['poly', 'roots', 'polyint', 'polyder', 'polyadd',
+ 'polysub', 'polymul', 'polydiv', 'polyval', 'poly1d',
+ 'polyfit']
+
+import re
+import numeric as NX
+
+from type_check import isscalar
+from twodim_base import diag, vander
+from shape_base import hstack, atleast_1d
+from function_base import trim_zeros, sort_complex
+eigvals = None
+lstsq = None
+
+def get_linalg_funcs():
+ "Look for linear algebra functions in scipy"
+ global eigvals, lstsq
+ from scipy.corelinalg import eigvals, lstsq
+ return
+
+def _eigvals(arg):
+ "Return the eigenvalues of the argument"
+ try:
+ return eigvals(arg)
+ except TypeError:
+ get_linalg_funcs()
+ return eigvals(arg)
+
+def _lstsq(X, y):
+ "Do least squares on the arguments"
+ try:
+ return lstsq(X, y)
+ except TypeError:
+ get_linalg_funcs()
+ return lstsq(X, y)
+
+def poly(seq_of_zeros):
+ """ Return a sequence representing a polynomial given a sequence of roots.
+
+ If the input is a matrix, return the characteristic polynomial.
+
+ Example:
+
+ >>> b = roots([1,3,1,5,6])
+ >>> poly(b)
+ array([1., 3., 1., 5., 6.])
+ """
+ seq_of_zeros = atleast_1d(seq_of_zeros)
+ sh = seq_of_zeros.shape
+ if len(sh) == 2 and sh[0] == sh[1]:
+ seq_of_zeros = _eigvals(seq_of_zeros)
+ elif len(sh) ==1:
+ pass
+ else:
+ raise ValueError, "input must be 1d or square 2d array."
+
+ if len(seq_of_zeros) == 0:
+ return 1.0
+
+ a = [1]
+ for k in range(len(seq_of_zeros)):
+ a = NX.convolve(a, [1, -seq_of_zeros[k]], mode='full')
+
+ if issubclass(a.dtype, NX.complexfloating):
+ # if complex roots are all complex conjugates, the roots are real.
+ roots = NX.asarray(seq_of_zeros, complex)
+ pos_roots = sort_complex(NX.compress(roots.imag > 0, roots))
+ neg_roots = NX.conjugate(sort_complex(
+ NX.compress(roots.imag < 0,roots)))
+ if (len(pos_roots) == len(neg_roots) and
+ NX.alltrue(neg_roots == pos_roots)):
+ a = a.real.copy()
+
+ return a
+
+def roots(p):
+ """ Return the roots of the polynomial coefficients in p.
+
+ The values in the rank-1 array p are coefficients of a polynomial.
+ If the length of p is n+1 then the polynomial is
+ p[0] * x**n + p[1] * x**(n-1) + ... + p[n-1]*x + p[n]
+ """
+ # If input is scalar, this makes it an array
+ p = atleast_1d(p)
+ if len(p.shape) != 1:
+ raise ValueError,"Input must be a rank-1 array."
+
+ # find non-zero array entries
+ non_zero = NX.nonzero(NX.ravel(p))
+
+ # find the number of trailing zeros -- this is the number of roots at 0.
+ trailing_zeros = len(p) - non_zero[-1] - 1
+
+ # strip leading and trailing zeros
+ p = p[int(non_zero[0]):int(non_zero[-1])+1]
+
+ # casting: if incoming array isn't floating point, make it floating point.
+ if not issubclass(p.dtype, (NX.floating, NX.complexfloating)):
+ p = p.astype(float)
+
+ N = len(p)
+ if N > 1:
+ # build companion matrix and find its eigenvalues (the roots)
+ A = diag(NX.ones((N-2,), p.dtype), -1)
+ A[0, :] = -p[1:] / p[0]
+ roots = _eigvals(A)
+ else:
+ return NX.array([])
+
+ # tack any zeros onto the back of the array
+ roots = hstack((roots, NX.zeros(trailing_zeros, roots.dtype)))
+ return roots
+
+def polyint(p, m=1, k=None):
+ """Return the mth analytical integral of the polynomial p.
+
+ If k is None, then zero-valued constants of integration are used.
+ otherwise, k should be a list of length m (or a scalar if m=1) to
+ represent the constants of integration to use for each integration
+ (starting with k[0])
+ """
+ m = int(m)
+ if m < 0:
+ raise ValueError, "Order of integral must be positive (see polyder)"
+ if k is None:
+ k = NX.zeros(m, float)
+ k = atleast_1d(k)
+ if len(k) == 1 and m > 1:
+ k = k[0]*NX.ones(m, float)
+ if len(k) < m:
+ raise ValueError, \
+ "k must be a scalar or a rank-1 array of length 1 or >m."
+ if m == 0:
+ return p
+ else:
+ truepoly = isinstance(p, poly1d)
+ p = NX.asarray(p)
+ y = NX.zeros(len(p)+1, float)
+ y[:-1] = p*1.0/NX.arange(len(p), 0, -1)
+ y[-1] = k[0]
+ val = polyint(y, m-1, k=k[1:])
+ if truepoly:
+ val = poly1d(val)
+ return val
+
+def polyder(p, m=1):
+ """Return the mth derivative of the polynomial p.
+ """
+ m = int(m)
+ truepoly = isinstance(p, poly1d)
+ p = NX.asarray(p)
+ n = len(p)-1
+ y = p[:-1] * NX.arange(n, 0, -1)
+ if m < 0:
+ raise ValueError, "Order of derivative must be positive (see polyint)"
+ if m == 0:
+ return p
+ else:
+ val = polyder(y, m-1)
+ if truepoly:
+ val = poly1d(val)
+ return val
+
+def polyfit(x, y, N):
+ """
+
+ Do a best fit polynomial of order N of y to x. Return value is a
+ vector of polynomial coefficients [pk ... p1 p0]. Eg, for N=2
+
+ p2*x0^2 + p1*x0 + p0 = y1
+ p2*x1^2 + p1*x1 + p0 = y1
+ p2*x2^2 + p1*x2 + p0 = y2
+ .....
+ p2*xk^2 + p1*xk + p0 = yk
+
+
+ Method: if X is a the Vandermonde Matrix computed from x (see
+ http://mathworld.wolfram.com/VandermondeMatrix.html), then the
+ polynomial least squares solution is given by the 'p' in
+
+ X*p = y
+
+ where X is a len(x) x N+1 matrix, p is a N+1 length vector, and y
+ is a len(x) x 1 vector
+
+ This equation can be solved as
+
+ p = (XT*X)^-1 * XT * y
+
+ where XT is the transpose of X and -1 denotes the inverse.
+
+ For more info, see
+ http://mathworld.wolfram.com/LeastSquaresFittingPolynomial.html,
+ but note that the k's and n's in the superscripts and subscripts
+ on that page. The linear algebra is correct, however.
+
+ See also polyval
+
+ """
+ x = NX.asarray(x)+0.
+ y = NX.asarray(y)+0.
+ y = NX.reshape(y, (len(y), 1))
+ X = vander(x, N+1)
+ c, resids, rank, s = _lstsq(X, y)
+ c.shape = (N+1,)
+ return c
+
+
+
+def polyval(p, x):
+ """Evaluate the polynomial p at x. If x is a polynomial then composition.
+
+ Description:
+
+ If p is of length N, this function returns the value:
+ p[0]*(x**N-1) + p[1]*(x**N-2) + ... + p[N-2]*x + p[N-1]
+
+ x can be a sequence and p(x) will be returned for all elements of x.
+ or x can be another polynomial and the composite polynomial p(x) will be
+ returned.
+
+ Notice: This can produce inaccurate results for polynomials with
+ significant variability. Use carefully.
+ """
+ p = NX.asarray(p)
+ if isinstance(x, poly1d):
+ y = 0
+ else:
+ x = NX.asarray(x)
+ y = NX.zeros_like(x)
+ for i in range(len(p)):
+ y = x * y + p[i]
+ return y
+
+def polyadd(a1, a2):
+ """Adds two polynomials represented as sequences
+ """
+ truepoly = (isinstance(a1, poly1d) or isinstance(a2, poly1d))
+ a1 = atleast_1d(a1)
+ a2 = atleast_1d(a2)
+ diff = len(a2) - len(a1)
+ if diff == 0:
+ return a1 + a2
+ elif diff > 0:
+ zr = NX.zeros(diff, a1.dtype)
+ val = NX.concatenate((zr, a1)) + a2
+ else:
+ zr = NX.zeros(abs(diff), a2.dtype)
+ val = a1 + NX.concatenate((zr, a2))
+ if truepoly:
+ val = poly1d(val)
+ return val
+
+def polysub(a1, a2):
+ """Subtracts two polynomials represented as sequences
+ """
+ truepoly = (isinstance(a1, poly1d) or isinstance(a2, poly1d))
+ a1 = atleast_1d(a1)
+ a2 = atleast_1d(a2)
+ diff = len(a2) - len(a1)
+ if diff == 0:
+ return a1 - a2
+ elif diff > 0:
+ zr = NX.zeros(diff, a1)
+ val = NX.concatenate((zr, a1)) - a2
+ else:
+ zr = NX.zeros(abs(diff), a2)
+ val = a1 - NX.concatenate((zr, a2))
+ if truepoly:
+ val = poly1d(val)
+ return val
+
+
+def polymul(a1, a2):
+ """Multiplies two polynomials represented as sequences.
+ """
+ truepoly = (isinstance(a1, poly1d) or isinstance(a2, poly1d))
+ val = NX.convolve(a1, a2)
+ if truepoly:
+ val = poly1d(val)
+ return val
+
+
+def deconvolve(signal, divisor):
+ """Deconvolves divisor out of signal. Requires scipy.signal library
+ """
+ import scipy.signal
+ num = atleast_1d(signal)
+ den = atleast_1d(divisor)
+ N = len(num)
+ D = len(den)
+ if D > N:
+ quot = [];
+ rem = num;
+ else:
+ input = NX.ones(N-D+1, float)
+ input[1:] = 0
+ quot = scipy.signal.lfilter(num, den, input)
+ rem = num - NX.convolve(den, quot, mode='full')
+ return quot, rem
+
+def polydiv(u, v):
+ """Computes q and r polynomials so that u(s) = q(s)*v(s) + r(s)
+ and deg r < deg v.
+ """
+ truepoly = (isinstance(u, poly1d) or isinstance(u, poly1d))
+ u = atleast_1d(u)
+ v = atleast_1d(v)
+ m = len(u) - 1
+ n = len(v) - 1
+ scale = 1. / v[0]
+ q = NX.zeros((m-n+1,), float)
+ r = u.copy()
+ for k in range(0, m-n+1):
+ d = scale * r[k]
+ q[k] = d
+ r[k:k+n+1] -= d*v
+ while NX.allclose(r[0], 0, rtol=1e-14) and (r.shape[-1] > 1):
+ r = r[1:]
+ if truepoly:
+ q = poly1d(q)
+ r = poly1d(r)
+ return q, r
+
+_poly_mat = re.compile(r"[*][*]([0-9]*)")
+def _raise_power(astr, wrap=70):
+ n = 0
+ line1 = ''
+ line2 = ''
+ output = ' '
+ while 1:
+ mat = _poly_mat.search(astr, n)
+ if mat is None:
+ break
+ span = mat.span()
+ power = mat.groups()[0]
+ partstr = astr[n:span[0]]
+ n = span[1]
+ toadd2 = partstr + ' '*(len(power)-1)
+ toadd1 = ' '*(len(partstr)-1) + power
+ if ((len(line2)+len(toadd2) > wrap) or \
+ (len(line1)+len(toadd1) > wrap)):
+ output += line1 + "\n" + line2 + "\n "
+ line1 = toadd1
+ line2 = toadd2
+ else:
+ line2 += partstr + ' '*(len(power)-1)
+ line1 += ' '*(len(partstr)-1) + power
+ output += line1 + "\n" + line2
+ return output + astr[n:]
+
+
+class poly1d(object):
+ """A one-dimensional polynomial class.
+
+ p = poly1d([1,2,3]) constructs the polynomial x**2 + 2 x + 3
+
+ p(0.5) evaluates the polynomial at the location
+ p.r is a list of roots
+ p.c is the coefficient array [1,2,3]
+ p.order is the polynomial order (after leading zeros in p.c are removed)
+ p[k] is the coefficient on the kth power of x (backwards from
+ sequencing the coefficient array.
+
+ polynomials can be added, substracted, multplied and divided (returns
+ quotient and remainder).
+ asarray(p) will also give the coefficient array, so polynomials can
+ be used in all functions that accept arrays.
+
+ p = poly1d([1,2,3], variable='lambda') will use lambda in the
+ string representation of p.
+ """
+ def __init__(self, c_or_r, r=0, variable=None):
+ if isinstance(c_or_r, poly1d):
+ for key in c_or_r.__dict__.keys():
+ self.__dict__[key] = c_or_r.__dict__[key]
+ if variable is not None:
+ self.__dict__['variable'] = variable
+ return
+ if r:
+ c_or_r = poly(c_or_r)
+ c_or_r = atleast_1d(c_or_r)
+ if len(c_or_r.shape) > 1:
+ raise ValueError, "Polynomial must be 1d only."
+ c_or_r = trim_zeros(c_or_r, trim='f')
+ if len(c_or_r) == 0:
+ c_or_r = NX.array([0.])
+ self.__dict__['coeffs'] = c_or_r
+ self.__dict__['order'] = len(c_or_r) - 1
+ if variable is None:
+ variable = 'x'
+ self.__dict__['variable'] = variable
+
+ def __array__(self, t=None):
+ if t:
+ return NX.asarray(self.coeffs, t)
+ else:
+ return NX.asarray(self.coeffs)
+
+ def __repr__(self):
+ vals = repr(self.coeffs)
+ vals = vals[6:-1]
+ return "poly1d(%s)" % vals
+
+ def __len__(self):
+ return self.order
+
+ def __str__(self):
+ N = self.order
+ thestr = "0"
+ var = self.variable
+ for k in range(len(self.coeffs)):
+ coefstr ='%.4g' % abs(self.coeffs[k])
+ if coefstr[-4:] == '0000':
+ coefstr = coefstr[:-5]
+ power = (N-k)
+ if power == 0:
+ if coefstr != '0':
+ newstr = '%s' % (coefstr,)
+ else:
+ if k == 0:
+ newstr = '0'
+ else:
+ newstr = ''
+ elif power == 1:
+ if coefstr == '0':
+ newstr = ''
+ elif coefstr == 'b':
+ newstr = var
+ else:
+ newstr = '%s %s' % (coefstr, var)
+ else:
+ if coefstr == '0':
+ newstr = ''
+ elif coefstr == 'b':
+ newstr = '%s**%d' % (var, power,)
+ else:
+ newstr = '%s %s**%d' % (coefstr, var, power)
+
+ if k > 0:
+ if newstr != '':
+ if self.coeffs[k] < 0:
+ thestr = "%s - %s" % (thestr, newstr)
+ else:
+ thestr = "%s + %s" % (thestr, newstr)
+ elif (k == 0) and (newstr != '') and (self.coeffs[k] < 0):
+ thestr = "-%s" % (newstr,)
+ else:
+ thestr = newstr
+ return _raise_power(thestr)
+
+
+ def __call__(self, val):
+ return polyval(self.coeffs, val)
+
+ def __mul__(self, other):
+ if isscalar(other):
+ return poly1d(self.coeffs * other)
+ else:
+ other = poly1d(other)
+ return poly1d(polymul(self.coeffs, other.coeffs))
+
+ def __rmul__(self, other):
+ if isscalar(other):
+ return poly1d(other * self.coeffs)
+ else:
+ other = poly1d(other)
+ return poly1d(polymul(self.coeffs, other.coeffs))
+
+ def __add__(self, other):
+ other = poly1d(other)
+ return poly1d(polyadd(self.coeffs, other.coeffs))
+
+ def __radd__(self, other):
+ other = poly1d(other)
+ return poly1d(polyadd(self.coeffs, other.coeffs))
+
+ def __pow__(self, val):
+ if not isscalar(val) or int(val) != val or val < 0:
+ raise ValueError, "Power to non-negative integers only."
+ res = [1]
+ for k in range(val):
+ res = polymul(self.coeffs, res)
+ return poly1d(res)
+
+ def __sub__(self, other):
+ other = poly1d(other)
+ return poly1d(polysub(self.coeffs, other.coeffs))
+
+ def __rsub__(self, other):
+ other = poly1d(other)
+ return poly1d(polysub(other.coeffs, self.coeffs))
+
+ def __div__(self, other):
+ if isscalar(other):
+ return poly1d(self.coeffs/other)
+ else:
+ other = poly1d(other)
+ return polydiv(self, other)
+
+ def __rdiv__(self, other):
+ if isscalar(other):
+ return poly1d(other/self.coeffs)
+ else:
+ other = poly1d(other)
+ return polydiv(other, self)
+
+ def __setattr__(self, key, val):
+ raise ValueError, "Attributes cannot be changed this way."
+
+ def __getattr__(self, key):
+ if key in ['r', 'roots']:
+ return roots(self.coeffs)
+ elif key in ['c','coef','coefficients']:
+ return self.coeffs
+ elif key in ['o']:
+ return self.order
+ else:
+ return self.__dict__[key]
+
+ def __getitem__(self, val):
+ ind = self.order - val
+ if val > self.order:
+ return 0
+ if val < 0:
+ return 0
+ return self.coeffs[ind]
+
+ def __setitem__(self, key, val):
+ ind = self.order - key
+ if key < 0:
+ raise ValueError, "Does not support negative powers."
+ if key > self.order:
+ zr = NX.zeros(key-self.order, self.coeffs.dtype)
+ self.__dict__['coeffs'] = NX.concatenate((zr, self.coeffs))
+ self.__dict__['order'] = key
+ ind = 0
+ self.__dict__['coeffs'][ind] = val
+ return
+
+ def integ(self, m=1, k=0):
+ """Return the mth analytical integral of this polynomial.
+ See the documentation for polyint.
+ """
+ return poly1d(polyint(self.coeffs, m=m, k=k))
+
+ def deriv(self, m=1):
+ """Return the mth derivative of this polynomial.
+ """
+ return poly1d(polyder(self.coeffs, m=m))
diff --git a/numpy/lib/scimath.py b/numpy/lib/scimath.py
new file mode 100644
index 000000000..4202fa640
--- /dev/null
+++ b/numpy/lib/scimath.py
@@ -0,0 +1,77 @@
+"""
+Wrapper functions to more user-friendly calling of certain math functions
+whose output data-type is different than the input data-type in certain domains of the input.
+"""
+
+__all__ = ['sqrt', 'log', 'log2', 'logn','log10', 'power', 'arccos',
+ 'arcsin', 'arctanh']
+
+import numeric as nx
+from numeric import *
+
+from type_check import isreal, asscalar
+
+__all__.extend([key for key in dir(nx.umath)
+ if key[0] != '_' and key not in __all__])
+
+_ln2 = log(2.0)
+
+def _tocomplex(arr):
+ if isinstance(arr.dtype, (nx.single, nx.byte, nx.short, nx.ubyte,
+ nx.ushort)):
+ return arr.astype(nx.csingle)
+ else:
+ return arr.astype(nx.cdouble)
+
+def _fix_real_lt_zero(x):
+ x = asarray(x)
+ if any(isreal(x) & (x<0)):
+ x = _tocomplex(x)
+ return asscalar(x)
+
+def _fix_real_abs_gt_1(x):
+ x = asarray(x)
+ if any(isreal(x) & (abs(x)>1)):
+ x = _tocomplex(x)
+ return x
+
+def sqrt(x):
+ x = _fix_real_lt_zero(x)
+ return nx.sqrt(x)
+
+def log(x):
+ x = _fix_real_lt_zero(x)
+ return nx.log(x)
+
+def log10(x):
+ x = _fix_real_lt_zero(x)
+ return nx.log10(x)
+
+def logn(n, x):
+ """ Take log base n of x.
+ """
+ x = _fix_real_lt_zero(x)
+ n = _fix_real_lt_zero(n)
+ return log(x)/log(n)
+
+def log2(x):
+ """ Take log base 2 of x.
+ """
+ x = _fix_real_lt_zero(x)
+ return log(x)/_ln2
+
+def power(x, p):
+ x = _fix_real_lt_zero(x)
+ return nx.power(x, p)
+
+def arccos(x):
+ x = _fix_real_abs_gt_1(x)
+ return arccos(x)
+
+def arcsin(x):
+ x = _fix_real_abs_gt_1(x)
+ return arcsin(x)
+
+def arctanh(x):
+ x = _fix_real_abs_gt_1(x)
+ return arctanh(x)
diff --git a/numpy/lib/shape_base.py b/numpy/lib/shape_base.py
new file mode 100644
index 000000000..8d66b41d1
--- /dev/null
+++ b/numpy/lib/shape_base.py
@@ -0,0 +1,539 @@
+__all__ = ['atleast_1d','atleast_2d','atleast_3d','vstack','hstack',
+ 'column_stack','dstack','array_split','split','hsplit',
+ 'vsplit','dsplit','apply_over_axes','expand_dims',
+ 'apply_along_axis']
+
+import numeric as _nx
+from numeric import *
+from oldnumeric import product
+
+from type_check import isscalar
+
+def apply_along_axis(func1d,axis,arr,*args):
+ """ Execute func1d(arr[i],*args) where func1d takes 1-D arrays
+ and arr is an N-d array. i varies so as to apply the function
+ along the given axis for each 1-d subarray in arr.
+ """
+ arr = asarray(arr)
+ nd = arr.ndim
+ if axis < 0:
+ axis += nd
+ if (axis >= nd):
+ raise ValueError("axis must be less than arr.ndim; axis=%d, rank=%d."
+ % (axis,nd))
+ ind = [0]*(nd-1)
+ i = zeros(nd,'O')
+ indlist = range(nd)
+ indlist.remove(axis)
+ i[axis] = slice(None,None)
+ outshape = asarray(arr.shape).take(indlist)
+ i.put(ind, indlist)
+ res = func1d(arr[tuple(i.tolist())],*args)
+ # if res is a number, then we have a smaller output array
+ if isscalar(res):
+ outarr = zeros(outshape,asarray(res).dtypechar)
+ outarr[ind] = res
+ Ntot = product(outshape)
+ k = 1
+ while k < Ntot:
+ # increment the index
+ ind[-1] += 1
+ n = -1
+ while (ind[n] >= outshape[n]) and (n > (1-nd)):
+ ind[n-1] += 1
+ ind[n] = 0
+ n -= 1
+ i.put(ind,indlist)
+ res = func1d(arr[tuple(i.tolist())],*args)
+ outarr[ind] = res
+ k += 1
+ return outarr
+ else:
+ Ntot = product(outshape)
+ holdshape = outshape
+ outshape = list(arr.shape)
+ outshape[axis] = len(res)
+ outarr = zeros(outshape,asarray(res).dtypechar)
+ outarr[tuple(i.tolist())] = res
+ k = 1
+ while k < Ntot:
+ # increment the index
+ ind[-1] += 1
+ n = -1
+ while (ind[n] >= holdshape[n]) and (n > (1-nd)):
+ ind[n-1] += 1
+ ind[n] = 0
+ n -= 1
+ i.put(ind, indlist)
+ res = func1d(arr[tuple(i.tolist())],*args)
+ outarr[tuple(i.tolist())] = res
+ k += 1
+ return outarr
+
+
+def apply_over_axes(func, a, axes):
+ """Apply a function repeatedly over multiple axes, keeping the same shape
+ for the resulting array.
+
+ func is called as res = func(a, axis). The result is assumed
+ to be either the same shape as a or have one less dimension.
+ This call is repeated for each axis in the axes sequence.
+ """
+ val = asarray(a)
+ N = a.ndim
+ if array(axes).ndim == 0:
+ axes = (axes,)
+ for axis in axes:
+ if axis < 0: axis = N + axis
+ args = (val, axis)
+ res = func(*args)
+ if res.ndim == val.ndim:
+ val = res
+ else:
+ res = expand_dims(res,axis)
+ if res.ndim == val.ndim:
+ val = res
+ else:
+ raise ValueError, "function is not returning"\
+ " an array of correct shape"
+ return val
+
+def expand_dims(a, axis):
+ """Expand the shape of a by including newaxis before given axis.
+ """
+ a = asarray(a)
+ shape = a.shape
+ if axis < 0:
+ axis = axis + len(shape) + 1
+ return a.reshape(shape[:axis] + (1,) + shape[axis:])
+
+
+def atleast_1d(*arys):
+ """ Force a sequence of arrays to each be at least 1D.
+
+ Description:
+ Force an array to be at least 1D. If an array is 0D, the
+ array is converted to a single row of values. Otherwise,
+ the array is unaltered.
+ Arguments:
+ *arys -- arrays to be converted to 1 or more dimensional array.
+ Returns:
+ input array converted to at least 1D array.
+ """
+ res = []
+ for ary in arys:
+ ary = asarray(ary)
+ if len(ary.shape) == 0:
+ ary = ary.reshape(1)
+ res.append(ary)
+ if len(res) == 1:
+ return res[0]
+ else:
+ return res
+
+def atleast_2d(*arys):
+ """ Force a sequence of arrays to each be at least 2D.
+
+ Description:
+ Force an array to each be at least 2D. If the array
+ is 0D or 1D, the array is converted to a single
+ row of values. Otherwise, the array is unaltered.
+ Arguments:
+ arys -- arrays to be converted to 2 or more dimensional array.
+ Returns:
+ input array converted to at least 2D array.
+ """
+ res = []
+ for ary in arys:
+ ary = asarray(ary)
+ if len(ary.shape) == 0:
+ result = ary.reshape(1,1)
+ elif len(ary.shape) == 1:
+ result = ary[newaxis,:]
+ else:
+ result = ary
+ res.append(result)
+ if len(res) == 1:
+ return res[0]
+ else:
+ return res
+
+def atleast_3d(*arys):
+ """ Force a sequence of arrays to each be at least 3D.
+
+ Description:
+ Force an array each be at least 3D. If the array is 0D or 1D,
+ the array is converted to a single 1xNx1 array of values where
+ N is the orginal length of the array. If the array is 2D, the
+ array is converted to a single MxNx1 array of values where MxN
+ is the orginal shape of the array. Otherwise, the array is
+ unaltered.
+ Arguments:
+ arys -- arrays to be converted to 3 or more dimensional array.
+ Returns:
+ input array converted to at least 3D array.
+ """
+ res = []
+ for ary in arys:
+ ary = asarray(ary)
+ if len(ary.shape) == 0:
+ result = ary.reshape(1,1,1)
+ elif len(ary.shape) == 1:
+ result = ary[newaxis,:,newaxis]
+ elif len(ary.shape) == 2:
+ result = ary[:,:,newaxis]
+ else:
+ result = ary
+ res.append(result)
+ if len(res) == 1:
+ return res[0]
+ else:
+ return res
+
+
+def vstack(tup):
+ """ Stack arrays in sequence vertically (row wise)
+
+ Description:
+ Take a sequence of arrays and stack them veritcally
+ to make a single array. All arrays in the sequence
+ must have the same shape along all but the first axis.
+ vstack will rebuild arrays divided by vsplit.
+ Arguments:
+ tup -- sequence of arrays. All arrays must have the same
+ shape.
+ Examples:
+ >>> import scipy
+ >>> a = array((1,2,3))
+ >>> b = array((2,3,4))
+ >>> scipy.vstack((a,b))
+ array([[1, 2, 3],
+ [2, 3, 4]])
+ >>> a = array([[1],[2],[3]])
+ >>> b = array([[2],[3],[4]])
+ >>> scipy.vstack((a,b))
+ array([[1],
+ [2],
+ [3],
+ [2],
+ [3],
+ [4]])
+
+ """
+ return _nx.concatenate(map(atleast_2d,tup),0)
+
+def hstack(tup):
+ """ Stack arrays in sequence horizontally (column wise)
+
+ Description:
+ Take a sequence of arrays and stack them horizontally
+ to make a single array. All arrays in the sequence
+ must have the same shape along all but the second axis.
+ hstack will rebuild arrays divided by hsplit.
+ Arguments:
+ tup -- sequence of arrays. All arrays must have the same
+ shape.
+ Examples:
+ >>> import scipy
+ >>> a = array((1,2,3))
+ >>> b = array((2,3,4))
+ >>> scipy.hstack((a,b))
+ array([1, 2, 3, 2, 3, 4])
+ >>> a = array([[1],[2],[3]])
+ >>> b = array([[2],[3],[4]])
+ >>> scipy.hstack((a,b))
+ array([[1, 2],
+ [2, 3],
+ [3, 4]])
+
+ """
+ return _nx.concatenate(map(atleast_1d,tup),1)
+
+def column_stack(tup):
+ """ Stack 1D arrays as columns into a 2D array
+
+ Description:
+ Take a sequence of 1D arrays and stack them as columns
+ to make a single 2D array. All arrays in the sequence
+ must have the same length.
+ Arguments:
+ tup -- sequence of 1D arrays. All arrays must have the same
+ length.
+ Examples:
+ >>> import scipy
+ >>> a = array((1,2,3))
+ >>> b = array((2,3,4))
+ >>> scipy.column_stack((a,b))
+ array([[1, 2],
+ [2, 3],
+ [3, 4]])
+
+ """
+ arrays = map(_nx.transpose,map(atleast_2d,tup))
+ return _nx.concatenate(arrays,1)
+
+def dstack(tup):
+ """ Stack arrays in sequence depth wise (along third dimension)
+
+ Description:
+ Take a sequence of arrays and stack them along the third axis.
+ All arrays in the sequence must have the same shape along all
+ but the third axis. This is a simple way to stack 2D arrays
+ (images) into a single 3D array for processing.
+ dstack will rebuild arrays divided by dsplit.
+ Arguments:
+ tup -- sequence of arrays. All arrays must have the same
+ shape.
+ Examples:
+ >>> import scipy
+ >>> a = array((1,2,3))
+ >>> b = array((2,3,4))
+ >>> scipy.dstack((a,b))
+ array([ [[1, 2],
+ [2, 3],
+ [3, 4]]])
+ >>> a = array([[1],[2],[3]])
+ >>> b = array([[2],[3],[4]])
+ >>> scipy.dstack((a,b))
+ array([[ [1, 2]],
+ [ [2, 3]],
+ [ [3, 4]]])
+ """
+ return _nx.concatenate(map(atleast_3d,tup),2)
+
+def _replace_zero_by_x_arrays(sub_arys):
+ for i in range(len(sub_arys)):
+ if len(_nx.shape(sub_arys[i])) == 0:
+ sub_arys[i] = _nx.array([])
+ elif _nx.sometrue(_nx.equal(_nx.shape(sub_arys[i]),0)):
+ sub_arys[i] = _nx.array([])
+ return sub_arys
+
+def array_split(ary,indices_or_sections,axis = 0):
+ """ Divide an array into a list of sub-arrays.
+
+ Description:
+ Divide ary into a list of sub-arrays along the
+ specified axis. If indices_or_sections is an integer,
+ ary is divided into that many equally sized arrays.
+ If it is impossible to make an equal split, each of the
+ leading arrays in the list have one additional member. If
+ indices_or_sections is a list of sorted integers, its
+ entries define the indexes where ary is split.
+
+ Arguments:
+ ary -- N-D array.
+ Array to be divided into sub-arrays.
+ indices_or_sections -- integer or 1D array.
+ If integer, defines the number of (close to) equal sized
+ sub-arrays. If it is a 1D array of sorted indices, it
+ defines the indexes at which ary is divided. Any empty
+ list results in a single sub-array equal to the original
+ array.
+ axis -- integer. default=0.
+ Specifies the axis along which to split ary.
+ Caveats:
+ Currently, the default for axis is 0. This
+ means a 2D array is divided into multiple groups
+ of rows. This seems like the appropriate default, but
+ we've agreed most other functions should default to
+ axis=-1. Perhaps we should use axis=-1 for consistency.
+ However, we could also make the argument that SciPy
+ works on "rows" by default. sum() sums up rows of
+ values. split() will split data into rows. Opinions?
+ """
+ try:
+ Ntotal = ary.shape[axis]
+ except AttributeError:
+ Ntotal = len(ary)
+ try: # handle scalar case.
+ Nsections = len(indices_or_sections) + 1
+ div_points = [0] + list(indices_or_sections) + [Ntotal]
+ except TypeError: #indices_or_sections is a scalar, not an array.
+ Nsections = int(indices_or_sections)
+ if Nsections <= 0:
+ raise ValueError, 'number sections must be larger than 0.'
+ Neach_section,extras = divmod(Ntotal,Nsections)
+ section_sizes = [0] + \
+ extras * [Neach_section+1] + \
+ (Nsections-extras) * [Neach_section]
+ div_points = _nx.array(section_sizes).cumsum()
+
+ sub_arys = []
+ sary = _nx.swapaxes(ary,axis,0)
+ for i in range(Nsections):
+ st = div_points[i]; end = div_points[i+1]
+ sub_arys.append(_nx.swapaxes(sary[st:end],axis,0))
+
+ # there is a wierd issue with array slicing that allows
+ # 0x10 arrays and other such things. The following cluge is needed
+ # to get around this issue.
+ sub_arys = _replace_zero_by_x_arrays(sub_arys)
+ # end cluge.
+
+ return sub_arys
+
+def split(ary,indices_or_sections,axis=0):
+ """ Divide an array into a list of sub-arrays.
+
+ Description:
+ Divide ary into a list of sub-arrays along the
+ specified axis. If indices_or_sections is an integer,
+ ary is divided into that many equally sized arrays.
+ If it is impossible to make an equal split, an error is
+ raised. This is the only way this function differs from
+ the array_split() function. If indices_or_sections is a
+ list of sorted integers, its entries define the indexes
+ where ary is split.
+
+ Arguments:
+ ary -- N-D array.
+ Array to be divided into sub-arrays.
+ indices_or_sections -- integer or 1D array.
+ If integer, defines the number of (close to) equal sized
+ sub-arrays. If it is a 1D array of sorted indices, it
+ defines the indexes at which ary is divided. Any empty
+ list results in a single sub-array equal to the original
+ array.
+ axis -- integer. default=0.
+ Specifies the axis along which to split ary.
+ Caveats:
+ Currently, the default for axis is 0. This
+ means a 2D array is divided into multiple groups
+ of rows. This seems like the appropriate default, but
+ we've agreed most other functions should default to
+ axis=-1. Perhaps we should use axis=-1 for consistency.
+ However, we could also make the argument that SciPy
+ works on "rows" by default. sum() sums up rows of
+ values. split() will split data into rows. Opinions?
+ """
+ try: len(indices_or_sections)
+ except TypeError:
+ sections = indices_or_sections
+ N = ary.shape[axis]
+ if N % sections:
+ raise ValueError, 'array split does not result in an equal division'
+ res = array_split(ary,indices_or_sections,axis)
+ return res
+
+def hsplit(ary,indices_or_sections):
+ """ Split ary into multiple columns of sub-arrays
+
+ Description:
+ Split a single array into multiple sub arrays. The array is
+ divided into groups of columns. If indices_or_sections is
+ an integer, ary is divided into that many equally sized sub arrays.
+ If it is impossible to make the sub-arrays equally sized, the
+ operation throws a ValueError exception. See array_split and
+ split for other options on indices_or_sections.
+ Arguments:
+ ary -- N-D array.
+ Array to be divided into sub-arrays.
+ indices_or_sections -- integer or 1D array.
+ If integer, defines the number of (close to) equal sized
+ sub-arrays. If it is a 1D array of sorted indices, it
+ defines the indexes at which ary is divided. Any empty
+ list results in a single sub-array equal to the original
+ array.
+ Returns:
+ sequence of sub-arrays. The returned arrays have the same
+ number of dimensions as the input array.
+ Related:
+ hstack, split, array_split, vsplit, dsplit.
+ Examples:
+ >>> import scipy
+ >>> a= array((1,2,3,4))
+ >>> scipy.hsplit(a,2)
+ [array([1, 2]), array([3, 4])]
+ >>> a = array([[1,2,3,4],[1,2,3,4]])
+ [array([[1, 2],
+ [1, 2]]), array([[3, 4],
+ [3, 4]])]
+
+ """
+ if len(_nx.shape(ary)) == 0:
+ raise ValueError, 'hsplit only works on arrays of 1 or more dimensions'
+ if len(ary.shape) > 1:
+ return split(ary,indices_or_sections,1)
+ else:
+ return split(ary,indices_or_sections,0)
+
+def vsplit(ary,indices_or_sections):
+ """ Split ary into multiple rows of sub-arrays
+
+ Description:
+ Split a single array into multiple sub arrays. The array is
+ divided into groups of rows. If indices_or_sections is
+ an integer, ary is divided into that many equally sized sub arrays.
+ If it is impossible to make the sub-arrays equally sized, the
+ operation throws a ValueError exception. See array_split and
+ split for other options on indices_or_sections.
+ Arguments:
+ ary -- N-D array.
+ Array to be divided into sub-arrays.
+ indices_or_sections -- integer or 1D array.
+ If integer, defines the number of (close to) equal sized
+ sub-arrays. If it is a 1D array of sorted indices, it
+ defines the indexes at which ary is divided. Any empty
+ list results in a single sub-array equal to the original
+ array.
+ Returns:
+ sequence of sub-arrays. The returned arrays have the same
+ number of dimensions as the input array.
+ Caveats:
+ How should we handle 1D arrays here? I am currently raising
+ an error when I encounter them. Any better approach?
+
+ Should we reduce the returned array to their minium dimensions
+ by getting rid of any dimensions that are 1?
+ Related:
+ vstack, split, array_split, hsplit, dsplit.
+ Examples:
+ import scipy
+ >>> a = array([[1,2,3,4],
+ ... [1,2,3,4]])
+ >>> scipy.vsplit(a)
+ [array([ [1, 2, 3, 4]]), array([ [1, 2, 3, 4]])]
+
+ """
+ if len(_nx.shape(ary)) < 2:
+ raise ValueError, 'vsplit only works on arrays of 2 or more dimensions'
+ return split(ary,indices_or_sections,0)
+
+def dsplit(ary,indices_or_sections):
+ """ Split ary into multiple sub-arrays along the 3rd axis (depth)
+
+ Description:
+ Split a single array into multiple sub arrays. The array is
+ divided into groups along the 3rd axis. If indices_or_sections is
+ an integer, ary is divided into that many equally sized sub arrays.
+ If it is impossible to make the sub-arrays equally sized, the
+ operation throws a ValueError exception. See array_split and
+ split for other options on indices_or_sections.
+ Arguments:
+ ary -- N-D array.
+ Array to be divided into sub-arrays.
+ indices_or_sections -- integer or 1D array.
+ If integer, defines the number of (close to) equal sized
+ sub-arrays. If it is a 1D array of sorted indices, it
+ defines the indexes at which ary is divided. Any empty
+ list results in a single sub-array equal to the original
+ array.
+ Returns:
+ sequence of sub-arrays. The returned arrays have the same
+ number of dimensions as the input array.
+ Caveats:
+ See vsplit caveats.
+ Related:
+ dstack, split, array_split, hsplit, vsplit.
+ Examples:
+ >>> a = array([[[1,2,3,4],[1,2,3,4]]])
+ [array([ [[1, 2],
+ [1, 2]]]), array([ [[3, 4],
+ [3, 4]]])]
+
+ """
+ if len(_nx.shape(ary)) < 3:
+ raise ValueError, 'vsplit only works on arrays of 3 or more dimensions'
+ return split(ary,indices_or_sections,2)
+
diff --git a/numpy/lib/src/_compiled_base.c b/numpy/lib/src/_compiled_base.c
new file mode 100644
index 000000000..3ce3743d7
--- /dev/null
+++ b/numpy/lib/src/_compiled_base.c
@@ -0,0 +1,453 @@
+#include "Python.h"
+#include "structmember.h"
+#include "scipy/arrayobject.h"
+
+static PyObject *ErrorObject;
+#define Py_Try(BOOLEAN) {if (!(BOOLEAN)) goto fail;}
+#define Py_Assert(BOOLEAN,MESS) {if (!(BOOLEAN)) { \
+ PyErr_SetString(ErrorObject, (MESS)); \
+ goto fail;} \
+ }
+
+static intp
+incr_slot_ (double x, double *bins, intp lbins)
+{
+ intp i ;
+ for ( i = 0 ; i < lbins ; i ++ )
+ if ( x < bins [i] )
+ return i ;
+ return lbins ;
+}
+
+static intp
+decr_slot_ (double x, double * bins, intp lbins)
+{
+ intp i ;
+ for ( i = lbins - 1 ; i >= 0; i -- )
+ if (x < bins [i])
+ return i + 1 ;
+ return 0 ;
+}
+
+static int
+monotonic_ (double * a, int lena)
+{
+ int i;
+ if (a [0] <= a [1]) /* possibly monotonic increasing */
+ {
+ for (i = 1 ; i < lena - 1; i ++)
+ if (a [i] > a [i + 1]) return 0 ;
+ return 1 ;
+ }
+ else /* possibly monotonic decreasing */
+ {
+ for (i = 1 ; i < lena - 1; i ++)
+ if (a [i] < a [i + 1]) return 0 ;
+ return -1 ;
+ }
+}
+
+
+
+static intp
+mxx (intp *i , intp len)
+{
+ /* find the index of the maximum element of an integer array */
+ intp mx = 0, max = i[0] ;
+ intp j ;
+ for ( j = 1 ; j < len; j ++ )
+ if ( i [j] > max )
+ {max = i [j] ;
+ mx = j ;}
+ return mx;
+}
+
+static intp
+mnx (intp *i , intp len)
+{
+ /* find the index of the minimum element of an integer array */
+ intp mn = 0, min = i [0] ;
+ intp j ;
+ for ( j = 1 ; j < len; j ++ )
+ if ( i [j] < min )
+ {min = i [j] ;
+ mn = j ;}
+ return mn;
+}
+
+
+static PyObject *
+arr_bincount(PyObject *self, PyObject *args, PyObject *kwds)
+{
+ /* histogram accepts one or two arguments. The first is an array
+ * of non-negative integers and the second, if present, is an
+ * array of weights, which must be promotable to double.
+ * Call these arguments list and weight. Both must be one-
+ * dimensional. len (weight) == len(list)
+ * If weight is not present:
+ * histogram (list) [i] is the number of occurrences of i in list.
+ * If weight is present:
+ * histogram (list, weight) [i] is the sum of all weight [j]
+ * where list [j] == i. */
+ /* self is not used */
+ PyArray_Descr *type;
+ PyObject *list = NULL, *weight=Py_None ;
+ PyObject *lst=NULL, *ans=NULL, *wts=NULL;
+ intp *numbers, *ians, len , mxi, mni, ans_size;
+ int i;
+ double *weights , *dans;
+ static char *kwlist[] = {"list", "weights", NULL};
+
+
+ Py_Try(PyArg_ParseTupleAndKeywords(args, kwds, "O|O", kwlist,
+ &list, &weight));
+ Py_Try(lst = PyArray_ContiguousFromAny(list, PyArray_INTP, 1, 1));
+ len = PyArray_SIZE(lst);
+ numbers = (intp *) PyArray_DATA(lst);
+ mxi = mxx (numbers, len) ;
+ mni = mnx (numbers, len) ;
+ Py_Assert(numbers[mni] >= 0,
+ "irst argument of bincount must be non-negative");
+ ans_size = numbers [mxi] + 1 ;
+ type = PyArray_DescrFromType(PyArray_INTP);
+ if (weight == Py_None) {
+ Py_Try(ans = PyArray_Zeros(1, &ans_size, type, 0));
+ ians = (intp *)(PyArray_DATA(ans));
+ for (i = 0 ; i < len ; i++)
+ ians [numbers [i]] += 1 ;
+ Py_DECREF(lst);
+ }
+ else {
+ Py_Try(wts = PyArray_ContiguousFromAny(weight,
+ PyArray_DOUBLE, 1, 1));
+ weights = (double *)PyArray_DATA (wts);
+ Py_Assert(PyArray_SIZE(wts) == len, "bincount: length of weights " \
+ "does not match that of list");
+ type = PyArray_DescrFromType(PyArray_DOUBLE);
+ Py_Try(ans = PyArray_Zeros(1, &ans_size, type, 0));
+ dans = (double *)PyArray_DATA (ans);
+ for (i = 0 ; i < len ; i++) {
+ dans[numbers[i]] += weights[i];
+ }
+ Py_DECREF(lst);
+ Py_DECREF(wts);
+ }
+ return ans;
+
+ fail:
+ Py_XDECREF(lst);
+ Py_XDECREF(wts);
+ Py_XDECREF(ans);
+ return NULL;
+}
+
+
+static PyObject *
+arr_digitize(PyObject *self, PyObject *args, PyObject *kwds)
+{
+ /* digitize (x, bins) returns an array of python integers the same
+ length of x. The values i returned are such that
+ bins [i - 1] <= x < bins [i] if bins is monotonically increasing,
+ or bins [i - 1] > x >= bins [i] if bins is monotonically decreasing.
+ Beyond the bounds of bins, returns either i = 0 or i = len (bins)
+ as appropriate. */
+ /* self is not used */
+ PyObject *ox, *obins ;
+ PyObject *ax=NULL, *abins=NULL, *aret=NULL;
+ double *dx, *dbins ;
+ intp lbins, lx ; /* lengths */
+ intp *iret;
+ int m, i ;
+ static char *kwlist[] = {"x", "bins", NULL};
+ PyArray_Descr *type;
+
+ Py_Try(PyArg_ParseTupleAndKeywords(args, kwds, "OO", kwlist,
+ &ox, &obins));
+
+ type = PyArray_DescrFromType(PyArray_DOUBLE);
+ Py_Try(ax=PyArray_FromAny(ox, type, 1, 1, CARRAY_FLAGS));
+ Py_Try(abins = PyArray_FromAny(obins, type, 1, 1, CARRAY_FLAGS));
+
+ lx = PyArray_SIZE(ax);
+ dx = (double *)PyArray_DATA(ax);
+ lbins = PyArray_SIZE(abins);
+ dbins = (double *)PyArray_DATA(abins);
+ Py_Try(aret = PyArray_SimpleNew(1, &lx, PyArray_INTP));
+ iret = (intp *)PyArray_DATA(aret);
+
+ Py_Assert(lx > 0 && lbins > 0,
+ "x and bins both must have non-zero length");
+
+ if (lbins == 1) {
+ for (i=0 ; i<lx ; i++)
+ if (dx [i] >= dbins[0])
+ iret[i] = 1;
+ else
+ iret[i] = 0;
+ }
+ else {
+ m = monotonic_ (dbins, lbins) ;
+ if ( m == -1 ) {
+ for ( i = 0 ; i < lx ; i ++ )
+ iret [i] = decr_slot_ (dx [i], dbins, lbins) ;
+ }
+ else if ( m == 1 ) {
+ for ( i = 0 ; i < lx ; i ++ )
+ iret [i] = incr_slot_ ((float)dx [i], dbins, lbins) ;
+ }
+ else Py_Assert(0, "bins must be montonically increasing or decreasing");
+ }
+
+ Py_DECREF(ax);
+ Py_DECREF(abins);
+ return aret;
+
+ fail:
+ Py_XDECREF(ax);
+ Py_XDECREF(abins);
+ Py_XDECREF(aret);
+ return NULL;
+}
+
+
+
+static char arr_insert__doc__[] = "Insert vals sequentially into equivalent 1-d positions indicated by mask.";
+
+static PyObject *
+arr_insert(PyObject *self, PyObject *args, PyObject *kwdict)
+{
+ /* Returns input array with values inserted sequentially into places
+ indicated by the mask
+ */
+ PyObject *mask=NULL, *vals=NULL;
+ PyArrayObject *ainput=NULL, *amask=NULL, *avals=NULL,
+ *tmp=NULL;
+ int numvals, totmask, sameshape;
+ char *input_data, *mptr, *vptr, *zero=NULL;
+ int melsize, delsize, copied, nd;
+ intp *instrides, *inshape;
+ int mindx, rem_indx, indx, i, k, objarray;
+
+ static char *kwlist[] = {"input","mask","vals",NULL};
+
+ if (!PyArg_ParseTupleAndKeywords(args, kwdict, "O&OO", kwlist,
+ PyArray_Converter, &ainput,
+ &mask, &vals))
+ goto fail;
+
+ amask = (PyArrayObject *) PyArray_FROM_OF(mask, CARRAY_FLAGS);
+ if (amask == NULL) goto fail;
+ /* Cast an object array */
+ if (amask->descr->type_num == PyArray_OBJECT) {
+ tmp = (PyArrayObject *)PyArray_Cast(amask, PyArray_INTP);
+ if (tmp == NULL) goto fail;
+ Py_DECREF(amask);
+ amask = tmp;
+ }
+
+ sameshape = 1;
+ if (amask->nd == ainput->nd) {
+ for (k=0; k < amask->nd; k++)
+ if (amask->dimensions[k] != ainput->dimensions[k])
+ sameshape = 0;
+ }
+ else { /* Test to see if amask is 1d */
+ if (amask->nd != 1) sameshape = 0;
+ else if ((PyArray_SIZE(ainput)) != PyArray_SIZE(amask)) sameshape = 0;
+ }
+ if (!sameshape) {
+ PyErr_SetString(PyExc_TypeError,
+ "mask array must be 1-d or same shape as input array");
+ goto fail;
+ }
+
+ avals = (PyArrayObject *)PyArray_FromObject(vals, ainput->descr->type_num, 0, 1);
+ if (avals == NULL) goto fail;
+
+ numvals = PyArray_SIZE(avals);
+ nd = ainput->nd;
+ input_data = ainput->data;
+ mptr = amask->data;
+ melsize = amask->descr->elsize;
+ vptr = avals->data;
+ delsize = avals->descr->elsize;
+ zero = PyArray_Zero(amask);
+ if (zero == NULL)
+ goto fail;
+ objarray = (ainput->descr->type_num == PyArray_OBJECT);
+
+ /* Handle zero-dimensional case separately */
+ if (nd == 0) {
+ if (memcmp(mptr,zero,melsize) != 0) {
+ /* Copy value element over to input array */
+ memcpy(input_data,vptr,delsize);
+ if (objarray) Py_INCREF(*((PyObject **)vptr));
+ }
+ Py_DECREF(amask);
+ Py_DECREF(avals);
+ PyDataMem_FREE(zero);
+ Py_INCREF(Py_None);
+ return Py_None;
+ }
+
+ /* Walk through mask array, when non-zero is encountered
+ copy next value in the vals array to the input array.
+ If we get through the value array, repeat it as necessary.
+ */
+ totmask = (int) PyArray_SIZE(amask);
+ copied = 0;
+ instrides = ainput->strides;
+ inshape = ainput->dimensions;
+ for (mindx = 0; mindx < totmask; mindx++) {
+ if (memcmp(mptr,zero,melsize) != 0) {
+ /* compute indx into input array
+ */
+ rem_indx = mindx;
+ indx = 0;
+ for(i=nd-1; i > 0; --i) {
+ indx += (rem_indx % inshape[i]) * instrides[i];
+ rem_indx /= inshape[i];
+ }
+ indx += rem_indx * instrides[0];
+ /* fprintf(stderr, "mindx = %d, indx=%d\n", mindx, indx); */
+ /* Copy value element over to input array */
+ memcpy(input_data+indx,vptr,delsize);
+ if (objarray) Py_INCREF(*((PyObject **)vptr));
+ vptr += delsize;
+ copied += 1;
+ /* If we move past value data. Reset */
+ if (copied >= numvals) vptr = avals->data;
+ }
+ mptr += melsize;
+ }
+
+ Py_DECREF(amask);
+ Py_DECREF(avals);
+ PyDataMem_FREE(zero);
+ Py_DECREF(ainput);
+ Py_INCREF(Py_None);
+ return Py_None;
+
+ fail:
+ PyDataMem_FREE(zero);
+ Py_XDECREF(ainput);
+ Py_XDECREF(amask);
+ Py_XDECREF(avals);
+ return NULL;
+}
+
+
+static PyTypeObject *PyMemberDescr_TypePtr=NULL;
+static PyTypeObject *PyGetSetDescr_TypePtr=NULL;
+
+/* Can only be called if doc is currently NULL
+*/
+static PyObject *
+arr_add_docstring(PyObject *dummy, PyObject *args)
+{
+ PyObject *obj;
+ PyObject *str;
+ char *docstr;
+ static char *msg = "already has a docstring";
+
+ if (!PyArg_ParseTuple(args, "OO!", &obj, &PyString_Type, &str))
+ return NULL;
+
+ docstr = PyString_AS_STRING(str);
+
+#define _TESTDOC1(typebase) (obj->ob_type == &Py##typebase##_Type)
+#define _TESTDOC2(typebase) (obj->ob_type == Py##typebase##_TypePtr)
+#define _ADDDOC(typebase, doc, name) { \
+ Py##typebase##Object *new = (Py##typebase##Object *)obj; \
+ if (!(doc)) { \
+ doc = docstr; \
+ } \
+ else { \
+ PyErr_Format(PyExc_RuntimeError, \
+ "%s method %s",name, msg); \
+ return NULL; \
+ } \
+ }
+
+ if _TESTDOC1(CFunction)
+ _ADDDOC(CFunction, new->m_ml->ml_doc, new->m_ml->ml_name)
+ else if _TESTDOC1(Type)
+ _ADDDOC(Type, new->tp_doc, new->tp_name)
+ else if _TESTDOC2(MemberDescr)
+ _ADDDOC(MemberDescr, new->d_member->doc, new->d_member->name)
+ else if _TESTDOC2(GetSetDescr)
+ _ADDDOC(GetSetDescr, new->d_getset->doc, new->d_getset->name)
+ else {
+ PyErr_SetString(PyExc_TypeError,
+ "Cannot set a docstring for that object");
+ return NULL;
+ }
+
+#undef _TESTDOC1
+#undef _TESTDOC2
+#undef _ADDDOC
+
+ Py_INCREF(str);
+ Py_INCREF(Py_None);
+ return Py_None;
+}
+
+static struct PyMethodDef methods[] = {
+ {"_insert", (PyCFunction)arr_insert, METH_VARARGS | METH_KEYWORDS,
+ arr_insert__doc__},
+ {"bincount", (PyCFunction)arr_bincount,
+ METH_VARARGS | METH_KEYWORDS, NULL},
+ {"digitize", (PyCFunction)arr_digitize, METH_VARARGS | METH_KEYWORDS,
+ NULL},
+ {"add_docstring", (PyCFunction)arr_add_docstring, METH_VARARGS,
+ NULL},
+ {NULL, NULL} /* sentinel */
+};
+
+static void
+define_types(void)
+{
+ PyObject *tp_dict;
+ PyObject *myobj;
+
+ tp_dict = PyArrayDescr_Type.tp_dict;
+ /* Get "subdescr" */
+ myobj = PyDict_GetItemString(tp_dict, "fields");
+ if (myobj == NULL) return;
+ PyGetSetDescr_TypePtr = myobj->ob_type;
+ myobj = PyDict_GetItemString(tp_dict, "alignment");
+ if (myobj == NULL) return;
+ PyMemberDescr_TypePtr = myobj->ob_type;
+ return;
+}
+
+/* Initialization function for the module (*must* be called initArray) */
+
+DL_EXPORT(void) init_compiled_base(void) {
+ PyObject *m, *d, *s;
+
+ /* Create the module and add the functions */
+ m = Py_InitModule("scipy.base._compiled_base", methods);
+
+ /* Import the array and ufunc objects */
+ import_array();
+
+ /* Add some symbolic constants to the module */
+ d = PyModule_GetDict(m);
+
+ s = PyString_FromString("0.5");
+ PyDict_SetItemString(d, "__version__", s);
+ Py_DECREF(s);
+
+ ErrorObject = PyString_FromString("scipy.base._compiled_base.error");
+ PyDict_SetItemString(d, "error", ErrorObject);
+ Py_DECREF(ErrorObject);
+
+
+ /* define PyGetSetDescr_Type and PyMemberDescr_Type */
+ define_types();
+
+ /* Check for errors */
+ if (PyErr_Occurred())
+ Py_FatalError("can't initialize module _compiled_base");
+}
diff --git a/numpy/lib/test_shape_base.py b/numpy/lib/test_shape_base.py
new file mode 100644
index 000000000..005868e96
--- /dev/null
+++ b/numpy/lib/test_shape_base.py
@@ -0,0 +1,364 @@
+
+from scipy.testing import *
+set_package_path()
+import scipy.base;
+from scipy.base import *
+restore_path()
+
+class test_apply_along_axis(ScipyTestCase):
+ def check_simple(self):
+ a = ones((20,10),'d')
+ assert_array_equal(apply_along_axis(len,0,a),len(a)*ones(shape(a)[1]))
+ def check_simple101(self,level=11):
+ # This test causes segmentation fault (Numeric 23.3,23.6,Python 2.3.4)
+ # when enabled and shape(a)[1]>100. See Issue 202.
+ a = ones((10,101),'d')
+ assert_array_equal(apply_along_axis(len,0,a),len(a)*ones(shape(a)[1]))
+
+class test_array_split(ScipyTestCase):
+ def check_integer_0_split(self):
+ a = arange(10)
+ try:
+ res = array_split(a,0)
+ assert(0) # it should have thrown a value error
+ except ValueError:
+ pass
+ def check_integer_split(self):
+ a = arange(10)
+ res = array_split(a,1)
+ desired = [arange(10)]
+ compare_results(res,desired)
+
+ res = array_split(a,2)
+ desired = [arange(5),arange(5,10)]
+ compare_results(res,desired)
+
+ res = array_split(a,3)
+ desired = [arange(4),arange(4,7),arange(7,10)]
+ compare_results(res,desired)
+
+ res = array_split(a,4)
+ desired = [arange(3),arange(3,6),arange(6,8),arange(8,10)]
+ compare_results(res,desired)
+
+ res = array_split(a,5)
+ desired = [arange(2),arange(2,4),arange(4,6),arange(6,8),arange(8,10)]
+ compare_results(res,desired)
+
+ res = array_split(a,6)
+ desired = [arange(2),arange(2,4),arange(4,6),arange(6,8),arange(8,9),
+ arange(9,10)]
+ compare_results(res,desired)
+
+ res = array_split(a,7)
+ desired = [arange(2),arange(2,4),arange(4,6),arange(6,7),arange(7,8),
+ arange(8,9), arange(9,10)]
+ compare_results(res,desired)
+
+ res = array_split(a,8)
+ desired = [arange(2),arange(2,4),arange(4,5),arange(5,6),arange(6,7),
+ arange(7,8), arange(8,9), arange(9,10)]
+ compare_results(res,desired)
+
+ res = array_split(a,9)
+ desired = [arange(2),arange(2,3),arange(3,4),arange(4,5),arange(5,6),
+ arange(6,7), arange(7,8), arange(8,9), arange(9,10)]
+ compare_results(res,desired)
+
+ res = array_split(a,10)
+ desired = [arange(1),arange(1,2),arange(2,3),arange(3,4),
+ arange(4,5),arange(5,6), arange(6,7), arange(7,8),
+ arange(8,9), arange(9,10)]
+ compare_results(res,desired)
+
+ res = array_split(a,11)
+ desired = [arange(1),arange(1,2),arange(2,3),arange(3,4),
+ arange(4,5),arange(5,6), arange(6,7), arange(7,8),
+ arange(8,9), arange(9,10),array([])]
+ compare_results(res,desired)
+ def check_integer_split_2D_rows(self):
+ a = array([arange(10),arange(10)])
+ res = array_split(a,3,axis=0)
+ desired = [array([arange(10)]),array([arange(10)]),array([])]
+ compare_results(res,desired)
+ def check_integer_split_2D_cols(self):
+ a = array([arange(10),arange(10)])
+ res = array_split(a,3,axis=-1)
+ desired = [array([arange(4),arange(4)]),
+ array([arange(4,7),arange(4,7)]),
+ array([arange(7,10),arange(7,10)])]
+ compare_results(res,desired)
+ def check_integer_split_2D_default(self):
+ """ This will fail if we change default axis
+ """
+ a = array([arange(10),arange(10)])
+ res = array_split(a,3)
+ desired = [array([arange(10)]),array([arange(10)]),array([])]
+ compare_results(res,desired)
+ #perhaps should check higher dimensions
+
+ def check_index_split_simple(self):
+ a = arange(10)
+ indices = [1,5,7]
+ res = array_split(a,indices,axis=-1)
+ desired = [arange(0,1),arange(1,5),arange(5,7),arange(7,10)]
+ compare_results(res,desired)
+
+ def check_index_split_low_bound(self):
+ a = arange(10)
+ indices = [0,5,7]
+ res = array_split(a,indices,axis=-1)
+ desired = [array([]),arange(0,5),arange(5,7),arange(7,10)]
+ compare_results(res,desired)
+ def check_index_split_high_bound(self):
+ a = arange(10)
+ indices = [0,5,7,10,12]
+ res = array_split(a,indices,axis=-1)
+ desired = [array([]),arange(0,5),arange(5,7),arange(7,10),
+ array([]),array([])]
+ compare_results(res,desired)
+
+class test_split(ScipyTestCase):
+ """* This function is essentially the same as array_split,
+ except that it test if splitting will result in an
+ equal split. Only test for this case.
+ *"""
+ def check_equal_split(self):
+ a = arange(10)
+ res = split(a,2)
+ desired = [arange(5),arange(5,10)]
+ compare_results(res,desired)
+
+ def check_unequal_split(self):
+ a = arange(10)
+ try:
+ res = split(a,3)
+ assert(0) # should raise an error
+ except ValueError:
+ pass
+
+class test_atleast_1d(ScipyTestCase):
+ def check_0D_array(self):
+ a = array(1); b = array(2);
+ res=map(atleast_1d,[a,b])
+ desired = [array([1]),array([2])]
+ assert_array_equal(res,desired)
+ def check_1D_array(self):
+ a = array([1,2]); b = array([2,3]);
+ res=map(atleast_1d,[a,b])
+ desired = [array([1,2]),array([2,3])]
+ assert_array_equal(res,desired)
+ def check_2D_array(self):
+ a = array([[1,2],[1,2]]); b = array([[2,3],[2,3]]);
+ res=map(atleast_1d,[a,b])
+ desired = [a,b]
+ assert_array_equal(res,desired)
+ def check_3D_array(self):
+ a = array([[1,2],[1,2]]); b = array([[2,3],[2,3]]);
+ a = array([a,a]);b = array([b,b]);
+ res=map(atleast_1d,[a,b])
+ desired = [a,b]
+ assert_array_equal(res,desired)
+ def check_r1array(self):
+ """ Test to make sure equivalent Travis O's r1array function
+ """
+ assert(atleast_1d(3).shape == (1,))
+ assert(atleast_1d(3j).shape == (1,))
+ assert(atleast_1d(3L).shape == (1,))
+ assert(atleast_1d(3.0).shape == (1,))
+ assert(atleast_1d([[2,3],[4,5]]).shape == (2,2))
+
+class test_atleast_2d(ScipyTestCase):
+ def check_0D_array(self):
+ a = array(1); b = array(2);
+ res=map(atleast_2d,[a,b])
+ desired = [array([[1]]),array([[2]])]
+ assert_array_equal(res,desired)
+ def check_1D_array(self):
+ a = array([1,2]); b = array([2,3]);
+ res=map(atleast_2d,[a,b])
+ desired = [array([[1,2]]),array([[2,3]])]
+ assert_array_equal(res,desired)
+ def check_2D_array(self):
+ a = array([[1,2],[1,2]]); b = array([[2,3],[2,3]]);
+ res=map(atleast_2d,[a,b])
+ desired = [a,b]
+ assert_array_equal(res,desired)
+ def check_3D_array(self):
+ a = array([[1,2],[1,2]]); b = array([[2,3],[2,3]]);
+ a = array([a,a]);b = array([b,b]);
+ res=map(atleast_2d,[a,b])
+ desired = [a,b]
+ assert_array_equal(res,desired)
+ def check_r2array(self):
+ """ Test to make sure equivalent Travis O's r2array function
+ """
+ assert(atleast_2d(3).shape == (1,1))
+ assert(atleast_2d([3j,1]).shape == (1,2))
+ assert(atleast_2d([[[3,1],[4,5]],[[3,5],[1,2]]]).shape == (2,2,2))
+
+class test_atleast_3d(ScipyTestCase):
+ def check_0D_array(self):
+ a = array(1); b = array(2);
+ res=map(atleast_3d,[a,b])
+ desired = [array([[[1]]]),array([[[2]]])]
+ assert_array_equal(res,desired)
+ def check_1D_array(self):
+ a = array([1,2]); b = array([2,3]);
+ res=map(atleast_3d,[a,b])
+ desired = [array([[[1],[2]]]),array([[[2],[3]]])]
+ assert_array_equal(res,desired)
+ def check_2D_array(self):
+ a = array([[1,2],[1,2]]); b = array([[2,3],[2,3]]);
+ res=map(atleast_3d,[a,b])
+ desired = [a[:,:,NewAxis],b[:,:,NewAxis]]
+ assert_array_equal(res,desired)
+ def check_3D_array(self):
+ a = array([[1,2],[1,2]]); b = array([[2,3],[2,3]]);
+ a = array([a,a]);b = array([b,b]);
+ res=map(atleast_3d,[a,b])
+ desired = [a,b]
+ assert_array_equal(res,desired)
+
+class test_hstack(ScipyTestCase):
+ def check_0D_array(self):
+ a = array(1); b = array(2);
+ res=hstack([a,b])
+ desired = array([1,2])
+ assert_array_equal(res,desired)
+ def check_1D_array(self):
+ a = array([1]); b = array([2]);
+ res=hstack([a,b])
+ desired = array([1,2])
+ assert_array_equal(res,desired)
+ def check_2D_array(self):
+ a = array([[1],[2]]); b = array([[1],[2]]);
+ res=hstack([a,b])
+ desired = array([[1,1],[2,2]])
+ assert_array_equal(res,desired)
+
+class test_vstack(ScipyTestCase):
+ def check_0D_array(self):
+ a = array(1); b = array(2);
+ res=vstack([a,b])
+ desired = array([[1],[2]])
+ assert_array_equal(res,desired)
+ def check_1D_array(self):
+ a = array([1]); b = array([2]);
+ res=vstack([a,b])
+ desired = array([[1],[2]])
+ assert_array_equal(res,desired)
+ def check_2D_array(self):
+ a = array([[1],[2]]); b = array([[1],[2]]);
+ res=vstack([a,b])
+ desired = array([[1],[2],[1],[2]])
+ assert_array_equal(res,desired)
+ def check_2D_array2(self):
+ a = array([1,2]); b = array([1,2]);
+ res=vstack([a,b])
+ desired = array([[1,2],[1,2]])
+ assert_array_equal(res,desired)
+
+class test_dstack(ScipyTestCase):
+ def check_0D_array(self):
+ a = array(1); b = array(2);
+ res=dstack([a,b])
+ desired = array([[[1,2]]])
+ assert_array_equal(res,desired)
+ def check_1D_array(self):
+ a = array([1]); b = array([2]);
+ res=dstack([a,b])
+ desired = array([[[1,2]]])
+ assert_array_equal(res,desired)
+ def check_2D_array(self):
+ a = array([[1],[2]]); b = array([[1],[2]]);
+ res=dstack([a,b])
+ desired = array([[[1,1]],[[2,2,]]])
+ assert_array_equal(res,desired)
+ def check_2D_array2(self):
+ a = array([1,2]); b = array([1,2]);
+ res=dstack([a,b])
+ desired = array([[[1,1],[2,2]]])
+ assert_array_equal(res,desired)
+
+""" array_split has more comprehensive test of splitting.
+ only do simple test on hsplit, vsplit, and dsplit
+"""
+class test_hsplit(ScipyTestCase):
+ """ only testing for integer splits.
+ """
+ def check_0D_array(self):
+ a= array(1)
+ try:
+ hsplit(a,2)
+ assert(0)
+ except ValueError:
+ pass
+ def check_1D_array(self):
+ a= array([1,2,3,4])
+ res = hsplit(a,2)
+ desired = [array([1,2]),array([3,4])]
+ compare_results(res,desired)
+ def check_2D_array(self):
+ a= array([[1,2,3,4],
+ [1,2,3,4]])
+ res = hsplit(a,2)
+ desired = [array([[1,2],[1,2]]),array([[3,4],[3,4]])]
+ compare_results(res,desired)
+
+class test_vsplit(ScipyTestCase):
+ """ only testing for integer splits.
+ """
+ def check_1D_array(self):
+ a= array([1,2,3,4])
+ try:
+ vsplit(a,2)
+ assert(0)
+ except ValueError:
+ pass
+ def check_2D_array(self):
+ a= array([[1,2,3,4],
+ [1,2,3,4]])
+ res = vsplit(a,2)
+ desired = [array([[1,2,3,4]]),array([[1,2,3,4]])]
+ compare_results(res,desired)
+
+class test_dsplit(ScipyTestCase):
+ """ only testing for integer splits.
+ """
+ def check_2D_array(self):
+ a= array([[1,2,3,4],
+ [1,2,3,4]])
+ try:
+ dsplit(a,2)
+ assert(0)
+ except ValueError:
+ pass
+ def check_3D_array(self):
+ a= array([[[1,2,3,4],
+ [1,2,3,4]],
+ [[1,2,3,4],
+ [1,2,3,4]]])
+ res = dsplit(a,2)
+ desired = [array([[[1,2],[1,2]],[[1,2],[1,2]]]),
+ array([[[3,4],[3,4]],[[3,4],[3,4]]])]
+ compare_results(res,desired)
+
+class test_squeeze(ScipyTestCase):
+ def check_basic(self):
+ a = rand(20,10,10,1,1)
+ b = rand(20,1,10,1,20)
+ c = rand(1,1,20,10)
+ assert_array_equal(squeeze(a),reshape(a,(20,10,10)))
+ assert_array_equal(squeeze(b),reshape(b,(20,10,20)))
+ assert_array_equal(squeeze(c),reshape(c,(20,10)))
+
+# Utility
+
+def compare_results(res,desired):
+ for i in range(len(desired)):
+ assert_array_equal(res[i],desired[i])
+
+
+if __name__ == "__main__":
+ ScipyTest().run()
diff --git a/numpy/lib/tests/test_function_base.py b/numpy/lib/tests/test_function_base.py
new file mode 100644
index 000000000..fafd75eef
--- /dev/null
+++ b/numpy/lib/tests/test_function_base.py
@@ -0,0 +1,338 @@
+
+import sys
+
+from scipy.testing import *
+set_package_path()
+import scipy.base;reload(scipy.base)
+from scipy.base import *
+del sys.path[0]
+
+class test_any(ScipyTestCase):
+ def check_basic(self):
+ y1 = [0,0,1,0]
+ y2 = [0,0,0,0]
+ y3 = [1,0,1,0]
+ assert(any(y1))
+ assert(any(y3))
+ assert(not any(y2))
+
+ def check_nd(self):
+ y1 = [[0,0,0],[0,1,0],[1,1,0]]
+ assert(any(y1))
+ assert_array_equal(sometrue(y1),[1,1,0])
+ assert_array_equal(sometrue(y1,axis=1),[0,1,1])
+
+class test_all(ScipyTestCase):
+ def check_basic(self):
+ y1 = [0,1,1,0]
+ y2 = [0,0,0,0]
+ y3 = [1,1,1,1]
+ assert(not all(y1))
+ assert(all(y3))
+ assert(not all(y2))
+ assert(all(~array(y2)))
+
+ def check_nd(self):
+ y1 = [[0,0,1],[0,1,1],[1,1,1]]
+ assert(not all(y1))
+ assert_array_equal(alltrue(y1),[0,0,1])
+ assert_array_equal(alltrue(y1,axis=1),[0,0,1])
+
+class test_average(ScipyTestCase):
+ def check_basic(self):
+ y1 = array([1,2,3])
+ assert(average(y1) == 2.)
+ y2 = array([1.,2.,3.])
+ assert(average(y2) == 2.)
+ y3 = [0.,0.,0.]
+ assert(average(y3) == 0.)
+
+ y4 = ones((4,4))
+ y4[0,1] = 0
+ y4[1,0] = 2
+ assert_array_equal(y4.mean(0), average(y4, 0))
+ assert_array_equal(y4.mean(1), average(y4, 1))
+
+ y5 = rand(5,5)
+ assert_array_equal(y5.mean(0), average(y5, 0))
+ assert_array_equal(y5.mean(1), average(y5, 1))
+
+class test_logspace(ScipyTestCase):
+ def check_basic(self):
+ y = logspace(0,6)
+ assert(len(y)==50)
+ y = logspace(0,6,num=100)
+ assert(y[-1] == 10**6)
+ y = logspace(0,6,endpoint=0)
+ assert(y[-1] < 10**6)
+ y = logspace(0,6,num=7)
+ assert_array_equal(y,[1,10,100,1e3,1e4,1e5,1e6])
+
+class test_linspace(ScipyTestCase):
+ def check_basic(self):
+ y = linspace(0,10)
+ assert(len(y)==50)
+ y = linspace(2,10,num=100)
+ assert(y[-1] == 10)
+ y = linspace(2,10,endpoint=0)
+ assert(y[-1] < 10)
+ y,st = linspace(2,10,retstep=1)
+ assert_almost_equal(st,8/49.0)
+ assert_array_almost_equal(y,mgrid[2:10:50j],13)
+
+ def check_corner(self):
+ y = list(linspace(0,1,1))
+ assert y == [0.0], y
+ y = list(linspace(0,1,2.5))
+ assert y == [0.0, 1.0]
+
+class test_amax(ScipyTestCase):
+ def check_basic(self):
+ a = [3,4,5,10,-3,-5,6.0]
+ assert_equal(amax(a),10.0)
+ b = [[3,6.0, 9.0],
+ [4,10.0,5.0],
+ [8,3.0,2.0]]
+ assert_equal(amax(b,axis=0),[8.0,10.0,9.0])
+ assert_equal(amax(b,axis=1),[9.0,10.0,8.0])
+
+class test_amin(ScipyTestCase):
+ def check_basic(self):
+ a = [3,4,5,10,-3,-5,6.0]
+ assert_equal(amin(a),-5.0)
+ b = [[3,6.0, 9.0],
+ [4,10.0,5.0],
+ [8,3.0,2.0]]
+ assert_equal(amin(b,axis=0),[3.0,3.0,2.0])
+ assert_equal(amin(b,axis=1),[3.0,4.0,2.0])
+
+class test_ptp(ScipyTestCase):
+ def check_basic(self):
+ a = [3,4,5,10,-3,-5,6.0]
+ assert_equal(ptp(a),15.0)
+ b = [[3,6.0, 9.0],
+ [4,10.0,5.0],
+ [8,3.0,2.0]]
+ assert_equal(ptp(b,axis=0),[5.0,7.0,7.0])
+ assert_equal(ptp(b,axis=-1),[6.0,6.0,6.0])
+
+class test_cumsum(ScipyTestCase):
+ def check_basic(self):
+ ba = [1,2,10,11,6,5,4]
+ ba2 = [[1,2,3,4],[5,6,7,9],[10,3,4,5]]
+ for ctype in [int8,uint8,int16,uint16,int32,uint32,
+ float32,float64,complex64,complex128]:
+ a = array(ba,ctype)
+ a2 = array(ba2,ctype)
+ assert_array_equal(cumsum(a), array([1,3,13,24,30,35,39],ctype))
+ assert_array_equal(cumsum(a2,axis=0), array([[1,2,3,4],[6,8,10,13],
+ [16,11,14,18]],ctype))
+ assert_array_equal(cumsum(a2,axis=1),
+ array([[1,3,6,10],
+ [5,11,18,27],
+ [10,13,17,22]],ctype))
+
+class test_prod(ScipyTestCase):
+ def check_basic(self):
+ ba = [1,2,10,11,6,5,4]
+ ba2 = [[1,2,3,4],[5,6,7,9],[10,3,4,5]]
+ for ctype in [int16,uint16,int32,uint32,
+ float32,float64,complex64,complex128]:
+ a = array(ba,ctype)
+ a2 = array(ba2,ctype)
+ if ctype in ['1', 'b']:
+ self.failUnlessRaises(ArithmeticError, prod, a)
+ self.failUnlessRaises(ArithmeticError, prod, a2, 1)
+ self.failUnlessRaises(ArithmeticError, prod, a)
+ else:
+ assert_equal(prod(a),26400)
+ assert_array_equal(prod(a2,axis=0),
+ array([50,36,84,180],ctype))
+ assert_array_equal(prod(a2,axis=-1),array([24, 1890, 600],ctype))
+
+class test_cumprod(ScipyTestCase):
+ def check_basic(self):
+ ba = [1,2,10,11,6,5,4]
+ ba2 = [[1,2,3,4],[5,6,7,9],[10,3,4,5]]
+ for ctype in [int16,uint16,int32,uint32,
+ float32,float64,complex64,complex128]:
+ a = array(ba,ctype)
+ a2 = array(ba2,ctype)
+ if ctype in ['1', 'b']:
+ self.failUnlessRaises(ArithmeticError, cumprod, a)
+ self.failUnlessRaises(ArithmeticError, cumprod, a2, 1)
+ self.failUnlessRaises(ArithmeticError, cumprod, a)
+ else:
+ assert_array_equal(cumprod(a,axis=-1),
+ array([1, 2, 20, 220,
+ 1320, 6600, 26400],ctype))
+ assert_array_equal(cumprod(a2,axis=0),
+ array([[ 1, 2, 3, 4],
+ [ 5, 12, 21, 36],
+ [50, 36, 84, 180]],ctype))
+ assert_array_equal(cumprod(a2,axis=-1),
+ array([[ 1, 2, 6, 24],
+ [ 5, 30, 210, 1890],
+ [10, 30, 120, 600]],ctype))
+
+class test_diff(ScipyTestCase):
+ def check_basic(self):
+ x = [1,4,6,7,12]
+ out = array([3,2,1,5])
+ out2 = array([-1,-1,4])
+ out3 = array([0,5])
+ assert_array_equal(diff(x),out)
+ assert_array_equal(diff(x,n=2),out2)
+ assert_array_equal(diff(x,n=3),out3)
+
+ def check_nd(self):
+ x = 20*rand(10,20,30)
+ out1 = x[:,:,1:] - x[:,:,:-1]
+ out2 = out1[:,:,1:] - out1[:,:,:-1]
+ out3 = x[1:,:,:] - x[:-1,:,:]
+ out4 = out3[1:,:,:] - out3[:-1,:,:]
+ assert_array_equal(diff(x),out1)
+ assert_array_equal(diff(x,n=2),out2)
+ assert_array_equal(diff(x,axis=0),out3)
+ assert_array_equal(diff(x,n=2,axis=0),out4)
+
+class test_angle(ScipyTestCase):
+ def check_basic(self):
+ x = [1+3j,sqrt(2)/2.0+1j*sqrt(2)/2,1,1j,-1,-1j,1-3j,-1+3j]
+ y = angle(x)
+ yo = [arctan(3.0/1.0),arctan(1.0),0,pi/2,pi,-pi/2.0,
+ -arctan(3.0/1.0),pi-arctan(3.0/1.0)]
+ z = angle(x,deg=1)
+ zo = array(yo)*180/pi
+ assert_array_almost_equal(y,yo,11)
+ assert_array_almost_equal(z,zo,11)
+
+class test_trim_zeros(ScipyTestCase):
+ """ only testing for integer splits.
+ """
+ def check_basic(self):
+ a= array([0,0,1,2,3,4,0])
+ res = trim_zeros(a)
+ assert_array_equal(res,array([1,2,3,4]))
+ def check_leading_skip(self):
+ a= array([0,0,1,0,2,3,4,0])
+ res = trim_zeros(a)
+ assert_array_equal(res,array([1,0,2,3,4]))
+ def check_trailing_skip(self):
+ a= array([0,0,1,0,2,3,0,4,0])
+ res = trim_zeros(a)
+ assert_array_equal(res,array([1,0,2,3,0,4]))
+
+
+class test_extins(ScipyTestCase):
+ def check_basic(self):
+ a = array([1,3,2,1,2,3,3])
+ b = extract(a>1,a)
+ assert_array_equal(b,[3,2,2,3,3])
+ def check_insert(self):
+ a = array([1,4,3,2,5,8,7])
+ insert(a,[0,1,0,1,0,1,0],[2,4,6])
+ assert_array_equal(a,[1,2,3,4,5,6,7])
+ def check_both(self):
+ a = rand(10)
+ mask = a > 0.5
+ ac = a.copy()
+ c = extract(mask, a)
+ insert(a,mask,0)
+ insert(a,mask,c)
+ assert_array_equal(a,ac)
+
+class test_vectorize(ScipyTestCase):
+ def check_simple(self):
+ def addsubtract(a,b):
+ if a > b:
+ return a - b
+ else:
+ return a + b
+ f = vectorize(addsubtract)
+ r = f([0,3,6,9],[1,3,5,7])
+ assert_array_equal(r,[1,6,1,2])
+ def check_scalar(self):
+ def addsubtract(a,b):
+ if a > b:
+ return a - b
+ else:
+ return a + b
+ f = vectorize(addsubtract)
+ r = f([0,3,6,9],5)
+ assert_array_equal(r,[5,8,1,4])
+
+
+
+class test_unwrap(ScipyTestCase):
+ def check_simple(self):
+ #check that unwrap removes jumps greather that 2*pi
+ assert_array_equal(unwrap([1,1+2*pi]),[1,1])
+ #check that unwrap maintans continuity
+ assert(all(diff(unwrap(rand(10)*100))<pi))
+
+
+class test_filterwindows(ScipyTestCase):
+ def check_hanning(self):
+ #check symmetry
+ w=hanning(10)
+ assert_array_almost_equal(w,flipud(w),7)
+ #check known value
+ assert_almost_equal(sum(w),4.500,4)
+
+ def check_hamming(self):
+ #check symmetry
+ w=hamming(10)
+ assert_array_almost_equal(w,flipud(w),7)
+ #check known value
+ assert_almost_equal(sum(w),4.9400,4)
+
+ def check_bartlett(self):
+ #check symmetry
+ w=bartlett(10)
+ assert_array_almost_equal(w,flipud(w),7)
+ #check known value
+ assert_almost_equal(sum(w),4.4444,4)
+
+ def check_blackman(self):
+ #check symmetry
+ w=blackman(10)
+ assert_array_almost_equal(w,flipud(w),7)
+ #check known value
+ assert_almost_equal(sum(w),3.7800,4)
+
+
+class test_trapz(ScipyTestCase):
+ def check_simple(self):
+ r=trapz(exp(-1.0/2*(arange(-10,10,.1))**2)/sqrt(2*pi),dx=0.1)
+ #check integral of normal equals 1
+ assert_almost_equal(sum(r),1,7)
+
+class test_sinc(ScipyTestCase):
+ def check_simple(self):
+ assert(sinc(0)==1)
+ w=sinc(linspace(-1,1,100))
+ #check symmetry
+ assert_array_almost_equal(w,flipud(w),7)
+
+class test_histogram(ScipyTestCase):
+ def check_simple(self):
+ n=100
+ v=rand(n)
+ (a,b)=histogram(v)
+ #check if the sum of the bins equals the number of samples
+ assert(sum(a)==n)
+ #check that the bin counts are evenly spaced when the data is from a linear function
+ (a,b)=histogram(linspace(0,10,100))
+ assert(all(a==10))
+
+
+
+
+
+def compare_results(res,desired):
+ for i in range(len(desired)):
+ assert_array_equal(res[i],desired[i])
+
+if __name__ == "__main__":
+ ScipyTest('scipy.base.function_base').run()
diff --git a/numpy/lib/tests/test_getlimits.py b/numpy/lib/tests/test_getlimits.py
new file mode 100644
index 000000000..99a6f5160
--- /dev/null
+++ b/numpy/lib/tests/test_getlimits.py
@@ -0,0 +1,38 @@
+""" Test functions for limits module.
+"""
+
+from scipy.testing import *
+set_package_path()
+import scipy.base;reload(scipy.base)
+from scipy.base.getlimits import finfo
+from scipy import single,double,longdouble
+restore_path()
+
+##################################################
+
+class test_python_float(ScipyTestCase):
+ def check_singleton(self):
+ ftype = finfo(float)
+ ftype2 = finfo(float)
+ assert_equal(id(ftype),id(ftype2))
+
+class test_single(ScipyTestCase):
+ def check_singleton(self):
+ ftype = finfo(single)
+ ftype2 = finfo(single)
+ assert_equal(id(ftype),id(ftype2))
+
+class test_double(ScipyTestCase):
+ def check_singleton(self):
+ ftype = finfo(double)
+ ftype2 = finfo(double)
+ assert_equal(id(ftype),id(ftype2))
+
+class test_longdouble(ScipyTestCase):
+ def check_singleton(self,level=2):
+ ftype = finfo(longdouble)
+ ftype2 = finfo(longdouble)
+ assert_equal(id(ftype),id(ftype2))
+
+if __name__ == "__main__":
+ ScipyTest().run()
diff --git a/numpy/lib/tests/test_index_tricks.py b/numpy/lib/tests/test_index_tricks.py
new file mode 100644
index 000000000..96e9dff84
--- /dev/null
+++ b/numpy/lib/tests/test_index_tricks.py
@@ -0,0 +1,53 @@
+
+from scipy.testing import *
+set_package_path()
+import scipy.base;reload(scipy.base)
+from scipy.base import *
+restore_path()
+
+class test_grid(ScipyTestCase):
+ def check_basic(self):
+ a = mgrid[-1:1:10j]
+ b = mgrid[-1:1:0.1]
+ assert(a.shape == (10,))
+ assert(b.shape == (20,))
+ assert(a[0] == -1)
+ assert_almost_equal(a[-1],1)
+ assert(b[0] == -1)
+ assert_almost_equal(b[1]-b[0],0.1,11)
+ assert_almost_equal(b[-1],b[0]+19*0.1,11)
+ assert_almost_equal(a[1]-a[0],2.0/9.0,11)
+
+ def check_nd(self):
+ c = mgrid[-1:1:10j,-2:2:10j]
+ d = mgrid[-1:1:0.1,-2:2:0.2]
+ assert(c.shape == (2,10,10))
+ assert(d.shape == (2,20,20))
+ assert_array_equal(c[0][0,:],-ones(10,'d'))
+ assert_array_equal(c[1][:,0],-2*ones(10,'d'))
+ assert_array_almost_equal(c[0][-1,:],ones(10,'d'),11)
+ assert_array_almost_equal(c[1][:,-1],2*ones(10,'d'),11)
+ assert_array_almost_equal(d[0,1,:]-d[0,0,:], 0.1*ones(20,'d'),11)
+ assert_array_almost_equal(d[1,:,1]-d[1,:,0], 0.2*ones(20,'d'),11)
+
+class test_concatenator(ScipyTestCase):
+ def check_1d(self):
+ assert_array_equal(r_[1,2,3,4,5,6],array([1,2,3,4,5,6]))
+ b = ones(5)
+ c = r_[b,0,0,b]
+ assert_array_equal(c,[1,1,1,1,1,0,0,1,1,1,1,1])
+
+ def check_2d(self):
+ b = rand(5,5)
+ c = rand(5,5)
+ d = r_[b,c,'1'] # append columns
+ assert(d.shape == (5,10))
+ assert_array_equal(d[:,:5],b)
+ assert_array_equal(d[:,5:],c)
+ d = r_[b,c]
+ assert(d.shape == (10,5))
+ assert_array_equal(d[:5,:],b)
+ assert_array_equal(d[5:,:],c)
+
+if __name__ == "__main__":
+ ScipyTest().run()
diff --git a/numpy/lib/tests/test_polynomial.py b/numpy/lib/tests/test_polynomial.py
new file mode 100644
index 000000000..51d4b5707
--- /dev/null
+++ b/numpy/lib/tests/test_polynomial.py
@@ -0,0 +1,83 @@
+"""
+>>> import scipy.base as nx
+>>> from scipy.base.polynomial import poly1d, polydiv
+
+>>> p = poly1d([1.,2,3])
+>>> p
+poly1d([ 1., 2., 3.])
+>>> print p
+ 2
+1 x + 2 x + 3
+>>> q = poly1d([3.,2,1])
+>>> q
+poly1d([ 3., 2., 1.])
+>>> print q
+ 2
+3 x + 2 x + 1
+
+>>> p(0)
+3.0
+>>> p(5)
+38.0
+>>> q(0)
+1.0
+>>> q(5)
+86.0
+
+>>> p * q
+poly1d([ 3., 8., 14., 8., 3.])
+>>> p / q
+(poly1d([ 0.33333333]), poly1d([ 1.33333333, 2.66666667]))
+>>> p + q
+poly1d([ 4., 4., 4.])
+>>> p - q
+poly1d([-2., 0., 2.])
+>>> p ** 4
+poly1d([ 1., 8., 36., 104., 214., 312., 324., 216., 81.])
+
+>>> p(q)
+poly1d([ 9., 12., 16., 8., 6.])
+>>> q(p)
+poly1d([ 3., 12., 32., 40., 34.])
+
+>>> nx.asarray(p)
+array([ 1., 2., 3.])
+>>> len(p)
+2
+
+>>> p[0], p[1], p[2], p[3]
+(3.0, 2.0, 1.0, 0)
+
+>>> p.integ()
+poly1d([ 0.33333333, 1. , 3. , 0. ])
+>>> p.integ(1)
+poly1d([ 0.33333333, 1. , 3. , 0. ])
+>>> p.integ(5)
+poly1d([ 0.00039683, 0.00277778, 0.025 , 0. , 0. ,
+ 0. , 0. , 0. ])
+>>> p.deriv()
+poly1d([ 2., 2.])
+>>> p.deriv(2)
+poly1d([ 2.])
+
+>>> q = poly1d([1.,2,3], variable='y')
+>>> print q
+ 2
+1 y + 2 y + 3
+>>> q = poly1d([1.,2,3], variable='lambda')
+>>> print q
+ 2
+1 lambda + 2 lambda + 3
+
+>>> polydiv(poly1d([1,0,-1]), poly1d([1,1]))
+(poly1d([ 1., -1.]), poly1d([ 0.]))
+"""
+
+from scipy.testing import *
+
+import doctest
+def test_suite(level=1):
+ return doctest.DocTestSuite()
+
+if __name__ == "__main__":
+ ScipyTest().run()
diff --git a/numpy/lib/tests/test_twodim_base.py b/numpy/lib/tests/test_twodim_base.py
new file mode 100644
index 000000000..b061d4a5d
--- /dev/null
+++ b/numpy/lib/tests/test_twodim_base.py
@@ -0,0 +1,134 @@
+""" Test functions for matrix module
+
+"""
+
+from scipy.testing import *
+set_package_path()
+import scipy.base;reload(scipy.base)
+from scipy.base import *
+restore_path()
+
+##################################################
+
+
+def get_mat(n):
+ data = arange(n)
+ data = add.outer(data,data)
+ return data
+
+class test_eye(ScipyTestCase):
+ def check_basic(self):
+ assert_equal(eye(4),array([[1,0,0,0],
+ [0,1,0,0],
+ [0,0,1,0],
+ [0,0,0,1]]))
+ assert_equal(eye(4,dtype='f'),array([[1,0,0,0],
+ [0,1,0,0],
+ [0,0,1,0],
+ [0,0,0,1]],'f'))
+ def check_diag(self):
+ assert_equal(eye(4,k=1),array([[0,1,0,0],
+ [0,0,1,0],
+ [0,0,0,1],
+ [0,0,0,0]]))
+ assert_equal(eye(4,k=-1),array([[0,0,0,0],
+ [1,0,0,0],
+ [0,1,0,0],
+ [0,0,1,0]]))
+ def check_2d(self):
+ assert_equal(eye(4,3),array([[1,0,0],
+ [0,1,0],
+ [0,0,1],
+ [0,0,0]]))
+ assert_equal(eye(3,4),array([[1,0,0,0],
+ [0,1,0,0],
+ [0,0,1,0]]))
+ def check_diag2d(self):
+ assert_equal(eye(3,4,k=2),array([[0,0,1,0],
+ [0,0,0,1],
+ [0,0,0,0]]))
+ assert_equal(eye(4,3,k=-2),array([[0,0,0],
+ [0,0,0],
+ [1,0,0],
+ [0,1,0]]))
+
+class test_diag(ScipyTestCase):
+ def check_vector(self):
+ vals = (100*arange(5)).astype('l')
+ b = zeros((5,5))
+ for k in range(5):
+ b[k,k] = vals[k]
+ assert_equal(diag(vals),b)
+ b = zeros((7,7))
+ c = b.copy()
+ for k in range(5):
+ b[k,k+2] = vals[k]
+ c[k+2,k] = vals[k]
+ assert_equal(diag(vals,k=2), b)
+ assert_equal(diag(vals,k=-2), c)
+
+ def check_matrix(self):
+ vals = (100*get_mat(5)+1).astype('l')
+ b = zeros((5,))
+ for k in range(5):
+ b[k] = vals[k,k]
+ assert_equal(diag(vals),b)
+ b = b*0
+ for k in range(3):
+ b[k] = vals[k,k+2]
+ assert_equal(diag(vals,2),b[:3])
+ for k in range(3):
+ b[k] = vals[k+2,k]
+ assert_equal(diag(vals,-2),b[:3])
+
+class test_fliplr(ScipyTestCase):
+ def check_basic(self):
+ self.failUnlessRaises(ValueError, fliplr, ones(4))
+ a = get_mat(4)
+ b = a[:,::-1]
+ assert_equal(fliplr(a),b)
+ a = [[0,1,2],
+ [3,4,5]]
+ b = [[2,1,0],
+ [5,4,3]]
+ assert_equal(fliplr(a),b)
+
+class test_flipud(ScipyTestCase):
+ def check_basic(self):
+ a = get_mat(4)
+ b = a[::-1,:]
+ assert_equal(flipud(a),b)
+ a = [[0,1,2],
+ [3,4,5]]
+ b = [[3,4,5],
+ [0,1,2]]
+ assert_equal(flipud(a),b)
+
+class test_rot90(ScipyTestCase):
+ def check_basic(self):
+ self.failUnlessRaises(ValueError, rot90, ones(4))
+
+ a = [[0,1,2],
+ [3,4,5]]
+ b1 = [[2,5],
+ [1,4],
+ [0,3]]
+ b2 = [[5,4,3],
+ [2,1,0]]
+ b3 = [[3,0],
+ [4,1],
+ [5,2]]
+ b4 = [[0,1,2],
+ [3,4,5]]
+
+ for k in range(-3,13,4):
+ assert_equal(rot90(a,k=k),b1)
+ for k in range(-2,13,4):
+ assert_equal(rot90(a,k=k),b2)
+ for k in range(-1,13,4):
+ assert_equal(rot90(a,k=k),b3)
+ for k in range(0,13,4):
+ assert_equal(rot90(a,k=k),b4)
+
+if __name__ == "__main__":
+ ScipyTest().run()
diff --git a/numpy/lib/twodim_base.py b/numpy/lib/twodim_base.py
new file mode 100644
index 000000000..b21532ea6
--- /dev/null
+++ b/numpy/lib/twodim_base.py
@@ -0,0 +1,123 @@
+""" Basic functions for manipulating 2d arrays
+
+"""
+
+__all__ = ['diag','eye','fliplr','flipud','rot90','tri','triu','tril',
+ 'vander']
+
+from numeric import *
+import sys
+
+def fliplr(m):
+ """ returns an array m with the rows preserved and columns flipped
+ in the left/right direction. Works on the first two dimensions of m.
+ """
+ m = asarray(m)
+ if m.ndim < 2:
+ raise ValueError, "Input must be >= 2-d."
+ return m[:, ::-1]
+
+def flipud(m):
+ """ returns an array with the columns preserved and rows flipped in
+ the up/down direction. Works on the first dimension of m.
+ """
+ m = asarray(m)
+ if m.ndim < 1:
+ raise ValueError, "Input must be >= 1-d."
+ return m[::-1]
+
+def rot90(m, k=1):
+ """ returns the array found by rotating m by k*90
+ degrees in the counterclockwise direction. Works on the first two
+ dimensions of m.
+ """
+ m = asarray(m)
+ if m.ndim < 2:
+ raise ValueError, "Input must >= 2-d."
+ k = k % 4
+ if k == 0: return m
+ elif k == 1: return fliplr(m).transpose()
+ elif k == 2: return fliplr(flipud(m))
+ else: return fliplr(m.transpose()) # k==3
+
+def eye(N, M=None, k=0, dtype=int_):
+ """ eye returns a N-by-M 2-d array where the k-th diagonal is all ones,
+ and everything else is zeros.
+ """
+ if M is None: M = N
+ m = equal(subtract.outer(arange(N), arange(M)),-k)
+ return m.astype(dtype)
+
+def diag(v, k=0):
+ """ returns the k-th diagonal if v is a array or returns a array
+ with v as the k-th diagonal if v is a vector.
+ """
+ v = asarray(v)
+ s = v.shape
+ if len(s)==1:
+ n = s[0]+abs(k)
+ res = zeros((n,n), v.dtype)
+ if (k>=0):
+ i = arange(0,n-k)
+ fi = i+k+i*n
+ else:
+ i = arange(0,n+k)
+ fi = i+(i-k)*n
+ res.flat[fi] = v
+ return res
+ elif len(s)==2:
+ N1,N2 = s
+ if k >= 0:
+ M = min(N1,N2-k)
+ i = arange(0,M)
+ fi = i+k+i*N2
+ else:
+ M = min(N1+k,N2)
+ i = arange(0,M)
+ fi = i + (i-k)*N2
+ return v.flat[fi]
+ else:
+ raise ValueError, "Input must be 1- or 2-d."
+
+
+def tri(N, M=None, k=0, dtype=int_):
+ """ returns a N-by-M array where all the diagonals starting from
+ lower left corner up to the k-th are all ones.
+ """
+ if M is None: M = N
+ m = greater_equal(subtract.outer(arange(N), arange(M)),-k)
+ return m.astype(dtype)
+
+def tril(m, k=0):
+ """ returns the elements on and below the k-th diagonal of m. k=0 is the
+ main diagonal, k > 0 is above and k < 0 is below the main diagonal.
+ """
+ m = asarray(m)
+ out = multiply(tri(m.shape[0], m.shape[1], k=k, dtype=m.dtype),m)
+ return out
+
+def triu(m, k=0):
+ """ returns the elements on and above the k-th diagonal of m. k=0 is the
+ main diagonal, k > 0 is above and k < 0 is below the main diagonal.
+ """
+ m = asarray(m)
+ out = multiply((1-tri(m.shape[0], m.shape[1], k-1, m.dtype)),m)
+ return out
+
+
+# borrowed from John Hunter and matplotlib
+def vander(x, N=None):
+ """
+ X = vander(x,N=None)
+
+ The Vandermonde matrix of vector x. The i-th column of X is the
+ the i-th power of x. N is the maximum power to compute; if N is
+ None it defaults to len(x).
+
+ """
+ x = asarray(x)
+ if N is None: N=len(x)
+ X = ones( (len(x),N), x.dtypechar)
+ for i in range(N-1):
+ X[:,i] = x**(N-i-1)
+ return X
diff --git a/numpy/lib/type_check.py b/numpy/lib/type_check.py
new file mode 100644
index 000000000..4c802ca86
--- /dev/null
+++ b/numpy/lib/type_check.py
@@ -0,0 +1,180 @@
+## Automatically adapted for scipy Sep 19, 2005 by convertcode.py
+
+__all__ = ['iscomplexobj','isrealobj','imag','iscomplex',
+ 'isscalar',
+ 'isreal','nan_to_num','real','real_if_close',
+ 'typename','asfarray','mintypecode','asscalar',
+ 'common_type']
+
+import numeric as _nx
+from numeric import ndarray, asarray, array, isinf, isnan, isfinite, signbit, \
+ ufunc, ScalarType, obj2dtype
+from ufunclike import isneginf, isposinf
+import umath
+
+_typecodes_by_elsize = 'GDFgdfQqLlIiHhBb?'
+
+def mintypecode(typechars,typeset='GDFgdf',default='d'):
+ """ Return a minimum data type character from typeset that
+ handles all typechars given
+
+ The returned type character must be the smallest size such that
+ an array of the returned type can handle the data from an array of
+ type t for each t in typechars (or if typechars is an array,
+ then its dtypechar).
+
+ If the typechars does not intersect with the typeset, then default
+ is returned.
+
+ If t in typechars is not a string then t=asarray(t).dtypechar is
+ applied.
+ """
+ typecodes = [(type(t) is type('') and t) or asarray(t).dtypechar\
+ for t in typechars]
+ intersection = [t for t in typecodes if t in typeset]
+ if not intersection:
+ return default
+ if 'F' in intersection and 'd' in intersection:
+ return 'D'
+ l = []
+ for t in intersection:
+ i = _typecodes_by_elsize.index(t)
+ l.append((i,t))
+ l.sort()
+ return l[0][1]
+
+def asfarray(a, dtype=_nx.float_):
+ """asfarray(a,dtype=None) returns a as a float array."""
+ dtype = _nx.obj2dtype(dtype)
+ if not issubclass(dtype, _nx.inexact):
+ dtype = _nx.float_
+ a = asarray(a,dtype=dtype)
+ return a
+
+def isscalar(num):
+ if isinstance(num, _nx.generic):
+ return True
+ else:
+ return type(num) in ScalarType
+
+def real(val):
+ return asarray(val).real
+
+def imag(val):
+ return asarray(val).imag
+
+def iscomplex(x):
+ return imag(x) != _nx.zeros_like(x)
+
+def isreal(x):
+ return imag(x) == _nx.zeros_like(x)
+
+def iscomplexobj(x):
+ return issubclass( asarray(x).dtype, _nx.complexfloating)
+
+def isrealobj(x):
+ return not issubclass( asarray(x).dtype, _nx.complexfloating)
+
+#-----------------------------------------------------------------------------
+
+def _getmaxmin(t):
+ import getlimits
+ f = getlimits.finfo(t)
+ return f.max, f.min
+
+def nan_to_num(x):
+ # mapping:
+ # NaN -> 0
+ # Inf -> limits.double_max
+ # -Inf -> limits.double_min
+ try:
+ t = x.dtype
+ except AttributeError:
+ t = obj2dtype(type(x))
+ if issubclass(t, _nx.complexfloating):
+ y = nan_to_num(x.real) + 1j * nan_to_num(x.imag)
+ elif issubclass(t, _nx.integer):
+ y = array(x)
+ else:
+ y = array(x)
+ if not y.shape:
+ y = array([x])
+ scalar = True
+ else:
+ scalar = False
+ are_inf = isposinf(y)
+ are_neg_inf = isneginf(y)
+ are_nan = isnan(y)
+ maxf, minf = _getmaxmin(y.dtype)
+ y[are_nan] = 0
+ y[are_inf] = maxf
+ y[are_neg_inf] = minf
+ if scalar:
+ y = y[0]
+ return y
+
+#-----------------------------------------------------------------------------
+
+def real_if_close(a,tol=100):
+ a = asarray(a)
+ if a.dtypechar not in 'FDG':
+ return a
+ if tol > 1:
+ import getlimits
+ f = getlimits.finfo(a.dtype)
+ tol = f.eps * tol
+ if _nx.allclose(a.imag, 0, atol=tol):
+ a = a.real
+ return a
+
+
+def asscalar(a):
+ return a.item()
+
+#-----------------------------------------------------------------------------
+
+_namefromtype = {'S1' : 'character',
+ '?' : 'bool',
+ 'b' : 'signed char',
+ 'B' : 'unsigned char',
+ 'h' : 'short',
+ 'H' : 'unsigned short',
+ 'i' : 'integer',
+ 'I' : 'unsigned integer',
+ 'l' : 'long integer',
+ 'L' : 'unsigned long integer',
+ 'q' : 'long long integer',
+ 'Q' : 'unsigned long long integer',
+ 'f' : 'single precision',
+ 'd' : 'double precision',
+ 'g' : 'long precision',
+ 'F' : 'complex single precision',
+ 'D' : 'complex double precision',
+ 'G' : 'complex long double precision',
+ 'S' : 'string',
+ 'U' : 'unicode',
+ 'V' : 'void',
+ 'O' : 'object'
+ }
+
+def typename(char):
+ """Return an english description for the given data type character.
+ """
+ return _namefromtype[char]
+
+#-----------------------------------------------------------------------------
+
+#determine the "minimum common type code" for a group of arrays.
+array_kind = {'i':0, 'l': 0, 'f': 0, 'd': 0, 'g':0, 'F': 1, 'D': 1, 'G':1}
+array_precision = {'i': 1, 'l': 1,
+ 'f': 0, 'd': 1, 'g':2,
+ 'F': 0, 'D': 1, 'G':2}
+array_type = [['f', 'd', 'g'], ['F', 'D', 'G']]
+def common_type(*arrays):
+ kind = 0
+ precision = 0
+ for a in arrays:
+ t = a.dtypechar
+ kind = max(kind, array_kind[t])
+ precision = max(precision, array_precision[t])
+ return array_type[kind][precision]
diff --git a/numpy/lib/ufunclike.py b/numpy/lib/ufunclike.py
new file mode 100644
index 000000000..7e8d44c7d
--- /dev/null
+++ b/numpy/lib/ufunclike.py
@@ -0,0 +1,77 @@
+"""
+Module of functions that are like ufuncs in acting on arrays and optionally
+storing results in an output array.
+"""
+__all__ = ['fix', 'isneginf', 'isposinf', 'sign', 'log2']
+
+import numeric as nx
+from numeric import asarray, empty, empty_like, isinf, signbit, zeros
+import umath
+
+def fix(x, y=None):
+ """ Round x to nearest integer towards zero.
+ """
+ x = asarray(x)
+ if y is None:
+ y = nx.floor(x)
+ else:
+ nx.floor(x, y)
+ if x.ndim == 0:
+ if (x<0):
+ y += 1
+ else:
+ y[x<0] = y[x<0]+1
+ return y
+
+def isposinf(x, y=None):
+ """Return a boolean array y with y[i] True for x[i] = +Inf.
+
+ If y is an array, the result replaces the contents of y.
+ """
+ if y is None:
+ y = empty(x.shape, dtype=nx.bool_)
+ umath.logical_and(isinf(x), ~signbit(x), y)
+ return y
+
+def isneginf(x, y=None):
+ """Return a boolean array y with y[i] True for x[i] = -Inf.
+
+ If y is an array, the result replaces the contents of y.
+ """
+ if y is None:
+ y = empty(x.shape, dtype=nx.bool_)
+ umath.logical_and(isinf(x), signbit(x), y)
+ return y
+
+def sign(x, y=None):
+ """sign(x) gives an array with shape of x with elexents defined by sign
+ function: where x is less than 0 return -1, where x greater than 0, a=1,
+ elsewhere a=0.
+ """
+ x = asarray(x)
+ if y is None:
+ y = zeros(x.shape, dtype=nx.int_)
+ if x.ndim == 0:
+ if x < 0:
+ y -= 1
+ elif x > 0:
+ y += 1
+ else:
+ y[x<0] = -1
+ y[x>0] = 1
+ return y
+
+_log2 = umath.log(2)
+def log2(x, y=None):
+ """Returns the base 2 logarithm of x
+
+ If y is an array, the result replaces the contents of y.
+ """
+ x = asarray(x)
+ if y is None:
+ y = umath.log(x)
+ else:
+ umath.log(x, y)
+ y /= _log2
+ return y
+
diff --git a/numpy/lib/utils.py b/numpy/lib/utils.py
new file mode 100644
index 000000000..19fb18d4d
--- /dev/null
+++ b/numpy/lib/utils.py
@@ -0,0 +1,28 @@
+from numerictypes import obj2dtype
+
+__all__ = ['issubclass_', 'get_scipy_include', 'issubdtype']
+
+def issubclass_(arg1, arg2):
+ try:
+ return issubclass(arg1, arg2)
+ except TypeError:
+ return False
+
+def issubdtype(arg1, arg2):
+ return issubclass(obj2dtype(arg1), obj2dtype(arg2))
+
+def get_scipy_include():
+ """Return the directory in the package that contains the scipy/*.h header
+ files.
+
+ Extension modules that need to compile against scipy.base should use this
+ function to locate the appropriate include directory. Using distutils:
+
+ import scipy
+ Extension('extension_name', ...
+ include_dirs=[scipy.get_scipy_include()])
+ """
+ from scipy.distutils.misc_util import get_scipy_include_dirs
+ include_dirs = get_scipy_include_dirs()
+ assert len(include_dirs)==1,`include_dirs`
+ return include_dirs[0]
View Lorry Controller Status