[Scipy-svn] r4154 - in trunk/scipy: cluster cluster/tests fftpack integrate integrate/tests io io/arff io/arff/tests io/matlab io/matlab/tests io/tests lib lib/blas lib/lapack lib/lapack/tests linalg linalg/benchmarks linalg/tests ndimage ndimage/tests optimize sandbox/numexpr sparse sparse/linalg sparse/linalg/dsolve sparse/linalg/dsolve/tests sparse/linalg/dsolve/umfpack/tests sparse/linalg/eigen sparse/linalg/eigen/arpack sparse/linalg/eigen/arpack/tests sparse/linalg/eigen/lobpcg sparse/linalg/eigen/lobpcg/tests sparse/linalg/isolve sparse/linalg/isolve/tests sparse/linalg/tests sparse/sparsetools sparse/tests splinalg stats stats/models stats/models/tests stats/tests testing testing/examples weave weave/tests
scipy-svn@scip...
scipy-svn@scip...
Sun Apr 20 07:16:11 CDT 2008
Author: jarrod.millman
Date: 2008-04-20 07:15:19 -0500 (Sun, 20 Apr 2008)
New Revision: 4154
Modified:
trunk/scipy/cluster/hierarchy.py
trunk/scipy/cluster/tests/test_hierarchy.py
trunk/scipy/fftpack/setupscons.py
trunk/scipy/integrate/ode.py
trunk/scipy/integrate/tests/test_integrate.py
trunk/scipy/io/__init__.py
trunk/scipy/io/arff/arffread.py
trunk/scipy/io/arff/tests/test_data.py
trunk/scipy/io/data_store.py
trunk/scipy/io/matlab/mio4.py
trunk/scipy/io/matlab/tests/test_mio.py
trunk/scipy/io/mmio.py
trunk/scipy/io/npfile.py
trunk/scipy/io/pickler.py
trunk/scipy/io/tests/test_mmio.py
trunk/scipy/io/tests/test_recaster.py
trunk/scipy/lib/blas/scons_support.py
trunk/scipy/lib/lapack/scons_support.py
trunk/scipy/lib/lapack/tests/test_lapack.py
trunk/scipy/lib/setupscons.py
trunk/scipy/linalg/basic.py
trunk/scipy/linalg/benchmarks/bench_decom.py
trunk/scipy/linalg/blas.py
trunk/scipy/linalg/decomp.py
trunk/scipy/linalg/iterative.py
trunk/scipy/linalg/matfuncs.py
trunk/scipy/linalg/scons_support.py
trunk/scipy/linalg/tests/test_decomp.py
trunk/scipy/ndimage/_registration.py
trunk/scipy/ndimage/_segmenter.py
trunk/scipy/ndimage/tests/test_segment.py
trunk/scipy/optimize/slsqp.py
trunk/scipy/sandbox/numexpr/__init__.py
trunk/scipy/sparse/base.py
trunk/scipy/sparse/bsr.py
trunk/scipy/sparse/compressed.py
trunk/scipy/sparse/construct.py
trunk/scipy/sparse/coo.py
trunk/scipy/sparse/csc.py
trunk/scipy/sparse/csr.py
trunk/scipy/sparse/data.py
trunk/scipy/sparse/dia.py
trunk/scipy/sparse/dok.py
trunk/scipy/sparse/info.py
trunk/scipy/sparse/lil.py
trunk/scipy/sparse/linalg/dsolve/linsolve.py
trunk/scipy/sparse/linalg/dsolve/setupscons.py
trunk/scipy/sparse/linalg/dsolve/tests/test_linsolve.py
trunk/scipy/sparse/linalg/dsolve/umfpack/tests/test_umfpack.py
trunk/scipy/sparse/linalg/eigen/arpack/arpack.py
trunk/scipy/sparse/linalg/eigen/arpack/tests/test_arpack.py
trunk/scipy/sparse/linalg/eigen/lobpcg/info.py
trunk/scipy/sparse/linalg/eigen/lobpcg/lobpcg.py
trunk/scipy/sparse/linalg/eigen/lobpcg/tests/test_lobpcg.py
trunk/scipy/sparse/linalg/eigen/setup.py
trunk/scipy/sparse/linalg/eigen/setupscons.py
trunk/scipy/sparse/linalg/isolve/__init__.py
trunk/scipy/sparse/linalg/isolve/iterative.py
trunk/scipy/sparse/linalg/isolve/minres.py
trunk/scipy/sparse/linalg/isolve/tests/test_iterative.py
trunk/scipy/sparse/linalg/isolve/utils.py
trunk/scipy/sparse/linalg/setup.py
trunk/scipy/sparse/linalg/setupscons.py
trunk/scipy/sparse/linalg/tests/test_interface.py
trunk/scipy/sparse/setupscons.py
trunk/scipy/sparse/sparsetools/__init__.py
trunk/scipy/sparse/sparsetools/bsr.py
trunk/scipy/sparse/sparsetools/coo.py
trunk/scipy/sparse/sparsetools/csc.py
trunk/scipy/sparse/sparsetools/csr.py
trunk/scipy/sparse/sparsetools/dia.py
trunk/scipy/sparse/sparsetools/setupscons.py
trunk/scipy/sparse/spfuncs.py
trunk/scipy/sparse/sputils.py
trunk/scipy/sparse/tests/bench_sparse.py
trunk/scipy/sparse/tests/test_base.py
trunk/scipy/sparse/tests/test_construct.py
trunk/scipy/sparse/tests/test_spfuncs.py
trunk/scipy/sparse/tests/test_sputils.py
trunk/scipy/splinalg/__init__.py
trunk/scipy/stats/mmorestats.py
trunk/scipy/stats/models/contrast.py
trunk/scipy/stats/models/formula.py
trunk/scipy/stats/models/tests/test_formula.py
trunk/scipy/stats/mstats.py
trunk/scipy/stats/stats.py
trunk/scipy/stats/tests/test_mmorestats.py
trunk/scipy/stats/tests/test_mstats.py
trunk/scipy/testing/decorators.py
trunk/scipy/testing/examples/test_foo.py
trunk/scipy/testing/nosetester.py
trunk/scipy/testing/nulltester.py
trunk/scipy/testing/pkgtester.py
trunk/scipy/testing/utils.py
trunk/scipy/weave/size_check.py
trunk/scipy/weave/tests/test_blitz_tools.py
trunk/scipy/weave/tests/test_c_spec.py
Log:
ran reindent
Modified: trunk/scipy/cluster/hierarchy.py
===================================================================
--- trunk/scipy/cluster/hierarchy.py 2008-04-20 08:27:24 UTC (rev 4153)
+++ trunk/scipy/cluster/hierarchy.py 2008-04-20 12:15:19 UTC (rev 4154)
@@ -106,7 +106,7 @@
The Wolfram Research, Inc. http://reference.wolfram.com/...
...mathematica/HierarchicalClustering/tutorial/...
HierarchicalClustering.html. Accessed October 1, 2007.
-
+
[3] Gower, JC and Ross, GJS. "Minimum Spanning Trees and Single Linkage
Cluster Analysis." Applied Statistics. 18(1): pp. 54--64. 1969.
@@ -209,7 +209,7 @@
return numpy.float64(a)
else:
return a
-
+
def _copy_arrays_if_base_present(T):
"""
Accepts a tuple of arrays T. Copies the array T[i] if its base array
@@ -219,8 +219,8 @@
"""
l = [_copy_array_if_base_present(a) for a in T]
return l
-
+
def copying():
""" Displays the license for this package."""
print _copyingtxt
@@ -286,7 +286,7 @@
def centroid(y):
"""
Z = centroid(y)
-
+
Performs centroid/UPGMC linkage on the condensed distance matrix Z.
See linkage for more information on the return structure and
algorithm.
@@ -297,7 +297,7 @@
Performs centroid/UPGMC linkage on the observation matrix X using
Euclidean distance as the distance metric. See linkage for more
- information on the return structure and algorithm.
+ information on the return structure and algorithm.
"""
return linkage(y, method='centroid', metric='euclidean')
@@ -314,7 +314,7 @@
Performs median/WPGMC linkage on the observation matrix X using
Euclidean distance as the distance metric. See linkage for more
- information on the return structure and algorithm.
+ information on the return structure and algorithm.
(a condensed alias for linkage)
"""
@@ -323,7 +323,7 @@
def ward(y):
"""
Z = ward(y)
-
+
Performs Ward's linkage on the condensed distance matrix Z. See
linkage for more information on the return structure and algorithm.
@@ -331,13 +331,13 @@
Performs Ward's linkage on the observation matrix X using Euclidean
distance as the distance metric. See linkage for more information
- on the return structure and algorithm.
+ on the return structure and algorithm.
(a condensed alias for linkage)
"""
return linkage(y, method='ward', metric='euclidean')
-
+
def linkage(y, method='single', metric='euclidean'):
""" Z = linkage(y, method)
@@ -367,11 +367,11 @@
A distance matrix is maintained at each iteration. The
d[i,j] entry corresponds to the distance between cluster
i and j in the original forest.
-
+
At each iteration, the algorithm must update the distance
matrix to reflect the distance of the newly formed cluster
u with the remaining clusters in the forest.
-
+
Suppose there are |u| original observations u[0], ..., u[|u|-1]
in cluster u and |v| original objects v[0], ..., v[|v|-1]
in cluster v. Recall s and t are combined to form cluster
@@ -380,7 +380,7 @@
The following are methods for calculating the distance between
the newly formed cluster u and each v.
-
+
* method='single' assigns dist(u,v) = MIN(dist(u[i],v[j])
for all points i in cluster u and j in cluster v.
@@ -402,7 +402,7 @@
* method='weighted' assigns
dist(u,v) = (dist(s,v) + dist(t,v))/2
-
+
where cluster u was formed with cluster s and t and v
is a remaining cluster in the forest. (also called WPGMA)
@@ -422,11 +422,11 @@
the Euclidean distance between the centroid of u and the
centroid of a remaining cluster v in the forest.
(also called UPGMC)
-
+
* method='median' assigns dist(s,t) as above. When two clusters
s and t are combined into a new cluster u, the average of
centroids s and t give the new centroid u. (also called WPGMC)
-
+
* method='ward' uses the Ward variance minimization algorithm.
The new entry dist(u, v) is computed as follows,
@@ -452,7 +452,7 @@
if type(y) != _array_type:
raise TypeError("Argument 'y' must be a numpy array.")
-
+
s = y.shape
if len(s) == 1:
is_valid_y(y, throw=True, name='y')
@@ -517,7 +517,7 @@
def getId(self):
"""
i = nd.getId()
-
+
Returns the id number of the node nd. For 0 <= i < n, i
corresponds to original observation i. For n <= i < 2n - 1,
i corresponds to non-singleton cluster formed at iteration i-n.
@@ -560,24 +560,24 @@
def preOrder(self, func=(lambda x: x.id)):
"""
vlst = preOrder(func)
-
+
Performs preorder traversal without recursive function calls.
When a leaf node is first encountered, func is called with the
leaf node as its argument, and its result is appended to the
list vlst.
-
+
For example, the statement
-
+
ids = root.preOrder(lambda x: x.id)
-
+
returns a list of the node ids corresponding to the leaf
nodes of the tree as they appear from left to right.
"""
-
+
# Do a preorder traversal, caching the result. To avoid having to do
# recursion, we'll store the previous index we've visited in a vector.
n = self.count
-
+
curNode = [None] * (2 * n)
lvisited = numpy.zeros((2 * n,), dtype='bool')
rvisited = numpy.zeros((2 * n,), dtype='bool')
@@ -603,7 +603,7 @@
# node already, go up in the tree.
else:
k = k - 1
-
+
return preorder
_cnode_bare = cnode(0)
@@ -612,11 +612,11 @@
def totree(Z, rd=False):
"""
r = totree(Z)
-
+
Converts a hierarchical clustering encoded in the matrix Z
(by linkage) into an easy-to-use tree object. The reference r
to the root cnode object is returned.
-
+
Each cnode object has a left, right, dist, id, and count
attribute. The left and right attributes point to cnode
objects that were combined to generate the cluster. If
@@ -638,7 +638,7 @@
"""
is_valid_linkage(Z, throw=True, name='Z')
-
+
# The number of original objects is equal to the number of rows minus
# 1.
n = Z.shape[0] + 1
@@ -690,7 +690,7 @@
... = squareform(...)
Converts a vector-form distance vector to a square-form distance
- matrix, and vice-versa.
+ matrix, and vice-versa.
v = squareform(X)
@@ -719,7 +719,7 @@
ignored any way so they do not disrupt the squareform
transformation.
"""
-
+
if type(X) is not _array_type:
raise TypeError('The parameter passed must be an array.')
@@ -739,7 +739,7 @@
# Check that v is of valid dimensions.
if d * (d - 1) / 2 != int(s[0]):
raise ValueError('Incompatible vector size. It must be a binomial coefficient n choose 2 for some integer n >= 2.')
-
+
# Allocate memory for the distance matrix.
M = numpy.zeros((d, d), 'double')
@@ -766,7 +766,7 @@
# One-side of the dimensions is set here.
d = s[0]
-
+
# Create a vector.
v = numpy.zeros(((d * (d - 1) / 2),), 'double')
@@ -785,7 +785,7 @@
def minkowski(u, v, p):
"""
d = minkowski(u, v, p)
-
+
Returns the Minkowski distance between two vectors u and v,
||u-v||_p = (\sum {|u_i - v_i|^p})^(1/p).
@@ -797,7 +797,7 @@
def euclidean(u, v):
"""
d = euclidean(u, v)
-
+
Computes the Euclidean distance between two n-vectors u and v, ||u-v||_2
"""
q=numpy.matrix(u-v)
@@ -825,7 +825,7 @@
def correlation(u, v):
"""
d = correlation(u, v)
-
+
Computes the correlation distance between two n-vectors u and v,
1 - (u - n|u|_1)(v - n|v|_1)^T
@@ -846,7 +846,7 @@
def hamming(u, v):
"""
d = hamming(u, v)
-
+
Computes the Hamming distance between two n-vectors u and v,
which is simply the proportion of disagreeing components in u
and v. If u and v are boolean vectors, the hamming distance is
@@ -904,7 +904,7 @@
def seuclidean(u, v, V):
"""
d = seuclidean(u, v, V)
-
+
Returns the standardized Euclidean distance between two
n-vectors u and v. V is a m-dimensional vector of component
variances. It is usually computed among a larger collection vectors.
@@ -925,7 +925,7 @@
def mahalanobis(u, v, VI):
"""
d = mahalanobis(u, v, VI)
-
+
Computes the Mahalanobis distance between two n-vectors u and v,
(u-v)VI(u-v)^T
where VI is the inverse covariance matrix.
@@ -937,7 +937,7 @@
def chebyshev(u, v):
"""
d = chebyshev(u, v)
-
+
Computes the Chebyshev distance between two n-vectors u and v,
\max {|u_i-v_i|}.
"""
@@ -946,7 +946,7 @@
def braycurtis(u, v):
"""
d = braycurtis(u, v)
-
+
Computes the Bray-Curtis distance between two n-vectors u and v,
\sum{|u_i-v_i|} / \sum{|u_i+v_i|}.
"""
@@ -1018,7 +1018,7 @@
def dice(u, v):
"""
d = dice(u, v)
-
+
Computes the Dice dissimilarity between two boolean n-vectors
u and v, which is
@@ -1039,7 +1039,7 @@
def rogerstanimoto(u, v):
"""
d = rogerstanimoto(u, v)
-
+
Computes the Rogers-Tanimoto dissimilarity between two boolean
n-vectors u and v,
@@ -1062,7 +1062,7 @@
def russellrao(u, v):
"""
d = russellrao(u, v)
-
+
Computes the Russell-Rao dissimilarity between two boolean n-vectors
u and v, (n - c_{TT}) / n where c_{ij} is the number of occurrences
of u[k] == i and v[k] == j for k < n.
@@ -1121,7 +1121,7 @@
def pdist(X, metric='euclidean', p=2, V=None, VI=None):
""" Y = pdist(X, method='euclidean', p=2)
-
+
Computes the distance between m original observations in
n-dimensional space. Returns a condensed distance matrix Y.
For each i and j (i<j), the metric dist(u=X[i], v=X[j]) is
@@ -1163,7 +1163,7 @@
6. Y = pdist(X, 'cosine')
Computes the cosine distance between vectors u and v,
-
+
1 - uv^T
-----------
|u|_2 |v|_2
@@ -1270,7 +1270,7 @@
boolean vectors. (see sokalsneath function documentation)
21. Y = pdist(X, f)
-
+
Computes the distance between all pairs of vectors in X
using the user supplied 2-arity function f. For example,
Euclidean distance between the vectors could be computed
@@ -1297,7 +1297,7 @@
# using the distance metric Y but with a more succint,
# verifiable, but less efficient implementation.
-
+
if type(X) is not _array_type:
raise TypeError('The parameter passed must be an array.')
@@ -1367,7 +1367,7 @@
else:
raise TypeError('Invalid input array value type %s for jaccard.' % str(X.dtype))
elif mstr in set(['chebychev', 'chebyshev', 'cheby', 'cheb', 'ch']):
- _hierarchy_wrap.pdist_chebyshev_wrap(X, dm)
+ _hierarchy_wrap.pdist_chebyshev_wrap(X, dm)
elif mstr in set(['minkowski', 'mi', 'm']):
_hierarchy_wrap.pdist_minkowski_wrap(X, dm, p)
elif mstr in set(['seuclidean', 'se', 's']):
@@ -1511,7 +1511,7 @@
(c, d) = cophenet(Z, Y, [])
Also returns the cophenetic distance matrix in condensed form.
-
+
"""
nargs = len(args)
@@ -1535,7 +1535,7 @@
Y = args[1]
Ys = Y.shape
is_valid_y(Y, throw=True, name='Y')
-
+
z = zz.mean()
y = Y.mean()
Yy = Y - y
@@ -1555,7 +1555,7 @@
def inconsistent(Z, d=2):
"""
R = inconsistent(Z, d=2)
-
+
Calculates statistics on links up to d levels below each
non-singleton cluster defined in the (n-1)x4 linkage matrix Z.
@@ -1587,17 +1587,17 @@
_hierarchy_wrap.inconsistent_wrap(Z, R, int(n), int(d));
return R
-
+
def from_mlab_linkage(Z):
"""
Z2 = from_mlab_linkage(Z)
-
+
Converts a linkage matrix Z generated by MATLAB(TM) to a new linkage
matrix Z2 compatible with this module. The conversion does two
things:
* the indices are converted from 1..N to 0..(N-1) form, and
-
+
* a fourth column Z[:,3] is added where Z[i,3] is equal to
the number of original observations (leaves) in the non-singleton
cluster i.
@@ -1625,13 +1625,13 @@
1..N indexing.
"""
is_valid_linkage(Z, throw=True, name='Z')
-
+
return numpy.hstack([Z[:,0:2] + 1, Z[:,2]])
def is_monotonic(Z):
"""
is_monotonic(Z)
-
+
Returns True if the linkage Z is monotonic. The linkage is monotonic
if for every cluster s and t joined, the distance between them is
no less than the distance between any previously joined clusters.
@@ -1644,7 +1644,7 @@
def is_valid_im(R, warning=False, throw=False, name=None):
"""
is_valid_im(R)
-
+
Returns True if the inconsistency matrix passed is valid. It must
be a n by 4 numpy array of doubles. The standard deviations R[:,1]
must be nonnegative. The link counts R[:,2] must be positive and
@@ -1805,7 +1805,7 @@
def is_valid_dm(D, t=0.0):
"""
is_valid_dm(D)
-
+
Returns True if the variable D passed is a valid distance matrix.
Distance matrices must be 2-dimensional numpy arrays containing
doubles. They must have a zero-diagonal, and they must be symmetric.
@@ -1889,7 +1889,7 @@
def numobs_dm(D):
"""
numobs_dm(D)
-
+
Returns the number of original observations that correspond to a
square, non-condensed distance matrix D.
"""
@@ -1899,7 +1899,7 @@
def numobs_y(Y):
"""
numobs_y(Y)
-
+
Returns the number of original observations that correspond to a
condensed distance matrix Y.
"""
@@ -1910,7 +1910,7 @@
def Z_y_correspond(Z, Y):
"""
yesno = Z_y_correspond(Z, Y)
-
+
Returns True if a linkage matrix Z and condensed distance matrix
Y could possibly correspond to one another. They must have the same
number of original observations. This function is useful as a sanity
@@ -1931,7 +1931,7 @@
original observation i belongs.
The criterion parameter can be any of the following values,
-
+
* 'inconsistent': If a cluster node and all its decendents have an
inconsistent value less than or equal to c then all its leaf
descendents belong to the same flat cluster. When no non-singleton
@@ -1968,19 +1968,19 @@
cluster node c when monocrit[i] <= r for all cluster indices i below
and including c. r is minimized such that no more than t flat clusters
are formed. monocrit must be monotonic.
-
+
For example, to minimize the threshold t on maximum inconsistency
values so that no more than 3 flat clusters are formed, do:
MI = maxinconsts(Z, R)
cluster(Z, t=3, criterion='maxclust_monocrit', monocrit=MI)
-
+
"""
is_valid_linkage(Z, throw=True, name='Z')
n = Z.shape[0] + 1
T = numpy.zeros((n,), dtype='int32')
-
+
# Since the C code does not support striding using strides.
# The dimensions are used instead.
[Z] = _copy_arrays_if_base_present([Z])
@@ -2033,12 +2033,12 @@
specified.
Named parameters are described below.
-
+
criterion: specifies the criterion for forming flat clusters.
Valid values are 'inconsistent', 'distance', or
'maxclust' cluster formation algorithms. See
cluster for descriptions.
-
+
method: the linkage method to use. See linkage for
descriptions.
@@ -2046,7 +2046,7 @@
distances. See pdist for descriptions and
linkage to verify compatibility with the linkage
method.
-
+
t: the cut-off threshold for the cluster function or
the maximum number of clusters (criterion='maxclust').
@@ -2083,15 +2083,15 @@
_hierarchy_wrap.prelist_wrap(Z, ML, int(n))
return ML
-# Let's do a conditional import. If matplotlib is not available,
+# Let's do a conditional import. If matplotlib is not available,
try:
-
+
import matplotlib
import matplotlib.pylab
import matplotlib.patches
#import matplotlib.collections
_mpl = True
-
+
# Maps number of leaves to text size.
#
# p <= 20, size="12"
@@ -2182,7 +2182,7 @@
if leaf_font_size:
matplotlib.pylab.setp(lbls, 'size', leaf_font_size)
else:
- matplotlib.pylab.setp(lbls, 'size', float(_get_tick_text_size(p)))
+ matplotlib.pylab.setp(lbls, 'size', float(_get_tick_text_size(p)))
axis.xaxis.set_ticks_position('top')
# Make the tick marks invisible because they cover up the links
for line in axis.get_xticklines():
@@ -2272,7 +2272,7 @@
e.set_clip_box(axis.bbox)
e.set_alpha(0.5)
e.set_facecolor('k')
-
+
#matplotlib.pylab.plot(xs, ys, 'go', markeredgecolor='k', markersize=3)
#matplotlib.pylab.plot(ys, xs, 'go', markeredgecolor='k', markersize=3)
@@ -2377,7 +2377,7 @@
* 'top': plots the root at the top, and plot descendent
links going downwards. (default).
-
+
* 'bottom': plots the root at the bottom, and plot descendent
links going upwards.
@@ -2396,7 +2396,7 @@
When labels=None, the index of the original observation is used
used.
-
+
R = dendrogram(..., count_sort)
When plotting a cluster node and its directly descendent links,
@@ -2405,7 +2405,7 @@
of count_sort are:
* False: nothing is done.
-
+
* 'ascending'/True: the child with the minimum number of
original objects in its cluster is plotted first.
@@ -2464,7 +2464,7 @@
When a callable function is passed, leaf_label_func is passed
cluster index k, and returns a string with the label for the
leaf.
-
+
Indices k < n correspond to original observations while indices
k >= n correspond to non-singleton clusters.
@@ -2618,7 +2618,7 @@
if show_leaf_counts:
ivl.append("(" + str(int(Z[i-n, 3])) + ")")
else:
- ivl.append("")
+ ivl.append("")
def _append_contraction_marks(Z, iv, i, n, contraction_marks):
_append_contraction_marks_sub(Z, iv, Z[i-n, 0], n, contraction_marks)
@@ -2629,8 +2629,8 @@
contraction_marks.append((iv, Z[i-n, 2]))
_append_contraction_marks_sub(Z, iv, Z[i-n, 0], n, contraction_marks)
_append_contraction_marks_sub(Z, iv, Z[i-n, 1], n, contraction_marks)
-
+
def _dendrogram_calculate_info(Z, p, truncate_mode, \
colorthreshold=scipy.inf, get_leaves=True, \
orientation='top', labels=None, \
@@ -2649,7 +2649,7 @@
variable value to plot the left-most leaf node below the root node i
(if orientation='top', this would be the left-most x value where the
plotting of this root node i and its descendents should begin).
-
+
ivl is a list to store the labels of the leaf nodes. The leaf_label_func
is called whenever ivl != None, labels == None, and
leaf_label_func != None. When ivl != None and labels != None, the
@@ -2668,7 +2668,7 @@
* left is the independent variable coordinate of the center of the
the U of the subtree
-
+
* w is the amount of space used for the subtree (in independent
variable units)
@@ -2676,7 +2676,7 @@
* md is the max(Z[*,2]) for all nodes * below and including
the target node.
-
+
"""
if n == 0:
raise ValueError("Invalid singleton cluster count n.")
@@ -2712,7 +2712,7 @@
elif truncate_mode in ('mlab',):
pass
-
+
# Otherwise, only truncate if we have a leaf node.
#
# If the truncate_mode is mlab, the linkage has been modified
@@ -2900,7 +2900,7 @@
else:
d[T1[i]] = T2[i]
return True
-
+
def maxdists(Z):
"""
MD = maxdists(Z)
@@ -2910,12 +2910,12 @@
and including the node with index i. More specifically,
MD[i] = Z[Q(i)-n, 2].max() where Q(i) is the set of all node indices
below and including node i.
-
+
Note that when Z[:,2] is monotonic, Z[:,2] and MD should not differ.
See linkage for more information on this issue.
"""
is_valid_linkage(Z, throw=True, name='Z')
-
+
n = Z.shape[0] + 1
MD = numpy.zeros((n-1,))
[Z] = _copy_arrays_if_base_present([Z])
@@ -2933,7 +2933,7 @@
"""
is_valid_linkage(Z, throw=True, name='Z')
is_valid_im(R, throw=True, name='R')
-
+
n = Z.shape[0] + 1
MI = numpy.zeros((n-1,))
[Z, R] = _copy_arrays_if_base_present([Z, R])
@@ -2989,7 +2989,7 @@
is_valid_linkage(Z, throw=True, name='Z')
if len(T) != Z.shape[0] + 1:
raise ValueError('Mismatch: len(T)!=Z.shape[0] + 1.')
-
+
Cl = numpy.unique(T)
kk = len(Cl)
L = numpy.zeros((kk,), dtype='int32')
Modified: trunk/scipy/cluster/tests/test_hierarchy.py
===================================================================
--- trunk/scipy/cluster/tests/test_hierarchy.py 2008-04-20 08:27:24 UTC (rev 4153)
+++ trunk/scipy/cluster/tests/test_hierarchy.py 2008-04-20 12:15:19 UTC (rev 4154)
@@ -368,7 +368,7 @@
Y_test1 = pdist(X, 'minkowski', 3.2)
#print "minkowski", numpy.abs(Y_test1 - Y_right).max()
self.failUnless(within_tol(Y_test1, Y_right, eps))
-
+
def test_pdist_minkowski_random_nonC(self):
"Tests pdist(X, 'test_minkowski') [the non-C implementation] on random data."
eps = 1e-05
@@ -388,7 +388,7 @@
Y_test1 = pdist(X, 'minkowski', 3.2)
#print "minkowski-iris-3.2", numpy.abs(Y_test1 - Y_right).max()
self.failUnless(within_tol(Y_test1, Y_right, eps))
-
+
def test_pdist_minkowski_iris_nonC(self):
"Tests pdist(X, 'test_minkowski') [the non-C implementation] on iris data."
eps = 1e-07
@@ -408,7 +408,7 @@
Y_test1 = pdist(X, 'minkowski', 5.8)
#print "minkowski-iris-5.8", numpy.abs(Y_test1 - Y_right).max()
self.failUnless(within_tol(Y_test1, Y_right, eps))
-
+
def test_pdist_minkowski_iris_nonC(self):
"Tests pdist(X, 'test_minkowski') [the non-C implementation] on iris data."
eps = 1e-07
Modified: trunk/scipy/fftpack/setupscons.py
===================================================================
--- trunk/scipy/fftpack/setupscons.py 2008-04-20 08:27:24 UTC (rev 4153)
+++ trunk/scipy/fftpack/setupscons.py 2008-04-20 12:15:19 UTC (rev 4154)
@@ -10,7 +10,7 @@
config.add_sconscript('SConstruct')
config.add_data_dir('tests')
-
+
return config
if __name__ == '__main__':
Modified: trunk/scipy/integrate/ode.py
===================================================================
--- trunk/scipy/integrate/ode.py 2008-04-20 08:27:24 UTC (rev 4153)
+++ trunk/scipy/integrate/ode.py 2008-04-20 12:15:19 UTC (rev 4154)
@@ -8,7 +8,7 @@
d y(t)[i]
--------- = f(t,y(t))[i],
d t
-
+
y(t=0)[i] = y0[i],
where::
@@ -20,7 +20,7 @@
A generic interface class to numeric integrators. It has the following
methods::
-
+
integrator = ode(f,jac=None)
integrator = integrator.set_integrator(name,**params)
integrator = integrator.set_initial_value(y0,t0=0.0)
@@ -108,22 +108,22 @@
# To wrap cvode to Python, one must write extension module by
# hand. Its interface is too much 'advanced C' that using f2py
# would be too complicated (or impossible).
-#
+#
# How to define a new integrator:
# ===============================
-#
+#
# class myodeint(IntegratorBase):
-#
+#
# runner = <odeint function> or None
-#
+#
# def __init__(self,...): # required
# <initialize>
-#
+#
# def reset(self,n,has_jac): # optional
# # n - the size of the problem (number of equations)
# # has_jac - whether user has supplied its own routine for Jacobian
# <allocate memory,initialize further>
-#
+#
# def run(self,f,jac,y0,t0,t1,f_params,jac_params): # required
# # this method is called to integrate from t=t0 to t=t1
# # with initial condition y0. f and jac are user-supplied functions
@@ -134,11 +134,11 @@
# if <calculation was unsuccesful>:
# self.success = 0
# return t1,y1
-#
+#
# # In addition, one can define step() and run_relax() methods (they
# # take the same arguments as run()) if the integrator can support
# # these features (see IntegratorBase doc strings).
-#
+#
# if myodeint.runner:
# IntegratorBase.integrator_classes.append(myodeint)
@@ -158,7 +158,7 @@
A generic interface class to numeric integrators.
See also
---------
+--------
odeint : an integrator with a simpler interface based on lsoda from ODEPACK
quad : for finding the area under a curve
@@ -533,7 +533,7 @@
rwork[5] = self.max_step
rwork[6] = self.min_step
self.rwork = rwork
-
+
iwork = zeros((liw,), int32)
if self.ml is not None:
iwork[0] = self.ml
@@ -543,7 +543,7 @@
iwork[5] = self.nsteps
iwork[6] = 2 # mxhnil
self.iwork = iwork
-
+
self.call_args = [self.rtol,self.atol,1,1,
self.zwork,self.rwork,self.iwork,mf]
self.success = 1
Modified: trunk/scipy/integrate/tests/test_integrate.py
===================================================================
--- trunk/scipy/integrate/tests/test_integrate.py 2008-04-20 08:27:24 UTC (rev 4153)
+++ trunk/scipy/integrate/tests/test_integrate.py 2008-04-20 12:15:19 UTC (rev 4154)
@@ -35,7 +35,7 @@
Check integrate.ode
"""
def _do_problem(self, problem, integrator, method='adams'):
-
+
# ode has callback arguments in different order than odeint
f = lambda t, z: problem.f(z, t)
jac = None
Modified: trunk/scipy/io/__init__.py
===================================================================
--- trunk/scipy/io/__init__.py 2008-04-20 08:27:24 UTC (rev 4153)
+++ trunk/scipy/io/__init__.py 2008-04-20 12:15:19 UTC (rev 4154)
@@ -16,17 +16,17 @@
convert_objectarray
fread = deprecate_with_doc("""
-scipy.io.fread is can be replaced with raw reading capabilities of NumPy
-including fromfile as well as memory-mapping capabilities.
+scipy.io.fread is can be replaced with raw reading capabilities of NumPy
+including fromfile as well as memory-mapping capabilities.
""")(fread)
fwrite = deprecate_with_doc("""
scipy.io.fwrite can be replaced with raw writing capabilities of
NumPy. Also, remember that files can be directly memory-mapped into NumPy
-arrays which is often a better way of reading especially large files.
+arrays which is often a better way of reading especially large files.
Look at the tofile methods as well as save and savez for writing arrays into
-easily transported files of data.
+easily transported files of data.
""")(fwrite)
bswap = deprecate_with_doc("""
@@ -54,7 +54,7 @@
unpackbits = deprecate_with_doc("""
The functionality of scipy.io.unpackbits is now available in numpy.unpackbits
The calling convention is different however as the 2-d case is no longer
-specialized.
+specialized.
Thus, the scipy.unpackbits behavior must be simulated using numpy.unpackbits.
Modified: trunk/scipy/io/arff/arffread.py
===================================================================
--- trunk/scipy/io/arff/arffread.py 2008-04-20 08:27:24 UTC (rev 4153)
+++ trunk/scipy/io/arff/arffread.py 2008-04-20 12:15:19 UTC (rev 4154)
@@ -12,7 +12,7 @@
__all__ = ['MetaData', 'loadarff', 'ArffError', 'ParseArffError']
-# An Arff file is basically two parts:
+# An Arff file is basically two parts:
# - header
# - data
#
@@ -42,7 +42,7 @@
r_comattrval = re.compile(r"'(..+)'\s+(..+$)")
# To get attributes name enclosed with '', possibly spread accross multilines
r_mcomattrval = re.compile(r"'([..\n]+)'\s+(..+$)")
-# To get normal attributes
+# To get normal attributes
r_wcomattrval = re.compile(r"(\S+)\s+(..+$)")
#-------------------------
@@ -61,7 +61,7 @@
# An attribute is defined as @attribute name value
def parse_type(attrtype):
"""Given an arff attribute value (meta data), returns its type.
-
+
Expect the value to be a name."""
uattribute = attrtype.lower().strip()
if uattribute[0] == '{':
@@ -83,7 +83,7 @@
def get_nominal(attribute):
"""If attribute is nominal, returns a list of the values"""
return attribute.split(',')
-
+
def read_data_list(ofile):
"""Read each line of the iterable and put it in a list."""
data = [ofile.next()]
@@ -105,9 +105,9 @@
def maxnomlen(atrv):
"""Given a string contening a nominal type definition, returns the string
len of the biggest component.
-
+
A nominal type is defined as seomthing framed between brace ({}).
-
+
Example: maxnomlen("{floup, bouga, fl, ratata}") returns 6 (the size of
ratata, the longest nominal value)."""
nomtp = get_nom_val(atrv)
@@ -115,10 +115,10 @@
def get_nom_val(atrv):
"""Given a string contening a nominal type, returns a tuple of the possible
- values.
-
+ values.
+
A nominal type is defined as something framed between brace ({}).
-
+
Example: get_nom_val("{floup, bouga, fl, ratata}") returns ("floup",
"bouga", "fl", "ratata")."""
r_nominal = re.compile('{(..+)}')
@@ -130,7 +130,7 @@
def go_data(ofile):
"""Skip header.
-
+
the first next() call of the returned iterator will be the @data line"""
return itertools.dropwhile(lambda x : not r_datameta.match(x), ofile)
@@ -139,18 +139,18 @@
#----------------
def tokenize_attribute(iterable, attribute):
"""Parse a raw string in header (eg starts by @attribute).
-
+
Given a raw string attribute, try to get the name and type of the
attribute. Constraints:
- The first line must start with @attribute (case insensitive, and
space like characters begore @attribute are allowed)
- - Works also if the attribute is spread on multilines.
+ - Works also if the attribute is spread on multilines.
- Works if empty lines or comments are in between
-
+
:Parameters:
attribute : str
- the attribute string.
-
+ the attribute string.
+
:Returns:
name : str
name of the attribute
@@ -205,7 +205,7 @@
else:
raise ValueError("Cannot parse attribute names spread over multi "\
"lines yet")
-
+
def tokenize_single_comma(val):
# XXX we match twice the same string (here and at the caller level). It is
# stupid, but it is easier for now...
@@ -299,7 +299,7 @@
class MetaData:
"""Small container to keep useful informations on a ARFF dataset.
-
+
Knows about attributes names and types.
:Example:
@@ -318,7 +318,7 @@
Also maintains the list of attributes in order, i.e. doing for i in
meta, where meta is an instance of MetaData, will return the different
attribute names in the order they were defined.
-
+
"""
def __init__(self, rel, attr):
self.name = rel
@@ -343,7 +343,7 @@
msg += ", range is %s" % str(self._attributes[i][1])
msg += '\n'
return msg
-
+
def __iter__(self):
return iter(self._attrnames)
@@ -386,7 +386,7 @@
"""
ofile = open(filename)
- # Parse the header file
+ # Parse the header file
try:
rel, attr = read_header(ofile)
except ValueError, e:
@@ -459,9 +459,9 @@
def generator(row_iter, delim = ','):
# TODO: this is where we are spending times (~80%). I think things
- # could be made more efficiently:
+ # could be made more efficiently:
# - We could for example "compile" the function, because some values
- # do not change here.
+ # do not change here.
# - The function to convert a line to dtyped values could also be
# generated on the fly from a string and be executed instead of
# looping.
Modified: trunk/scipy/io/arff/tests/test_data.py
===================================================================
--- trunk/scipy/io/arff/tests/test_data.py 2008-04-20 08:27:24 UTC (rev 4153)
+++ trunk/scipy/io/arff/tests/test_data.py 2008-04-20 12:15:19 UTC (rev 4154)
@@ -11,8 +11,8 @@
test4 = os.path.join(data_path, 'test4.arff')
test5 = os.path.join(data_path, 'test5.arff')
-expect4_data = [(0.1, 0.2, 0.3, 0.4, 'class1'),
- (-0.1, -0.2, -0.3, -0.4, 'class2'),
+expect4_data = [(0.1, 0.2, 0.3, 0.4, 'class1'),
+ (-0.1, -0.2, -0.3, -0.4, 'class2'),
(1, 2, 3, 4, 'class3')]
missing = os.path.join(data_path, 'missing.arff')
Modified: trunk/scipy/io/data_store.py
===================================================================
--- trunk/scipy/io/data_store.py 2008-04-20 08:27:24 UTC (rev 4153)
+++ trunk/scipy/io/data_store.py 2008-04-20 12:15:19 UTC (rev 4154)
@@ -15,7 +15,7 @@
1
"""
-__all__ = ['save_as_module',
+__all__ = ['save_as_module',
# The rest of these are all deprecated
'save', 'create_module',
'create_shelf', 'load']
Modified: trunk/scipy/io/matlab/mio4.py
===================================================================
--- trunk/scipy/io/matlab/mio4.py 2008-04-20 08:27:24 UTC (rev 4153)
+++ trunk/scipy/io/matlab/mio4.py 2008-04-20 12:15:19 UTC (rev 4154)
@@ -168,13 +168,13 @@
dims = res[-1,0:2]
I = N.ascontiguousarray(tmp[:,0],dtype='intc') #fixes byte order also
J = N.ascontiguousarray(tmp[:,1],dtype='intc')
- I -= 1 # for 1-based indexing
+ I -= 1 # for 1-based indexing
J -= 1
if res.shape[1] == 3:
V = N.ascontiguousarray(tmp[:,2],dtype='float')
else:
V = N.ascontiguousarray(tmp[:,2],dtype='complex')
- V.imag = tmp[:,3]
+ V.imag = tmp[:,3]
if have_sparse:
return scipy.sparse.coo_matrix((V,(I,J)), dims)
return (dims, I, J, V)
Modified: trunk/scipy/io/matlab/tests/test_mio.py
===================================================================
--- trunk/scipy/io/matlab/tests/test_mio.py 2008-04-20 08:27:24 UTC (rev 4153)
+++ trunk/scipy/io/matlab/tests/test_mio.py 2008-04-20 12:15:19 UTC (rev 4154)
@@ -238,5 +238,3 @@
expected = case['expected']
format = case in case_table4 and '4' or '5'
yield _make_rt_check_case, name, expected, format
-
-
Modified: trunk/scipy/io/mmio.py
===================================================================
--- trunk/scipy/io/mmio.py 2008-04-20 08:27:24 UTC (rev 4153)
+++ trunk/scipy/io/mmio.py 2008-04-20 12:15:19 UTC (rev 4154)
@@ -124,7 +124,7 @@
SYMMETRY_SYMMETRIC = 'symmetric'
SYMMETRY_SKEW_SYMMETRIC = 'skew-symmetric'
SYMMETRY_HERMITIAN = 'hermitian'
- SYMMETRY_VALUES = ( SYMMETRY_GENERAL, SYMMETRY_SYMMETRIC,
+ SYMMETRY_VALUES = ( SYMMETRY_GENERAL, SYMMETRY_SYMMETRIC,
SYMMETRY_SKEW_SYMMETRIC, SYMMETRY_HERMITIAN)
@classmethod
@@ -217,7 +217,7 @@
stream = bz2.BZ2File(filespec, 'r')
else:
stream = open(filespec, mode)
-
+
# open for writing
else:
if filespec[-4:] != '.mtx':
@@ -257,7 +257,7 @@
@staticmethod
def _field_template(field, precision):
return {
- MMFile.FIELD_REAL: '%%.%ie\n' % precision,
+ MMFile.FIELD_REAL: '%%.%ie\n' % precision,
MMFile.FIELD_INTEGER: '%i\n',
MMFile.FIELD_COMPLEX: '%%.%ie %%.%ie\n' % (precision,precision)
}.get(field, None)
@@ -296,7 +296,7 @@
attrs = self.__class__.__slots__
public_attrs = [attr[1:] for attr in attrs]
invalid_keys = set(kwargs.keys()) - set(public_attrs)
-
+
if invalid_keys:
raise ValueError, \
'found %s invalid keyword arguments, please only use %s' % \
@@ -395,10 +395,10 @@
except:
# fallback - fromfile fails for some file-like objects
flat_data = fromstring(stream.read(), sep=' ')
-
+
# TODO use iterator (e.g. xreadlines) to avoid reading
# the whole file into memory
-
+
if is_pattern:
flat_data = flat_data.reshape(-1,2)
I = ascontiguousarray(flat_data[:,0], dtype='intc')
Modified: trunk/scipy/io/npfile.py
===================================================================
--- trunk/scipy/io/npfile.py 2008-04-20 08:27:24 UTC (rev 4153)
+++ trunk/scipy/io/npfile.py 2008-04-20 12:15:19 UTC (rev 4154)
@@ -225,8 +225,8 @@
npfile = N.deprecate_with_doc("""
You can achieve the same effect as using npfile, using ndarray.tofile
-and numpy.fromfile.
+and numpy.fromfile.
-Even better you can use memory-mapped arrays and data-types to map out a
+Even better you can use memory-mapped arrays and data-types to map out a
file format for direct manipulation in NumPy.
""")(npfile)
Modified: trunk/scipy/io/pickler.py
===================================================================
--- trunk/scipy/io/pickler.py 2008-04-20 08:27:24 UTC (rev 4153)
+++ trunk/scipy/io/pickler.py 2008-04-20 12:15:19 UTC (rev 4154)
@@ -25,7 +25,7 @@
fid.close()
@deprecate_with_doc("""
-Just use cPickle.load or numpy.load.
+Just use cPickle.load or numpy.load.
""")
def objload(file, allglobals):
"""Load a previously pickled dictionary and insert into given dictionary.
Modified: trunk/scipy/io/tests/test_mmio.py
===================================================================
--- trunk/scipy/io/tests/test_mmio.py 2008-04-20 08:27:24 UTC (rev 4153)
+++ trunk/scipy/io/tests/test_mmio.py 2008-04-20 12:15:19 UTC (rev 4154)
@@ -150,37 +150,37 @@
_skew_example = '''\
%%MatrixMarket matrix coordinate real skew-symmetric
5 5 7
- 1 1 1.0
- 2 2 10.5
- 4 2 250.5
- 3 3 1.5e-2
- 4 4 -2.8e2
- 5 5 12.
- 5 4 0
+ 1 1 1.0
+ 2 2 10.5
+ 4 2 250.5
+ 3 3 1.5e-2
+ 4 4 -2.8e2
+ 5 5 12.
+ 5 4 0
'''
_symmetric_example = '''\
%%MatrixMarket matrix coordinate real symmetric
5 5 7
- 1 1 1.0
- 2 2 10.5
- 4 2 250.5
- 3 3 1.5e-2
- 4 4 -2.8e2
- 5 5 12.
- 5 4 8
+ 1 1 1.0
+ 2 2 10.5
+ 4 2 250.5
+ 3 3 1.5e-2
+ 4 4 -2.8e2
+ 5 5 12.
+ 5 4 8
'''
_symmetric_pattern_example = '''\
%%MatrixMarket matrix coordinate pattern symmetric
5 5 7
- 1 1
- 2 2
- 4 2
- 3 3
- 4 4
- 5 5
- 5 4
+ 1 1
+ 2 2
+ 4 2
+ 3 3
+ 4 4
+ 5 5
+ 5 4
'''
class TestMMIOCoordinate(TestCase):
@@ -213,7 +213,7 @@
[0, 0, 0, 33.32j, 12]]
b = mmread(fn).todense()
assert_array_almost_equal(a,b)
-
+
def test_read_skew(self):
"""read a skew-symmetric matrix"""
fn = mktemp()
@@ -228,7 +228,7 @@
[0, 0, 0, 0, 12]]
b = mmread(fn).todense()
assert_array_almost_equal(a,b)
-
+
def test_read_symmetric(self):
"""read a symmetric matrix"""
fn = mktemp()
Modified: trunk/scipy/io/tests/test_recaster.py
===================================================================
--- trunk/scipy/io/tests/test_recaster.py 2008-04-20 08:27:24 UTC (rev 4153)
+++ trunk/scipy/io/tests/test_recaster.py 2008-04-20 12:15:19 UTC (rev 4154)
@@ -171,7 +171,6 @@
dtt = arr.dtype.type
assert dtt is outp, \
'Expected %s from %s, got %s' % (outp, inp, dtt)
-
+
if __name__ == "__main__":
nose.run(argv=['', __file__])
-
Modified: trunk/scipy/lib/blas/scons_support.py
===================================================================
--- trunk/scipy/lib/blas/scons_support.py 2008-04-20 08:27:24 UTC (rev 4153)
+++ trunk/scipy/lib/blas/scons_support.py 2008-04-20 12:15:19 UTC (rev 4154)
@@ -15,7 +15,7 @@
# XXX handle skip names
name = splitext(pbasename(target_name))[0]
#generate_interface(name, source_name, target_name)
-
+
f = open(target_name, 'w')
f.write('python module '+name+'\n')
f.write('usercode void empty_module(void) {}\n')
@@ -27,4 +27,3 @@
f.close()
return 0
-
Modified: trunk/scipy/lib/lapack/scons_support.py
===================================================================
--- trunk/scipy/lib/lapack/scons_support.py 2008-04-20 08:27:24 UTC (rev 4153)
+++ trunk/scipy/lib/lapack/scons_support.py 2008-04-20 12:15:19 UTC (rev 4154)
@@ -15,7 +15,7 @@
# XXX handle skip names
name = splitext(pbasename(target_name))[0]
#generate_interface(name, source_name, target_name)
-
+
f = open(target_name, 'w')
f.write('python module '+name+'\n')
f.write('usercode void empty_module(void) {}\n')
@@ -27,4 +27,3 @@
f.close()
return 0
-
Modified: trunk/scipy/lib/lapack/tests/test_lapack.py
===================================================================
--- trunk/scipy/lib/lapack/tests/test_lapack.py 2008-04-20 08:27:24 UTC (rev 4153)
+++ trunk/scipy/lib/lapack/tests/test_lapack.py 2008-04-20 12:15:19 UTC (rev 4154)
@@ -11,7 +11,7 @@
method names. There are no subclasses of TestCase. Thus nose will
pick up nothing but the final test_all_lapack generator function.
This does the work of collecting the test methods and checking if they
-can be run (see the isrunnable method).
+can be run (see the isrunnable method).
'''
import os
@@ -137,4 +137,3 @@
methods += [getattr(o, n) for n in dir(o) if o.isrunnable(n) is True]
for method in methods:
yield (method, )
-
Modified: trunk/scipy/lib/setupscons.py
===================================================================
--- trunk/scipy/lib/setupscons.py 2008-04-20 08:27:24 UTC (rev 4153)
+++ trunk/scipy/lib/setupscons.py 2008-04-20 12:15:19 UTC (rev 4154)
@@ -3,7 +3,7 @@
def configuration(parent_package='',top_path=None):
from numpy.distutils.misc_util import Configuration
- config = Configuration('lib',parent_package,top_path,
+ config = Configuration('lib',parent_package,top_path,
setup_name = 'setupscons.py')
config.add_subpackage('blas')
config.add_subpackage('lapack')
Modified: trunk/scipy/linalg/basic.py
===================================================================
--- trunk/scipy/linalg/basic.py 2008-04-20 08:27:24 UTC (rev 4153)
+++ trunk/scipy/linalg/basic.py 2008-04-20 12:15:19 UTC (rev 4154)
@@ -26,7 +26,7 @@
def lu_solve((lu, piv), b, trans=0, overwrite_b=0):
"""Solve an equation system, a x = b, given the LU factorization of a
-
+
Parameters
----------
(lu, piv)
@@ -52,7 +52,7 @@
See also
--------
lu_factor : LU factorize a matrix
-
+
"""
b1 = asarray_chkfinite(b)
overwrite_b = overwrite_b or (b1 is not b and not hasattr(b,'__array__'))
@@ -83,7 +83,7 @@
See also
--------
cho_factor : Cholesky factorization of a matrix
-
+
"""
b1 = asarray_chkfinite(b)
overwrite_b = overwrite_b or (b1 is not b and not hasattr(b,'__array__'))
@@ -114,14 +114,14 @@
Allow overwriting data in a (may enhance performance)
overwrite_b : boolean
Allow overwriting data in b (may enhance performance)
-
+
Returns
-------
x : array, shape (M,) or (M, N) depending on b
Solution to the system a x = b
Raises LinAlgError if a is singular
-
+
"""
a1, b1 = map(asarray_chkfinite,(a,b))
if len(a1.shape) != 2 or a1.shape[0] != a1.shape[1]:
@@ -164,8 +164,8 @@
* a01 a12 a23 a34 a45
a00 a11 a22 a33 a44 a55
- a10 a21 a32 a43 a54 *
- a20 a31 a42 a53 * *
+ a10 a21 a32 a43 a54 *
+ a20 a31 a42 a53 * *
Parameters
----------
@@ -184,7 +184,7 @@
-------
x : array, shape (M,) or (M, K)
The solution to the system a x = b
-
+
"""
a1, b1 = map(asarray_chkfinite,(ab,b))
overwrite_b = overwrite_b or (b1 is not b and not hasattr(b,'__array__'))
@@ -218,7 +218,7 @@
* * a02 a13 a24 a35
* a01 a12 a23 a34 a45
a00 a11 a22 a33 a44 a55
-
+
lower form:
a00 a11 a22 a33 a44 a55
a10 a21 a32 a43 a54 *
@@ -245,7 +245,7 @@
Cholesky factorization of a, in the same banded format as ab
x : array, shape (M,) or (M, K)
The solution to the system a x = b
-
+
"""
ab, b = map(asarray_chkfinite,(ab,b))
@@ -263,25 +263,25 @@
def cholesky_banded(ab, overwrite_ab=0, lower=0):
"""Cholesky decompose a banded Hermitian positive-definite matrix
-
+
The matrix a is stored in ab either in lower diagonal or upper
diagonal ordered form:
-
+
ab[u + i - j, j] == a[i,j] (if upper form; i <= j)
ab[ i - j, j] == a[i,j] (if lower form; i >= j)
-
+
Example of ab (shape of a is (6,6), u=2)::
-
+
upper form:
* * a02 a13 a24 a35
* a01 a12 a23 a34 a45
a00 a11 a22 a33 a44 a55
-
+
lower form:
a00 a11 a22 a33 a44 a55
a10 a21 a32 a43 a54 *
a20 a31 a42 a53 * *
-
+
Parameters
----------
ab : array, shape (M, u + 1)
@@ -290,12 +290,12 @@
Discard data in ab (may enhance performance)
lower : boolean
Is the matrix in the lower form. (Default is upper form)
-
+
Returns
-------
c : array, shape (M, u+1)
Cholesky factorization of a, in the same banded format as ab
-
+
"""
ab = asarray_chkfinite(ab)
@@ -315,12 +315,12 @@
# matrix inversion
def inv(a, overwrite_a=0):
"""Compute the inverse of a matrix.
-
+
Parameters
----------
a : array-like, shape (M, M)
Matrix to be inverted
-
+
Returns
-------
ainv : array-like, shape (M, M)
@@ -409,7 +409,7 @@
-2 smallest singular value as below
other - sum(abs(x)**ord)**(1./ord)
===== ============================ ==========================
-
+
Returns
-------
n : float
@@ -420,7 +420,7 @@
For values ord < 0, the result is, strictly speaking, not a
mathematical 'norm', but it may still be useful for numerical
purposes.
-
+
"""
x = asarray_chkfinite(x)
if ord is None: # check the default case first and handle it immediately
@@ -472,7 +472,7 @@
-------
det : float or complex
Determinant of a
-
+
Notes
-----
The determinant is computed via LU factorization, LAPACK routine z/dgetrf.
@@ -491,9 +491,9 @@
def lstsq(a, b, cond=None, overwrite_a=0, overwrite_b=0):
"""Compute least-squares solution to equation :m:`a x = b`
-
+
Compute a vector x such that the 2-norm :m:`|b - a x|` is minimised.
-
+
Parameters
----------
a : array, shape (M, N)
@@ -506,7 +506,7 @@
Discard data in a (may enhance performance)
overwrite_b : boolean
Discard data in b (may enhance performance)
-
+
Returns
-------
x : array, shape (N,) or (N, K) depending on shape of b
@@ -519,9 +519,9 @@
Effective rank of matrix a
s : array, shape (min(M,N),)
Singular values of a. The condition number of a is abs(s[0]/s[-1]).
-
+
Raises LinAlgError if computation does not converge
-
+
"""
a1, b1 = map(asarray_chkfinite,(a,b))
if len(a1.shape) != 2:
@@ -562,10 +562,10 @@
def pinv(a, cond=None, rcond=None):
"""Compute the (Moore-Penrose) pseudo-inverse of a matrix.
-
+
Calculate a generalized inverse of a matrix using a least-squares
solver.
-
+
Parameters
----------
a : array, shape (M, N)
@@ -574,11 +574,11 @@
Cutoff for 'small' singular values in the least-squares solver.
Singular values smaller than rcond*largest_singular_value are
considered zero.
-
+
Returns
-------
B : array, shape (N, M)
-
+
Raises LinAlgError if computation does not converge
Examples
@@ -590,7 +590,7 @@
True
>>> allclose(B, dot(B, dot(a, B)))
True
-
+
"""
a = asarray_chkfinite(a)
b = numpy.identity(a.shape[0], dtype=a.dtype)
@@ -605,11 +605,11 @@
_array_precision = {'f': 0, 'd': 1, 'F': 0, 'D': 1}
def pinv2(a, cond=None, rcond=None):
"""Compute the (Moore-Penrose) pseudo-inverse of a matrix.
-
+
Calculate a generalized inverse of a matrix using its
singular-value decomposition and including all 'large' singular
values.
-
+
Parameters
----------
a : array, shape (M, N)
@@ -620,11 +620,11 @@
considered zero.
If None or -1, suitable machine precision is used.
-
+
Returns
-------
B : array, shape (N, M)
-
+
Raises LinAlgError if SVD computation does not converge
Examples
@@ -636,7 +636,7 @@
True
>>> allclose(B, dot(B, dot(a, B)))
True
-
+
"""
a = asarray_chkfinite(a)
u, s, vh = decomp.svd(a)
@@ -689,7 +689,7 @@
array([[0, 0, 0, 0, 0],
[1, 0, 0, 0, 0],
[1, 1, 0, 0, 0]])
-
+
"""
if M is None: M = N
if type(M) == type('d'):
@@ -705,7 +705,7 @@
def tril(m, k=0):
"""Construct a copy of a matrix with elements above the k-th diagonal zeroed.
-
+
Parameters
----------
m : array
@@ -717,7 +717,7 @@
Returns
-------
A : array, shape m.shape, dtype m.dtype
-
+
Examples
--------
>>> from scipy.linalg import tril
@@ -726,7 +726,7 @@
[ 4, 0, 0],
[ 7, 8, 0],
[10, 11, 12]])
-
+
"""
svsp = getattr(m,'spacesaver',lambda:0)()
m = asarray(m)
@@ -736,7 +736,7 @@
def triu(m, k=0):
"""Construct a copy of a matrix with elements below the k-th diagonal zeroed.
-
+
Parameters
----------
m : array
@@ -748,7 +748,7 @@
Returns
-------
A : array, shape m.shape, dtype m.dtype
-
+
Examples
--------
>>> from scipy.linalg import tril
@@ -757,7 +757,7 @@
[ 4, 5, 6],
[ 0, 8, 9],
[ 0, 0, 12]])
-
+
"""
svsp = getattr(m,'spacesaver',lambda:0)()
m = asarray(m)
@@ -767,23 +767,23 @@
def toeplitz(c,r=None):
"""Construct a Toeplitz matrix.
-
+
The Toepliz matrix has constant diagonals, c as its first column,
and r as its first row (if not given, r == c is assumed).
-
+
Parameters
----------
c : array
First column of the matrix
r : array
First row of the matrix. If None, r == c is assumed.
-
+
Returns
-------
A : array, shape (len(c), len(r))
Constructed Toeplitz matrix.
dtype is the same as (c[0] + r[0]).dtype
-
+
Examples
--------
>>> from scipy.linalg import toeplitz
@@ -791,11 +791,11 @@
array([[1, 4, 5, 6],
[2, 1, 4, 5],
[3, 2, 1, 4]])
-
+
See also
--------
hankel : Hankel matrix
-
+
"""
isscalar = numpy.isscalar
if isscalar(c) or isscalar(r):
@@ -819,23 +819,23 @@
def hankel(c,r=None):
"""Construct a Hankel matrix.
-
+
The Hankel matrix has constant anti-diagonals, c as its first column,
and r as its last row (if not given, r == 0 os assumed).
-
+
Parameters
----------
c : array
First column of the matrix
r : array
Last row of the matrix. If None, r == 0 is assumed.
-
+
Returns
-------
A : array, shape (len(c), len(r))
Constructed Hankel matrix.
dtype is the same as (c[0] + r[0]).dtype
-
+
Examples
--------
>>> from scipy.linalg import hankel
@@ -844,11 +844,11 @@
[2, 3, 4, 7, 7],
[3, 4, 7, 7, 8],
[4, 7, 7, 8, 9]])
-
+
See also
--------
toeplitz : Toeplitz matrix
-
+
"""
isscalar = numpy.isscalar
if isscalar(c) or isscalar(r):
@@ -889,14 +889,14 @@
-------
A : array, shape (M*P, N*Q)
Kronecker product of a and b
-
+
Examples
--------
>>> from scipy import kron, array
>>> kron(array([[1,2],[3,4]]), array([[1,1,1]]))
array([[1, 1, 1, 2, 2, 2],
[3, 3, 3, 4, 4, 4]])
-
+
"""
if not a.flags['CONTIGUOUS']:
a = reshape(a, a.shape)
Modified: trunk/scipy/linalg/benchmarks/bench_decom.py
===================================================================
--- trunk/scipy/linalg/benchmarks/bench_decom.py 2008-04-20 08:27:24 UTC (rev 4153)
+++ trunk/scipy/linalg/benchmarks/bench_decom.py 2008-04-20 12:15:19 UTC (rev 4154)
@@ -20,16 +20,15 @@
print ' | contiguous '#'| non-contiguous '
print '----------------------------------------------'
print ' size | scipy '#'| core | scipy | core '
-
+
for size,repeat in [(20,150),(100,7),(200,2)]:
repeat *= 1
print '%5s' % size,
sys.stdout.flush()
-
+
a = random([size,size])
-
+
print '| %6.2f ' % measure('eigvals(a)',repeat),
sys.stdout.flush()
-
- print ' (secs for %s calls)' % (repeat)
+ print ' (secs for %s calls)' % (repeat)
Modified: trunk/scipy/linalg/blas.py
===================================================================
--- trunk/scipy/linalg/blas.py 2008-04-20 08:27:24 UTC (rev 4153)
+++ trunk/scipy/linalg/blas.py 2008-04-20 12:15:19 UTC (rev 4154)
@@ -29,7 +29,7 @@
"""Return available BLAS function objects with names.
arrays are used to determine the optimal prefix of
BLAS routines.
-
+
"""
ordering = []
for i in range(len(arrays)):
Modified: trunk/scipy/linalg/decomp.py
===================================================================
--- trunk/scipy/linalg/decomp.py 2008-04-20 08:27:24 UTC (rev 4153)
+++ trunk/scipy/linalg/decomp.py 2008-04-20 12:15:19 UTC (rev 4154)
@@ -101,12 +101,12 @@
"""Solve an ordinary or generalized eigenvalue problem of a square matrix.
Find eigenvalues w and right or left eigenvectors of a general matrix::
-
+
a vr[:,i] = w[i] b vr[:,i]
a.H vl[:,i] = w[i].conj() b.H vl[:,i]
-
+
where .H is the Hermitean conjugation.
-
+
Parameters
----------
a : array, shape (M, M)
@@ -119,12 +119,12 @@
Whether to calculate and return left eigenvectors
right : boolean
Whether to calculate and return right eigenvectors
-
+
overwrite_a : boolean
Whether to overwrite data in a (may improve performance)
overwrite_b : boolean
Whether to overwrite data in b (may improve performance)
-
+
Returns
-------
w : double or complex array, shape (M,)
@@ -134,18 +134,18 @@
vl : double or complex array, shape (M, M)
The normalized left eigenvector corresponding to the eigenvalue w[i]
is the column v[:,i].
-
+
(if right == True)
vr : double or complex array, shape (M, M)
The normalized right eigenvector corresponding to the eigenvalue w[i]
is the column vr[:,i].
-
+
Raises LinAlgError if eigenvalue computation does not converge
See Also
--------
eigh : eigenvalues and right eigenvectors for symmetric/Hermitian arrays
-
+
"""
a1 = asarray_chkfinite(a)
if len(a1.shape) != 2 or a1.shape[0] != a1.shape[1]:
@@ -207,10 +207,10 @@
"""Solve the eigenvalue problem for a Hermitian or real symmetric matrix.
Find eigenvalues w and optionally right eigenvectors v of a::
-
+
a v[:,i] = w[i] v[:,i]
v.H v = identity
-
+
Parameters
----------
a : array, shape (M, M)
@@ -224,24 +224,24 @@
(Default: both are calculated)
overwrite_a : boolean
Whether data in a is overwritten (may improve performance).
-
+
Returns
-------
w : double array, shape (M,)
The eigenvalues, in ascending order, each repeated according to its
multiplicity.
-
+
(if eigvals_only == False)
v : double or complex double array, shape (M, M)
The normalized eigenvector corresponding to the eigenvalue w[i] is
the column v[:,i].
-
+
Raises LinAlgError if eigenvalue computation does not converge
See Also
--------
eig : eigenvalues and right eigenvectors for non-symmetric arrays
-
+
"""
if eigvals_only or overwrite_a:
a1 = asarray_chkfinite(a)
@@ -299,10 +299,10 @@
"""Solve real symmetric or complex hermetian band matrix eigenvalue problem.
Find eigenvalues w and optionally right eigenvectors v of a::
-
+
a v[:,i] = w[i] v[:,i]
v.H v = identity
-
+
The matrix a is stored in ab either in lower diagonal or upper
diagonal ordered form:
@@ -315,7 +315,7 @@
* * a02 a13 a24 a35
* a01 a12 a23 a34 a45
a00 a11 a22 a33 a44 a55
-
+
lower form:
a00 a11 a22 a33 a44 a55
a10 a21 a32 a43 a54 *
@@ -349,7 +349,7 @@
max_ev : integer
For select=='v', maximum number of eigenvalues expected.
For other values of select, has no meaning.
-
+
In doubt, leave this parameter untouched.
Returns
@@ -361,9 +361,9 @@
v : double or complex double array, shape (M, M)
The normalized eigenvector corresponding to the eigenvalue w[i] is
the column v[:,i].
-
+
Raises LinAlgError if eigenvalue computation does not converge
-
+
"""
if eigvals_only or overwrite_a_band:
a1 = asarray_chkfinite(a_band)
@@ -446,7 +446,7 @@
"""Compute eigenvalues from an ordinary or generalized eigenvalue problem.
Find eigenvalues of a general matrix::
-
+
a vr[:,i] = w[i] b vr[:,i]
Parameters
@@ -465,9 +465,9 @@
w : double or complex array, shape (M,)
The eigenvalues, each repeated according to its multiplicity,
but not in any specific order.
-
+
Raises LinAlgError if eigenvalue computation does not converge
-
+
See Also
--------
eigvalsh : eigenvalues of symmetric or Hemitiean arrays
@@ -481,10 +481,10 @@
"""Solve the eigenvalue problem for a Hermitian or real symmetric matrix.
Find eigenvalues w of a::
-
+
a v[:,i] = w[i] v[:,i]
v.H v = identity
-
+
Parameters
----------
a : array, shape (M, M)
@@ -495,13 +495,13 @@
triangle of a. (Default: lower)
overwrite_a : boolean
Whether data in a is overwritten (may improve performance).
-
+
Returns
-------
w : double array, shape (M,)
The eigenvalues, in ascending order, each repeated according to its
multiplicity.
-
+
Raises LinAlgError if eigenvalue computation does not converge
See Also
@@ -509,7 +509,7 @@
eigvals : eigenvalues of general arrays
eigh : eigenvalues and right eigenvectors for symmetric/Hermitian arrays
eig : eigenvalues and right eigenvectors for non-symmetric arrays
-
+
"""
return eigh(a,lower=lower,eigvals_only=1,overwrite_a=overwrite_a)
@@ -518,10 +518,10 @@
"""Solve real symmetric or complex hermetian band matrix eigenvalue problem.
Find eigenvalues w of a::
-
+
a v[:,i] = w[i] v[:,i]
v.H v = identity
-
+
The matrix a is stored in ab either in lower diagonal or upper
diagonal ordered form:
@@ -534,7 +534,7 @@
* * a02 a13 a24 a35
* a01 a12 a23 a34 a45
a00 a11 a22 a33 a44 a55
-
+
lower form:
a00 a11 a22 a33 a44 a55
a10 a21 a32 a43 a54 *
@@ -570,14 +570,14 @@
multiplicity.
Raises LinAlgError if eigenvalue computation does not converge
-
+
See Also
--------
eig_banded : eigenvalues and right eigenvectors for symmetric/Hermitian band matrices
eigvals : eigenvalues of general arrays
eigh : eigenvalues and right eigenvectors for symmetric/Hermitian arrays
eig : eigenvalues and right eigenvectors for non-symmetric arrays
-
+
"""
return eig_banded(a_band,lower=lower,eigvals_only=1,
overwrite_a_band=overwrite_a_band, select=select,
@@ -585,14 +585,14 @@
def lu_factor(a, overwrite_a=0):
"""Compute pivoted LU decomposition of a matrix.
-
+
The decomposition is::
A = P L U
where P is a permutation matrix, L lower triangular with unit
diagonal elements, and U upper triangular.
-
+
Parameters
----------
a : array, shape (M, M)
@@ -612,11 +612,11 @@
See also
--------
lu_solve : solve an equation system using the LU factorization of a matrix
-
+
Notes
-----
This is a wrapper to the *GETRF routines from LAPACK.
-
+
"""
a1 = asarray(a)
if len(a1.shape) != 2 or (a1.shape[0] != a1.shape[1]):
@@ -632,7 +632,7 @@
def lu_solve(a_lu_pivots,b):
"""Solve an equation system, a x = b, given the LU factorization of a
-
+
Parameters
----------
(lu, piv)
@@ -676,7 +676,7 @@
where P is a permutation matrix, L lower triangular with unit
diagonal elements, and U upper triangular.
-
+
Parameters
----------
a : array, shape (M, N)
@@ -696,18 +696,18 @@
K = min(M, N)
u : array, shape (K, N)
Upper triangular or trapezoidal matrix
-
+
(If permute_l == True)
pl : array, shape (M, K)
Permuted L matrix.
K = min(M, N)
u : array, shape (K, N)
Upper triangular or trapezoidal matrix
-
+
Notes
-----
This is a LU factorization routine written for Scipy.
-
+
"""
a1 = asarray_chkfinite(a)
if len(a1.shape) != 2:
@@ -728,7 +728,7 @@
an 1d-array s of singular values (real, non-negative) such that
a == U S Vh if S is an suitably shaped matrix of zeros whose
main diagonal is s.
-
+
Parameters
----------
a : array, shape (M, N)
@@ -740,7 +740,7 @@
Whether to compute also U, Vh in addition to s (Default: true)
overwrite_a : boolean
Whether data in a is overwritten (may improve performance)
-
+
Returns
-------
U: array, shape (M,M) or (M,K) depending on full_matrices
@@ -759,14 +759,14 @@
>>> U, s, Vh = linalg.svd(a)
>>> U.shape, Vh.shape, s.shape
((9, 9), (6, 6), (6,))
-
+
>>> U, s, Vh = linalg.svd(a, full_matrices=False)
>>> U.shape, Vh.shape, s.shape
((9, 6), (6, 6), (6,))
>>> S = linalg.diagsvd(s, 6, 6)
>>> allclose(a, dot(U, dot(S, Vh)))
True
-
+
>>> s2 = linalg.svd(a, compute_uv=False)
>>> allclose(s, s2)
True
@@ -775,7 +775,7 @@
--------
svdvals : return singular values of a matrix
diagsvd : return the Sigma matrix, given the vector s
-
+
"""
# A hack until full_matrices == 0 support is fixed here.
if full_matrices == 0:
@@ -810,7 +810,7 @@
Matrix to decompose
overwrite_a : boolean
Whether data in a is overwritten (may improve performance)
-
+
Returns
-------
s: array, shape (K,)
@@ -822,7 +822,7 @@
--------
svd : return the full singular value decomposition of a matrix
diagsvd : return the Sigma matrix, given the vector s
-
+
"""
return svd(a,compute_uv=0,overwrite_a=overwrite_a)
@@ -841,7 +841,7 @@
-------
S : array, shape (M, N)
The S-matrix in the singular value decomposition
-
+
"""
part = diag(s)
typ = part.dtype.char
@@ -855,10 +855,10 @@
def cholesky(a,lower=0,overwrite_a=0):
"""Compute the Cholesky decomposition of a matrix.
-
+
Returns the Cholesky decomposition, :lm:`A = L L^*` or :lm:`A = U^* U`
of a Hermitian positive-definite matrix :lm:`A`.
-
+
Parameters
----------
a : array, shape (M, M)
@@ -868,14 +868,14 @@
(Default: upper-triangular)
overwrite_a : boolean
Whether to overwrite data in a (may improve performance)
-
+
Returns
-------
B : array, shape (M, M)
Upper- or lower-triangular Cholesky factor of A
-
+
Raises LinAlgError if decomposition fails
-
+
Examples
--------
>>> from scipy import array, linalg, dot
@@ -887,7 +887,7 @@
>>> dot(L, L.T.conj())
array([[ 1.+0.j, 0.-2.j],
[ 0.+2.j, 5.+0.j]])
-
+
"""
a1 = asarray_chkfinite(a)
if len(a1.shape) != 2 or a1.shape[0] != a1.shape[1]:
@@ -902,12 +902,12 @@
def cho_factor(a, lower=0, overwrite_a=0):
"""Compute the Cholesky decomposition of a matrix, to use in cho_solve
-
+
Returns the Cholesky decomposition, :lm:`A = L L^*` or :lm:`A = U^* U`
of a Hermitian positive-definite matrix :lm:`A`.
The return value can be directly used as the first parameter to cho_solve.
-
+
Parameters
----------
a : array, shape (M, M)
@@ -917,16 +917,16 @@
(Default: upper-triangular)
overwrite_a : boolean
Whether to overwrite data in a (may improve performance)
-
+
Returns
-------
c : array, shape (M, M)
Upper- or lower-triangular Cholesky factor of A
lower : array, shape (M, M)
Flag indicating whether the factor is lower or upper triangular
-
+
Raises LinAlgError if decomposition fails
-
+
"""
a1 = asarray_chkfinite(a)
if len(a1.shape) != 2 or a1.shape[0] != a1.shape[1]:
@@ -955,7 +955,7 @@
The return value from cho_factor can be used.
b : array
Right-hand side of the equation system
-
+
First input is a tuple (LorU, lower) which is the output to cho_factor.
Second input is the right-hand side.
@@ -963,7 +963,7 @@
-------
x : array
Solution to the equation system
-
+
"""
c, lower = clow
c = asarray_chkfinite(c)
@@ -1000,7 +1000,7 @@
mode : {'qr', 'r'}
Determines what information is to be returned: either both Q and R
or only R.
-
+
Returns
-------
(if mode == 'qr')
@@ -1029,11 +1029,11 @@
>>> r2 = linalg.qr(a, mode='r')
>>> allclose(r, r2)
-
+
>>> q3, r3 = linalg.qr(a, econ=True)
>>> q3.shape, r3.shape
((9, 6), (6, 6))
-
+
"""
a1 = asarray_chkfinite(a)
if len(a1.shape) != 2:
@@ -1108,7 +1108,7 @@
lwork : integer
Work array size, lwork >= a.shape[1]. If None or -1, an optimal size
is computed.
-
+
Returns
-------
Q : double or complex array, shape (M, M)
@@ -1116,7 +1116,7 @@
Size K = min(M, N)
Raises LinAlgError if decomposition fails
-
+
"""
a1 = asarray_chkfinite(a)
if len(a1.shape) != 2:
@@ -1163,15 +1163,15 @@
Work array size, lwork >= a.shape[1]. If None or -1, an optimal size
is computed.
econ : boolean
-
+
Returns
-------
R : double array, shape (M, N) or (K, N) for econ==True
Size K = min(M, N)
Q : double or complex array, shape (M, M) or (M, K) for econ==True
-
+
Raises LinAlgError if decomposition fails
-
+
"""
# TODO: implement support for non-square and complex arrays
a1 = asarray_chkfinite(a)
@@ -1211,16 +1211,16 @@
def schur(a,output='real',lwork=None,overwrite_a=0):
"""Compute Schur decomposition of a matrix.
-
+
The Schur decomposition is
-
+
A = Z T Z^H
-
+
where Z is unitary and T is either upper-triangular, or for real
Schur decomposition (output='real'), quasi-upper triangular. In
the quasi-triangular form, 2x2 blocks describing complex-valued
eigenvalue pairs may extrude from the diagonal.
-
+
Parameters
----------
a : array, shape (M, M)
@@ -1243,7 +1243,7 @@
See also
--------
rsf2csf : Convert real Schur form to complex Schur form
-
+
"""
if not output in ['real','complex','r','c']:
raise ValueError, "argument must be 'real', or 'complex'"
@@ -1305,28 +1305,28 @@
def rsf2csf(T, Z):
"""Convert real Schur form to complex Schur form.
-
+
Convert a quasi-diagonal real-valued Schur form to the upper triangular
complex-valued Schur form.
-
+
Parameters
----------
T : array, shape (M, M)
Real Schur form of the original matrix
Z : array, shape (M, M)
Schur transformation matrix
-
+
Returns
-------
T : array, shape (M, M)
Complex Schur form of the original matrix
Z : array, shape (M, M)
Schur transformation matrix corresponding to the complex form
-
+
See also
--------
schur : Schur decompose a matrix
-
+
"""
Z,T = map(asarray_chkfinite, (Z,T))
if len(Z.shape) !=2 or Z.shape[0] != Z.shape[1]:
@@ -1366,21 +1366,21 @@
def orth(A):
"""Construct an orthonormal basis for the range of A using SVD
-
+
Parameters
----------
A : array, shape (M, N)
-
+
Returns
-------
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
-
+
"""
u,s,vh = svd(A)
M,N = A.shape
@@ -1391,14 +1391,14 @@
def hessenberg(a,calc_q=0,overwrite_a=0):
"""Compute Hessenberg form of a matrix.
-
+
The Hessenberg decomposition is
-
+
A = Q H Q^H
-
+
where Q is unitary/orthogonal and H has only zero elements below the first
subdiagonal.
-
+
Parameters
----------
a : array, shape (M,M)
Modified: trunk/scipy/linalg/iterative.py
===================================================================
--- trunk/scipy/linalg/iterative.py 2008-04-20 08:27:24 UTC (rev 4153)
+++ trunk/scipy/linalg/iterative.py 2008-04-20 12:15:19 UTC (rev 4154)
@@ -7,7 +7,7 @@
for name in __all__:
oldfn = getattr(isolve, name)
- oldname='scipy.linalg.' + name
+ oldname='scipy.linalg.' + name
newname='scipy.sparse.linalg.' + name
newfn = deprecate(oldfn, oldname=oldname, newname=newname)
exec(name + ' = newfn')
Modified: trunk/scipy/linalg/matfuncs.py
===================================================================
--- trunk/scipy/linalg/matfuncs.py 2008-04-20 08:27:24 UTC (rev 4153)
+++ trunk/scipy/linalg/matfuncs.py 2008-04-20 12:15:19 UTC (rev 4154)
@@ -104,7 +104,7 @@
Matrix to be exponentiated
q : integer
Order of the Taylor series
-
+
Returns
-------
expA : array, shape(M,M)
@@ -154,7 +154,7 @@
Parameters
----------
A : array, shape(M,M)
-
+
Returns
-------
cosA : array, shape(M,M)
@@ -275,10 +275,10 @@
def funm(A,func,disp=1):
"""Evaluate a matrix function specified by a callable.
-
+
Returns the value of matrix-valued function f at A. The function f
is an extension of the scalar-valued function func to matrices.
-
+
Parameters
----------
A : array, shape(M,M)
@@ -350,9 +350,9 @@
def logm(A,disp=1):
"""Compute matrix logarithm.
-
+
The matrix logarithm is the inverse of expm: expm(logm(A)) == A
-
+
Parameters
----------
A : array, shape(M,M)
@@ -401,9 +401,9 @@
def signm(a,disp=1):
"""Matrix sign function.
-
+
Extension of the scalar sign(x) to matrices.
-
+
Parameters
----------
A : array, shape(M,M)
@@ -411,7 +411,7 @@
disp : boolean
Print warning if error in the result is estimated large
instead of returning estimated error. (Default: True)
-
+
Returns
-------
sgnA : array, shape(M,M)
@@ -420,7 +420,7 @@
(if disp == False)
errest : float
1-norm of the estimated error, ||err||_1 / ||A||_1
-
+
Examples
--------
>>> from scipy.linalg import signm, eigvals
@@ -429,7 +429,7 @@
array([ 4.12488542+0.j, -0.76155718+0.j, 0.63667176+0.j])
>>> eigvals(signm(a))
array([-1.+0.j, 1.+0.j, 1.+0.j])
-
+
"""
def rounded_sign(x):
rx = real(x)
@@ -478,7 +478,7 @@
def sqrtm(A,disp=1):
"""Matrix square root.
-
+
Parameters
----------
A : array, shape(M,M)
@@ -486,7 +486,7 @@
disp : boolean
Print warning if error in the result is estimated large
instead of returning estimated error. (Default: True)
-
+
Returns
-------
sgnA : array, shape(M,M)
@@ -499,7 +499,7 @@
Notes
-----
Uses algorithm by Nicholas J. Higham
-
+
"""
A = asarray(A)
if len(A.shape)!=2:
Modified: trunk/scipy/linalg/scons_support.py
===================================================================
--- trunk/scipy/linalg/scons_support.py 2008-04-20 08:27:24 UTC (rev 4153)
+++ trunk/scipy/linalg/scons_support.py 2008-04-20 12:15:19 UTC (rev 4154)
@@ -28,7 +28,7 @@
# XXX handle skip names
name = splitext(pbasename(target_name))[0]
generate_interface(name, source_name, target_name)
-
+
f = open(target_name, 'w')
f.write('python module '+name+'\n')
f.write('usercode void empty_module(void) {}\n')
@@ -40,4 +40,3 @@
f.close()
return 0
-
Modified: trunk/scipy/linalg/tests/test_decomp.py
===================================================================
--- trunk/scipy/linalg/tests/test_decomp.py 2008-04-20 08:27:24 UTC (rev 4153)
+++ trunk/scipy/linalg/tests/test_decomp.py 2008-04-20 12:15:19 UTC (rev 4154)
@@ -11,7 +11,7 @@
Run tests if scipy is installed:
python -c 'import scipy;scipy.linalg.test()'
Run tests if linalg is not installed:
- python tests/test_decomp.py
+ python tests/test_decomp.py
"""
import sys
Modified: trunk/scipy/ndimage/_registration.py
===================================================================
--- trunk/scipy/ndimage/_registration.py 2008-04-20 08:27:24 UTC (rev 4153)
+++ trunk/scipy/ndimage/_registration.py 2008-04-20 12:15:19 UTC (rev 4154)
@@ -30,951 +30,951 @@
#
def resize_image(imageG, imageF_mat):
- """
- zoom_image = resize_image(source_image, reference_image[mat])
+ """
+ zoom_image = resize_image(source_image, reference_image[mat])
- Fractional resample source_image to reference_imagesize. The
- resample is implemented with 3D cubic spline. The reference
- image [mat] is the 4x4 voxel-to-physical conversion matrix.
-
- Parameters
- ----------
+ Fractional resample source_image to reference_imagesize. The
+ resample is implemented with 3D cubic spline. The reference
+ image [mat] is the 4x4 voxel-to-physical conversion matrix.
+
+ Parameters
+ ----------
- imageG : {dictionary}
- imageG is the source image to be resized. it is a dictionary with
- the data as an ndarray in the ['data'] component.
+ imageG : {dictionary}
+ imageG is the source image to be resized. it is a dictionary with
+ the data as an ndarray in the ['data'] component.
- reference_image[mat] : {ndarray}
- refernce_image is the image whose sampling dimensions the source
- image is to be remapped to. [mat] refers to the component
- of the image dictionary, reference_image['mat'] that is the
- sampling dimensions.
+ reference_image[mat] : {ndarray}
+ refernce_image is the image whose sampling dimensions the source
+ image is to be remapped to. [mat] refers to the component
+ of the image dictionary, reference_image['mat'] that is the
+ sampling dimensions.
- Returns
- -------
- zoom_image : {dictionary}
+ Returns
+ -------
+ zoom_image : {dictionary}
- Examples
- --------
+ Examples
+ --------
- >>> import _registration as reg
- >>> measures, imageF_anat, fmri_series = reg.demo_MRI_coregistration()
+ >>> import _registration as reg
+ >>> measures, imageF_anat, fmri_series = reg.demo_MRI_coregistration()
- >>> resampled_fmri = reg.resize_image(fmri_series[10], imageF_anat['mat'])
+ >>> resampled_fmri = reg.resize_image(fmri_series[10], imageF_anat['mat'])
- image 10 in the fmri_series is resampled to imageF_anat coordinates
+ image 10 in the fmri_series is resampled to imageF_anat coordinates
- """
+ """
- Z = NP.zeros(3, dtype=NP.float64);
- # get the zoom
- Z[0] = imageG['mat'][0][0] / imageF_mat[0][0]
- Z[1] = imageG['mat'][1][1] / imageF_mat[1][1]
- Z[2] = imageG['mat'][2][2] / imageF_mat[2][2]
+ Z = NP.zeros(3, dtype=NP.float64);
+ # get the zoom
+ Z[0] = imageG['mat'][0][0] / imageF_mat[0][0]
+ Z[1] = imageG['mat'][1][1] / imageF_mat[1][1]
+ Z[2] = imageG['mat'][2][2] / imageF_mat[2][2]
- # new volume dimensions (rounded)
- D = NP.zeros(3, dtype=NP.int32);
- D[0] = int(float(imageG['dim'][0])*Z[0]+0.5)
- D[1] = int(float(imageG['dim'][1])*Z[1]+0.5)
- D[2] = int(float(imageG['dim'][2])*Z[2]+0.5)
+ # new volume dimensions (rounded)
+ D = NP.zeros(3, dtype=NP.int32);
+ D[0] = int(float(imageG['dim'][0])*Z[0]+0.5)
+ D[1] = int(float(imageG['dim'][1])*Z[1]+0.5)
+ D[2] = int(float(imageG['dim'][2])*Z[2]+0.5)
- M = NP.eye(4, dtype=NP.float64);
- # for the test data, set the xyz voxel sizes for fMRI volume
- M[0][0] = imageG['mat'][0][0]/Z[0]
- M[1][1] = imageG['mat'][1][1]/Z[1]
- M[2][2] = imageG['mat'][2][2]/Z[2]
+ M = NP.eye(4, dtype=NP.float64);
+ # for the test data, set the xyz voxel sizes for fMRI volume
+ M[0][0] = imageG['mat'][0][0]/Z[0]
+ M[1][1] = imageG['mat'][1][1]/Z[1]
+ M[2][2] = imageG['mat'][2][2]/Z[2]
- image = NP.zeros(D[2]*D[1]*D[0], dtype=NP.uint8).reshape(D[2], D[0], D[1])
- mode = 2
- scale = 0
- R.register_volume_resample(imageG['data'], image, Z, scale, mode)
- F = NP.zeros(3, dtype=NP.float64);
- zoom_image = {'data' : image, 'mat' : M, 'dim' : D, 'fwhm' : F}
+ image = NP.zeros(D[2]*D[1]*D[0], dtype=NP.uint8).reshape(D[2], D[0], D[1])
+ mode = 2
+ scale = 0
+ R.register_volume_resample(imageG['data'], image, Z, scale, mode)
+ F = NP.zeros(3, dtype=NP.float64);
+ zoom_image = {'data' : image, 'mat' : M, 'dim' : D, 'fwhm' : F}
- return zoom_image
+ return zoom_image
def remap_image(image, parm_vector, resample='linear'):
- """
- remaped_image = remap_image(image, parm_vector, resample='linear')
+ """
+ remaped_image = remap_image(image, parm_vector, resample='linear')
- rotates and translates image using the 3 angles and 3 translations in the 6-dim
- parm_vector. The mapping is stored in the 4x4 M_inverse matrix from the get_inverse_mapping
- method.
+ rotates and translates image using the 3 angles and 3 translations in the 6-dim
+ parm_vector. The mapping is stored in the 4x4 M_inverse matrix from the get_inverse_mapping
+ method.
- Parameters
- ----------
- image : {dictionary}
- image is the source image to be remapped. it is a dictionary with
- the data as an ndarray in the ['data'] component.
+ Parameters
+ ----------
+ image : {dictionary}
+ image is the source image to be remapped. it is a dictionary with
+ the data as an ndarray in the ['data'] component.
- parm_vector : {ndarray}
- parm_vector is the 6-dimensional vector (3 angles, 3 translations)
- generated from the registration.
+ parm_vector : {ndarray}
+ parm_vector is the 6-dimensional vector (3 angles, 3 translations)
+ generated from the registration.
- resample : {'linear', 'cubic'}, optional
+ resample : {'linear', 'cubic'}, optional
- Returns
- -------
- remaped_image : {dictionary}
+ Returns
+ -------
+ remaped_image : {dictionary}
- Examples
- --------
- image = fmri_series[i]
- x[0:6] = measures[i]['align_rotate'][0:6]
- # overwrite the fMRI volume with the aligned volume
- fmri_series[i] = remap_image(image, x, resample='cubic')
+ Examples
+ --------
+ image = fmri_series[i]
+ x[0:6] = measures[i]['align_rotate'][0:6]
+ # overwrite the fMRI volume with the aligned volume
+ fmri_series[i] = remap_image(image, x, resample='cubic')
- """
+ """
- #
- # remap imageG to coordinates of imageF (creates imageG')
- # use the 6 dim parm_vector (3 angles, 3 translations) to remap
- #
- M_inverse = get_inverse_mappings(parm_vector)
- (layers, rows, cols) = image['data'].shape
- # allocate the zero image
- remaped_image = NP.zeros(layers*rows*cols, dtype=NP.uint8).reshape(layers, rows, cols)
- remaped_image = {'data' : remaped_image, 'mat' : image['mat'],
- 'dim' : image['dim'], 'fwhm' : image['fwhm']}
- imdata = build_structs()
+ #
+ # remap imageG to coordinates of imageF (creates imageG')
+ # use the 6 dim parm_vector (3 angles, 3 translations) to remap
+ #
+ M_inverse = get_inverse_mappings(parm_vector)
+ (layers, rows, cols) = image['data'].shape
+ # allocate the zero image
+ remaped_image = NP.zeros(layers*rows*cols, dtype=NP.uint8).reshape(layers, rows, cols)
+ remaped_image = {'data' : remaped_image, 'mat' : image['mat'],
+ 'dim' : image['dim'], 'fwhm' : image['fwhm']}
+ imdata = build_structs()
- if resample == 'linear':
- # trilinear interpolation mapping.
- R.register_linear_resample(image['data'], remaped_image['data'], M_inverse, imdata['step'])
- elif resample == 'cubic':
- # tricubic convolve interpolation mapping.
- R.register_cubic_resample(image['data'], remaped_image['data'], M_inverse, imdata['step'])
+ if resample == 'linear':
+ # trilinear interpolation mapping.
+ R.register_linear_resample(image['data'], remaped_image['data'], M_inverse, imdata['step'])
+ elif resample == 'cubic':
+ # tricubic convolve interpolation mapping.
+ R.register_cubic_resample(image['data'], remaped_image['data'], M_inverse, imdata['step'])
- return remaped_image
+ return remaped_image
def get_inverse_mappings(parm_vector):
- """
- M_inverse = get_inverse_mappings(parm_vector)
+ """
+ M_inverse = get_inverse_mappings(parm_vector)
- takes the 6-dimensional rotation/translation vector and builds the inverse
- 4x4 mapping matrix M_inverse that will map imageG to imageF orientation
+ takes the 6-dimensional rotation/translation vector and builds the inverse
+ 4x4 mapping matrix M_inverse that will map imageG to imageF orientation
- Parameters
- ----------
- parm_vector : {nd_array}
+ Parameters
+ ----------
+ parm_vector : {nd_array}
- Returns
- -------
- M_inverse : {nd_array}
+ Returns
+ -------
+ M_inverse : {nd_array}
- Examples
- --------
+ Examples
+ --------
- >>> import numpy as NP
- >>> import _registration as reg
- >>> array = NP.zeros(6, dtype=float)
- >>> M = reg.get_inverse_mappings(array)
- >>> M
+ >>> import numpy as NP
+ >>> import _registration as reg
+ >>> array = NP.zeros(6, dtype=float)
+ >>> M = reg.get_inverse_mappings(array)
+ >>> M
- array([
- [ 1., 0., 0., 0.],
- [ 0., 1., 0., 0.],
- [ 0., 0., 1., 0.],
- [ 0., 0., 0., 1.]])
+ array([
+ [ 1., 0., 0., 0.],
+ [ 0., 1., 0., 0.],
+ [ 0., 0., 1., 0.],
+ [ 0., 0., 0., 1.]])
- """
- # get the inverse mapping to rotate the G matrix to F space following registration
- imdata = build_structs()
- # inverse angles and translations
- imdata['parms'][0] = -parm_vector[0]
- imdata['parms'][1] = -parm_vector[1]
- imdata['parms'][2] = -parm_vector[2]
- imdata['parms'][3] = -parm_vector[3]
- imdata['parms'][4] = -parm_vector[4]
- imdata['parms'][5] = -parm_vector[5]
- M_inverse = build_rotate_matrix(imdata['parms'])
- return M_inverse
+ """
+ # get the inverse mapping to rotate the G matrix to F space following registration
+ imdata = build_structs()
+ # inverse angles and translations
+ imdata['parms'][0] = -parm_vector[0]
+ imdata['parms'][1] = -parm_vector[1]
+ imdata['parms'][2] = -parm_vector[2]
+ imdata['parms'][3] = -parm_vector[3]
+ imdata['parms'][4] = -parm_vector[4]
+ imdata['parms'][5] = -parm_vector[5]
+ M_inverse = build_rotate_matrix(imdata['parms'])
+ return M_inverse
def python_coreg(image1, image2, imdata, ftype=1, smimage=0, lite=0, smhist=0,
- method='nmi', opt_method='powell'):
- """
- parm_vector = python_coreg(image1, image2, imdata, ftype=1, smimage=0, lite=0,
- smhist=0, method='nmi', opt_method='powell'):
+ method='nmi', opt_method='powell'):
+ """
+ parm_vector = python_coreg(image1, image2, imdata, ftype=1, smimage=0, lite=0,
+ smhist=0, method='nmi', opt_method='powell'):
- takes two images and the image data descriptor (imdata) and determines the optimal
- alignment of the two images (measured by mutual information or cross correlation)
- using optimization search of 3 angle and 3 translation parameters. The optimization
- uses either the Powell or Conjugate Gradient methods in the scipy optimization
- package. The optimal parameter is returned.
+ takes two images and the image data descriptor (imdata) and determines the optimal
+ alignment of the two images (measured by mutual information or cross correlation)
+ using optimization search of 3 angle and 3 translation parameters. The optimization
+ uses either the Powell or Conjugate Gradient methods in the scipy optimization
+ package. The optimal parameter is returned.
- Parameters
- ----------
- image1 : {dictionary}
- image1 is the source image to be remapped during the registration.
- it is a dictionary with the data as an ndarray in the ['data'] component.
- image2 : {dictionary}
- image2 is the reference image that image1 gets mapped to.
- imdata : {dictionary}
- image sampling and optimization information.
- ftype : {0, 1}, optional
- flag for type of low pass filter. 0 is Gauss-Spline
- 1 is pure Gauss. Sigma determined from volume sampling info.
- smimage : {0, 1}, optional
- flag for volume 3D low pass filtering of image 2.
- 0 for no filter, 1 for do filter.
- lite : {0, 1}, optional
- lite of 1 is to jitter both images during resampling. 0
- is to not jitter. jittering is for non-aliased volumes.
- smhist: {0, 1}, optional
- flag for joint histogram low pass filtering. 0 for no filter,
- 1 for do filter.
- method: {'nmi', 'mi', 'ncc', 'ecc'}, optional
- flag for type of registration metric. nmi is normalized mutual
- information; mi is mutual information; ecc is entropy cross
- correlation; ncc is normalized cross correlation.
- opt_method: {'powell', 'hybrid'}, optional
- registration is two pass. Pass 1 is low res to get close to alignment
- and pass 2 starts at the pass 1 optimal alignment. In powell pass 1 and
- 2 are powell, in hybrid pass 2 is conjugate gradient.
+ Parameters
+ ----------
+ image1 : {dictionary}
+ image1 is the source image to be remapped during the registration.
+ it is a dictionary with the data as an ndarray in the ['data'] component.
+ image2 : {dictionary}
+ image2 is the reference image that image1 gets mapped to.
+ imdata : {dictionary}
+ image sampling and optimization information.
+ ftype : {0, 1}, optional
+ flag for type of low pass filter. 0 is Gauss-Spline
+ 1 is pure Gauss. Sigma determined from volume sampling info.
+ smimage : {0, 1}, optional
+ flag for volume 3D low pass filtering of image 2.
+ 0 for no filter, 1 for do filter.
+ lite : {0, 1}, optional
+ lite of 1 is to jitter both images during resampling. 0
+ is to not jitter. jittering is for non-aliased volumes.
+ smhist: {0, 1}, optional
+ flag for joint histogram low pass filtering. 0 for no filter,
+ 1 for do filter.
+ method: {'nmi', 'mi', 'ncc', 'ecc'}, optional
+ flag for type of registration metric. nmi is normalized mutual
+ information; mi is mutual information; ecc is entropy cross
+ correlation; ncc is normalized cross correlation.
+ opt_method: {'powell', 'hybrid'}, optional
+ registration is two pass. Pass 1 is low res to get close to alignment
+ and pass 2 starts at the pass 1 optimal alignment. In powell pass 1 and
+ 2 are powell, in hybrid pass 2 is conjugate gradient.
- Returns
- -------
- parm_vector : {nd_array}
- this is the optimal alignment (6-dim) array with 3 angles and
- 3 translations.
+ Returns
+ -------
+ parm_vector : {nd_array}
+ this is the optimal alignment (6-dim) array with 3 angles and
+ 3 translations.
- Examples
- --------
+ Examples
+ --------
- >>> import numpy as NP
- >>> import _registration as reg
+ >>> import numpy as NP
+ >>> import _registration as reg
- >>> image1, image2, imdata = reg.demo_MRI_volume_align()
- >>> parm_vector = python_coreg(image1, image2, imdata)
+ >>> image1, image2, imdata = reg.demo_MRI_volume_align()
+ >>> parm_vector = python_coreg(image1, image2, imdata)
- """
- start = time.time()
- # smooth of the images
- if smimage:
- image_F_xyz2 = filter_image_3D(image2['data'], image2['fwhm'], ftype)
- image2['data'] = image_F_xyz2
- parm_vector = multires_registration(image1, image2, imdata, lite, smhist, method, opt_method)
- stop = time.time()
- print 'Total Optimizer Time is ', (stop-start)
- return parm_vector
+ """
+ start = time.time()
+ # smooth of the images
+ if smimage:
+ image_F_xyz2 = filter_image_3D(image2['data'], image2['fwhm'], ftype)
+ image2['data'] = image_F_xyz2
+ parm_vector = multires_registration(image1, image2, imdata, lite, smhist, method, opt_method)
+ stop = time.time()
+ print 'Total Optimizer Time is ', (stop-start)
+ return parm_vector
def multires_registration(image1, image2, imdata, lite, smhist, method, opt_method):
- """
- x = multires_registration(image1, image2, imdata, lite, smhist, method, opt_method)
+ """
+ x = multires_registration(image1, image2, imdata, lite, smhist, method, opt_method)
- to be called by python_coreg() which optionally does 3D image filtering and
- provies timing for registration.
+ to be called by python_coreg() which optionally does 3D image filtering and
+ provies timing for registration.
- Parameters
- ----------
+ Parameters
+ ----------
- image1 : {dictionary}
- image1 is the source image to be remapped during the registration.
- it is a dictionary with the data as an ndarray in the ['data'] component.
- image2 : {dictionary}
- image2 is the reference image that image1 gets mapped to.
- imdata : {dictionary}
- image sampling and optimization information.
- lite : {integer}
- lite of 1 is to jitter both images during resampling. 0
- is to not jitter. jittering is for non-aliased volumes.
- smhist: {integer}
- flag for joint histogram low pass filtering. 0 for no filter,
- 1 for do filter.
- method: {'nmi', 'mi', 'ncc', 'ecc'}
- flag for type of registration metric. nmi is normalized mutual
- information; mi is mutual information; ecc is entropy cross
- correlation; ncc is normalized cross correlation.
- opt_method: {'powell', 'hybrid'}
- registration is two pass. Pass 1 is low res to get close to alignment
- and pass 2 starts at the pass 1 optimal alignment. In powell pass 1 and
- 2 are powell, in hybrid pass 2 is conjugate gradient.
+ image1 : {dictionary}
+ image1 is the source image to be remapped during the registration.
+ it is a dictionary with the data as an ndarray in the ['data'] component.
+ image2 : {dictionary}
+ image2 is the reference image that image1 gets mapped to.
+ imdata : {dictionary}
+ image sampling and optimization information.
+ lite : {integer}
+ lite of 1 is to jitter both images during resampling. 0
+ is to not jitter. jittering is for non-aliased volumes.
+ smhist: {integer}
+ flag for joint histogram low pass filtering. 0 for no filter,
+ 1 for do filter.
+ method: {'nmi', 'mi', 'ncc', 'ecc'}
+ flag for type of registration metric. nmi is normalized mutual
+ information; mi is mutual information; ecc is entropy cross
+ correlation; ncc is normalized cross correlation.
+ opt_method: {'powell', 'hybrid'}
+ registration is two pass. Pass 1 is low res to get close to alignment
+ and pass 2 starts at the pass 1 optimal alignment. In powell pass 1 and
+ 2 are powell, in hybrid pass 2 is conjugate gradient.
- Returns
- -------
- x : {nd_array}
- this is the optimal alignment (6-dim) array with 3 angles and
- 3 translations.
+ Returns
+ -------
+ x : {nd_array}
+ this is the optimal alignment (6-dim) array with 3 angles and
+ 3 translations.
- Examples
- --------
+ Examples
+ --------
- (calling this from python_coreg which optionally filters image2)
- >>> import numpy as NP
- >>> import _registration as reg
- >>> image1, image2, imdata = reg.demo_MRI_volume_align()
- >>> parm_vector = python_coreg(image1, image2, imdata)
+ (calling this from python_coreg which optionally filters image2)
+ >>> import numpy as NP
+ >>> import _registration as reg
+ >>> image1, image2, imdata = reg.demo_MRI_volume_align()
+ >>> parm_vector = python_coreg(image1, image2, imdata)
- """
- ret_histo=0
- # zero out the start parameter; but this may be set to large values
- # if the head is out of range and well off the optimal alignment skirt
- imdata['parms'][0:5] = 0.0
- # make the step a scalar to can put in a multi-res loop
- loop = range(imdata['sample'].size)
- x = imdata['parms']
- for i in loop:
- step = imdata['sample'][i]
- imdata['step'][:] = step
- optfunc_args = (image1, image2, imdata['step'], imdata['fwhm'], lite, smhist,
- method, ret_histo)
- p_args = (optfunc_args,)
- if opt_method=='powell':
- print 'POWELL multi-res registration step size ', step
- print 'vector ', x
- x = OPT.fmin_powell(optimize_function, x, args=p_args,
- callback=callback_powell)
- elif opt_method=='cg':
- print 'CG multi-res registration step size ', step
- print 'vector ', x
- x = OPT.fmin_cg(optimize_function, x, args=p_args, callback=callback_cg)
- elif opt_method=='hybrid':
- if i==0:
- print 'Hybrid POWELL multi-res registration step size ', step
- print 'vector ', x
- lite = 0
- optfunc_args = (image1, image2, imdata['step'], imdata['fwhm'], lite, smhist,
- method, ret_histo)
- p_args = (optfunc_args,)
- x = OPT.fmin_powell(optimize_function, x, args=p_args, callback=callback_powell)
- elif i==1:
- print 'Hybrid CG multi-res registration step size ', step
- print 'vector ', x
- lite = 1
- optfunc_args = (image1, image2, imdata['step'], imdata['fwhm'], lite,
- smhist, method, ret_histo)
- p_args = (optfunc_args,)
- x = OPT.fmin_cg(optimize_function, x, args=p_args, callback=callback_cg)
+ """
+ ret_histo=0
+ # zero out the start parameter; but this may be set to large values
+ # if the head is out of range and well off the optimal alignment skirt
+ imdata['parms'][0:5] = 0.0
+ # make the step a scalar to can put in a multi-res loop
+ loop = range(imdata['sample'].size)
+ x = imdata['parms']
+ for i in loop:
+ step = imdata['sample'][i]
+ imdata['step'][:] = step
+ optfunc_args = (image1, image2, imdata['step'], imdata['fwhm'], lite, smhist,
+ method, ret_histo)
+ p_args = (optfunc_args,)
+ if opt_method=='powell':
+ print 'POWELL multi-res registration step size ', step
+ print 'vector ', x
+ x = OPT.fmin_powell(optimize_function, x, args=p_args,
+ callback=callback_powell)
+ elif opt_method=='cg':
+ print 'CG multi-res registration step size ', step
+ print 'vector ', x
+ x = OPT.fmin_cg(optimize_function, x, args=p_args, callback=callback_cg)
+ elif opt_method=='hybrid':
+ if i==0:
+ print 'Hybrid POWELL multi-res registration step size ', step
+ print 'vector ', x
+ lite = 0
+ optfunc_args = (image1, image2, imdata['step'], imdata['fwhm'], lite, smhist,
+ method, ret_histo)
+ p_args = (optfunc_args,)
+ x = OPT.fmin_powell(optimize_function, x, args=p_args, callback=callback_powell)
+ elif i==1:
+ print 'Hybrid CG multi-res registration step size ', step
+ print 'vector ', x
+ lite = 1
+ optfunc_args = (image1, image2, imdata['step'], imdata['fwhm'], lite,
+ smhist, method, ret_histo)
+ p_args = (optfunc_args,)
+ x = OPT.fmin_cg(optimize_function, x, args=p_args, callback=callback_cg)
- return x
+ return x
def callback_powell(x):
- """
- called by optimize.powell only. prints current parameter vector.
- """
- print 'Parameter Vector from Powell: - '
- print x
- return
+ """
+ called by optimize.powell only. prints current parameter vector.
+ """
+ print 'Parameter Vector from Powell: - '
+ print x
+ return
def callback_cg(x):
- """
- called by optimize.cg only. prints current parameter vector.
- """
- print 'Parameter Vector from Conjugate Gradient: - '
- print x
- return
+ """
+ called by optimize.cg only. prints current parameter vector.
+ """
+ print 'Parameter Vector from Conjugate Gradient: - '
+ print x
+ return
def smooth_kernel(fwhm, x, ktype=1):
- """
- kernel = smooth_kernel(fwhm, x, ktype=1)
+ """
+ kernel = smooth_kernel(fwhm, x, ktype=1)
- smooth_kernel creates filter kernel based on image sampling parameters.
- provide domain of kernel and sampling parameters.
+ smooth_kernel creates filter kernel based on image sampling parameters.
+ provide domain of kernel and sampling parameters.
- Parameters
- ----------
- fwhm : {int}
- used for kernel width
- x : {nd_array}
- domain of kernel
- ktype: {1, 2}, optional
- kernel type. 1 is Gauss convoled with spline, 2 is Gauss
+ Parameters
+ ----------
+ fwhm : {int}
+ used for kernel width
+ x : {nd_array}
+ domain of kernel
+ ktype: {1, 2}, optional
+ kernel type. 1 is Gauss convoled with spline, 2 is Gauss
- Returns
- -------
- kernel : {nd_array}
+ Returns
+ -------
+ kernel : {nd_array}
- Examples
- --------
+ Examples
+ --------
- >>> import numpy as NP
- >>> import _registration as reg
- >>> fwhm = 3
- >>> ftype = 2
- >>> p = NP.ceil(2*fwhm).astype(int)
- >>> x = NP.array(range(-p, p+1))
- >>> kernel = reg.smooth_kernel(fwhm, x, ktype=ftype)
- >>> kernel
+ >>> import numpy as NP
+ >>> import _registration as reg
+ >>> fwhm = 3
+ >>> ftype = 2
+ >>> p = NP.ceil(2*fwhm).astype(int)
+ >>> x = NP.array(range(-p, p+1))
+ >>> kernel = reg.smooth_kernel(fwhm, x, ktype=ftype)
+ >>> kernel
- array([
- 4.77853772e-06, 1.41575516e-04, 2.26516955e-03,
- 1.95718875e-02, 9.13238336e-02, 2.30120330e-01,
- 3.13144850e-01, 2.30120330e-01, 9.13238336e-02,
- 1.95718875e-02, 2.26516955e-03, 1.41575516e-04,
- 4.77853772e-06])
+ array([
+ 4.77853772e-06, 1.41575516e-04, 2.26516955e-03,
+ 1.95718875e-02, 9.13238336e-02, 2.30120330e-01,
+ 3.13144850e-01, 2.30120330e-01, 9.13238336e-02,
+ 1.95718875e-02, 2.26516955e-03, 1.41575516e-04,
+ 4.77853772e-06])
- """
- eps = 0.00001
- s = NP.square((fwhm/math.sqrt(8.0*math.log(2.0)))) + eps
- if ktype==1:
- # from SPM: Gauss kernel convolved with 1st degree B spline
- w1 = 0.5 * math.sqrt(2.0/s)
- w2 = -0.5 / s
- w3 = math.sqrt((s*math.pi) /2.0)
- kernel = 0.5*(SP.erf(w1*(x+1))*(x+1) + SP.erf(w1*(x-1))*(x-1) - 2.0*SP.erf(w1*(x))*(x) +
- w3*(NP.exp(w2*NP.square(x+1))) + NP.exp(w2*(NP.square(x-1))) - 2.0*NP.exp(w2*NP.square(x)))
- kernel[kernel<0] = 0
- kernel = kernel / kernel.sum()
- else:
- # Gauss kernel
- kernel = (1.0/math.sqrt(2.0*math.pi*s)) * NP.exp(-NP.square(x)/(2.0*s))
- kernel = kernel / kernel.sum()
+ """
+ eps = 0.00001
+ s = NP.square((fwhm/math.sqrt(8.0*math.log(2.0)))) + eps
+ if ktype==1:
+ # from SPM: Gauss kernel convolved with 1st degree B spline
+ w1 = 0.5 * math.sqrt(2.0/s)
+ w2 = -0.5 / s
+ w3 = math.sqrt((s*math.pi) /2.0)
+ kernel = 0.5*(SP.erf(w1*(x+1))*(x+1) + SP.erf(w1*(x-1))*(x-1) - 2.0*SP.erf(w1*(x))*(x) +
+ w3*(NP.exp(w2*NP.square(x+1))) + NP.exp(w2*(NP.square(x-1))) - 2.0*NP.exp(w2*NP.square(x)))
+ kernel[kernel<0] = 0
+ kernel = kernel / kernel.sum()
+ else:
+ # Gauss kernel
+ kernel = (1.0/math.sqrt(2.0*math.pi*s)) * NP.exp(-NP.square(x)/(2.0*s))
+ kernel = kernel / kernel.sum()
- return kernel
+ return kernel
def filter_image_3D(imageRaw, fwhm, ftype=2):
- """
- image_F_xyz = filter_image_3D(imageRaw, fwhm, ftype=2):
- does 3D separable digital filtering using scipy.ndimage.correlate1d
+ """
+ image_F_xyz = filter_image_3D(imageRaw, fwhm, ftype=2):
+ does 3D separable digital filtering using scipy.ndimage.correlate1d
- Parameters
- ----------
- imageRaw : {nd_array}
- the unfiltered 3D volume image
- fwhm : {int}
- used for kernel width
- ktype: {1, 2}, optional
- kernel type. 1 is Gauss convoled with spline, 2 is Gauss
+ Parameters
+ ----------
+ imageRaw : {nd_array}
+ the unfiltered 3D volume image
+ fwhm : {int}
+ used for kernel width
+ ktype: {1, 2}, optional
+ kernel type. 1 is Gauss convoled with spline, 2 is Gauss
- Returns
- -------
- image_F_xyz : {nd_array}
- 3D filtered volume image
+ Returns
+ -------
+ image_F_xyz : {nd_array}
+ 3D filtered volume image
- Examples
- --------
+ Examples
+ --------
- >>> import _registration as reg
- >>> image1, image2, imdata = reg.demo_MRI_volume_align()
- >>> ftype = 1
- >>> image_Filter_xyz = filter_image_3D(image1['data'], image1['fwhm'], ftype)
- >>> image1['data'] = image_Filter_xyz
- """
+ >>> import _registration as reg
+ >>> image1, image2, imdata = reg.demo_MRI_volume_align()
+ >>> ftype = 1
+ >>> image_Filter_xyz = filter_image_3D(image1['data'], image1['fwhm'], ftype)
+ >>> image1['data'] = image_Filter_xyz
+ """
- p = NP.ceil(2*fwhm[0]).astype(int)
- x = NP.array(range(-p, p+1))
- kernel_x = smooth_kernel(fwhm[0], x, ktype=ftype)
- p = NP.ceil(2*fwhm[1]).astype(int)
- x = NP.array(range(-p, p+1))
- kernel_y = smooth_kernel(fwhm[1], x, ktype=ftype)
- p = NP.ceil(2*fwhm[2]).astype(int)
- x = NP.array(range(-p, p+1))
- kernel_z = smooth_kernel(fwhm[2], x, ktype=ftype)
- output=None
- # 3D filter in 3 1D separable stages
- axis = 0
- image_F_x = NDI.correlate1d(imageRaw, kernel_x, axis, output)
- axis = 1
- image_F_xy = NDI.correlate1d(image_F_x, kernel_y, axis, output)
- axis = 2
- image_F_xyz = NDI.correlate1d(image_F_xy, kernel_z, axis, output)
- return image_F_xyz
+ p = NP.ceil(2*fwhm[0]).astype(int)
+ x = NP.array(range(-p, p+1))
+ kernel_x = smooth_kernel(fwhm[0], x, ktype=ftype)
+ p = NP.ceil(2*fwhm[1]).astype(int)
+ x = NP.array(range(-p, p+1))
+ kernel_y = smooth_kernel(fwhm[1], x, ktype=ftype)
+ p = NP.ceil(2*fwhm[2]).astype(int)
+ x = NP.array(range(-p, p+1))
+ kernel_z = smooth_kernel(fwhm[2], x, ktype=ftype)
+ output=None
+ # 3D filter in 3 1D separable stages
+ axis = 0
+ image_F_x = NDI.correlate1d(imageRaw, kernel_x, axis, output)
+ axis = 1
+ image_F_xy = NDI.correlate1d(image_F_x, kernel_y, axis, output)
+ axis = 2
+ image_F_xyz = NDI.correlate1d(image_F_xy, kernel_z, axis, output)
+ return image_F_xyz
def build_fwhm(M, S):
- """
- fwhm = build_fwhm(M, S)
+ """
+ fwhm = build_fwhm(M, S)
- builds the low pass filter kernel sigma value from the image pixel sampling
+ builds the low pass filter kernel sigma value from the image pixel sampling
- Parameters
- ----------
- M : {nd_array}
- input 4x4 voxel to physical map matrix (called 'MAT')
+ Parameters
+ ----------
+ M : {nd_array}
+ input 4x4 voxel to physical map matrix (called 'MAT')
- S : {nd_array}
- 1x3 sample increment array. should be = (1, 1, 1)
+ S : {nd_array}
+ 1x3 sample increment array. should be = (1, 1, 1)
- Returns
- -------
- fwhm : {nd_array}
- the 3D Gaussian kernel width
+ Returns
+ -------
+ fwhm : {nd_array}
+ the 3D Gaussian kernel width
- Examples
- --------
+ Examples
+ --------
- >>> import numpy as NP
- >>> import _registration as reg
- >>> anat_desc = reg.load_anatMRI_desc()
- >>> image1 = reg.load_volume(anat_desc, imagename='ANAT1_V0001.img')
- >>> imdata = reg.build_structs()
- >>> image1['fwhm'] = reg.build_fwhm(image1['mat'], imdata['step'])
+ >>> import numpy as NP
+ >>> import _registration as reg
+ >>> anat_desc = reg.load_anatMRI_desc()
+ >>> image1 = reg.load_volume(anat_desc, imagename='ANAT1_V0001.img')
+ >>> imdata = reg.build_structs()
+ >>> image1['fwhm'] = reg.build_fwhm(image1['mat'], imdata['step'])
- """
- view_3x3 = NP.square(M[0:3, 0:3])
- # sum the elements inn the first row
- vxg = NP.sqrt(view_3x3.sum(axis=0))
- # assumes that sampling is the same for xyz
- size = NP.array([1,1,1])*S[0]
- x = NP.square(size) - NP.square(vxg)
- # clip
- x[x<0] = 0
- fwhm = NP.sqrt(x) / vxg
- # pathology when stepsize = 1 for MAT equal to the identity matrix
- fwhm[fwhm==0] = 1
- # return the 3D Gaussian kernel width (xyz)
- return fwhm
+ """
+ view_3x3 = NP.square(M[0:3, 0:3])
+ # sum the elements inn the first row
+ vxg = NP.sqrt(view_3x3.sum(axis=0))
+ # assumes that sampling is the same for xyz
+ size = NP.array([1,1,1])*S[0]
+ x = NP.square(size) - NP.square(vxg)
+ # clip
+ x[x<0] = 0
+ fwhm = NP.sqrt(x) / vxg
+ # pathology when stepsize = 1 for MAT equal to the identity matrix
+ fwhm[fwhm==0] = 1
+ # return the 3D Gaussian kernel width (xyz)
+ return fwhm
def optimize_function(x, optfunc_args):
- """
- cost = optimize_function(x, optfunc_args) --- OR ---
- cost, joint_histogram = optimize_function(x, optfunc_args)
+ """
+ cost = optimize_function(x, optfunc_args) --- OR ---
+ cost, joint_histogram = optimize_function(x, optfunc_args)
- computes the alignment between 2 volumes using cross correlation or mutual
- information metrics. In both the 8 bit joint histogram of the 2 images is
- computed. The 8 bit scaling is done using an integrated histogram method and
- is called prior to this.
+ computes the alignment between 2 volumes using cross correlation or mutual
+ information metrics. In both the 8 bit joint histogram of the 2 images is
+ computed. The 8 bit scaling is done using an integrated histogram method and
+ is called prior to this.
- Parameters
- ----------
- x : {nd_array}
- this is the current (6-dim) array with 3 angles and 3 translations.
+ Parameters
+ ----------
+ x : {nd_array}
+ this is the current (6-dim) array with 3 angles and 3 translations.
- optfunc_args : {tuple}
- this is a tuple of 8 elements that is formed in the scipy.optimize powell
- and cg (conjugate gradient) functions. this is why the elements need to be
- a tuple. The elements of optfunc_args are:
+ optfunc_args : {tuple}
+ this is a tuple of 8 elements that is formed in the scipy.optimize powell
+ and cg (conjugate gradient) functions. this is why the elements need to be
+ a tuple. The elements of optfunc_args are:
- image_F : {dictionary}
- image_F is the source image to be remapped during the registration.
- it is a dictionary with the data as an ndarray in the ['data'] component.
- image_G : {dictionary}
- image_G is the reference image that image_F gets mapped to.
- sample_vector : {nd_array}
- sample in x,y,x. should be (1,1,1)
- fwhm : {nd_array}
- Gaussian sigma
- do_lite : {0, 1}
- lite of 1 is to jitter both images during resampling.
- 0 is to not jitter. jittering is for non-aliased volumes.
- smooth : {0, 1}
- flag for joint histogram low pass filtering. 0 for no filter,
- 1 for do filter.
- method : {'nmi', 'mi', 'ncc', 'ecc', 'mse'}
- flag for type of registration metric. nmi is normalized mutual
- information; mi is mutual information; ecc is entropy cross
- correlation; ncc is normalized cross correlation. mse is mean
- square error. with mse there is no joint histogram.
- ret_histo : {0, 1}
- if 0 return is: cost
- if 0 return is: cost, joint_histogram
+ image_F : {dictionary}
+ image_F is the source image to be remapped during the registration.
+ it is a dictionary with the data as an ndarray in the ['data'] component.
+ image_G : {dictionary}
+ image_G is the reference image that image_F gets mapped to.
+ sample_vector : {nd_array}
+ sample in x,y,x. should be (1,1,1)
+ fwhm : {nd_array}
+ Gaussian sigma
+ do_lite : {0, 1}
+ lite of 1 is to jitter both images during resampling.
+ 0 is to not jitter. jittering is for non-aliased volumes.
+ smooth : {0, 1}
+ flag for joint histogram low pass filtering. 0 for no filter,
+ 1 for do filter.
+ method : {'nmi', 'mi', 'ncc', 'ecc', 'mse'}
+ flag for type of registration metric. nmi is normalized mutual
+ information; mi is mutual information; ecc is entropy cross
+ correlation; ncc is normalized cross correlation. mse is mean
+ square error. with mse there is no joint histogram.
+ ret_histo : {0, 1}
+ if 0 return is: cost
+ if 0 return is: cost, joint_histogram
- Returns
- -------
- cost : {float}
- the negative of one of the mutual information metrics
- or negative cross correlation. use negative as the optimization
- is minimization.
+ Returns
+ -------
+ cost : {float}
+ the negative of one of the mutual information metrics
+ or negative cross correlation. use negative as the optimization
+ is minimization.
- --- OR --- (if ret_histo = 1)
+ --- OR --- (if ret_histo = 1)
- cost : {float}
- the negative of one of the mutual information metrics
- or negative cross correlation. use negative as the optimization
- is minimization.
+ cost : {float}
+ the negative of one of the mutual information metrics
+ or negative cross correlation. use negative as the optimization
+ is minimization.
- joint_histogram : {nd_array}
- the 2D (256x256) joint histogram of the two volumes
+ joint_histogram : {nd_array}
+ the 2D (256x256) joint histogram of the two volumes
- Examples
- --------
+ Examples
+ --------
- >>> import numpy as NP
- >>> import _registration as reg
- >>> anat_desc = reg.load_anatMRI_desc()
- >>> image1 = reg.load_volume(anat_desc, imagename='ANAT1_V0001.img')
- >>> image2 = reg.load_volume(anat_desc, imagename='ANAT1_V0001.img')
- >>> imdata = reg.build_structs()
- >>> image1['fwhm'] = reg.build_fwhm(image1['mat'], imdata['step'])
- >>> image2['fwhm'] = reg.build_fwhm(image2['mat'], imdata['step'])
- >>> method = 'ncc'
- >>> lite = 1
- >>> smhist = 0
- >>> ret_histo = 1
- >>> optfunc_args = (image1, image2, imdata['step'], imdata['fwhm'], lite, smhist, method, ret_histo)
- >>> x = NP.zeros(6, dtype=NP.float64)
- >>> return cost, joint_histogram = reg.optimize_function(x, optfunc_args)
+ >>> import numpy as NP
+ >>> import _registration as reg
+ >>> anat_desc = reg.load_anatMRI_desc()
+ >>> image1 = reg.load_volume(anat_desc, imagename='ANAT1_V0001.img')
+ >>> image2 = reg.load_volume(anat_desc, imagename='ANAT1_V0001.img')
+ >>> imdata = reg.build_structs()
+ >>> image1['fwhm'] = reg.build_fwhm(image1['mat'], imdata['step'])
+ >>> image2['fwhm'] = reg.build_fwhm(image2['mat'], imdata['step'])
+ >>> method = 'ncc'
+ >>> lite = 1
+ >>> smhist = 0
+ >>> ret_histo = 1
+ >>> optfunc_args = (image1, image2, imdata['step'], imdata['fwhm'], lite, smhist, method, ret_histo)
+ >>> x = NP.zeros(6, dtype=NP.float64)
+ >>> return cost, joint_histogram = reg.optimize_function(x, optfunc_args)
- """
+ """
- image_F = optfunc_args[0]
- image_G = optfunc_args[1]
- sample_vector = optfunc_args[2]
- fwhm = optfunc_args[3]
- do_lite = optfunc_args[4]
- smooth = optfunc_args[5]
- method = optfunc_args[6]
- ret_histo = optfunc_args[7]
+ image_F = optfunc_args[0]
+ image_G = optfunc_args[1]
+ sample_vector = optfunc_args[2]
+ fwhm = optfunc_args[3]
+ do_lite = optfunc_args[4]
+ smooth = optfunc_args[5]
+ method = optfunc_args[6]
+ ret_histo = optfunc_args[7]
- rot_matrix = build_rotate_matrix(x)
- cost = 0.0
- epsilon = 2.2e-16
- # image_G is base image
- # image_F is the to-be-rotated image
- # rot_matrix is the 4x4 constructed (current angles and translates) transform matrix
- # sample_vector is the subsample vector for x-y-z
+ rot_matrix = build_rotate_matrix(x)
+ cost = 0.0
+ epsilon = 2.2e-16
+ # image_G is base image
+ # image_F is the to-be-rotated image
+ # rot_matrix is the 4x4 constructed (current angles and translates) transform matrix
+ # sample_vector is the subsample vector for x-y-z
- F_inv = NP.linalg.inv(image_F['mat'])
- composite = NP.dot(F_inv, image_G['mat'])
- composite = NP.dot(composite, rot_matrix)
+ F_inv = NP.linalg.inv(image_F['mat'])
+ composite = NP.dot(F_inv, image_G['mat'])
+ composite = NP.dot(composite, rot_matrix)
- if method == 'mse':
- #
- #