[Scipysvn] r5216  in trunk/doc: . source/tutorial
scipysvn@scip...
scipysvn@scip...
Wed Dec 3 17:24:57 CST 2008
Author: ptvirtan
Date: 20081203 17:24:44 0600 (Wed, 03 Dec 2008)
New Revision: 5216
Modified:
trunk/doc/Makefile
trunk/doc/source/tutorial/index.rst
Log:
tutorial: remove, update or comment out parts that are outofdate. Remove the updating warning
Modified: trunk/doc/Makefile
===================================================================
 trunk/doc/Makefile 20081203 23:23:53 UTC (rev 5215)
+++ trunk/doc/Makefile 20081203 23:24:44 UTC (rev 5216)
@@ 40,7 +40,7 @@
mkdir p build/dist
cp r build/html build/dist/reference
touch build/dist/index.html
 perl pi e 's#^\s*(<li><a href=".*?">SciPy.*?Reference Guide.*?»</li>)$$#<li><a href="/">Numpy and Scipy Documentation</a> »</li> $$1#;' build/dist/*.html build/dist/*/*.html build/dist/*/*/*.html
+ perl pi e 's#^\s*(<li><a href=".*?">SciPy.*?Reference Guide.*?»</li>)\s*$$#<li><a href="/">Numpy and Scipy Documentation</a> »</li> $$1#;' build/dist/*.html build/dist/*/*.html build/dist/*/*/*.html
(cd build/html && zip 9qr ../dist/scipyhtml.zip .)
cp build/latex/*.pdf build/dist
cd build/dist && tar czf ../dist.tar.gz *
Modified: trunk/doc/source/tutorial/index.rst
===================================================================
 trunk/doc/source/tutorial/index.rst 20081203 23:23:53 UTC (rev 5215)
+++ trunk/doc/source/tutorial/index.rst 20081203 23:24:44 UTC (rev 5216)
@@ 6,12 +6,7 @@
.. currentmodule:: scipy
.. warning::
 This is adapted from the old Scipy tutorial, and it is being revised.
 Parts of it are still outofdate.


Introduction
============
@@ 52,8 +47,13 @@
General Help

Python provides the facility of documentation strings. The functions
and classes available in SciPy use this method for online
+Scipy and Numpy have HTML and PDF versions of their documentation
+available at http://docs.scipy.org/, which currently details nearly
+all available functionality. However, this documentation is still
+workinprogress, and some parts may be incomplete or sparse.
+
+Python also provides the facility of documentation strings. The
+functions and classes available in SciPy use this method for online
documentation. There are two methods for reading these messages and
getting help. Python provides the command :func:`help` in the pydoc
module. Entering this command with no arguments (i.e. ``>>> help`` )
@@ 68,9 +68,10 @@
is also available under the command scipy.info. The signature and
documentation string for the object passed to the help command are
printed to standard output (or to a writeable object passed as the
third argument). The second keyword argument of "scipy.info "defines the maximum width of the line for printing. If a module is
passed as the argument to help than a list of the functions and
classes defined in that module is printed. For example:
+third argument). The second keyword argument of "scipy.info" defines
+the maximum width of the line for printing. If a module is passed as
+the argument to help than a list of the functions and classes defined
+in that module is printed. For example:
.. literalinclude:: examples/11
@@ 86,113 +87,91 @@

SciPy is organized into subpackages covering different scientific
computing domains. Many common functions which several subpackages
rely are in :mod:`numpy`. This allows for the possibility of
separately distributing the subpackages of scipy as long as Numpy is
provided as well.
+computing domains. These are summarized in the following table:
Two other packages are installed at the higherlevel: scipy_distutils
and weave. These two packages while distributed with main scipy
package could see use independently of scipy and so are treated as
separate packages and described elsewhere.
+.. currentmodule:: scipy
The remaining subpackages are summarized in the following table (a *
denotes an optional subpackage that requires additional libraries to
function or is not available on all platforms).
+================== =====================================================================
+Subpackage Description
+================== =====================================================================
+:mod:`cluster` Clustering algorithms
+:mod:`constants` Physical and mathematical constants
+:mod:`fftpack` Fast Fourier Transform routines
+:mod:`integrate` Integration and ordinary differential equation solvers
+:mod:`interpolate` Interpolation and smoothing splines
+:mod:`io` Input and Output
+:mod:`linalg` Linear algebra
+:mod:`maxentropy` Maximum entropy methods
+:mod:`ndimage` Ndimensional image processing
+:mod:`odr` Orthogonal distance regression
+:mod:`optimize` Optimization and rootfinding routines
+:mod:`signal` Signal processing
+:mod:`sparse` Sparse matrices and associated routines
+:mod:`spatial` Spatial data structures and algorithms
+:mod:`special` Special functions
+:mod:`stats` Statistical distributions and functions
+:mod:`weave` C/C++ integration
+================== =====================================================================

=========== =====================================================================
Subpackage Description
=========== =====================================================================
cluster Clustering algorithms
cow Cluster of Workstations code for parallel programming
fftpack FFT based on fftpack  default
fftw* FFT based on fftw  requires FFTW libraries (is this still needed?)
ga Genetic algorithms
gplt* Plotting  requires gnuplot
integrate Integration
interpolate Interpolation
io Input and Output
linalg Linear algebra
optimize Optimization and rootfinding routines
plt* Plotting  requires wxPython
signal Signal processing
special Special functions
stats Statistical distributions and functions
xplt Plotting with gist
=========== =====================================================================



Because of their ubiquitousness, some of the functions in these
subpackages are also made available in the scipy namespace to ease
their use in interactive sessions and programs. In addition, many
convenience functions are located in :mod:`numpy` and the in the
+basic array functions from :mod:`numpy` are also available at the
toplevel of the :mod:`scipy` package. Before looking at the
subpackages individually, we will first look at some of these common
functions.
Basic functions in Numpy and toplevel scipy
============================================
+Basic functions in Numpy (and toplevel scipy)
+==============================================
+.. currentmodule:: numpy
+
Interaction with Numpy

To begin with, all of the Numpy functions have been subsumed into
the scipy namespace so that all of those functions are available
+To begin with, all of the Numpy functions have been subsumed into the
+:mod:`scipy` namespace so that all of those functions are available
without additionally importing Numpy. In addition, the universal
functions (addition, subtraction, division) have been altered to not
raise exceptions if floatingpoint errors are encountered,
instead NaN's and Inf's are returned in the arrays. To assist in
detection of these events new universal functions (isnan, isfinite,
isinf) have been added.
+raise exceptions if floatingpoint errors are encountered; instead,
+NaN's and Inf's are returned in the arrays. To assist in detection of
+these events, several functions (:func:`isnan`, :func:`isfinite`,
+:func:`isinf`) are available.
In addition, the comparision operators have been changed to allow
comparisons and logical operations of complex numbers (only the real
part is compared). Also, with the new universal functions in SciPy,
the logical operations (except logical_XXX functions) all return
arrays of unsigned bytes (8bits per element instead of the old
32bits, or even 64bits) per element [#]_ .

In an effort to get a consistency for default arguments, some of the
default arguments have changed from Numpy. The idea is for you to use
scipy as a base package instead of Numpy anyway.

Finally, some of the basic functions like log, sqrt, and inverse trig
functions have been modified to return complex numbers instead of
NaN's where appropriate (*i.e.* ``scipy.sqrt(1)`` returns ``1j``).
.. [#] Be careful when treating logical expressions as integers as the 8bit
 integers may silently overflow at 256.


Toplevel scipy routines

The purpose of the top level of scipy is to collect generalpurpose
routines that the other subpackages can use and to provide a simple
replacement for Numpy. Anytime you might think to import Numpy,
you can import scipy_base instead and remove yourself from direct
dependence on Numpy. These routines are divided into several files
for organizational purposes, but they are all available under the
scipy_base namespace (and the scipy namespace). There are routines for
type handling and type checking, shape and matrix manipulation,
polynomial processing, and other useful functions. Rather than giving
a detailed description of each of these functions (which is available
using the :func:`help`, :func:`info` and :func:`source` commands),
this tutorial will discuss some of the more useful commands which
require a little introduction to use to their full potential.
+replacement for Numpy. Anytime you might think to import Numpy, you
+can import scipy instead and remove yourself from direct dependence on
+Numpy. These routines are divided into several files for
+organizational purposes, but they are all available under the numpy
+namespace (and the scipy namespace). There are routines for type
+handling and type checking, shape and matrix manipulation, polynomial
+processing, and other useful functions. Rather than giving a detailed
+description of each of these functions (which is available in the
+Numpy Reference Guide or by using the :func:`help`, :func:`info` and
+:func:`source` commands), this tutorial will discuss some of the more
+useful commands which require a little introduction to use to their
+full potential.
Type handling
^^^^^^^^^^^^^
Note the difference between :func:`iscomplex` ( :func:`isreal` ) and :func:`iscomplexobj` ( :func:`isrealobj` ). The former command is array based and returns byte arrays of ones
and zeros providing the result of the elementwise test. The latter
command is object based and returns a scalar describing the result of
the test on the entire object.
+Note the difference between :func:`iscomplex` (:func:`isreal`) and
+:func:`iscomplexobj` (:func:`isrealobj`). The former command is
+array based and returns byte arrays of ones and zeros providing the
+result of the elementwise test. The latter command is object based
+and returns a scalar describing the result of the test on the entire
+object.
Often it is required to get just the real and/or imaginary part of a
complex number. While complex numbers and arrays have attributes that
@@ 213,9 +192,8 @@
the use of the :obj:`cast` dictionary. The dictionary is keyed by the
type it is desired to cast to and the dictionary stores functions to
perform the casting. Thus, ``>>> a = cast['f'](d)`` returns an array
of float32 from d. This function is also useful as an easy way to get
a scalar of a certain type: ``>>> fpi = cast['f'](pi)`` although this
should not be needed if you use alter_numeric().
+of :class:`float32` from *d*. This function is also useful as an easy
+way to get a scalar of a certain type: ``>>> fpi = cast['f'](pi)``.
Index Tricks
@@ 223,7 +201,8 @@
Thre are some class instances that make special use of the slicing
functionality to provide efficient means for array construction. This
part will discuss the operation of :obj:`mgrid` , :obj:`ogrid` , :obj:`r_` , and :obj:`c_` for quickly constructing arrays.
+part will discuss the operation of :obj:`mgrid` , :obj:`ogrid` ,
+:obj:`r_` , and :obj:`c_` for quickly constructing arrays.
One familiar with Matlab may complain that it is difficult to
construct arrays from the interactive session with Python. Suppose,
@@ 310,19 +289,11 @@
available.
Matrix manipulations
^^^^^^^^^^^^^^^^^^^^

These are functions specifically suited for 2dimensional arrays that
were part of MLab in the Numpy distribution, but have been placed in
scipy_base for completeness so that users are not importing Numpy.


Polynomials
^^^^^^^^^^^
There are two (interchangeable) ways to deal with 1d polynomials in
SciPy. The first is to use the :class:`poly1d` class in Numpy. This
+SciPy. The first is to use the :class:`poly1d` class from Numpy. This
class accepts coefficients or polynomial roots to initialize a
polynomial. The polynomial object can then be manipulated in algebraic
expressions, integrated, differentiated, and evaluated. It even prints
@@ 353,9 +324,9 @@
Vectorizing functions (vectorize)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
One of the features that SciPy provides is a class :obj:`vectorize` to
+One of the features that NumPy provides is a class :obj:`vectorize` to
convert an ordinary Python function which accepts scalars and returns
scalars into a "vectorizedfunction "with the same broadcasting rules
+scalars into a "vectorizedfunction" with the same broadcasting rules
as other Numpy functions (*i.e.* the Universal functions, or
ufuncs). For example, suppose you have a Python function named
:obj:`addsubtract` defined as:
@@ 387,20 +358,21 @@
^^^^^^^^^^^^^^^^^^^^^^
There are several other functions in the scipy_base package including
most of the other functions that are also in MLab that comes with the
Numpy package. The reason for duplicating these functions is to
allow SciPy to potentially alter their original interface and make it
easier for users to know how to get access to functions
+most of the other functions that are also in the Numpy package. The
+reason for duplicating these functions is to allow SciPy to
+potentially alter their original interface and make it easier for
+users to know how to get access to functions
>>> from scipy import *
New functions which should be mentioned are :obj:`mod(x,y)` which can
+Functions which should be mentioned are :obj:`mod(x,y)` which can
replace ``x % y`` when it is desired that the result take the sign of
*y* instead of *x* . Also included is :obj:`fix` which always rounds
to the nearest integer towards zero. For doing phase processing, the
functions :func:`angle`, and :obj:`unwrap` are also useful. Also, the
:obj:`linspace` and :obj:`logspace` functions return equally spaced samples
in a linear or log scale. Finally, mention should be made of the new
+in a linear or log scale. Finally, it's useful to be aware of the indexing
+capabilities of Numpy.mention should be made of the new
function :obj:`select` which extends the functionality of :obj:`where` to
include multiple conditions and multiple choices. The calling
convention is ``select(condlist,choicelist,default=0).`` :obj:`select` is
@@ 435,16 +407,21 @@
returns a common image used in image processing: :obj:`lena`.
Finally, two functions are provided that are useful for approximating
derivatives of functions using discretedifferences. The function :obj:`central_diff_weights` returns weighting coefficients for an equallyspaced :math:`N` point approximation to the derivative of order :math:`o.` These weights must be multiplied by the function corresponding to
these points and the results added to obtain the derivative
approximation. This function is intended for use when only samples of
the function are avaiable. When the function is an object that can be
handed to a routine and evaluated, the function :obj:`derivative` can be used to automatically evaluate the object at the correct points
to obtain an Npoint approximation to the :math:`o^{\textrm{th}}` derivative at a given point.
+derivatives of functions using discretedifferences. The function
+:obj:`central_diff_weights` returns weighting coefficients for an
+equallyspaced :math:`N`point approximation to the derivative of
+order *o*. These weights must be multiplied by the function
+corresponding to these points and the results added to obtain the
+derivative approximation. This function is intended for use when only
+samples of the function are avaiable. When the function is an object
+that can be handed to a routine and evaluated, the function
+:obj:`derivative` can be used to automatically evaluate the object at
+the correct points to obtain an Npoint approximation to the *o*th
+derivative at a given point.
Special functions (special)
===========================
+Special functions (:mod:`scipy.special`)
+========================================
.. currentmodule:: scipy.special
@@ 459,15 +436,15 @@
broadcasting rules as other math functions in Numerical Python. Many
of these functions also accept complexnumbers as input. For a
complete list of the available functions with a oneline description
type ``>>> info(special).`` Each function also has it's own
+type ``>>> help(special).`` Each function also has it's own
documentation accessible using help. If you don't see a function you
need, consider writing it and contributing it to the library. You can
write the function in either C, Fortran, or Python. Look in the source
code of the library for examples of each of these kind of functions.
Integration (integrate)
=======================
+Integration (:mod:`scipy.integrate`)
+====================================
.. currentmodule:: scipy.integrate
@@ 478,8 +455,8 @@
.. literalinclude:: examples/41
General integration (integrate.quad)

+General integration (:func:`quad`)
+
The function :obj:`quad` is provided to integrate a function of one
variable between two points. The points can be :math:`\pm\infty`
@@ 506,7 +483,8 @@
>>> print abs(result[0]I)
1.03761443881e11
The first argument to quad is a "callable "Python object ( *i.e* a function, method, or class instance). Notice the use of a lambda
+The first argument to quad is a "callable" Python object (*i.e* a
+function, method, or class instance). Notice the use of a lambda
function in this case as the argument. The next two arguments are the
limits of integration. The return value is a tuple, with the first
element holding the estimated value of the integral and the second
@@ 610,42 +588,49 @@
:mod:`special.orthogonal` which can calculate the roots and quadrature
weights of a large variety of orthogonal polynomials (the polynomials
themselves are available as special functions returning instances of
the polynomial class  e.g. ``special.legendre`` ).
+the polynomial class  e.g. :obj:`special.legendre <scipy.special.legendre>`).
Integrating using samples

There are three functions for computing integrals given only samples: :obj:`trapz` , :obj:`simps` , and :obj:`romb` . The first two functions use NewtonCoates formulas of order 1 and 2
respectively to perform integration. These two functions can handle,
+There are three functions for computing integrals given only samples:
+:obj:`trapz` , :obj:`simps`, and :obj:`romb` . The first two
+functions use NewtonCoates formulas of order 1 and 2 respectively to
+perform integration. These two functions can handle,
nonequallyspaced samples. The trapezoidal rule approximates the
function as a straight line between adjacent points, while Simpson's
rule approximates the function between three adjacent points as a
parabola.
If the samples are equallyspaced and the number of samples available
is :math:`2^{k}+1` for some integer :math:`k` , then Romberg integration can be used to obtain highprecision
estimates of the integral using the available samples. Romberg
integration uses the trapezoid rule at stepsizes related by a power
of two and then performs Richardson extrapolation on these estimates
to approximate the integral with a higherdegree of accuracy. (A
different interface to Romberg integration useful when the function
can be provided is also available as :func:`integrate.romberg`).
+is :math:`2^{k}+1` for some integer :math:`k`, then Romberg
+integration can be used to obtain highprecision estimates of the
+integral using the available samples. Romberg integration uses the
+trapezoid rule at stepsizes related by a power of two and then
+performs Richardson extrapolation on these estimates to approximate
+the integral with a higherdegree of accuracy. (A different interface
+to Romberg integration useful when the function can be provided is
+also available as :func:`romberg`).
Ordinary differential equations (integrate.odeint)

+Ordinary differential equations (:func:`odeint`)
+
Integrating a set of ordinary differential equations (ODEs) given
initial conditions is another useful example. The function :obj:`odeint` is available in SciPy for integrating a firstorder vector
differential equation:
+initial conditions is another useful example. The function
+:obj:`odeint` is available in SciPy for integrating a firstorder
+vector differential equation:
.. math::
:nowrap:
\[ \frac{d\mathbf{y}}{dt}=\mathbf{f}\left(\mathbf{y},t\right),\]
given initial conditions :math:`\mathbf{y}\left(0\right)=y_{0},` where :math:`\mathbf{y}` is a length :math:`N` vector and :math:`\mathbf{f}` is a mapping from :math:`\mathcal{R}^{N}` to :math:`\mathcal{R}^{N}.` A higherorder ordinary differential equation can always be reduced to
+given initial conditions :math:`\mathbf{y}\left(0\right)=y_{0}`, where
+:math:`\mathbf{y}` is a length :math:`N` vector and :math:`\mathbf{f}`
+is a mapping from :math:`\mathcal{R}^{N}` to :math:`\mathcal{R}^{N}.`
+A higherorder ordinary differential equation can always be reduced to
a differential equation of this type by introducing intermediate
derivatives into the :math:`\mathbf{y}` vector.
@@ 665,9 +650,11 @@
\[ w=\textrm{Ai}\left(z\right),\]
which gives a means to check the integrator using :func:`special.airy`.
+which gives a means to check the integrator using :func:`special.airy <scipy.special.airy>`.
First, convert this ODE into standard form by setting :math:`\mathbf{y}=\left[\frac{dw}{dz},w\right]` and :math:`t=z.` Thus, the differential equation becomes
+First, convert this ODE into standard form by setting
+:math:`\mathbf{y}=\left[\frac{dw}{dz},w\right]` and :math:`t=z`. Thus,
+the differential equation becomes
.. math::
:nowrap:
@@ 681,9 +668,9 @@
\[ \mathbf{f}\left(\mathbf{y},t\right)=\mathbf{A}\left(t\right)\mathbf{y}.\]


As an interesting reminder, if :math:`\mathbf{A}\left(t\right)` commutes with :math:`\int_{0}^{t}\mathbf{A}\left(\tau\right)\, d\tau` under matrix multiplication, then this linear differential equation
+As an interesting reminder, if :math:`\mathbf{A}\left(t\right)`
+commutes with :math:`\int_{0}^{t}\mathbf{A}\left(\tau\right)\, d\tau`
+under matrix multiplication, then this linear differential equation
has an exact solution using the matrix exponential:
.. math::
@@ 739,6 +726,8 @@
Optimization (optimize)
=======================
+.. currentmodule:: scipy.optimize
+
There are several classical optimization algorithms provided by SciPy
in the :mod:`scipy.optimize` package. An overview of the module is
available using :func:`help` (or :func:`pydoc.help`):
@@ 746,17 +735,16 @@
.. literalinclude:: examples/51
The first four algorithms are unconstrained minimization algorithms
(fmin: NelderMead simplex, fmin_bfgs: BFGS, fmin_ncg: Newton
Conjugate Gradient, and leastsq: LevenburgMarquardt). The fourth
algorithm only works for functions of a single variable but allows
minimization over a specified interval. The last algorithm actually
finds the roots of a general function of possibly many variables. It
is included in the optimization package because at the (nonboundary)
extreme points of a function, the gradient is equal to zero.
+(:func:`fmin`: NelderMead simplex, :func:`fmin_bfgs`: BFGS,
+:func:`fmin_ncg`: Newton Conjugate Gradient, and :func:`leastsq`:
+LevenburgMarquardt). The last algorithm actually finds the roots of a
+general function of possibly many variables. It is included in the
+optimization package because at the (nonboundary) extreme points of a
+function, the gradient is equal to zero.
NelderMead Simplex algorithm (optimize.fmin)

+NelderMead Simplex algorithm (:func:`fmin`)
+
The simplex algorithm is probably the simplest way to minimize a
fairly wellbehaved function. The simplex algorithm requires only
@@ 789,11 +777,11 @@
[ 1. 1. 1. 1. 1.]
Another optimization algorithm that needs only function calls to find
the minimum is Powell's method available as :func:`optimize.fmin_powell`.
+the minimum is Powell's method available as :func:`fmin_powell`.
BroydenFletcherGoldfarbShanno algorithm (optimize.fmin_bfgs)

+BroydenFletcherGoldfarbShanno algorithm (:func:`fmin_bfgs`)
+
In order to converge more quickly to the solution, this routine uses
the gradient of the objective function. If the gradient is not given
@@ 849,8 +837,8 @@
[ 1. 1. 1. 1. 1.]
NewtonConjugateGradient (optimize.fmin_ncg)

+NewtonConjugateGradient (:func:`fmin_ncg`)
+
The method which requires the fewest function calls and is therefore
often the fastest method to minimize functions of many variables is
@@ 983,16 +971,20 @@
[ 1. 1. 1. 1. 1.]
Leastsquare fitting (minimize.leastsq)

+Leastsquare fitting (:func:`leastsq`)
+
All of the previouslyexplained minimization procedures can be used to
solve a leastsquares problem provided the appropriate objective
function is constructed. For example, suppose it is desired to fit a
set of data :math:`\left\{ \mathbf{x}_{i},\mathbf{y}_{i}\right\} ` to a known model, :math:`\mathbf{y}=\mathbf{f}\left(\mathbf{x},\mathbf{p}\right)` where :math:`\mathbf{p}` is a vector of parameters for the model that need to be found. A
common method for determining which parameter vector gives the best
fit to the data is to minimize the sum of squares of the residuals.
The residual is usually defined for each observed datapoint as
+set of data :math:`\left\{\mathbf{x}_{i}, \mathbf{y}_{i}\right\}`
+to a known model,
+:math:`\mathbf{y}=\mathbf{f}\left(\mathbf{x},\mathbf{p}\right)`
+where :math:`\mathbf{p}` is a vector of parameters for the model that
+need to be found. A common method for determining which parameter
+vector gives the best fit to the data is to minimize the sum of squares
+of the residuals. The residual is usually defined for each observed
+datapoint as
.. math::
:nowrap:
@@ 1009,8 +1001,12 @@
The :obj:`leastsq` algorithm performs this squaring and summing of the residuals
automatically. It takes as an input argument the vector function :math:`\mathbf{e}\left(\mathbf{p}\right)` and returns the value of :math:`\mathbf{p}` which minimizes :math:`J\left(\mathbf{p}\right)=\mathbf{e}^{T}\mathbf{e}` directly. The user is also encouraged to provide the Jacobian matrix
+The :obj:`leastsq` algorithm performs this squaring and summing of the
+residuals automatically. It takes as an input argument the vector
+function :math:`\mathbf{e}\left(\mathbf{p}\right)` and returns the
+value of :math:`\mathbf{p}` which minimizes
+:math:`J\left(\mathbf{p}\right)=\mathbf{e}^{T}\mathbf{e}`
+directly. The user is also encouraged to provide the Jacobian matrix
of the function (with derivatives down the columns or across the
rows). If the Jacobian is not provided, it is estimated.
@@ 1081,8 +1077,8 @@
been developed that can work faster.
Unconstrained minimization (optimize.brent)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+Unconstrained minimization (:func:`brent`)
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
There are actually two methods that can be used to minimize a scalar
function (:obj:`brent` and :func:`golden`), but :obj:`golden` is
@@ 1098,16 +1094,16 @@
for your function and result in an unexpected minimum being returned).
Bounded minimization (optimize.fminbound)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+Bounded minimization (:func:`fminbound`)
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Thus far all of the minimization routines described have been
unconstrained minimization routines. Very often, however, there are
constraints that can be placed on the solution space before
minimization occurs. The :obj:`fminbound` function is an example of a constrained minimization procedure that
provides a rudimentary interval constraint for scalar functions. The
interval constraint allows the minimization to occur only between two
fixed endpoints.
+minimization occurs. The :obj:`fminbound` function is an example of a
+constrained minimization procedure that provides a rudimentary
+interval constraint for scalar functions. The interval constraint
+allows the minimization to occur only between two fixed endpoints.
For example, to find the minimum of :math:`J_{1}\left(x\right)` near :math:`x=5` , :obj:`fminbound` can be called using the interval :math:`\left[4,7\right]` as a constraint. The result is :math:`x_{\textrm{min}}=5.3314` :
@@ 1125,9 +1121,11 @@
Sets of equations
^^^^^^^^^^^^^^^^^
To find the roots of a polynomial, the command :obj:`roots` is useful. To find a root of a set of nonlinear equations, the
command :obj:`optimize.fsolve` is needed. For example, the following example finds the roots of the
singlevariable transcendental equation
+To find the roots of a polynomial, the command :obj:`roots
+<scipy.roots>` is useful. To find a root of a set of nonlinear
+equations, the command :obj:`fsolve` is needed. For example, the
+following example finds the roots of the singlevariable
+transcendental equation
.. math::
:nowrap:
@@ 1168,8 +1166,9 @@
If one has a singlevariable equation, there are four different root
finder algorithms that can be tried. Each of these root finding
algorithms requires the endpoints of an interval where a root is
suspected (because the function changes signs). In general :obj:`brentq` is the best choice, but the other methods may be useful in certain
circumstances or for academic purposes.
+suspected (because the function changes signs). In general
+:obj:`brentq` is the best choice, but the other methods may be useful
+in certain circumstances or for academic purposes.
Fixedpoint solving
@@ 1178,13 +1177,20 @@
A problem closely related to finding the zeros of a function is the
problem of finding a fixedpoint of a function. A fixed point of a
function is the point at which evaluation of the function returns the
point: :math:`g\left(x\right)=x.` Clearly the fixed point of :math:`g` is the root of :math:`f\left(x\right)=g\left(x\right)x.` Equivalently, the root of :math:`f` is the fixed_point of :math:`g\left(x\right)=f\left(x\right)+x.` The routine :obj:`fixed_point` provides a simple iterative method using Aitkens sequence acceleration
to estimate the fixed point of :math:`g` given a starting point.
+point: :math:`g\left(x\right)=x.` Clearly the fixed point of :math:`g`
+is the root of :math:`f\left(x\right)=g\left(x\right)x.`
+Equivalently, the root of :math:`f` is the fixed_point of
+:math:`g\left(x\right)=f\left(x\right)+x.` The routine
+:obj:`fixed_point` provides a simple iterative method using Aitkens
+sequence acceleration to estimate the fixed point of :math:`g` given a
+starting point.
Interpolation (interpolate)
===========================
+Interpolation (:mod:`scipy.interpolate`)
+========================================
+.. currentmodule:: scipy.interpolate
+
There are two general interpolation facilities available in SciPy. The
first facility is an interpolation class which performs linear
1dimensional interpolation. The second facility is based on the
@@ 1192,19 +1198,19 @@
2dimensional (smoothed) cubicspline interpolation.
Linear 1d interpolation (interpolate.interp1d)

+Linear 1d interpolation (:class:`interp1d`)
+
The interp1d class in scipy.interpolate is a convenient method to
create a function based on fixed data points which can be evaluated
anywhere within the domain defined by the given data using linear
interpolation. An instance of this class is created by passing the 1d
vectors comprising the data. The instance of this class defines a
*__call__* method and can therefore by treated like a function which
interpolates between known data values to obtain unknown values (it
even has a docstring for help). Behavior at the boundary can be
specified at instantiation time. The following example demonstrates
it's use.
+:meth:`__call__ <interp1d.__call__>` method and can therefore by
+treated like a function which interpolates between known data values
+to obtain unknown values (it also has a docstring for help). Behavior
+at the boundary can be specified at instantiation time. The following
+example demonstrates it's use.
.. plot::
@@ 1232,7 +1238,7 @@
representation, there are two different was to represent a curve and
obtain (smoothing) spline coefficients: directly and parametrically.
The direct method finds the spline representation of a curve in a two
dimensional plane using the function :obj:`interpolate.splrep`. The
+dimensional plane using the function :obj:`splrep`. The
first two arguments are the only ones required, and these provide the
:math:`x` and :math:`y` components of the curve. The normal output is
a 3tuple, :math:`\left(t,c,k\right)` , containing the knotpoints,
@@ 1241,7 +1247,7 @@
with the input keyword, *k.*
For curves in :math:`N` dimensional space the function
:obj:`interpolate.splprep` allows defining the curve
+:obj:`splprep` allows defining the curve
parametrically. For this function only 1 input argument is
required. This input is a list of :math:`N` arrays representing the
curve in :math:`N` dimensional space. The length of each array is the
@@ 1261,12 +1267,12 @@
Once the spline representation of the data has been determined,
functions are available for evaluating the spline
(:func:`interpolate.splev`) and its derivatives
(:func:`interpolate.splev`, :func:`interpolate.splade`) at any point
+(:func:`splev`) and its derivatives
+(:func:`splev`, :func:`splade`) at any point
and the integral of the spline between any two points (
:func:`interpolate.splint`). In addition, for cubic splines ( :math:`k=3`
+:func:`splint`). In addition, for cubic splines ( :math:`k=3`
) with 8 or more knots, the roots of the spline can be estimated (
:func:`interpolate.sproot`). These functions are demonstrated in the
+:func:`sproot`). These functions are demonstrated in the
example that follows.
.. plot::
@@ 1339,31 +1345,31 @@
>>> plt.show()
Twodimensional spline representation (interpolate.bisplrep)

+Twodimensional spline representation (:func:`bisplrep`)
+
For (smooth) splinefitting to a two dimensional surface, the function
:obj:`interpolate.bisplrep` is available. This function takes as
required inputs the **1D** arrays *x*, *y*, and *z* which represent
points on the surface :math:`z=f\left(x,y\right).` The default output
is a list :math:`\left[tx,ty,c,kx,ky\right]` whose entries represent
+:func:`bisplrep` is available. This function takes as required inputs
+the **1D** arrays *x*, *y*, and *z* which represent points on the
+surface :math:`z=f\left(x,y\right).` The default output is a list
+:math:`\left[tx,ty,c,kx,ky\right]` whose entries represent
respectively, the components of the knot positions, the coefficients
of the spline, and the order of the spline in each coordinate. It is
convenient to hold this list in a single object, *tck,* so that it can
be passed easily to the function :obj:`interpolate.bisplev`. The
+be passed easily to the function :obj:`bisplev`. The
keyword, *s* , can be used to change the amount of smoothing performed
on the data while determining the appropriate spline. The default
value is :math:`s=m\sqrt{2m}` where :math:`m` is the number of data
points in the *x, y,* and *z* vectors. As a result, if no smoothing is
desired, then :math:`s=0` should be passed to
:obj:`interpolate.bisplrep` .
+:obj:`bisplrep` .
To evaluate the twodimensional spline and it's partial derivatives
(up to the order of the spline), the function
:obj:`interpolate.bisplev` is required. This function takes as the
+:obj:`bisplev` is required. This function takes as the
first two arguments **two 1D arrays** whose crossproduct specifies
the domain over which to evaluate the spline. The third argument is
the *tck* list returned from :obj:`interpolate.bisplrep`. If desired,
+the *tck* list returned from :obj:`bisplrep`. If desired,
the fourth and fifth arguments provide the orders of the partial
derivative in the :math:`x` and :math:`y` direction respectively.
@@ 1373,13 +1379,14 @@
processing toolbox contains more appropriate algorithms for finding
the spline representation of an image. The two dimensional
interpolation commands are intended for use when interpolating a two
dimensional function as shown in the example that follows (See also
Figure `4 <#fig2dinterp>`__ ). This example uses the :obj:`mgrid` command in SciPy which is useful for defining a "meshgrid "in many dimensions. (See also the :obj:`ogrid` command if the fullmesh is not needed). The number of output
arguments and the number of dimensions of each argument is determined
by the number of indexing objects passed in :obj:`mgrid`[].
+dimensional function as shown in the example that follows. This
+example uses the :obj:`mgrid <numpy.mgrid>` command in SciPy which is
+useful for defining a "meshgrid "in many dimensions. (See also the
+:obj:`ogrid <numpy.ogrid>` command if the fullmesh is not
+needed). The number of output arguments and the number of dimensions
+of each argument is determined by the number of indexing objects
+passed in :obj:`mgrid <numpy.mgrid>`[].
.. _`fig:2d_interp`:

.. plot::
>>> from numpy import *
@@ 1528,8 +1535,6 @@
:obj:`signal.convolve2d` which convolves arbitrary twodimensional
filters and allows for choosing mirrorsymmetric boundary conditions.
.. _`fig:lena_edge_spline`:

.. plot::
>>> from numpy import *
@@ 1602,7 +1607,11 @@
This equation can only be implemented directly if we limit the
sequences to finite support sequences that can be stored in a
computer, choose :math:`n=0` to be the starting point of both sequences, let :math:`K+1` be that value for which :math:`y\left[n\right]=0` for all :math:`n>K+1` and :math:`M+1` be that value for which :math:`x\left[n\right]=0` for all :math:`n>M+1` , then the discrete convolution expression is
+computer, choose :math:`n=0` to be the starting point of both
+sequences, let :math:`K+1` be that value for which
+:math:`y\left[n\right]=0` for all :math:`n>K+1` and :math:`M+1` be
+that value for which :math:`x\left[n\right]=0` for all :math:`n>M+1` ,
+then the discrete convolution expression is
.. math::
:nowrap:
@@ 1618,12 +1627,25 @@
Thus, the full discrete convolution of two finite sequences of lengths :math:`K+1` and :math:`M+1` respectively results in a finite sequence of length :math:`K+M+1=\left(K+1\right)+\left(M+1\right)1.`
One dimensional convolution is implemented in SciPy with the function ``signal.convolve`` . This function takes as inputs the signals :math:`x,` :math:`h` , and an optional flag and returns the signal :math:`y.` The optional flag allows for specification of which part of the output
signal to return. The default value of 'full' returns the entire
signal. If the flag has a value of 'same' then only the middle :math:`K` values are returned starting at :math:`y\left[\left\lfloor \frac{M1}{2}\right\rfloor \right]` so that the output has the same length as the largest input. If the
flag has a value of 'valid' then only the middle :math:`KM+1=\left(K+1\right)\left(M+1\right)+1` output values are returned where :math:`z` depends on all of the values of the smallest input from :math:`h\left[0\right]` to :math:`h\left[M\right].` In other words only the values :math:`y\left[M\right]` to :math:`y\left[K\right]` inclusive are returned.
+One dimensional convolution is implemented in SciPy with the function
+``signal.convolve`` . This function takes as inputs the signals
+:math:`x,` :math:`h` , and an optional flag and returns the signal
+:math:`y.` The optional flag allows for specification of which part of
+the output signal to return. The default value of 'full' returns the
+entire signal. If the flag has a value of 'same' then only the middle
+:math:`K` values are returned starting at :math:`y\left[\left\lfloor
+\frac{M1}{2}\right\rfloor \right]` so that the output has the same
+length as the largest input. If the flag has a value of 'valid' then
+only the middle :math:`KM+1=\left(K+1\right)\left(M+1\right)+1`
+output values are returned where :math:`z` depends on all of the
+values of the smallest input from :math:`h\left[0\right]` to
+:math:`h\left[M\right].` In other words only the values
+:math:`y\left[M\right]` to :math:`y\left[K\right]` inclusive are
+returned.
This same function ``signal.convolve`` can actually take :math:`N` dimensional arrays as inputs and will return the :math:`N` dimensional convolution of the two arrays. The same input flags are
+This same function ``signal.convolve`` can actually take :math:`N`
+dimensional arrays as inputs and will return the :math:`N`
+dimensional convolution of the two arrays. The same input flags are
available for that case as well.
Correlation is very similar to convolution except for the minus sign
@@ 1650,19 +1672,29 @@
The SciPy function ``signal.correlate`` implements this operation. Equivalent flags are available for this
operation to return the full :math:`K+M+1` length sequence ('full') or a sequence with the same size as the
largest sequence starting at :math:`w\left[K+\left\lfloor \frac{M1}{2}\right\rfloor \right]` ('same') or a sequence where the values depend on all the values of
the smallest sequence ('valid'). This final option returns the :math:`KM+1` values :math:`w\left[MK\right]` to :math:`w\left[0\right]` inclusive.
+The SciPy function ``signal.correlate`` implements this
+operation. Equivalent flags are available for this operation to return
+the full :math:`K+M+1` length sequence ('full') or a sequence with the
+same size as the largest sequence starting at
+:math:`w\left[K+\left\lfloor \frac{M1}{2}\right\rfloor \right]`
+('same') or a sequence where the values depend on all the values of
+the smallest sequence ('valid'). This final option returns the
+:math:`KM+1` values :math:`w\left[MK\right]` to
+:math:`w\left[0\right]` inclusive.
The function :obj:`signal.correlate` can also take arbitrary :math:`N` dimensional arrays as input and return the :math:`N` dimensional convolution of the two arrays on output.
+The function :obj:`signal.correlate` can also take arbitrary :math:`N`
+dimensional arrays as input and return the :math:`N` dimensional
+convolution of the two arrays on output.
When :math:`N=2,` :obj:`signal.correlate` and/or :obj:`signal.convolve` can be used to construct arbitrary image filters to perform actions
such as blurring, enhancing, and edgedetection for an image.
+When :math:`N=2,` :obj:`signal.correlate` and/or
+:obj:`signal.convolve` can be used to construct arbitrary image
+filters to perform actions such as blurring, enhancing, and
+edgedetection for an image.
Convolution is mainly used for filtering when one of the signals is
much smaller than the other ( :math:`K\gg M` ), otherwise linear filtering is more easily accomplished in the
frequency domain (see Fourier Transforms).
+much smaller than the other ( :math:`K\gg M` ), otherwise linear
+filtering is more easily accomplished in the frequency domain (see
+Fourier Transforms).
Differenceequation filtering
@@ 1676,9 +1708,15 @@
\[ \sum_{k=0}^{N}a_{k}y\left[nk\right]=\sum_{k=0}^{M}b_{k}x\left[nk\right]\]
where :math:`x\left[n\right]` is the input sequence and :math:`y\left[n\right]` is the output sequence. If we assume initial rest so that :math:`y\left[n\right]=0` for :math:`n<0` , then this kind of filter can be implemented using convolution.
However, the convolution filter sequence :math:`h\left[n\right]` could be infinite if :math:`a_{k}\neq0` for :math:`k\geq1.` In addition, this general class of linear filter allows initial
conditions to be placed on :math:`y\left[n\right]` for :math:`n<0` resulting in a filter that cannot be expressed using convolution.
+where :math:`x\left[n\right]` is the input sequence and
+:math:`y\left[n\right]` is the output sequence. If we assume initial
+rest so that :math:`y\left[n\right]=0` for :math:`n<0` , then this
+kind of filter can be implemented using convolution. However, the
+convolution filter sequence :math:`h\left[n\right]` could be infinite
+if :math:`a_{k}\neq0` for :math:`k\geq1.` In addition, this general
+class of linear filter allows initial conditions to be placed on
+:math:`y\left[n\right]` for :math:`n<0` resulting in a filter that
+cannot be expressed using convolution.
The difference equation filter can be thought of as finding :math:`y\left[n\right]` recursively in terms of it's previous values
@@ 1687,28 +1725,46 @@
\[ a_{0}y\left[n\right]=a_{1}y\left[n1\right]\cdotsa_{N}y\left[nN\right]+\cdots+b_{0}x\left[n\right]+\cdots+b_{M}x\left[nM\right].\]
Often :math:`a_{0}=1` is chosen for normalization. The implementation in SciPy of this
general difference equation filter is a little more complicated then
would be implied by the previous equation. It is implemented so that
only one signal needs to be delayed. The actual implementation
equations are (assuming :math:`a_{0}=1` ).
+Often :math:`a_{0}=1` is chosen for normalization. The implementation
+in SciPy of this general difference equation filter is a little more
+complicated then would be implied by the previous equation. It is
+implemented so that only one signal needs to be delayed. The actual
+implementation equations are (assuming :math:`a_{0}=1` ).
.. math::
:nowrap:
\begin{eqnarray*} y\left[n\right] & = & b_{0}x\left[n\right]+z_{0}\left[n1\right]\\ z_{0}\left[n\right] & = & b_{1}x\left[n\right]+z_{1}\left[n1\right]a_{1}y\left[n\right]\\ z_{1}\left[n\right] & = & b_{2}x\left[n\right]+z_{2}\left[n1\right]a_{2}y\left[n\right]\\ \vdots & \vdots & \vdots\\ z_{K2}\left[n\right] & = & b_{K1}x\left[n\right]+z_{K1}\left[n1\right]a_{K1}y\left[n\right]\\ z_{K1}\left[n\right] & = & b_{K}x\left[n\right]a_{K}y\left[n\right],\end{eqnarray*}
where :math:`K=\max\left(N,M\right).` Note that :math:`b_{K}=0` if :math:`K>M` and :math:`a_{K}=0` if :math:`K>N.` In this way, the output at time :math:`n` depends only on the input at time :math:`n` and the value of :math:`z_{0}` at the previous time. This can always be calculated as long as the :math:`K` values :math:`z_{0}\left[n1\right]\ldots z_{K1}\left[n1\right]` are computed and stored at each time step.
+where :math:`K=\max\left(N,M\right).` Note that :math:`b_{K}=0` if
+:math:`K>M` and :math:`a_{K}=0` if :math:`K>N.` In this way, the
+output at time :math:`n` depends only on the input at time :math:`n`
+and the value of :math:`z_{0}` at the previous time. This can always
+be calculated as long as the :math:`K` values
+:math:`z_{0}\left[n1\right]\ldots z_{K1}\left[n1\right]` are
+computed and stored at each time step.
The differenceequation filter is called using the command :obj:`signal.lfilter` in SciPy. This command takes as inputs the vector :math:`b,` the vector, :math:`a,` a signal :math:`x` and returns the vector :math:`y` (the same length as :math:`x` ) computed using the equation given above. If :math:`x` is :math:`N` dimensional, then the filter is computed along the axis provided. If,
desired, initial conditions providing the values of :math:`z_{0}\left[1\right]` to :math:`z_{K1}\left[1\right]` can be provided or else it will be assumed that they are all zero. If
initial conditions are provided, then the final conditions on the
intermediate variables are also returned. These could be used, for
example, to restart the calculation in the same state.
+The differenceequation filter is called using the command
+:obj:`signal.lfilter` in SciPy. This command takes as inputs the
+vector :math:`b,` the vector, :math:`a,` a signal :math:`x` and
+returns the vector :math:`y` (the same length as :math:`x` ) computed
+using the equation given above. If :math:`x` is :math:`N`
+dimensional, then the filter is computed along the axis provided. If,
+desired, initial conditions providing the values of
+:math:`z_{0}\left[1\right]` to :math:`z_{K1}\left[1\right]` can be
+provided or else it will be assumed that they are all zero. If initial
+conditions are provided, then the final conditions on the intermediate
+variables are also returned. These could be used, for example, to
+restart the calculation in the same state.
Sometimes it is more convenient to express the initial conditions in
terms of the signals :math:`x\left[n\right]` and :math:`y\left[n\right].` In other words, perhaps you have the values of :math:`x\left[M\right]` to :math:`x\left[1\right]` and the values of :math:`y\left[N\right]` to :math:`y\left[1\right]` and would like to determine what values of :math:`z_{m}\left[1\right]` should be delivered as initial conditions to the differenceequation
filter. It is not difficult to show that for :math:`0\leq m<K,`
+terms of the signals :math:`x\left[n\right]` and
+:math:`y\left[n\right].` In other words, perhaps you have the values
+of :math:`x\left[M\right]` to :math:`x\left[1\right]` and the values
+of :math:`y\left[N\right]` to :math:`y\left[1\right]` and would like
+to determine what values of :math:`z_{m}\left[1\right]` should be
+delivered as initial conditions to the differenceequation filter. It
+is not difficult to show that for :math:`0\leq m<K,`
.. math::
:nowrap:
@@ 1736,8 +1792,9 @@
neighborhood values. If there are an even number of elements in the
neighborhood, then the average of the middle two values is used as the
median. A general purpose median filter that works on Ndimensional
arrays is :obj:`signal.medfilt` . A specialized version that works only for twodimensional arrays is
available as :obj:`signal.medfilt2d` .
+arrays is :obj:`signal.medfilt` . A specialized version that works
+only for twodimensional arrays is available as
+:obj:`signal.medfilt2d` .
Order Filter
@@ 1773,16 +1830,22 @@
\[ y=\left\{ \begin{array}{cc} \frac{\sigma^{2}}{\sigma_{x}^{2}}m_{x}+\left(1\frac{\sigma^{2}}{\sigma_{x}^{2}}\right)x & \sigma_{x}^{2}\geq\sigma^{2},\\ m_{x} & \sigma_{x}^{2}<\sigma^{2}.\end{array}\right.\]
Where :math:`m_{x}` is the local estimate of the mean and :math:`\sigma_{x}^{2}` is the local estimate of the variance. The window for these estimates
is an optional input parameter (default is :math:`3\times3` ). The parameter :math:`\sigma^{2}` is a threshold noise parameter. If :math:`\sigma` is not given then it is estimated as the average of the local
variances.
+Where :math:`m_{x}` is the local estimate of the mean and
+:math:`\sigma_{x}^{2}` is the local estimate of the variance. The
+window for these estimates is an optional input parameter (default is
+:math:`3\times3` ). The parameter :math:`\sigma^{2}` is a threshold
+noise parameter. If :math:`\sigma` is not given then it is estimated
+as the average of the local variances.
Hilbert filter
""""""""""""""
The Hilbert transform constructs the complexvalued analytic signal
from a real signal. For example if :math:`x=\cos\omega n` then :math:`y=\textrm{hilbert}\left(x\right)` would return (except near the edges) :math:`y=\exp\left(j\omega n\right).` In the frequency domain, the hilbert transform performs
+from a real signal. For example if :math:`x=\cos\omega n` then
+:math:`y=\textrm{hilbert}\left(x\right)` would return (except near the
+edges) :math:`y=\exp\left(j\omega n\right).` In the frequency domain,
+the hilbert transform performs
.. math::
:nowrap:
@@ 1791,118 +1854,120 @@
where :math:`H` is 2 for positive frequencies, :math:`0` for negative frequencies and :math:`1` for zerofrequencies.
+.. XXX: TODO
+..
+.. Detrend
+.. """""""
+..
+.. Filter design
+.. 
+..
+..
+.. Finiteimpulse response design
+.. ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+..
+..
+.. Inifiniteimpulse response design
+.. ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+..
+..
+.. Analog filter frequency response
+.. ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+..
+..
+.. Digital filter frequency response
+.. ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+..
+..
+.. Linear TimeInvariant Systems
+.. 
+..
+..
+.. LTI Object
+.. ^^^^^^^^^^
+..
+..
+.. ContinuousTime Simulation
+.. ^^^^^^^^^^^^^^^^^^^^^^^^^^
+..
+..
+.. Step response
+.. ^^^^^^^^^^^^^
+..
+..
+.. Impulse response
+.. ^^^^^^^^^^^^^^^^
+..
+..
+.. Input/Output
+.. ============
+..
+..
+.. Binary
+.. 
+..
+..
+.. Arbitrary binary input and output (fopen)
+.. ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+..
+..
+.. Read and write Matlab .mat files
+.. ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+..
+..
+.. Saving workspace
+.. ^^^^^^^^^^^^^^^^
+..
+..
+.. Textfile
+.. 
+..
+..
+.. Read textfiles (read_array)
+.. ^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+..
+..
+.. Write a textfile (write_array)
+.. ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+..
+..
+.. Fourier Transforms
+.. ==================
+..
+..
+.. Onedimensional
+.. 
+..
+..
+.. Twodimensional
+.. 
+..
+..
+.. Ndimensional
+.. 
+..
+..
+.. Shifting
+.. 
+..
+..
+.. Sample frequencies
+.. 
+..
+..
+.. Hilbert transform
+.. 
+..
+..
+.. Tilbert transform
+.. 
+..
Detrend
"""""""


Filter design



Finiteimpulse response design
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^


Inifiniteimpulse response design
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^


Analog filter frequency response
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^


Digital filter frequency response
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^


Linear TimeInvariant Systems



LTI Object
^^^^^^^^^^


ContinuousTime Simulation
^^^^^^^^^^^^^^^^^^^^^^^^^^


Step response
^^^^^^^^^^^^^


Impulse response
^^^^^^^^^^^^^^^^


Input/Output
============


Binary



Arbitrary binary input and output (fopen)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^


Read and write Matlab .mat files
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^


Saving workspace
^^^^^^^^^^^^^^^^


Textfile



Read textfiles (read_array)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^


Write a textfile (write_array)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^


Fourier Transforms
==================


Onedimensional



Twodimensional



Ndimensional



Shifting



Sample frequencies



Hilbert transform



Tilbert transform



Linear Algebra
==============
+.. currentmodule: scipy
+
When SciPy is built using the optimized ATLAS LAPACK and BLAS
libraries, it has very fast linear algebra capabilities. If you dig
deep enough, all of the raw lapack and blas libraries are available
@@ 1921,13 +1986,13 @@

The matrix class is initialized with the SciPy command :obj:`mat`
which is just convenient shorthand for :class:`matrix`. If
you are going to be doing a lot of matrixmath, it is convenient to
convert arrays into matrices using this command. One convencience of
using the mat command is that you can enter twodimensional matrices
using MATLABlike syntax with commas or spaces separating columns and
semicolons separting rows as long as the matrix is placed in a string
passed to :obj:`mat` .
+which is just convenient shorthand for :class:`matrix
+<numpy.matrix>`. If you are going to be doing a lot of matrixmath, it
+is convenient to convert arrays into matrices using this command. One
+advantage of using the :func:`mat` command is that you can enter
+twodimensional matrices using MATLABlike syntax with commas or
+spaces separating columns and semicolons separting rows as long as the
+matrix is placed in a string passed to :obj:`mat` .
Basic routines
@@ 1937,9 +2002,13 @@
Finding Inverse
^^^^^^^^^^^^^^^
The inverse of a matrix :math:`\mathbf{A}` is the matrix :math:`\mathbf{B}` such that :math:`\mathbf{AB}=\mathbf{I}` where :math:`\mathbf{I}` is the identity matrix consisting of ones down the main diagonal.
Usually :math:`\mathbf{B}` is denoted :math:`\mathbf{B}=\mathbf{A}^{1}` . In SciPy, the matrix inverse of the Numpy array, A, is obtained
using :obj:`linalg.inv` ``(A)`` , or using ``A.I`` if ``A`` is a Matrix. For example, let
+The inverse of a matrix :math:`\mathbf{A}` is the matrix
+:math:`\mathbf{B}` such that :math:`\mathbf{AB}=\mathbf{I}` where
+:math:`\mathbf{I}` is the identity matrix consisting of ones down the
+main diagonal. Usually :math:`\mathbf{B}` is denoted
+:math:`\mathbf{B}=\mathbf{A}^{1}` . In SciPy, the matrix inverse of
+the Numpy array, A, is obtained using :obj:`linalg.inv` ``(A)`` , or
+using ``A.I`` if ``A`` is a Matrix. For example, let
.. math::
:nowrap:
@@ 2011,7 +2080,13 @@
Finding Determinant
^^^^^^^^^^^^^^^^^^^
The determinant of a square matrix :math:`\mathbf{A}` is often denoted :math:`\left\mathbf{A}\right` and is a quantity often used in linear algebra. Suppose :math:`a_{ij}` are the elements of the matrix :math:`\mathbf{A}` and let :math:`M_{ij}=\left\mathbf{A}_{ij}\right` be the determinant of the matrix left by removing the :math:`i^{\textrm{th}}` row and :math:`j^{\textrm{th}}` column from :math:`\mathbf{A}` . Then for any row :math:`i,`
+The determinant of a square matrix :math:`\mathbf{A}` is often denoted
+:math:`\left\mathbf{A}\right` and is a quantity often used in linear
+algebra. Suppose :math:`a_{ij}` are the elements of the matrix
+:math:`\mathbf{A}` and let :math:`M_{ij}=\left\mathbf{A}_{ij}\right`
+be the determinant of the matrix left by removing the
+:math:`i^{\textrm{th}}` row and :math:`j^{\textrm{th}}` column from
+:math:`\mathbf{A}` . Then for any row :math:`i,`
.. math::
:nowrap:
@@ 2046,9 +2121,10 @@
Matrix and vector norms can also be computed with SciPy. A wide range
of norm definitions are available using different parameters to the
order argument of :obj:`linalg.norm` . This function takes a rank1 (vectors) or a rank2 (matrices) array
and an optional order argument (default is 2). Based on these inputs a
vector or matrix norm of the requested order is computed.
+order argument of :obj:`linalg.norm` . This function takes a rank1
+(vectors) or a rank2 (matrices) array and an optional order argument
+(default is 2). Based on these inputs a vector or matrix norm of the
+requested order is computed.
For vector *x* , the order parameter can be any real number including
``inf`` or ``inf``. The computed norm is
@@ 2076,15 +2152,18 @@
Linear leastsquares problems occur in many branches of applied
mathematics. In this problem a set of linear scaling coefficients is
sought that allow a model to fit data. In particular it is assumed
that data :math:`y_{i}` is related to data :math:`\mathbf{x}_{i}` through a set of coefficients :math:`c_{j}` and model functions :math:`f_{j}\left(\mathbf{x}_{i}\right)` via the model
+that data :math:`y_{i}` is related to data :math:`\mathbf{x}_{i}`
+through a set of coefficients :math:`c_{j}` and model functions
+:math:`f_{j}\left(\mathbf{x}_{i}\right)` via the model
.. math::
:nowrap:
\[ y_{i}=\sum_{j}c_{j}f_{j}\left(\mathbf{x}_{i}\right)+\epsilon_{i}\]
where :math:`\epsilon_{i}` represents uncertainty in the data. The strategy of least squares is
to pick the coefficients :math:`c_{j}` to minimize
+where :math:`\epsilon_{i}` represents uncertainty in the data. The
+strategy of least squares is to pick the coefficients :math:`c_{j}` to
+minimize
.. math::
:nowrap:
@@ 2121,25 +2200,35 @@
\[ \mathbf{c}=\left(\mathbf{A}^{H}\mathbf{A}\right)^{1}\mathbf{A}^{H}\mathbf{y}=\mathbf{A}^{\dagger}\mathbf{y}\]
where :math:`\mathbf{A}^{\dagger}` is called the pseudoinverse of :math:`\mathbf{A}.` Notice that using this definition of :math:`\mathbf{A}` the model can be written
+where :math:`\mathbf{A}^{\dagger}` is called the pseudoinverse of
+:math:`\mathbf{A}.` Notice that using this definition of
+:math:`\mathbf{A}` the model can be written
.. math::
:nowrap:
\[ \mathbf{y}=\mathbf{Ac}+\boldsymbol{\epsilon}.\]
The command :obj:`linalg.lstsq` will solve the linear least squares problem for :math:`\mathbf{c}` given :math:`\mathbf{A}` and :math:`\mathbf{y}` . In addition :obj:`linalg.pinv` or :obj:`linalg.pinv2` (uses a different method based on singular value decomposition) will
find :math:`\mathbf{A}^{\dagger}` given :math:`\mathbf{A}.`
+The command :obj:`linalg.lstsq` will solve the linear least squares
+problem for :math:`\mathbf{c}` given :math:`\mathbf{A}` and
+:math:`\mathbf{y}` . In addition :obj:`linalg.pinv` or
+:obj:`linalg.pinv2` (uses a different method based on singular value
+decomposition) will find :math:`\mathbf{A}^{\dagger}` given
+:math:`\mathbf{A}.`
The following example and figure demonstrate the use of :obj:`linalg.lstsq` and :obj:`linalg.pinv` for solving a datafitting problem. The data shown below were
generated using the model:
+The following example and figure demonstrate the use of
+:obj:`linalg.lstsq` and :obj:`linalg.pinv` for solving a datafitting
+problem. The data shown below were generated using the model:
.. math::
:nowrap:
\[ y_{i}=c_{1}e^{x_{i}}+c_{2}x_{i}\]
where :math:`x_{i}=0.1i` for :math:`i=1\ldots10` , :math:`c_{1}=5` , and :math:`c_{2}=4.` Noise is added to :math:`y_{i}` and the coefficients :math:`c_{1}` and :math:`c_{2}` are estimated using linear least squares.
+where :math:`x_{i}=0.1i` for :math:`i=1\ldots10` , :math:`c_{1}=5` ,
+and :math:`c_{2}=4.` Noise is added to :math:`y_{i}` and the
+coefficients :math:`c_{1}` and :math:`c_{2}` are estimated using
+linear least squares.
.. plot::
@@ 2170,9 +2259,12 @@
Generalized inverse
^^^^^^^^^^^^^^^^^^^
The generalized inverse is calculated using the command :obj:`linalg.pinv` or :obj:`linalg.pinv2`. These two commands differ in how they compute the generalized inverse.
The first uses the linalg.lstsq algorithm while the second uses
singular value decomposition. Let :math:`\mathbf{A}` be an :math:`M\times N` matrix, then if :math:`M>N` the generalized inverse is
+The generalized inverse is calculated using the command
+:obj:`linalg.pinv` or :obj:`linalg.pinv2`. These two commands differ
+in how they compute the generalized inverse. The first uses the
+linalg.lstsq algorithm while the second uses singular value
+decomposition. Let :math:`\mathbf{A}` be an :math:`M\times N` matrix,
+then if :math:`M>N` the generalized inverse is
.. math::
:nowrap:
@@ 2208,26 +2300,27 @@
The eigenvalueeigenvector problem is one of the most commonly
employed linear algebra operations. In one popular form, the
eigenvalueeigenvector problem is to find for some square matrix :math:`\mathbf{A}` scalars :math:`\lambda` and corresponding vectors :math:`\mathbf{v}` such that
+eigenvalueeigenvector problem is to find for some square matrix
+:math:`\mathbf{A}` scalars :math:`\lambda` and corresponding vectors
+:math:`\mathbf{v}` such that
.. math::
:nowrap:
\[ \mathbf{Av}=\lambda\mathbf{v}.\]
For an :math:`N\times N` matrix, there are :math:`N` (not necessarily distinct) eigenvalues  roots of the
(characteristic) polynomial
+For an :math:`N\times N` matrix, there are :math:`N` (not necessarily
+distinct) eigenvalues  roots of the (characteristic) polynomial
.. math::
:nowrap:
\[ \left\mathbf{A}\lambda\mathbf{I}\right=0.\]
+The eigenvectors, :math:`\mathbf{v}` , are also sometimes called right
+eigenvectors to distinguish them from another set of left eigenvectors
+that satisfy

The eigenvectors, :math:`\mathbf{v}` , are also sometimes called right eigenvectors to distinguish them
from another set of left eigenvectors that satisfy

.. math::
:nowrap:
@@ 2240,7 +2333,10 @@
\[ \mathbf{A}^{H}\mathbf{v}_{L}=\lambda^{*}\mathbf{v}_{L}.\]
With it's default optional arguments, the command :obj:`linalg.eig` returns :math:`\lambda` and :math:`\mathbf{v}.` However, it can also return :math:`\mathbf{v}_{L}` and just :math:`\lambda` by itself ( :obj:`linalg.eigvals` returns just :math:`\lambda` as well).
+With it's default optional arguments, the command :obj:`linalg.eig`
+returns :math:`\lambda` and :math:`\mathbf{v}.` However, it can also
+return :math:`\mathbf{v}_{L}` and just :math:`\lambda` by itself (
+:obj:`linalg.eigvals` returns just :math:`\lambda` as well).
In addtion, :obj:`linalg.eig` can also solve the more general eigenvalue problem
@@ 2249,20 +2345,25 @@
\begin{eqnarray*} \mathbf{Av} & = & \lambda\mathbf{Bv}\\ \mathbf{A}^{H}\mathbf{v}_{L} & = & \lambda^{*}\mathbf{B}^{H}\mathbf{v}_{L}\end{eqnarray*}
for square matrices :math:`\mathbf{A}` and :math:`\mathbf{B}.` The standard eigenvalue problem is an example of the general
eigenvalue problem for :math:`\mathbf{B}=\mathbf{I}.` When a generalized eigenvalue problem can be solved, then it provides
a decomposition of :math:`\mathbf{A}` as
+for square matrices :math:`\mathbf{A}` and :math:`\mathbf{B}.` The
+standard eigenvalue problem is an example of the general eigenvalue
+problem for :math:`\mathbf{B}=\mathbf{I}.` When a generalized
+eigenvalue problem can be solved, then it provides a decomposition of
+:math:`\mathbf{A}` as
.. math::
:nowrap:
\[ \mathbf{A}=\mathbf{BV}\boldsymbol{\Lambda}\mathbf{V}^{1}\]
where :math:`\mathbf{V}` is the collection of eigenvectors into columns and :math:`\boldsymbol{\Lambda}` is a diagonal matrix of eigenvalues.
+where :math:`\mathbf{V}` is the collection of eigenvectors into
+columns and :math:`\boldsymbol{\Lambda}` is a diagonal matrix of
+eigenvalues.
By definition, eigenvectors are only defined up to a constant scale
factor. In SciPy, the scaling factor for the eigenvectors is chosen so
that :math:`\left\Vert \mathbf{v}\right\Vert ^{2}=\sum_{i}v_{i}^{2}=1.`
+that :math:`\left\Vert \mathbf{v}\right\Vert
+^{2}=\sum_{i}v_{i}^{2}=1.`
As an example, consider finding the eigenvalues and eigenvectors of
the matrix
@@ 2315,16 +2416,39 @@
^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Singular Value Decompostion (SVD) can be thought of as an extension of
the eigenvalue problem to matrices that are not square. Let :math:`\mathbf{A}` be an :math:`M\times N` matrix with :math:`M` and :math:`N` arbitrary. The matrices :math:`\mathbf{A}^{H}\mathbf{A}` and :math:`\mathbf{A}\mathbf{A}^{H}` are square hermitian matrices [#]_ of size :math:`N\times N` and :math:`M\times M` respectively. It is known that the eigenvalues of square hermitian
matrices are real and nonnegative. In addtion, there are at most :math:`\min\left(M,N\right)` identical nonzero eigenvalues of :math:`\mathbf{A}^{H}\mathbf{A}` and :math:`\mathbf{A}\mathbf{A}^{H}.` Define these positive eigenvalues as :math:`\sigma_{i}^{2}.` The squareroot of these are called singular values of :math:`\mathbf{A}.` The eigenvectors of :math:`\mathbf{A}^{H}\mathbf{A}` are collected by columns into an :math:`N\times N` unitary [#]_ matrix :math:`\mathbf{V}` while the eigenvectors of :math:`\mathbf{A}\mathbf{A}^{H}` are collected by columns in the unitary matrix :math:`\mathbf{U}` , the singular values are collected in an :math:`M\times N` zero matrix :math:`\mathbf{\boldsymbol{\Sigma}}` with main diagonal entries set to the singular values. Then
+the eigenvalue problem to matrices that are not square. Let
+:math:`\mathbf{A}` be an :math:`M\times N` matrix with :math:`M` and
+:math:`N` arbitrary. The matrices :math:`\mathbf{A}^{H}\mathbf{A}` and
+:math:`\mathbf{A}\mathbf{A}^{H}` are square hermitian matrices [#]_ of
+size :math:`N\times N` and :math:`M\times M` respectively. It is known
+that the eigenvalues of square hermitian matrices are real and
+nonnegative. In addtion, there are at most
+:math:`\min\left(M,N\right)` identical nonzero eigenvalues of
+:math:`\mathbf{A}^{H}\mathbf{A}` and :math:`\mathbf{A}\mathbf{A}^{H}.`
+Define these positive eigenvalues as :math:`\sigma_{i}^{2}.` The
+squareroot of these are called singular values of :math:`\mathbf{A}.`
+The eigenvectors of :math:`\mathbf{A}^{H}\mathbf{A}` are collected by
+columns into an :math:`N\times N` unitary [#]_ matrix
+:math:`\mathbf{V}` while the eigenvectors of
+:math:`\mathbf{A}\mathbf{A}^{H}` are collected by columns in the
+unitary matrix :math:`\mathbf{U}` , the singular values are collected
+in an :math:`M\times N` zero matrix
+:math:`\mathbf{\boldsymbol{\Sigma}}` with main diagonal entries set to
+the singular values. Then
.. math::
:nowrap:
\[ \mathbf{A=U}\boldsymbol{\Sigma}\mathbf{V}^{H}\]
is the singularvalue decomposition of :math:`\mathbf{A}.` Every matrix has a singular value decomposition. Sometimes, the
singular values are called the spectrum of :math:`\mathbf{A}.` The command :obj:`linalg.svd` will return :math:`\mathbf{U}` , :math:`\mathbf{V}^{H}` , and :math:`\sigma_{i}` as an array of the singular values. To obtain the matrix :math:`\mathbf{\Sigma}` use :obj:`linalg.diagsvd`. The following example illustrates the use of :obj:`linalg.svd` .
+is the singularvalue decomposition of :math:`\mathbf{A}.` Every
+matrix has a singular value decomposition. Sometimes, the singular
+values are called the spectrum of :math:`\mathbf{A}.` The command
+:obj:`linalg.svd` will return :math:`\mathbf{U}` ,
+:math:`\mathbf{V}^{H}` , and :math:`\sigma_{i}` as an array of the
+singular values. To obtain the matrix :math:`\mathbf{\Sigma}` use
+:obj:`linalg.diagsvd`. The following example illustrates the use of
+:obj:`linalg.svd` .
>>> A = mat('[1 3 2; 1 2 3]')
>>> M,N = A.shape
@@ 2349,7 +2473,7 @@
[[ 1. 3. 2.]
[ 1. 2. 3.]]
.. [#] A hermition matrix :math:`\mathbf{D}` satisfies :math:`\mathbf{D}^{H}=\mathbf{D}.`
+.. [#] A hermitian matrix :math:`\mathbf{D}` satisfies :math:`\mathbf{D}^{H}=\mathbf{D}.`
.. [#] A unitary matrix :math:`\mathbf{D}` satisfies :math:`\mathbf{D}^{H}\mathbf{D}=\mathbf{I}=\mathbf{D}\mathbf{D}^{H}` so that :math:`\mathbf{D}^{1}=\mathbf{D}^{H}.`
@@ 2364,8 +2488,12 @@
\[ \mathbf{A}=\mathbf{PLU}\]
where :math:`\mathbf{P}` is an :math:`M\times M` permutation matrix (a permutation of the rows of the identity matrix), :math:`\mathbf{L}` is in :math:`M\times K` lower triangular or trapezoidal matrix ( :math:`K=\min\left(M,N\right)` ) with unitdiagonal, and :math:`\mathbf{U}` is an upper triangular or trapezoidal matrix. The SciPy command for
this decomposition is :obj:`linalg.lu` .
+where :math:`\mathbf{P}` is an :math:`M\times M` permutation matrix (a
+permutation of the rows of the identity matrix), :math:`\mathbf{L}` is
+in :math:`M\times K` lower triangular or trapezoidal matrix (
+:math:`K=\min\left(M,N\right)` ) with unitdiagonal, and
+:math:`\mathbf{U}` is an upper triangular or trapezoidal matrix. The
+SciPy command for this decomposition is :obj:`linalg.lu` .
Such a decomposition is often useful for solving many simultaneous
equations where the lefthandside does not change but the right hand
@@ 2383,32 +2511,48 @@
\[ \mathbf{PLUx}_{i}=\mathbf{b}_{i}.\]
Because :math:`\mathbf{L}` is lowertriangular, the equation can be solved for :math:`\mathbf{U}\mathbf{x}_{i}` and finally :math:`\mathbf{x}_{i}` very rapidly using forward and backsubstitution. An initial time
spent factoring :math:`\mathbf{A}` allows for very rapid solution of similar systems of equations in the
+Because :math:`\mathbf{L}` is lowertriangular, the equation can be
+solved for :math:`\mathbf{U}\mathbf{x}_{i}` and finally
+:math:`\mathbf{x}_{i}` very rapidly using forward and
+backsubstitution. An initial time spent factoring :math:`\mathbf{A}`
+allows for very rapid solution of similar systems of equations in the
future. If the intent for performing LU decomposition is for solving
linear systems then the command :obj:`linalg.lu_factor` should be used followed by repeated applications of the command :obj:`linalg.lu_solve` to solve the system for each new righthandside.
+linear systems then the command :obj:`linalg.lu_factor` should be used
+followed by repeated applications of the command
+:obj:`linalg.lu_solve` to solve the system for each new
+righthandside.
Cholesky decomposition
^^^^^^^^^^^^^^^^^^^^^^
Cholesky decomposition is a special case of LU decomposition
applicable to Hermitian positive definite matrices. When :math:`\mathbf{A}=\mathbf{A}^{H}` and :math:`\mathbf{x}^{H}\mathbf{Ax}\geq0` for all :math:`\mathbf{x}` , then decompositions of :math:`\mathbf{A}` can be found so that
+applicable to Hermitian positive definite matrices. When
+:math:`\mathbf{A}=\mathbf{A}^{H}` and
+:math:`\mathbf{x}^{H}\mathbf{Ax}\geq0` for all :math:`\mathbf{x}` ,
+then decompositions of :math:`\mathbf{A}` can be found so that
.. math::
:nowrap:
\begin{eqnarray*} \mathbf{A} & = & \mathbf{U}^{H}\mathbf{U}\\ \mathbf{A} & = & \mathbf{L}\mathbf{L}^{H}\end{eqnarray*}
where :math:`\mathbf{L}` is lowertriangular and :math:`\mathbf{U}` is upper triangular. Notice that :math:`\mathbf{L}=\mathbf{U}^{H}.` The command :obj:`linagl.cholesky` computes the cholesky factorization. For using cholesky factorization
to solve systems of equations there are also :obj:`linalg.cho_factor` and :obj:`linalg.cho_solve` routines that work similarly to their LU decomposition counterparts.
+where :math:`\mathbf{L}` is lowertriangular and :math:`\mathbf{U}` is
+upper triangular. Notice that :math:`\mathbf{L}=\mathbf{U}^{H}.` The
+command :obj:`linagl.cholesky` computes the cholesky
+factorization. For using cholesky factorization to solve systems of
+equations there are also :obj:`linalg.cho_factor` and
+:obj:`linalg.cho_solve` routines that work similarly to their LU
+decomposition counterparts.
QR decomposition
^^^^^^^^^^^^^^^^
The QR decomposition (sometimes called a polar decomposition) works
for any :math:`M\times N` array and finds an :math:`M\times M` unitary matrix :math:`\mathbf{Q}` and an :math:`M\times N` uppertrapezoidal matrix :math:`\mathbf{R}` such that
+for any :math:`M\times N` array and finds an :math:`M\times M` unitary
+matrix :math:`\mathbf{Q}` and an :math:`M\times N` uppertrapezoidal
+matrix :math:`\mathbf{R}` such that
.. math::
:nowrap:
@@ 2422,26 +2566,37 @@
\[ \mathbf{A}=\mathbf{U}\boldsymbol{\Sigma}\mathbf{V}^{H}=\mathbf{QR}\]
implies that :math:`\mathbf{Q}=\mathbf{U}` and :math:`\mathbf{R}=\boldsymbol{\Sigma}\mathbf{V}^{H}.` Note, however, that in SciPy independent algorithms are used to find
QR and SVD decompositions. The command for QR decomposition is :obj:`linalg.qr` .
+implies that :math:`\mathbf{Q}=\mathbf{U}` and
+:math:`\mathbf{R}=\boldsymbol{\Sigma}\mathbf{V}^{H}.` Note, however,
+that in SciPy independent algorithms are used to find QR and SVD
+decompositions. The command for QR decomposition is :obj:`linalg.qr` .
Schur decomposition
^^^^^^^^^^^^^^^^^^^
For a square :math:`N\times N` matrix, :math:`\mathbf{A}` , the Schur decomposition finds (notnecessarily unique) matrices :math:`\mathbf{T}` and :math:`\mathbf{Z}` such that
+For a square :math:`N\times N` matrix, :math:`\mathbf{A}` , the Schur
+decomposition finds (notnecessarily unique) matrices
+:math:`\mathbf{T}` and :math:`\mathbf{Z}` such that
.. math::
:nowrap:
\[ \mathbf{A}=\mathbf{ZT}\mathbf{Z}^{H}\]
where :math:`\mathbf{Z}` is a unitary matrix and :math:`\mathbf{T}` is either uppertriangular or quasiupper triangular depending on
whether or not a real schur form or complex schur form is requested.
For a real schur form both :math:`\mathbf{T}` and :math:`\mathbf{Z}` are realvalued when :math:`\mathbf{A}` is realvalued. When :math:`\mathbf{A}` is a realvalued matrix the real schur form is only quasiupper
triangular because :math:`2\times2` blocks extrude from the main diagonal corresponding to any complex
valued eigenvalues. The command :obj:`linalg.schur` finds the Schur decomposition while the command :obj:`linalg.rsf2csf` converts :math:`\mathbf{T}` and :math:`\mathbf{Z}` from a real Schur form to a complex Schur form. The Schur form is
especially useful in calculating functions of matrices.
+where :math:`\mathbf{Z}` is a unitary matrix and :math:`\mathbf{T}` is
+either uppertriangular or quasiupper triangular depending on whether
+or not a real schur form or complex schur form is requested. For a
+real schur form both :math:`\mathbf{T}` and :math:`\mathbf{Z}` are
+realvalued when :math:`\mathbf{A}` is realvalued. When
+:math:`\mathbf{A}` is a realvalued matrix the real schur form is only
+quasiupper triangular because :math:`2\times2` blocks extrude from
+the main diagonal corresponding to any complex valued
+eigenvalues. The command :obj:`linalg.schur` finds the Schur
+decomposition while the command :obj:`linalg.rsf2csf` converts
+:math:`\mathbf{T}` and :math:`\mathbf{Z}` from a real Schur form to a
+complex Schur form. The Schur form is especially useful in calculating
+functions of matrices.
The following example illustrates the schur decomposition:
@@ 2537,7 +2692,8 @@
where the matrix exponential of the diagonal matrix :math:`\boldsymbol{\Lambda}` is just the exponential of its elements. This method is implemented in :obj:`linalg.expm2` .
The preferred method for implementing the matrix exponential is to use
scaling and a Padé approximation for :math:`e^{x}` . This algorithm is implemented as :obj:`linalg.expm` .
+scaling and a Padé approximation for :math:`e^{x}` . This algorithm is
+implemented as :obj:`linalg.expm` .
The inverse of the matrix exponential is the matrix logarithm defined
as the inverse of the matrix exponential.
@@ 2583,14 +2739,17 @@
Hyperbolic trigonometric functions
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
The hyperbolic trigonemetric functions :math:`\sinh` , :math:`\cosh` , and :math:`\tanh` can also be defined for matrices using the familiar definitions:
+The hyperbolic trigonemetric functions :math:`\sinh` , :math:`\cosh` ,
+and :math:`\tanh` can also be defined for matrices using the familiar
+definitions:
.. math::
:nowrap:
\begin{eqnarray*} \sinh\left(\mathbf{A}\right) & = & \frac{e^{\mathbf{A}}e^{\mathbf{A}}}{2}\\ \cosh\left(\mathbf{A}\right) & = & \frac{e^{\mathbf{A}}+e^{\mathbf{A}}}{2}\\ \tanh\left(\mathbf{A}\right) & = & \left[\cosh\left(\mathbf{A}\right)\right]^{1}\sinh\left(\mathbf{A}\right).\end{eqnarray*}
These matrix functions can be found using :obj:`linalg.sinhm` , :obj:`linalg.coshm` , and :obj:`linalg.tanhm` .
+These matrix functions can be found using :obj:`linalg.sinhm`,
+:obj:`linalg.coshm` , and :obj:`linalg.tanhm`.
Arbitrary function
@@ 2608,7 +2767,7 @@
>>> from scipy import special, random, linalg
>>> A = random.rand(3,3)
 >>> B = linalg.funm(A,lambda x: special.jv(0,real(x)))
+ >>> B = linalg.funm(A,lambda x: special.jv(0,x))
>>> print A
[[ 0.72578091 0.34105276 0.79570345]
[ 0.65767207 0.73855618 0.541453 ]
@@ 2624,7 +2783,10 @@
>>> print linalg.eigvals(B)
[ 0.27448286+0.j 0.98810383+0.j 0.99164854+0.j]
+Note how, by virtue of how matrix analytic functions are defined,
+the Bessel function has acted on the matrix eigenvalues.
+
Statistics
==========
@@ 2647,32 +2809,34 @@
files continuous.lyx and discrete.lyx in the stats subdirectories.
Interfacing with Python Imaging Library
=======================================

If you have the Python Imaging Library installed, SciPy provides some
convenient functions that make use of it's facilities particularly for
reading, writing, displaying, and rotating images. In SciPy an image
is always a two or threedimensional array. Grayscale, and colormap
images are always twodimensional arrays while RGB images are three
dimensional with the third dimension specifying the channel.

Commands available include

 fromimage  convert a PIL image to a Numpy array

 toimage  convert Numpy array to PIL image

 imsave  save Numpy array to an image file

 imread  read an image file into a Numpy array

 imrotate  rotate an image (a Numpy array) counterclockwise

 imresize  resize an image using the PIL

 imshow  external viewer of a Numpy array (using PIL)

 imfilter  fast, simple builtin filters provided by PIL

 radon  simple radon transform based on imrotate
+.. XXX: is this uptodate?
+..
+.. Interfacing with Python Imaging Library
+.. =======================================
+..
+.. If you have the Python Imaging Library installed, SciPy provides some
+.. convenient functions that make use of it's facilities particularly for
+.. reading, writing, displaying, and rotating images. In SciPy an image
+.. is always a two or threedimensional array. Grayscale, and colormap
+.. images are always twodimensional arrays while RGB images are three
+.. dimensional with the third dimension specifying the channel.
+..
+.. Commands available include
+..
+..  fromimage  convert a PIL image to a Numpy array
+..
+..  toimage  convert Numpy array to PIL image
+..
+..  imsave  save Numpy array to an image file
+..
+..  imread  read an image file into a Numpy array
+..
+..  imrotate  rotate an image (a Numpy array) counterclockwise
+..
+..  imresize  resize an image using the PIL
+..
+..  imshow  external viewer of a Numpy array (using PIL)
+..
+..  imfilter  fast, simple builtin filters provided by PIL
+..
+..  radon  simple radon transform based on imrotate
More information about the Scipysvn
mailing list