The following extensions to Numerical Python are also defined:
- solve_tridiag_linear_system:
returns the solution of a tridiagonal linear system
- wrap2callable:
tool for turning constants, discrete data, string formulas, function objects, or plain functions into an object that behaves as a function
- NumPy_array_iterator:
allows iterating over all array elements using a single, standard for loop (for value, index in iterator), has some additional features compared with numpy.ndenumerate
- asarray_cpwarn:
as asarray(a), but a warning or exception is issued if the array a is copied
- meshgrid:
extended version of numpy.meshgrid to 1D, 2D and 3D grids, with sparse or dense coordinate arrays and matrix or grid indexing
- ndgrid:
same as calling meshgrid with indexing=’ij’ (matrix indexing)
- float_eq:
operator == for float operands with tolerance, float_eq(a,b,tol) means abs(a-b) < tol works for both scalar and array arguments (similar functions for other operations exists: float_le, float_lt, float_ge, float_gt, float_ne)
- cut_noise:
set all small (noise) elements of an array to zero
- matrix_rank:
compute the rank of a matrix
- orth:
compute an orthonormal basis from a matrix (taken from scipy.linalg to avoid scipy dependence)
- null:
compute the null space of a matrix
- norm_L2, norm_l2, norm_L1, norm_l1, norm_inf:
discrete and continuous norms for multi-dimensional arrays viewed as vectors
- compute_historgram:
return x and y arrays of a histogram, given a vector of samples
- factorial:
compute the factorial n! by various methods (iterative, recursive, reduce, functional, scipy, etc)
- seq
seq(a,b,s, [type]) computes numbers from a up to and including b in steps of s and (default) type float_ sequence = seq (for backward compatibility)
- iseq:
as seq, but integer counters are computed (iseq is an alternative to range where the upper limit is included in the sequence - this can be important for direct mapping of indices between mathematics and Python code) isequence = iseq (for backward compatibility)
- arr:
simplified/unified interface to creating various types of NumPy arrays (see its doc string)
Apply the Gram-Schmidt orthogonalization algorithm to a set of vectors. vecs is a two-dimensional array where the vectors are stored row-wise, or vecs may be a list of vectors, where each vector can be a list or a one-dimensional array.
The argument tol is a tolerance for null vectors (the absolute value of all elements must be less than tol to have a null vector).
If normalize is True, the orthogonal vectors are normalized to form an orthonormal basis.
If remove_null_vectors is True, all null vectors are removed from the resulting basis.
If remove_noise is True, all elements whose absolute values are less than tol are set to zero.
An array basis is returned, where basis[i,:] (row_wise_storage is True) or basis[:,i] (row_wise_storage is False) is the i-th orthogonal vector in the basis.
This function handles null vectors, see Gram_Schmidt1 for a (faster) function that does not.
Apply the Gram-Schmidt orthogonalization algorithm to a set of vectors. vecs is a two-dimensional array where the vectors are stored row-wise, or vecs may be a list of vectors, where each vector can be a list or a one-dimensional array. An array basis is returned, where basis[i,:] (row_wise_storage is True) or basis[:,i] (row_wise_storage is False) is the i-th orthonormal vector in the basis.
This function does not handle null vectors, see Gram_Schmidt for a (slower) function that does.
Bases: _Pmw.Pmw_1_3.lib.PmwBlt.Vector
Copy a NumPy array to a BLT vector: # a: some NumPy array b = NumPy2BltVector(a) # b is BLT vector g = Pmw.Blt.Graph(someframe) # send b to g for plotting
Iterate over all elements in a NumPy array a. Return values: generator function and the code of this function. The numpy.ndenumerate iterator performs the same iteration over an array, but NumPy_array_iterator has some additional features (especially handy for coding finite difference stencils, see next paragraph).
The keyword arguments specify offsets in the start and stop value of the index in each dimension. Legal values are offset0_start, offset0_stop, offset1_start, offset1_stop, etc. Also offset_start and offset_stop are legal keyword arguments, these imply the same offset value for all dimensions.
Another keyword argument is no_value, which can be True or False. If the value is True, the iterator returns the indices as a tuple, otherwise (default) the iterator returns a two-tuple consisting of the value of the array and the corresponding indices (as a tuple).
Examples:
>>> q = linspace(1, 2*3*4, 2*3*4); q.shape = (2,3,4)
>>> it, code = NumPy_array_iterator(q)
>>> print code # generator function with 3 nested loops:
>>> type(it)
<type 'function'>
>>> for value, index in it(q):
... print 'a%s = %g' % (index, value)
...
a(0, 0, 0) = 1
a(0, 0, 1) = 2
a(0, 0, 2) = 3
a(0, 0, 3) = 4
a(0, 1, 0) = 5
a(0, 1, 1) = 6
a(0, 1, 2) = 7
a(0, 1, 3) = 8
a(0, 2, 0) = 9
a(0, 2, 1) = 10
a(0, 2, 2) = 11
a(0, 2, 3) = 12
a(1, 0, 0) = 13
a(1, 0, 1) = 14
a(1, 0, 2) = 15
a(1, 0, 3) = 16
a(1, 1, 0) = 17
a(1, 1, 1) = 18
a(1, 1, 2) = 19
a(1, 1, 3) = 20
a(1, 2, 0) = 21
a(1, 2, 1) = 22
a(1, 2, 2) = 23
a(1, 2, 3) = 24
Here is the version where only the indices and no the values are returned by the iterator:
>>> q = linspace(1, 1*3, 3); q.shape = (1,3)
>>> it, code = NumPy_array_iterator(q, no_value=True)
>>> print code
>>> for i,j in it(q):
... print i,j
0 0
0 1
0 2
Now let us try some offsets:
>>> it, code = NumPy_array_iterator(q, offset1_stop=1, offset_start=1)
>>> print code
>>> # note: the offsets appear in the xrange arguments
>>> for value, index in it(q):
... print 'a%s = %g' % (index, value)
...
a(1, 1, 1) = 18
a(1, 1, 2) = 19
a(1, 1, 3) = 20
Turn discrete data on a uniform grid into a callable function, i.e., equip the data with an interpolation function.
>>> x = linspace(0, 1, 11)
>>> y = 1+2*x
>>> f = WrapDiscreteData2Callable((x,y))
>>> # or just use the wrap2callable generic function:
>>> f = wrap2callable((x,y))
>>> f(0.5) # evaluate f(x) by interpolation
1.5
>>> f(0.5, 0.1) # discrete data with extra time prm: f(x,t)
1.5
Compact and flexible interface for creating NumPy arrays, including several consistency and error checks.
The array can be created in four ways:
- as zeros (just shape specified),
- as uniformly spaced coordinates in an interval [a,b]
- as a copy of or reference to (depending on copy=True,False resp.) a list, tuple, or NumPy array (provided as the data argument),
- from data in a file (for one- or two-dimensional real-valued arrays).
The function calls the underlying NumPy functions zeros, array and linspace (see the NumPy manual for the functionality of these functions). In case of data in a file, the first line determines the number of columns in the array. The file format is just rows and columns with numbers, no decorations (square brackets, commas, etc.) are allowed.
>>> arr((3,4))
array([[ 0., 0., 0., 0.],
[ 0., 0., 0., 0.],
[ 0., 0., 0., 0.]])
>>> arr(4, element_type=int) + 4 # integer array
array([4, 4, 4, 4])
>>> arr(3, interval=[0,2])
array([ 0., 1., 2.])
>>> somelist=[[0,1],[5,5]]
>>> a = arr(data=somelist)
>>> a # a has always float elements by default
array([[ 0., 1.],
[ 5., 5.]])
>>> a = arr(data=somelist, element_type=int)
>>> a
array([[0, 1],
[5, 5]])
>>> b = a + 1
>>> c = arr(data=b, copy=False) # let c share data with b
>>> b is c
True
>>> id(b) == id(c)
True
>>> # make a file with array data:
>>> f = open('tmp.dat', 'w')
>>> f.write(''' ... 1 3
... 2 6
... 3 12
... 3.5 20
... ''')
>>> f.close()
>>> # read array data from file:
>>> a = arr(file_='tmp.dat')
>>> a
array([[ 1. , 3. ],
[ 2. , 6. ],
[ 3. , 12. ],
[ 3.5, 20. ]])
Compute the factorial n! using long integers. Different implementations are available (see source code for the methods).
Here is an efficiency comparison of the methods (computing 80!): reduce | 1.00 lambda list comprehension | 1.70 lambda functional | 3.08 plain recursive | 5.83 lambda recursive | 21.73 scipy | 131.18
Extension of numpy.meshgrid to 1D, 2D and 3D problems, and also support of both “matrix” and “grid” numbering.
This extended version makes 1D/2D/3D coordinate arrays for vectorized evaluations of 1D/2D/3D scalar/vector fields over 1D/2D/3D grids, given one-dimensional coordinate arrays x, y, and/or, z.
>>> x=linspace(0,1,3) # coordinates along x axis
>>> y=linspace(0,1,2) # coordinates along y axis
>>> xv, yv = meshgrid(x,y) # extend x and y for a 2D xy grid
>>> xv
array([[ 0. , 0.5, 1. ],
[ 0. , 0.5, 1. ]])
>>> yv
array([[ 0., 0., 0.],
[ 1., 1., 1.]])
>>> xv, yv = meshgrid(x,y, sparse=True) # make sparse output arrays
>>> xv
array([[ 0. , 0.5, 1. ]])
>>> yv
array([[ 0.],
[ 1.]])
>>> # 2D slice of a 3D grid, with z=const:
>>> z=5
>>> xv, yv, zc = meshgrid(x,y,z)
>>> xv
array([[ 0. , 0.5, 1. ],
[ 0. , 0.5, 1. ]])
>>> yv
array([[ 0., 0., 0.],
[ 1., 1., 1.]])
>>> zc
5
>>> # 2D slice of a 3D grid, with x=const:
>>> meshgrid(2,y,x)
(2, array([[ 0., 1.],
[ 0., 1.],
[ 0., 1.]]), array([[ 0. , 0. ],
[ 0.5, 0.5],
[ 1. , 1. ]]))
>>> meshgrid(0,1,5, sparse=True) # just a 3D point
(0, 1, 5)
>>> meshgrid(y) # 1D grid; y is just returned
array([ 0., 1.])
>>> meshgrid(x,y, indexing='ij') # change to matrix indexing
(array([[ 0. , 0. ],
[ 0.5, 0.5],
[ 1. , 1. ]]), array([[ 0., 1.],
[ 0., 1.],
[ 0., 1.]]))
Why does SciTools has its own meshgrid function when NumPy has three similar functions, mgrid, ogrid, and meshgrid? The meshgrid function in NumPy is limited to two dimensions only, while the SciTools version can also work with 3D and 1D grids. In addition, the NumPy version of meshgrid has no option for generating sparse grids to conserve memory, like we have in SciTools by specifying the sparse argument:
The NumPy functions mgrid and ogrid does provide support for, respectively, full and sparse n-dimensional meshgrids, however, these functions uses slices to generate the meshgrids rather than one-dimensional coordinate arrays such as in Matlab. With slices, the user does not have the option to generate meshgrid with, e.g., irregular spacings, like:: >>> x = array([-1,-0.5,1,4,5], float) >>> y = array([0,-2,-5], float) >>> xv, yv = meshgrid(x, y, sparse=False) >>> xv array([[-1. , -0.5, 1. , 4. , 5. ],
[-1. , -0.5, 1. , 4. , 5. ], [-1. , -0.5, 1. , 4. , 5. ]])
>>> yv
array([[ 0., 0., 0., 0., 0.],
[-2., -2., -2., -2., -2.],
[-5., -5., -5., -5., -5.]])
>>>
In addition to the reasons mentioned above, the meshgrid function in NumPy supports only Cartesian indexing, i.e., x and y, not matrix indexing, i.e., rows and columns (mgrid and ogrid supports only matrix indexing). The meshgrid function in SciTools supports both indexing conventions through the indexing keyword argument. Giving the string ‘ij’ returns a meshgrid with matrix indexing, while ‘xy’ returns a meshgrid with Cartesian indexing. The difference is illustrated by the following code snippet:
nx = 10 ny = 15
x = linspace(-2,2,nx) y = linspace(-2,2,ny)
xv, yv = meshgrid(x, y, sparse=False, indexing=’ij’) for i in range(nx):
- for j in range(ny):
- # treat xv[i,j], yv[i,j]
xv, yv = meshgrid(x, y, sparse=False, indexing=’xy’) for i in range(nx):
- for j in range(ny):
- # treat xv[j,i], yv[j,i]
It is not entirely true that matrix indexing is not supported by the meshgrid function in NumPy because we can just switch the order of the first two input and output arguments:: >>> yv, xv = numpy.meshgrid(y, x) >>> # same as: >>> xv, yv = meshgrid(x, y, indexing=’ij’) However, we think it is clearer to have the logical “x, y” sequence on the left-hand side and instead adjust a keyword argument.
L1 norm of a multi-dimensional array u viewed as a vector (norm=sum(abs(u.ravel()))).
If u holds function values and the norm of u is supposed to approximate an integral (L1 norm) of the function, this (and not norm_l1) is the right norm function to use.
L2 norm of a multi-dimensional array u viewed as a vector (norm is norm_l2/n, where n is length of u (no of elements)).
If u holds function values and the norm of u is supposed to approximate an integral (L2 norm) of the function, this (and not norm_l2) is the right norm function to use.
Return the null space of a matrix A. If row_wise_storage is True, a two-dimensional array where the vectors that span the null space are stored as rows, otherwise they are stored as columns.
Code by Bastian Weber based on code by Robert Kern and Ryan Krauss.
(Plain copy from scipy.linalg.orth - this one here applies numpy.svd and avoids the need for having scipy installed.)
Construct an orthonormal basis for the range of A using SVD.
@param A: array, shape (M, N) @return:
Q : array, shape (M, K) Orthonormal basis for the range of A. K = effective rank of A, as determined by automatic cutoff
see also svd (singular value decomposition of a matrix in scipy.linalg)
Solve an n times n tridiagonal linear system of the form:
A[0,1]*x[0] + A[0,2]*x[1] = 0
A[1,0]*x[0] + A[1,1]*x[1] + A[1,2]*x[2] = 0
...
...
A[k,0]*x[k-1] + A[k,1]*x[k] + A[k,2]*x[k+1] = 0
...
A[n-2,0]*x[n-3] + A[n-2,1]*x[n-2] + A[n-2,2]*x[n-1] = 0
...
A[n-1,0]*x[n-2] + A[n-1,1]*x[n-1] = 0
The diagonal of the coefficent matrix is stored in A[:,1], the subdiagonal is stored in A[1:,0], and the superdiagonal is stored in A[:-1,2].
Allow constants, string formulas, discrete data points, user-defined functions and (callable) classes to be wrapped in a new callable function. That is, all the mentioned data structures can be used as a function, usually of space and/or time. (kwargs is used for string formulas)
>>> f1 = wrap2callable(2.0)
>>> f1(0.5)
2.0
>>> f2 = wrap2callable('1+2*x')
>>> f2(0.5)
2.0
>>> f3 = wrap2callable('1+2*t', independent_variable='t')
>>> f3(0.5)
2.0
>>> f4 = wrap2callable('a+b*t')
>>> f4(0.5)
Traceback (most recent call last):
...
NameError: name 'a' is not defined
>>> f4 = wrap2callable('a+b*t', independent_variable='t', a=1, b=2)
>>> f4(0.5)
2.0
>>> x = linspace(0, 1, 3); y=1+2*x
>>> f5 = wrap2callable((x,y))
>>> f5(0.5)
2.0
>>> def myfunc(x): return 1+2*x
>>> f6 = wrap2callable(myfunc)
>>> f6(0.5)
2.0
>>> f7 = wrap2callable(lambda x: 1+2*x)
>>> f7(0.5)
2.0
>>> class MyClass:
'Representation of a function f(x; a, b) =a + b*x'
def __init__(self, a=1, b=1):
self.a = a; self.b = b
def __call__(self, x):
return self.a + self.b*x
>>> myclass = MyClass(a=1, b=2)
>>> f8 = wrap2callable(myclass)
>>> f8(0.5)
2.0
>>> # 3D functions:
>>> f9 = wrap2callable('1+2*x+3*y+4*z', independent_variables=('x','y','z'))
>>> f9(0.5,1/3.,0.25)
4.0
>>> # discrete 3D data:
>>> y = linspace(0, 1, 3); z = linspace(-1, 0.5, 16)
>>> xv = reshape(x, (len(x),1,1))
>>> yv = reshape(y, (1,len(y),1))
>>> zv = reshape(z, (1,1,len(z)))
>>> def myfunc3(x,y,z): return 1+2*x+3*y+4*z
>>> values = myfunc3(xv, yv, zv)
>>> f10 = wrap2callable((x, y, z, values))
>>> f10(0.5, 1/3., 0.25)
4.0
One can also check what the object is wrapped as and do more specific operations, e.g.,
>>> f9.__class__.__name__
'StringFunction'
>>> str(f9) # look at function formula
'1+2*x+3*y+4*z'
>>> f8.__class__.__name__
'MyClass'
>>> f8.a, f8.b # access MyClass-specific data
(1, 2)
Troubleshooting regarding string functions: If you use a string formula with a NumPy array, you typically get error messages like:
TypeError: only rank-0 arrays can be converted to Python scalars.
You must then make the right import (numpy is recommended):
from Numeric/numarray/numpy/scitools.numpytools import *
in the calling code and supply the keyword argument:
globals=globals()
to wrap2callable. See also the documentation of class StringFunction for more information.