Relationship between SciPy and NumPy
From the SciPy Reference Guide:
... all of the Numpy functions have been subsumed into the
scipynamespace so that all of those functions are available without additionally importing Numpy.
The intention is for users not to have to know the distinction between the scipy and numpy namespaces, though apparently you've found an exception.
There is a short comment at the end of the introduction to SciPy documentation:
Another useful command is
source. When given a function written in Python as an argument, it prints out a listing of the source code for that function. This can be helpful in learning about an algorithm or understanding exactly what a function is doing with its arguments. Also don’t forget about the Python command dir which can be used to look at the namespace of a module or package.
I think this will allow someone with enough knowledge of all the packages involved to pick apart exactly what the differences are between some scipy and numpy functions (it didn't help me with the log10 question at all). I definitely don't have that knowledge but source does indicate that scipy.linalg.solve and numpy.linalg.solve interact with lapack in different ways;
Python 2.4.3 (#1, May 5 2011, 18:44:23)
[GCC 4.1.2 20080704 (Red Hat 4.1.2-50)] on linux2
>>> import scipy
>>> import scipy.linalg
>>> import numpy
>>> scipy.source(scipy.linalg.solve)
In file: /usr/lib64/python2.4/site-packages/scipy/linalg/basic.py
def solve(a, b, sym_pos=0, lower=0, overwrite_a=0, overwrite_b=0,
debug = 0):
""" solve(a, b, sym_pos=0, lower=0, overwrite_a=0, overwrite_b=0) -> x
Solve a linear system of equations a * x = b for x.
Inputs:
a -- An N x N matrix.
b -- An N x nrhs matrix or N vector.
sym_pos -- Assume a is symmetric and positive definite.
lower -- Assume a is lower triangular, otherwise upper one.
Only used if sym_pos is true.
overwrite_y - Discard data in y, where y is a or b.
Outputs:
x -- The solution to the system a * x = b
"""
a1, b1 = map(asarray_chkfinite,(a,b))
if len(a1.shape) != 2 or a1.shape[0] != a1.shape[1]:
raise ValueError, 'expected square matrix'
if a1.shape[0] != b1.shape[0]:
raise ValueError, 'incompatible dimensions'
overwrite_a = overwrite_a or (a1 is not a and not hasattr(a,'__array__'))
overwrite_b = overwrite_b or (b1 is not b and not hasattr(b,'__array__'))
if debug:
print 'solve:overwrite_a=',overwrite_a
print 'solve:overwrite_b=',overwrite_b
if sym_pos:
posv, = get_lapack_funcs(('posv',),(a1,b1))
c,x,info = posv(a1,b1,
lower = lower,
overwrite_a=overwrite_a,
overwrite_b=overwrite_b)
else:
gesv, = get_lapack_funcs(('gesv',),(a1,b1))
lu,piv,x,info = gesv(a1,b1,
overwrite_a=overwrite_a,
overwrite_b=overwrite_b)
if info==0:
return x
if info>0:
raise LinAlgError, "singular matrix"
raise ValueError,\
'illegal value in %-th argument of internal gesv|posv'%(-info)
>>> scipy.source(numpy.linalg.solve)
In file: /usr/lib64/python2.4/site-packages/numpy/linalg/linalg.py
def solve(a, b):
"""
Solve the equation ``a x = b`` for ``x``.
Parameters
----------
a : array_like, shape (M, M)
Input equation coefficients.
b : array_like, shape (M,)
Equation target values.
Returns
-------
x : array, shape (M,)
Raises
------
LinAlgError
If `a` is singular or not square.
Examples
--------
Solve the system of equations ``3 * x0 + x1 = 9`` and ``x0 + 2 * x1 = 8``:
>>> a = np.array([[3,1], [1,2]])
>>> b = np.array([9,8])
>>> x = np.linalg.solve(a, b)
>>> x
array([ 2., 3.])
Check that the solution is correct:
>>> (np.dot(a, x) == b).all()
True
"""
a, _ = _makearray(a)
b, wrap = _makearray(b)
one_eq = len(b.shape) == 1
if one_eq:
b = b[:, newaxis]
_assertRank2(a, b)
_assertSquareness(a)
n_eq = a.shape[0]
n_rhs = b.shape[1]
if n_eq != b.shape[0]:
raise LinAlgError, 'Incompatible dimensions'
t, result_t = _commonType(a, b)
# lapack_routine = _findLapackRoutine('gesv', t)
if isComplexType(t):
lapack_routine = lapack_lite.zgesv
else:
lapack_routine = lapack_lite.dgesv
a, b = _fastCopyAndTranspose(t, a, b)
pivots = zeros(n_eq, fortran_int)
results = lapack_routine(n_eq, n_rhs, a, n_eq, pivots, b, n_eq, 0)
if results['info'] > 0:
raise LinAlgError, 'Singular matrix'
if one_eq:
return wrap(b.ravel().astype(result_t))
else:
return wrap(b.transpose().astype(result_t))
This is also my first post so if I should change something here please let me know.
Last time I checked it, the scipy __init__ method executes a
from numpy import *
so that the whole numpy namespace is included into scipy when the scipy module is imported.
The log10 behavior you are describing is interesting, because both versions are coming from numpy. One is a ufunc, the other is a numpy.lib function. Why scipy is preferring the library function over the ufunc, I don't know off the top of my head.
EDIT: In fact, I can answer the log10 question. Looking in the scipy __init__ method I see this:
# Import numpy symbols to scipy name space
import numpy as _num
from numpy import oldnumeric
from numpy import *
from numpy.random import rand, randn
from numpy.fft import fft, ifft
from numpy.lib.scimath import *
The log10 function you get in scipy comes from numpy.lib.scimath. Looking at that code, it says:
"""
Wrapper functions to more user-friendly calling of certain math functions
whose output data-type is different than the input data-type in certain
domains of the input.
For example, for functions like log() with branch cuts, the versions in this
module provide the mathematically valid answers in the complex plane:
>>> import math
>>> from numpy.lib import scimath
>>> scimath.log(-math.exp(1)) == (1+1j*math.pi)
True
Similarly, sqrt(), other base logarithms, power() and trig functions are
correctly handled. See their respective docstrings for specific examples.
"""
It seems that module overlays the base numpy ufuncs for sqrt, log, log2, logn, log10, power, arccos, arcsin, and arctanh. That explains the behavior you are seeing. The underlying design reason why it is done like that is probably buried in a mailing list post somewhere.
It seems from the SciPy FAQ that some functions from NumPy are here for historical reasons while it should only be in SciPy:
What is the difference between NumPy and SciPy?
In an ideal world, NumPy would contain nothing but the array data type and the most basic operations: indexing, sorting, reshaping, basic elementwise functions, et cetera. All numerical code would reside in SciPy. However, one of NumPy’s important goals is compatibility, so NumPy tries to retain all features supported by either of its predecessors. Thus NumPy contains some linear algebra functions, even though these more properly belong in SciPy. In any case, SciPy contains more fully-featured versions of the linear algebra modules, as well as many other numerical algorithms. If you are doing scientific computing with python, you should probably install both NumPy and SciPy. Most new features belong in SciPy rather than NumPy.
That explains why scipy.linalg.solve offers some additional features over numpy.linalg.solve.
I did not see the answer of SethMMorton to the related question