summaryrefslogtreecommitdiff
path: root/function_base.py
diff options
context:
space:
mode:
Diffstat (limited to 'function_base.py')
-rw-r--r--function_base.py1492
1 files changed, 1492 insertions, 0 deletions
diff --git a/function_base.py b/function_base.py
new file mode 100644
index 000000000..711c609e8
--- /dev/null
+++ b/function_base.py
@@ -0,0 +1,1492 @@
+__docformat__ = "restructuredtext en"
+__all__ = ['logspace', 'linspace',
+ 'select', 'piecewise', 'trim_zeros',
+ 'copy', 'iterable',
+ 'diff', 'gradient', 'angle', 'unwrap', 'sort_complex', 'disp',
+ 'unique', 'extract', 'place', 'nansum', 'nanmax', 'nanargmax',
+ 'nanargmin', 'nanmin', 'vectorize', 'asarray_chkfinite', 'average',
+ 'histogram', 'histogramdd', 'bincount', 'digitize', 'cov',
+ 'corrcoef', 'msort', 'median', 'sinc', 'hamming', 'hanning',
+ 'bartlett', 'blackman', 'kaiser', 'trapz', 'i0', 'add_newdoc',
+ 'add_docstring', 'meshgrid', 'delete', 'insert', 'append',
+ 'interp'
+ ]
+
+import types
+import numpy.core.numeric as _nx
+from numpy.core.numeric import ones, zeros, arange, concatenate, array, \
+ asarray, asanyarray, empty, empty_like, asanyarray, ndarray, around
+from numpy.core.numeric import ScalarType, dot, where, newaxis, intp, \
+ integer, isscalar
+from numpy.core.umath import pi, multiply, add, arctan2, \
+ frompyfunc, isnan, cos, less_equal, sqrt, sin, mod, exp, log10
+from numpy.core.fromnumeric import ravel, nonzero, choose, sort
+from numpy.core.numerictypes import typecodes
+from numpy.lib.shape_base import atleast_1d, atleast_2d
+from numpy.lib.twodim_base import diag
+from _compiled_base import _insert, add_docstring
+from _compiled_base import digitize, bincount, interp
+from arraysetops import setdiff1d
+
+#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 (seq, step_value), where step_value used.
+
+ :Parameters:
+ start : {float}
+ The value the sequence starts at.
+ stop : {float}
+ The value the sequence stops at. If ``endpoint`` is false, then
+ this is not included in the sequence. Otherwise it is
+ guaranteed to be the last value.
+ num : {integer}
+ Number of samples to generate. Default is 50.
+ endpoint : {boolean}
+ If true, ``stop`` is the last sample. Otherwise, it is not
+ included. Default is true.
+ retstep : {boolean}
+ If true, return ``(samples, step)``, where ``step`` is the
+ spacing used in generating the samples.
+
+ :Returns:
+ samples : {array}
+ ``num`` equally spaced samples from the range [start, stop]
+ or [start, stop).
+ step : {float} (Only if ``retstep`` is true)
+ Size of spacing between samples.
+
+ :See Also:
+ `arange` : Similiar to linspace, however, when used with
+ a float endpoint, that endpoint may or may not be included.
+ `logspace`
+ """
+ num = int(num)
+ if num <= 0:
+ return array([], float)
+ if endpoint:
+ if num == 1:
+ return array([float(start)])
+ step = (stop-start)/float((num-1))
+ y = _nx.arange(0, num) * step + start
+ y[-1] = stop
+ 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 base**start to
+ base**stop. If endpoint=True, then last number is base**stop
+ """
+ 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):
+ """Compute the histogram from a set of data.
+
+ Parameters:
+
+ a : array
+ The data to histogram. n-D arrays will be flattened.
+
+ bins : int or sequence of floats
+ If an int, then the number of equal-width bins in the given range.
+ Otherwise, a sequence of the lower bound of each bin.
+
+ range : (float, float)
+ The lower and upper range of the bins. If not provided, then
+ (a.min(), a.max()) is used. Values outside of this range are
+ allocated to the closest bin.
+
+ normed : bool
+ If False, the result array will contain the number of samples in
+ each bin. If True, the result array is the value of the
+ probability *density* function at the bin normalized such that the
+ *integral* over the range is 1. Note that the sum of all of the
+ histogram values will not usually be 1; it is not a probability
+ *mass* function.
+
+ Returns:
+
+ hist : array
+ The values of the histogram. See `normed` for a description of the
+ possible semantics.
+
+ lower_edges : float array
+ The lower edges of each bin.
+
+ SeeAlso:
+
+ histogramdd
+
+ """
+ a = asarray(a).ravel()
+
+ if (range is not None):
+ mn, mx = range
+ if (mn > mx):
+ raise AttributeError, 'max must be larger than min in range parameter.'
+
+ 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)
+ else:
+ if(any(bins[1:]-bins[:-1] < 0)):
+ raise AttributeError, 'bins must increase monotonically.'
+
+ # best block size probably depends on processor cache size
+ block = 65536
+ n = sort(a[:block]).searchsorted(bins)
+ for i in xrange(block, a.size, block):
+ n += sort(a[i:i+block]).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 histogramdd(sample, bins=10, range=None, normed=False, weights=None):
+ """histogramdd(sample, bins=10, range=None, normed=False, weights=None)
+
+ Return the N-dimensional histogram of the sample.
+
+ Parameters:
+
+ sample : sequence or array
+ A sequence containing N arrays or an NxM array. Input data.
+
+ bins : sequence or scalar
+ A sequence of edge arrays, a sequence of bin counts, or a scalar
+ which is the bin count for all dimensions. Default is 10.
+
+ range : sequence
+ A sequence of lower and upper bin edges. Default is [min, max].
+
+ normed : boolean
+ If False, return the number of samples in each bin, if True,
+ returns the density.
+
+ weights : array
+ Array of weights. The weights are normed only if normed is True.
+ Should the sum of the weights not equal N, the total bin count will
+ not be equal to the number of samples.
+
+ Returns:
+
+ hist : array
+ Histogram array.
+
+ edges : list
+ List of arrays defining the lower bin edges.
+
+ SeeAlso:
+
+ histogram
+
+ Example
+
+ >>> x = random.randn(100,3)
+ >>> hist3d, edges = histogramdd(x, bins = (5, 6, 7))
+
+ """
+
+ try:
+ # Sample is an ND-array.
+ N, D = sample.shape
+ except (AttributeError, ValueError):
+ # Sample is a sequence of 1D arrays.
+ sample = atleast_2d(sample).T
+ N, D = sample.shape
+
+ nbin = empty(D, int)
+ edges = D*[None]
+ dedges = D*[None]
+ if weights is not None:
+ weights = asarray(weights)
+
+ try:
+ M = len(bins)
+ if M != D:
+ raise AttributeError, 'The dimension of bins must be a equal to the dimension of the sample x.'
+ except TypeError:
+ bins = D*[bins]
+
+ # Select range for each dimension
+ # Used only if number of bins is given.
+ if range is None:
+ smin = atleast_1d(array(sample.min(0), float))
+ smax = atleast_1d(array(sample.max(0), float))
+ else:
+ smin = zeros(D)
+ smax = zeros(D)
+ for i in arange(D):
+ smin[i], smax[i] = range[i]
+
+ # Make sure the bins have a finite width.
+ for i in arange(len(smin)):
+ if smin[i] == smax[i]:
+ smin[i] = smin[i] - .5
+ smax[i] = smax[i] + .5
+
+ # Create edge arrays
+ for i in arange(D):
+ if isscalar(bins[i]):
+ nbin[i] = bins[i] + 2 # +2 for outlier bins
+ edges[i] = linspace(smin[i], smax[i], nbin[i]-1)
+ else:
+ edges[i] = asarray(bins[i], float)
+ nbin[i] = len(edges[i])+1 # +1 for outlier bins
+ dedges[i] = diff(edges[i])
+
+ nbin = asarray(nbin)
+
+ # Compute the bin number each sample falls into.
+ Ncount = {}
+ for i in arange(D):
+ Ncount[i] = digitize(sample[:,i], edges[i])
+
+ # Using digitize, values that fall on an edge are put in the right bin.
+ # For the rightmost bin, we want values equal to the right
+ # edge to be counted in the last bin, and not as an outlier.
+ outliers = zeros(N, int)
+ for i in arange(D):
+ # Rounding precision
+ decimal = int(-log10(dedges[i].min())) +6
+ # Find which points are on the rightmost edge.
+ on_edge = where(around(sample[:,i], decimal) == around(edges[i][-1], decimal))[0]
+ # Shift these points one bin to the left.
+ Ncount[i][on_edge] -= 1
+
+ # Flattened histogram matrix (1D)
+ hist = zeros(nbin.prod(), float)
+
+ # Compute the sample indices in the flattened histogram matrix.
+ ni = nbin.argsort()
+ shape = []
+ xy = zeros(N, int)
+ for i in arange(0, D-1):
+ xy += Ncount[ni[i]] * nbin[ni[i+1:]].prod()
+ xy += Ncount[ni[-1]]
+
+ # Compute the number of repetitions in xy and assign it to the flattened histmat.
+ if len(xy) == 0:
+ return zeros(nbin-2, int), edges
+
+ flatcount = bincount(xy, weights)
+ a = arange(len(flatcount))
+ hist[a] = flatcount
+
+ # Shape into a proper matrix
+ hist = hist.reshape(sort(nbin))
+ for i in arange(nbin.size):
+ j = ni[i]
+ hist = hist.swapaxes(i,j)
+ ni[i],ni[j] = ni[j],ni[i]
+
+ # Remove outliers (indices 0 and -1 for each dimension).
+ core = D*[slice(1,-1)]
+ hist = hist[core]
+
+ # Normalize if normed is True
+ if normed:
+ s = hist.sum()
+ for i in arange(D):
+ shape = ones(D, int)
+ shape[i] = nbin[i]-2
+ hist = hist / dedges[i].reshape(shape)
+ hist /= s
+
+ return hist, edges
+
+
+def average(a, axis=None, 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) and to
+
+ a.sum(axis) / size(a, axis)
+
+ If weights are given, result is:
+ sum(a * weights,axis) / sum(weights,axis),
+ 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, ndarray):
+ 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.dtype.char in 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.
+ A constant may be used in funclist for a function that returns a
+ constant (e.g. val and lambda x: val are equivalent in a funclist).
+
+ 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
+ |--
+
+ """
+ x = asanyarray(x)
+ 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[k]
+ 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 in choicelist,
+ depending on the list of conditions.
+
+ :Parameters:
+ condlist : list of N boolean arrays of length M
+ The conditions C_0 through C_(N-1) which determine
+ from which vector the output elements are taken.
+ choicelist : list of N arrays of length M
+ Th vectors V_0 through V_(N-1), from which the output
+ elements are chosen.
+
+ :Returns:
+ output : 1-dimensional array of length M
+ The output at position m is the m-th element of the first
+ vector V_n for which C_n[m] is non-zero. Note that the
+ output depends on the order of conditions, since the
+ first satisfied condition is used.
+
+ Equivalent to:
+
+ output = []
+ for m in range(M):
+ output += [V[m] for V,C in zip(values,cond) if C[m]]
+ or [default]
+
+ """
+ n = len(condlist)
+ n2 = len(choicelist)
+ if n2 != n:
+ raise ValueError, "list of cases must be same length as list of conditions"
+ choicelist = [default] + choicelist
+ 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])
+ if type(S) in ScalarType:
+ S = S*ones(asarray(pfac).shape, type(S))
+ else:
+ S = S*ones(asarray(pfac).shape, S.dtype)
+ 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
+
+ outvals = []
+
+ # create slice objects --- initially all are [:, :, ..., :]
+ slice1 = [slice(None)]*N
+ slice2 = [slice(None)]*N
+ slice3 = [slice(None)]*N
+
+ otype = f.dtype.char
+ 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.dtype.char)
+ 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 ' + repr(n)
+ a = asanyarray(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]
+
+try:
+ add_docstring(digitize,
+r"""digitize(x,bins)
+
+Return the index of the bin to which each value of x belongs.
+
+Each index i returned is 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 the bins 0 or len(bins) is returned as appropriate.
+
+""")
+except RuntimeError:
+ pass
+
+try:
+ add_docstring(bincount,
+r"""bincount(x,weights=None)
+
+Return the number of occurrences of each value in x.
+
+x must be a list of non-negative integers. The output, b[i],
+represents the number of times that i is found in x. If weights
+is specified, every occurrence of i at a position p contributes
+weights[p] instead of 1.
+
+See also: histogram, digitize, unique.
+
+""")
+except RuntimeError:
+ pass
+
+try:
+ add_docstring(add_docstring,
+r"""docstring(obj, docstring)
+
+Add a docstring to a built-in obj if possible.
+If the obj already has a docstring raise a RuntimeError
+If this routine does not know how to add a docstring to the object
+raise a TypeError
+
+""")
+except RuntimeError:
+ pass
+
+try:
+ add_docstring(interp,
+r"""interp(x, xp, fp, left=None, right=None)
+
+Return the value of a piecewise-linear function at each value in x.
+
+The piecewise-linear function, f, is defined by the known data-points fp=f(xp).
+The xp points must be sorted in increasing order but this is not checked.
+
+For values of x < xp[0] return the value given by left. If left is None, then
+return fp[0].
+For values of x > xp[-1] return the value given by right. If right is None, then
+return fp[-1].
+"""
+ )
+except RuntimeError:
+ pass
+
+
+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.type, _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.type, _nx.complexfloating):
+ if b.dtype.char in 'bhBH':
+ return b.astype('F')
+ elif b.dtype.char == '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 numpy
+ >>> a = array((0, 0, 0, 1, 2, 3, 2, 1, 0))
+ >>> numpy.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]
+
+import sys
+if sys.hexversion < 0x2040000:
+ from sets import Set as set
+
+def unique(x):
+ """Return sorted unique items from an array or sequence.
+
+ Example:
+ >>> unique([5,2,4,0,4,4,2,2,1])
+ array([0, 1, 2, 4, 5])
+
+ """
+ try:
+ tmp = x.flatten()
+ if tmp.size == 0:
+ return tmp
+ tmp.sort()
+ idx = concatenate(([True],tmp[1:]!=tmp[:-1]))
+ return tmp[idx]
+ except AttributeError:
+ items = list(set(x))
+ items.sort()
+ return asarray(items)
+
+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))[0])
+
+def place(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=None):
+ """Sum the array over the given axis, treating NaNs as 0.
+ """
+ y = array(a,subok=True)
+ if not issubclass(y.dtype.type, _nx.integer):
+ y[isnan(a)] = 0
+ return y.sum(axis)
+
+def nanmin(a, axis=None):
+ """Find the minimium over the given axis, ignoring NaNs.
+ """
+ y = array(a,subok=True)
+ if not issubclass(y.dtype.type, _nx.integer):
+ y[isnan(a)] = _nx.inf
+ return y.min(axis)
+
+def nanargmin(a, axis=None):
+ """Find the indices of the minimium over the given axis ignoring NaNs.
+ """
+ y = array(a, subok=True)
+ if not issubclass(y.dtype.type, _nx.integer):
+ y[isnan(a)] = _nx.inf
+ return y.argmin(axis)
+
+def nanmax(a, axis=None):
+ """Find the maximum over the given axis ignoring NaNs.
+ """
+ y = array(a, subok=True)
+ if not issubclass(y.dtype.type, _nx.integer):
+ y[isnan(a)] = -_nx.inf
+ return y.max(axis)
+
+def nanargmax(a, axis=None):
+ """Find the maximum over the given axis ignoring NaNs.
+ """
+ y = array(a,subok=True)
+ if not issubclass(y.dtype.type, _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
+
+# return number of input arguments and
+# number of default arguments
+import re
+def _get_nargs(obj):
+ if not callable(obj):
+ raise TypeError, "Object is not callable."
+ if hasattr(obj,'func_code'):
+ fcode = obj.func_code
+ nargs = fcode.co_argcount
+ if obj.func_defaults is not None:
+ ndefaults = len(obj.func_defaults)
+ else:
+ ndefaults = 0
+ if isinstance(obj, types.MethodType):
+ nargs -= 1
+ return nargs, ndefaults
+ terr = re.compile(r'.*? takes exactly (?P<exargs>\d+) argument(s|) \((?P<gargs>\d+) given\)')
+ try:
+ obj()
+ return 0, 0
+ except TypeError, msg:
+ m = terr.match(str(msg))
+ if m:
+ nargs = int(m.group('exargs'))
+ ndefaults = int(m.group('gargs'))
+ if isinstance(obj, types.MethodType):
+ nargs -= 1
+ return nargs, ndefaults
+ raise ValueError, 'failed to determine the number of arguments for %s' % (obj)
+
+
+class vectorize(object):
+ """
+ vectorize(somefunction, otypes=None, doc=None)
+ Generalized Function class.
+
+ Description:
+
+ Define a vectorized function which takes nested sequence
+ of objects or numpy arrays as inputs and returns a
+ numpy 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 numpy.
+
+ Data-type of output of vectorized is determined by calling the function
+ with the first element of the input. This can be avoided by specifying
+ the otypes argument as either a string of typecode characters or a list
+ of data-types specifiers. There should be one data-type specifier for
+ each output.
+
+ 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):
+ self.thefunc = pyfunc
+ self.ufunc = None
+ nin, ndefault = _get_nargs(pyfunc)
+ if nin == 0 and ndefault == 0:
+ self.nin = None
+ self.nin_wo_defaults = None
+ else:
+ self.nin = nin
+ self.nin_wo_defaults = nin - ndefault
+ self.nout = None
+ if doc is None:
+ self.__doc__ = pyfunc.__doc__
+ else:
+ self.__doc__ = doc
+ if isinstance(otypes, types.StringType):
+ self.otypes = otypes
+ for char in self.otypes:
+ if char not in typecodes['All']:
+ raise ValueError, "invalid otype specified"
+ elif iterable(otypes):
+ self.otypes = ''.join([_nx.dtype(x).char for x in otypes])
+ else:
+ raise ValueError, "output types must be a string of typecode characters or a list of data-types"
+ 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 self.nin:
+ if (nargs > self.nin) or (nargs < self.nin_wo_defaults):
+ raise ValueError, "mismatch between python function inputs"\
+ " and received arguments"
+
+ # we need a new ufunc if this is being called with more arguments.
+ if (self.lastcallargs != nargs):
+ self.lastcallargs = nargs
+ self.ufunc = None
+ self.nout = None
+
+ 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]).dtype.char)
+ self.otypes = ''.join(otypes)
+
+ # Create ufunc if not already created
+ if (self.ufunc is None):
+ self.ufunc = frompyfunc(self.thefunc, nargs, self.nout)
+
+ # Convert to object arrays first
+ newargs = [array(arg,copy=False,subok=True,dtype=object) for arg in args]
+ if self.nout == 1:
+ _res = array(self.ufunc(*newargs),copy=False,
+ subok=True,dtype=self.otypes[0])
+ else:
+ _res = tuple([array(x,copy=False,subok=True,dtype=c) \
+ for x, c in zip(self.ufunc(*newargs), self.otypes)])
+ return _res
+
+def cov(m, y=None, rowvar=1, bias=0):
+ """Estimate the covariance matrix.
+
+ If m is a vector, return the variance. For matrices return the
+ covariance matrix.
+
+ If y is given it is treated as an additional (set of)
+ variable(s).
+
+ 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 non-zero (default), then each row is a variable with
+ observations in the columns, otherwise each column
+ is a variable and the observations are in the rows.
+ """
+
+ X = array(m, ndmin=2, dtype=float)
+ if X.shape[0] == 1:
+ rowvar = 1
+ if rowvar:
+ axis = 0
+ tup = (slice(None),newaxis)
+ else:
+ axis = 1
+ tup = (newaxis, slice(None))
+
+
+ if y is not None:
+ y = array(y, copy=False, ndmin=2, dtype=float)
+ X = concatenate((X,y),axis)
+
+ X -= X.mean(axis=1-axis)[tup]
+ if rowvar:
+ N = X.shape[1]
+ else:
+ N = X.shape[0]
+
+ if bias:
+ fact = N*1.0
+ else:
+ fact = N-1.0
+
+ if not rowvar:
+ return (dot(X.T, X.conj()) / fact).squeeze()
+ else:
+ return (dot(X, X.T.conj()) / fact).squeeze()
+
+def corrcoef(x, y=None, rowvar=1, bias=0):
+ """The correlation coefficients
+ """
+ c = cov(x, y, rowvar, bias)
+ try:
+ d = diag(c)
+ except ValueError: # scalar covariance
+ return 1
+ return c/sqrt(multiply.outer(d,d))
+
+def blackman(M):
+ """blackman(M) returns the M-point Blackman window.
+ """
+ if M < 1:
+ return array([])
+ if M == 1:
+ return ones(1, float)
+ 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.
+ """
+ if M < 1:
+ return array([])
+ if M == 1:
+ return ones(1, float)
+ 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.
+ """
+ if M < 1:
+ return array([])
+ if M == 1:
+ return ones(1, float)
+ 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.
+ """
+ if M < 1:
+ return array([])
+ if M == 1:
+ return ones(1,float)
+ n = arange(0,M)
+ return 0.54-0.46*cos(2.0*pi*n/(M-1))
+
+## Code from cephes for i0
+
+_i0A = [
+-4.41534164647933937950E-18,
+ 3.33079451882223809783E-17,
+-2.43127984654795469359E-16,
+ 1.71539128555513303061E-15,
+-1.16853328779934516808E-14,
+ 7.67618549860493561688E-14,
+-4.85644678311192946090E-13,
+ 2.95505266312963983461E-12,
+-1.72682629144155570723E-11,
+ 9.67580903537323691224E-11,
+-5.18979560163526290666E-10,
+ 2.65982372468238665035E-9,
+-1.30002500998624804212E-8,
+ 6.04699502254191894932E-8,
+-2.67079385394061173391E-7,
+ 1.11738753912010371815E-6,
+-4.41673835845875056359E-6,
+ 1.64484480707288970893E-5,
+-5.75419501008210370398E-5,
+ 1.88502885095841655729E-4,
+-5.76375574538582365885E-4,
+ 1.63947561694133579842E-3,
+-4.32430999505057594430E-3,
+ 1.05464603945949983183E-2,
+-2.37374148058994688156E-2,
+ 4.93052842396707084878E-2,
+-9.49010970480476444210E-2,
+ 1.71620901522208775349E-1,
+-3.04682672343198398683E-1,
+ 6.76795274409476084995E-1]
+
+_i0B = [
+-7.23318048787475395456E-18,
+-4.83050448594418207126E-18,
+ 4.46562142029675999901E-17,
+ 3.46122286769746109310E-17,
+-2.82762398051658348494E-16,
+-3.42548561967721913462E-16,
+ 1.77256013305652638360E-15,
+ 3.81168066935262242075E-15,
+-9.55484669882830764870E-15,
+-4.15056934728722208663E-14,
+ 1.54008621752140982691E-14,
+ 3.85277838274214270114E-13,
+ 7.18012445138366623367E-13,
+-1.79417853150680611778E-12,
+-1.32158118404477131188E-11,
+-3.14991652796324136454E-11,
+ 1.18891471078464383424E-11,
+ 4.94060238822496958910E-10,
+ 3.39623202570838634515E-9,
+ 2.26666899049817806459E-8,
+ 2.04891858946906374183E-7,
+ 2.89137052083475648297E-6,
+ 6.88975834691682398426E-5,
+ 3.36911647825569408990E-3,
+ 8.04490411014108831608E-1]
+
+def _chbevl(x, vals):
+ b0 = vals[0]
+ b1 = 0.0
+
+ for i in xrange(1,len(vals)):
+ b2 = b1
+ b1 = b0
+ b0 = x*b1 - b2 + vals[i]
+
+ return 0.5*(b0 - b2)
+
+def _i0_1(x):
+ return exp(x) * _chbevl(x/2.0-2, _i0A)
+
+def _i0_2(x):
+ return exp(x) * _chbevl(32.0/x - 2.0, _i0B) / sqrt(x)
+
+def i0(x):
+ x = atleast_1d(x).copy()
+ y = empty_like(x)
+ ind = (x<0)
+ x[ind] = -x[ind]
+ ind = (x<=8.0)
+ y[ind] = _i0_1(x[ind])
+ ind2 = ~ind
+ y[ind2] = _i0_2(x[ind2])
+ return y.squeeze()
+
+## End of cephes code for i0
+
+def kaiser(M,beta):
+ """kaiser(M, beta) returns a Kaiser window of length M with shape parameter
+ beta.
+ """
+ from numpy.dual 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,subok=True,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)
+ index = int(sorted.shape[0]/2)
+ if sorted.shape[0] % 2 == 1:
+ return sorted[index]
+ else:
+ 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)
+
+#always succeed
+def add_newdoc(place, obj, doc):
+ """Adds documentation to obj which is in module place.
+
+ If doc is a string add it to obj as a docstring
+
+ If doc is a tuple, then the first element is interpreted as
+ an attribute of obj and the second as the docstring
+ (method, docstring)
+
+ If doc is a list, then each element of the list should be a
+ sequence of length two --> [(method1, docstring1),
+ (method2, docstring2), ...]
+
+ This routine never raises an error.
+ """
+ try:
+ new = {}
+ exec 'from %s import %s' % (place, obj) in new
+ if isinstance(doc, str):
+ add_docstring(new[obj], doc.strip())
+ elif isinstance(doc, tuple):
+ add_docstring(getattr(new[obj], doc[0]), doc[1].strip())
+ elif isinstance(doc, list):
+ for val in doc:
+ add_docstring(getattr(new[obj], val[0]), val[1].strip())
+ except:
+ pass
+
+
+# From matplotlib
+def meshgrid(x,y):
+ """
+ For vectors x, y with lengths Nx=len(x) and Ny=len(y), return X, Y
+ where X and Y are (Ny, Nx) shaped arrays with the elements of x
+ and y repeated to fill the matrix
+
+ EG,
+
+ [X, Y] = meshgrid([1,2,3], [4,5,6,7])
+
+ X =
+ 1 2 3
+ 1 2 3
+ 1 2 3
+ 1 2 3
+
+
+ Y =
+ 4 4 4
+ 5 5 5
+ 6 6 6
+ 7 7 7
+ """
+ x = asarray(x)
+ y = asarray(y)
+ numRows, numCols = len(y), len(x) # yes, reversed
+ x = x.reshape(1,numCols)
+ X = x.repeat(numRows, axis=0)
+
+ y = y.reshape(numRows,1)
+ Y = y.repeat(numCols, axis=1)
+ return X, Y
+
+def delete(arr, obj, axis=None):
+ """Return a new array with sub-arrays along an axis deleted.
+
+ Return a new array with the sub-arrays (i.e. rows or columns)
+ deleted along the given axis as specified by obj
+
+ obj may be a slice_object (s_[3:5:2]) or an integer
+ or an array of integers indicated which sub-arrays to
+ remove.
+
+ If axis is None, then ravel the array first.
+
+ Example:
+ >>> arr = [[3,4,5],
+ ... [1,2,3],
+ ... [6,7,8]]
+
+ >>> delete(arr, 1, 1)
+ array([[3, 5],
+ [1, 3],
+ [6, 8]])
+ >>> delete(arr, 1, 0)
+ array([[3, 4, 5],
+ [6, 7, 8]])
+ """
+ wrap = None
+ if type(arr) is not ndarray:
+ try:
+ wrap = arr.__array_wrap__
+ except AttributeError:
+ pass
+
+
+ arr = asarray(arr)
+ ndim = arr.ndim
+ if axis is None:
+ if ndim != 1:
+ arr = arr.ravel()
+ ndim = arr.ndim;
+ axis = ndim-1;
+ if ndim == 0:
+ if wrap:
+ return wrap(arr)
+ else:
+ return arr.copy()
+ slobj = [slice(None)]*ndim
+ N = arr.shape[axis]
+ newshape = list(arr.shape)
+ if isinstance(obj, (int, long, integer)):
+ if (obj < 0): obj += N
+ if (obj < 0 or obj >=N):
+ raise ValueError, "invalid entry"
+ newshape[axis]-=1;
+ new = empty(newshape, arr.dtype, arr.flags.fnc)
+ slobj[axis] = slice(None, obj)
+ new[slobj] = arr[slobj]
+ slobj[axis] = slice(obj,None)
+ slobj2 = [slice(None)]*ndim
+ slobj2[axis] = slice(obj+1,None)
+ new[slobj] = arr[slobj2]
+ elif isinstance(obj, slice):
+ start, stop, step = obj.indices(N)
+ numtodel = len(xrange(start, stop, step))
+ if numtodel <= 0:
+ if wrap:
+ return wrap(new)
+ else:
+ return arr.copy()
+ newshape[axis] -= numtodel
+ new = empty(newshape, arr.dtype, arr.flags.fnc)
+ # copy initial chunk
+ if start == 0:
+ pass
+ else:
+ slobj[axis] = slice(None, start)
+ new[slobj] = arr[slobj]
+ # copy end chunck
+ if stop == N:
+ pass
+ else:
+ slobj[axis] = slice(stop-numtodel,None)
+ slobj2 = [slice(None)]*ndim
+ slobj2[axis] = slice(stop, None)
+ new[slobj] = arr[slobj2]
+ # copy middle pieces
+ if step == 1:
+ pass
+ else: # use array indexing.
+ obj = arange(start, stop, step, dtype=intp)
+ all = arange(start, stop, dtype=intp)
+ obj = setdiff1d(all, obj)
+ slobj[axis] = slice(start, stop-numtodel)
+ slobj2 = [slice(None)]*ndim
+ slobj2[axis] = obj
+ new[slobj] = arr[slobj2]
+ else: # default behavior
+ obj = array(obj, dtype=intp, copy=0, ndmin=1)
+ all = arange(N, dtype=intp)
+ obj = setdiff1d(all, obj)
+ slobj[axis] = obj
+ new = arr[slobj]
+ if wrap:
+ return wrap(new)
+ else:
+ return new
+
+def insert(arr, obj, values, axis=None):
+ """Return a new array with values inserted along the given axis
+ before the given indices
+
+ If axis is None, then ravel the array first.
+
+ The obj argument can be an integer, a slice, or a sequence of
+ integers.
+
+ Example:
+ >>> a = array([[1,2,3],
+ ... [4,5,6],
+ ... [7,8,9]])
+
+ >>> insert(a, [1,2], [[4],[5]], axis=0)
+ array([[1, 2, 3],
+ [4, 4, 4],
+ [4, 5, 6],
+ [5, 5, 5],
+ [7, 8, 9]])
+ """
+ wrap = None
+ if type(arr) is not ndarray:
+ try:
+ wrap = arr.__array_wrap__
+ except AttributeError:
+ pass
+
+ arr = asarray(arr)
+ ndim = arr.ndim
+ if axis is None:
+ if ndim != 1:
+ arr = arr.ravel()
+ ndim = arr.ndim
+ axis = ndim-1
+ if (ndim == 0):
+ arr = arr.copy()
+ arr[...] = values
+ if wrap:
+ return wrap(arr)
+ else:
+ return arr
+ slobj = [slice(None)]*ndim
+ N = arr.shape[axis]
+ newshape = list(arr.shape)
+ if isinstance(obj, (int, long, integer)):
+ if (obj < 0): obj += N
+ if obj < 0 or obj > N:
+ raise ValueError, "index (%d) out of range (0<=index<=%d) "\
+ "in dimension %d" % (obj, N, axis)
+ newshape[axis] += 1;
+ new = empty(newshape, arr.dtype, arr.flags.fnc)
+ slobj[axis] = slice(None, obj)
+ new[slobj] = arr[slobj]
+ slobj[axis] = obj
+ new[slobj] = values
+ slobj[axis] = slice(obj+1,None)
+ slobj2 = [slice(None)]*ndim
+ slobj2[axis] = slice(obj,None)
+ new[slobj] = arr[slobj2]
+ if wrap:
+ return wrap(new)
+ return new
+
+ elif isinstance(obj, slice):
+ # turn it into a range object
+ obj = arange(*obj.indices(N),**{'dtype':intp})
+
+ # get two sets of indices
+ # one is the indices which will hold the new stuff
+ # two is the indices where arr will be copied over
+
+ obj = asarray(obj, dtype=intp)
+ numnew = len(obj)
+ index1 = obj + arange(numnew)
+ index2 = setdiff1d(arange(numnew+N),index1)
+ newshape[axis] += numnew
+ new = empty(newshape, arr.dtype, arr.flags.fnc)
+ slobj2 = [slice(None)]*ndim
+ slobj[axis] = index1
+ slobj2[axis] = index2
+ new[slobj] = values
+ new[slobj2] = arr
+
+ if wrap:
+ return wrap(new)
+ return new
+
+def append(arr, values, axis=None):
+ """Append to the end of an array along axis (ravel first if None)
+ """
+ arr = asanyarray(arr)
+ if axis is None:
+ if arr.ndim != 1:
+ arr = arr.ravel()
+ values = ravel(values)
+ axis = arr.ndim-1
+ return concatenate((arr, values), axis=axis)
View Lorry Controller Status