[Scipy-svn] r4630 - in trunk/scipy/stats/models: . src tests
scipy-svn@scip...
scipy-svn@scip...
Fri Aug 8 16:35:32 CDT 2008
Author: chris.burns
Date: 2008-08-08 16:35:24 -0500 (Fri, 08 Aug 2008)
New Revision: 4630
Added:
trunk/scipy/stats/models/TODO.txt
trunk/scipy/stats/models/src/
trunk/scipy/stats/models/src/bspline_ext.c
trunk/scipy/stats/models/src/bspline_impl.c
Removed:
trunk/scipy/stats/models/_bspline.py
trunk/scipy/stats/models/bspline_module.py
trunk/scipy/stats/models/src/bspline_ext.c
trunk/scipy/stats/models/src/bspline_impl.c
Modified:
trunk/scipy/stats/models/bspline.py
trunk/scipy/stats/models/setup.py
trunk/scipy/stats/models/smoothers.py
trunk/scipy/stats/models/tests/test_bspline.py
trunk/scipy/stats/models/tests/test_formula.py
trunk/scipy/stats/models/tests/test_scale.py
Log:
Merge branch converting models.bspline weave code to a c extension. Suppress known test failures.
Copied: trunk/scipy/stats/models/TODO.txt (from rev 4629, branches/stats_models/TODO.txt)
Deleted: trunk/scipy/stats/models/_bspline.py
===================================================================
--- trunk/scipy/stats/models/_bspline.py 2008-08-08 21:10:46 UTC (rev 4629)
+++ trunk/scipy/stats/models/_bspline.py 2008-08-08 21:35:24 UTC (rev 4630)
@@ -1,657 +0,0 @@
-'''
-Bspines and smoothing splines.
-
-General references:
-
- Craven, P. and Wahba, G. (1978) "Smoothing noisy data with spline functions.
- Estimating the correct degree of smoothing by
- the method of generalized cross-validation."
- Numerische Mathematik, 31(4), 377-403.
-
- Hastie, Tibshirani and Friedman (2001). "The Elements of Statistical
- Learning." Springer-Verlag. 536 pages.
-
- Hutchison, M. and Hoog, F. "Smoothing noisy data with spline functions."
- Numerische Mathematik, 47(1), 99-106.
-'''
-
-import numpy as N
-import numpy.linalg as L
-
-from scipy.linalg import solveh_banded
-from scipy.optimize import golden
-from scipy.stats.models import _bspline
-
-
-# Issue warning regarding heavy development status of this module
-import warnings
-_msg = "The bspline code is technology preview and requires significant work\
-on the public API and documentation. The API will likely change in the future"
-warnings.warn(_msg, UserWarning)
-
-
-def _band2array(a, lower=0, symmetric=False, hermitian=False):
- """
- Take an upper or lower triangular banded matrix and return a
- numpy array.
-
- INPUTS:
- a -- a matrix in upper or lower triangular banded matrix
- lower -- is the matrix upper or lower triangular?
- symmetric -- if True, return the original result plus its transpose
- hermitian -- if True (and symmetric False), return the original
- result plus its conjugate transposed
-
- """
-
- n = a.shape[1]
- r = a.shape[0]
- _a = 0
-
- if not lower:
- for j in range(r):
- _b = N.diag(a[r-1-j],k=j)[j:(n+j),j:(n+j)]
- _a += _b
- if symmetric and j > 0: _a += _b.T
- elif hermitian and j > 0: _a += _b.conjugate().T
- else:
- for j in range(r):
- _b = N.diag(a[j],k=j)[0:n,0:n]
- _a += _b
- if symmetric and j > 0: _a += _b.T
- elif hermitian and j > 0: _a += _b.conjugate().T
- _a = _a.T
-
- return _a
-
-
-def _upper2lower(ub):
- """
- Convert upper triangular banded matrix to lower banded form.
-
- INPUTS:
- ub -- an upper triangular banded matrix
-
- OUTPUTS: lb
- lb -- a lower triangular banded matrix with same entries
- as ub
- """
-
- lb = N.zeros(ub.shape, ub.dtype)
- nrow, ncol = ub.shape
- for i in range(ub.shape[0]):
- lb[i,0:(ncol-i)] = ub[nrow-1-i,i:ncol]
- lb[i,(ncol-i):] = ub[nrow-1-i,0:i]
- return lb
-
-def _lower2upper(lb):
- """
- Convert lower triangular banded matrix to upper banded form.
-
- INPUTS:
- lb -- a lower triangular banded matrix
-
- OUTPUTS: ub
- ub -- an upper triangular banded matrix with same entries
- as lb
- """
-
- ub = N.zeros(lb.shape, lb.dtype)
- nrow, ncol = lb.shape
- for i in range(lb.shape[0]):
- ub[nrow-1-i,i:ncol] = lb[i,0:(ncol-i)]
- ub[nrow-1-i,0:i] = lb[i,(ncol-i):]
- return ub
-
-def _triangle2unit(tb, lower=0):
- """
- Take a banded triangular matrix and return its diagonal and the
- unit matrix: the banded triangular matrix with 1's on the diagonal,
- i.e. each row is divided by the corresponding entry on the diagonal.
-
- INPUTS:
- tb -- a lower triangular banded matrix
- lower -- if True, then tb is assumed to be lower triangular banded,
- in which case return value is also lower triangular banded.
-
- OUTPUTS: d, b
- d -- diagonal entries of tb
- b -- unit matrix: if lower is False, b is upper triangular
- banded and its rows of have been divided by d,
- else lower is True, b is lower triangular banded
- and its columns have been divieed by d.
-
- """
-
- if lower: d = tb[0].copy()
- else: d = tb[-1].copy()
-
- if lower: return d, (tb / d)
- else:
- l = _upper2lower(tb)
- return d, _lower2upper(l / d)
-
-def _trace_symbanded(a, b, lower=0):
- """
- Compute the trace(ab) for two upper or banded real symmetric matrices
- stored either in either upper or lower form.
-
- INPUTS:
- a, b -- two banded real symmetric matrices (either lower or upper)
- lower -- if True, a and b are assumed to be the lower half
-
-
- OUTPUTS: trace
- trace -- trace(ab)
-
- """
-
- if lower:
- t = _zero_triband(a * b, lower=1)
- return t[0].sum() + 2 * t[1:].sum()
- else:
- t = _zero_triband(a * b, lower=0)
- return t[-1].sum() + 2 * t[:-1].sum()
-
-
-def _zero_triband(a, lower=0):
- """
- Explicitly zero out unused elements of a real symmetric banded matrix.
-
- INPUTS:
- a -- a real symmetric banded matrix (either upper or lower hald)
- lower -- if True, a is assumed to be the lower half
-
- """
-
- nrow, ncol = a.shape
- if lower:
- for i in range(nrow): a[i,(ncol-i):] = 0.
- else:
- for i in range(nrow): a[i,0:i] = 0.
- return a
-
-
-class BSpline(object):
-
- '''
-
- Bsplines of a given order and specified knots.
-
- Implementation is based on description in Chapter 5 of
-
- Hastie, Tibshirani and Friedman (2001). "The Elements of Statistical
- Learning." Springer-Verlag. 536 pages.
-
-
- INPUTS:
- knots -- a sorted array of knots with knots[0] the lower boundary,
- knots[1] the upper boundary and knots[1:-1] the internal
- knots.
- order -- order of the Bspline, default is 4 which yields cubic
- splines
- M -- number of additional boundary knots, if None it defaults
- to order
- coef -- an optional array of real-valued coefficients for the Bspline
- of shape (knots.shape + 2 * (M - 1) - order,).
- x -- an optional set of x values at which to evaluate the
- Bspline to avoid extra evaluation in the __call__ method
-
- '''
-
- def __init__(self, knots, order=4, M=None, coef=None, x=None):
-
- knots = N.squeeze(N.unique(N.asarray(knots)))
-
- if knots.ndim != 1:
- raise ValueError, 'expecting 1d array for knots'
-
- self.m = order
- if M is None:
- M = self.m
- self.M = M
-
- self.tau = N.hstack([[knots[0]]*(self.M-1), knots, [knots[-1]]*(self.M-1)])
-
- self.K = knots.shape[0] - 2
- if coef is None:
- self.coef = N.zeros((self.K + 2 * self.M - self.m), N.float64)
- else:
- self.coef = N.squeeze(coef)
- if self.coef.shape != (self.K + 2 * self.M - self.m):
- raise ValueError, 'coefficients of Bspline have incorrect shape'
- if x is not None:
- self.x = x
-
- def _setx(self, x):
- self._x = x
- self._basisx = self.basis(self._x)
-
- def _getx(self):
- return self._x
-
- x = property(_getx, _setx)
-
- def __call__(self, *args):
- """
- Evaluate the BSpline at a given point, yielding
- a matrix B and return
-
- B * self.coef
-
-
- INPUTS:
- args -- optional arguments. If None, it returns self._basisx,
- the BSpline evaluated at the x values passed in __init__.
- Otherwise, return the BSpline evaluated at the
- first argument args[0].
-
- OUTPUTS: y
- y -- value of Bspline at specified x values
-
- BUGS:
- If self has no attribute x, an exception will be raised
- because self has no attribute _basisx.
-
- """
-
- if not args:
- b = self._basisx.T
- else:
- x = args[0]
- b = N.asarray(self.basis(x)).T
- return N.squeeze(N.dot(b, self.coef))
-
- def basis_element(self, x, i, d=0):
- """
- Evaluate a particular basis element of the BSpline,
- or its derivative.
-
- INPUTS:
- x -- x values at which to evaluate the basis element
- i -- which element of the BSpline to return
- d -- the order of derivative
-
- OUTPUTS: y
- y -- value of d-th derivative of the i-th basis element
- of the BSpline at specified x values
-
- """
-
- x = N.asarray(x, N.float64)
- _shape = x.shape
- if _shape == ():
- x.shape = (1,)
- x.shape = (N.product(_shape,axis=0),)
- if i < self.tau.shape[0] - 1:
- ## TODO: OWNDATA flags...
- v = _bspline.evaluate(x, self.tau, self.m, d, i, i+1)
- else:
- return N.zeros(x.shape, N.float64)
-
- if (i == self.tau.shape[0] - self.m):
- v = N.where(N.equal(x, self.tau[-1]), 1, v)
- v.shape = _shape
- return v
-
- def basis(self, x, d=0, lower=None, upper=None):
- """
- Evaluate the basis of the BSpline or its derivative.
- If lower or upper is specified, then only
- the [lower:upper] elements of the basis are returned.
-
- INPUTS:
- x -- x values at which to evaluate the basis element
- i -- which element of the BSpline to return
- d -- the order of derivative
- lower -- optional lower limit of the set of basis
- elements
- upper -- optional upper limit of the set of basis
- elements
-
- OUTPUTS: y
- y -- value of d-th derivative of the basis elements
- of the BSpline at specified x values
-
- """
- x = N.asarray(x)
- _shape = x.shape
- if _shape == ():
- x.shape = (1,)
- x.shape = (N.product(_shape,axis=0),)
-
- if upper is None:
- upper = self.tau.shape[0] - self.m
- if lower is None:
- lower = 0
- upper = min(upper, self.tau.shape[0] - self.m)
- lower = max(0, lower)
-
- d = N.asarray(d)
- if d.shape == ():
- v = _bspline.evaluate(x, self.tau, self.m, int(d), lower, upper)
- else:
- if d.shape[0] != 2:
- raise ValueError, "if d is not an integer, expecting a jx2 \
- array with first row indicating order \
- of derivative, second row coefficient in front."
-
- v = 0
- for i in range(d.shape[1]):
- v += d[1,i] * _bspline.evaluate(x, self.tau, self.m, d[0,i], lower, upper)
-
- v.shape = (upper-lower,) + _shape
- if upper == self.tau.shape[0] - self.m:
- v[-1] = N.where(N.equal(x, self.tau[-1]), 1, v[-1])
- return v
-
- def gram(self, d=0):
- """
- Compute Gram inner product matrix, storing it in lower
- triangular banded form.
-
- The (i,j) entry is
-
- G_ij = integral b_i^(d) b_j^(d)
-
- where b_i are the basis elements of the BSpline and (d) is the
- d-th derivative.
-
- If d is a matrix then, it is assumed to specify a differential
- operator as follows: the first row represents the order of derivative
- with the second row the coefficient corresponding to that order.
-
- For instance:
-
- [[2, 3],
- [3, 1]]
-
- represents 3 * f^(2) + 1 * f^(3).
-
- INPUTS:
- d -- which derivative to apply to each basis element,
- if d is a matrix, it is assumed to specify
- a differential operator as above
-
- OUTPUTS: gram
- gram -- the matrix of inner products of (derivatives)
- of the BSpline elements
-
- """
-
- d = N.squeeze(d)
- if N.asarray(d).shape == ():
- self.g = _bspline.gram(self.tau, self.m, int(d), int(d))
- else:
- d = N.asarray(d)
- if d.shape[0] != 2:
- raise ValueError, "if d is not an integer, expecting a jx2 \
- array with first row indicating order \
- of derivative, second row coefficient in front."
- if d.shape == (2,):
- d.shape = (2,1)
- self.g = 0
- for i in range(d.shape[1]):
- for j in range(d.shape[1]):
- self.g += d[1,i]* d[1,j] * _bspline.gram(self.tau, self.m, int(d[0,i]), int(d[0,j]))
- self.g = self.g.T
- self.d = d
- return N.nan_to_num(self.g)
-
-class SmoothingSpline(BSpline):
-
- penmax = 30.
- method = "target_df"
- target_df = 5
- default_pen = 1.0e-03
- optimize = True
-
- '''
- A smoothing spline, which can be used to smooth scatterplots, i.e.
- a list of (x,y) tuples.
-
- See fit method for more information.
-
- '''
-
- def fit(self, y, x=None, weights=None, pen=0.):
- """
- Fit the smoothing spline to a set of (x,y) pairs.
-
- INPUTS:
- y -- response variable
- x -- if None, uses self.x
- weights -- optional array of weights
- pen -- constant in front of Gram matrix
-
- OUTPUTS: None
- The smoothing spline is determined by self.coef,
- subsequent calls of __call__ will be the smoothing spline.
-
- ALGORITHM:
- Formally, this solves a minimization:
-
- fhat = ARGMIN_f SUM_i=1^n (y_i-f(x_i))^2 + pen * int f^(2)^2
-
- See Chapter 5 of
-
- Hastie, Tibshirani and Friedman (2001). "The Elements of Statistical
- Learning." Springer-Verlag. 536 pages.
-
- for more details.
-
- TODO:
- Should add arbitrary derivative penalty instead of just
- second derivative.
- """
-
- banded = True
-
- if x is None:
- x = self._x
- bt = self._basisx.copy()
- else:
- bt = self.basis(x)
-
- if pen == 0.: # can't use cholesky for singular matrices
- banded = False
-
- if x.shape != y.shape:
- raise ValueError, 'x and y shape do not agree, by default x are \
- the Bspline\'s internal knots'
-
- if pen >= self.penmax:
- pen = self.penmax
-
-
- if weights is not None:
- self.weights = weights
- else:
- self.weights = 1.
-
- _w = N.sqrt(self.weights)
- bt *= _w
-
- # throw out rows with zeros (this happens at boundary points!)
-
- mask = N.flatnonzero(1 - N.alltrue(N.equal(bt, 0), axis=0))
-
- bt = bt[:,mask]
- y = y[mask]
-
- self.df_total = y.shape[0]
-
- bty = N.squeeze(N.dot(bt, _w * y))
- self.N = y.shape[0]
-
- if not banded:
- self.btb = N.dot(bt, bt.T)
- _g = _band2array(self.g, lower=1, symmetric=True)
- self.coef, _, self.rank = L.lstsq(self.btb + pen*_g, bty)[0:3]
- self.rank = min(self.rank, self.btb.shape[0])
- del(_g)
- else:
- self.btb = N.zeros(self.g.shape, N.float64)
- nband, nbasis = self.g.shape
- for i in range(nbasis):
- for k in range(min(nband, nbasis-i)):
- self.btb[k,i] = (bt[i] * bt[i+k]).sum()
-
- bty.shape = (1,bty.shape[0])
- self.pen = pen
- self.chol, self.coef = solveh_banded(self.btb +
- pen*self.g,
- bty, lower=1)
-
- self.coef = N.squeeze(self.coef)
- self.resid = y * self.weights - N.dot(self.coef, bt)
- self.pen = pen
-
- del(bty); del(mask); del(bt)
-
- def smooth(self, y, x=None, weights=None):
-
- if self.method == "target_df":
- if hasattr(self, 'pen'):
- self.fit(y, x=x, weights=weights, pen=self.pen)
- else:
- self.fit_target_df(y, x=x, weights=weights, df=self.target_df)
- elif self.method == "optimize_gcv":
- self.fit_optimize_gcv(y, x=x, weights=weights)
-
-
- def gcv(self):
- """
- Generalized cross-validation score of current fit.
-
- Craven, P. and Wahba, G. "Smoothing noisy data with spline functions.
- Estimating the correct degree of smoothing by
- the method of generalized cross-validation."
- Numerische Mathematik, 31(4), 377-403.
- """
-
- norm_resid = (self.resid**2).sum()
- return norm_resid / (self.df_total - self.trace())
-
- def df_resid(self):
- """
- Residual degrees of freedom in the fit.
-
- self.N - self.trace()
-
- where self.N is the number of observations of last fit.
- """
-
- return self.N - self.trace()
-
- def df_fit(self):
- """
- How many degrees of freedom used in the fit?
-
- self.trace()
-
- """
- return self.trace()
-
- def trace(self):
- """
- Trace of the smoothing matrix S(pen)
-
- TODO: addin a reference to Wahba, and whoever else I used.
- """
-
- if self.pen > 0:
- _invband = _bspline.invband(self.chol.copy())
- tr = _trace_symbanded(_invband, self.btb, lower=1)
- return tr
- else:
- return self.rank
-
- def fit_target_df(self, y, x=None, df=None, weights=None, tol=1.0e-03,
- apen=0, bpen=1.0e-03):
-
- """
- Fit smoothing spline with approximately df degrees of freedom
- used in the fit, i.e. so that self.trace() is approximately df.
-
- Uses binary search strategy.
-
- In general, df must be greater than the dimension of the null space
- of the Gram inner product. For cubic smoothing splines, this means
- that df > 2.
-
- INPUTS:
- y -- response variable
- x -- if None, uses self.x
- df -- target degrees of freedom
- weights -- optional array of weights
- tol -- (relative) tolerance for convergence
- apen -- lower bound of penalty for binary search
- bpen -- upper bound of penalty for binary search
-
- OUTPUTS: None
- The smoothing spline is determined by self.coef,
- subsequent calls of __call__ will be the smoothing spline.
-
- """
-
- df = df or self.target_df
-
- olddf = y.shape[0] - self.m
-
- if hasattr(self, "pen"):
- self.fit(y, x=x, weights=weights, pen=self.pen)
- curdf = self.trace()
- if N.fabs(curdf - df) / df < tol:
- return
- if curdf > df:
- apen, bpen = self.pen, 2 * self.pen
- else:
- apen, bpen = 0., self.pen
-
- while True:
-
- curpen = 0.5 * (apen + bpen)
- self.fit(y, x=x, weights=weights, pen=curpen)
- curdf = self.trace()
- if curdf > df:
- apen, bpen = curpen, 2 * curpen
- else:
- apen, bpen = apen, curpen
- if apen >= self.penmax:
- raise ValueError, "penalty too large, try setting penmax \
- higher or decreasing df"
- if N.fabs(curdf - df) / df < tol:
- break
-
- def fit_optimize_gcv(self, y, x=None, weights=None, tol=1.0e-03,
- brack=(-100,20)):
- """
- Fit smoothing spline trying to optimize GCV.
-
- Try to find a bracketing interval for scipy.optimize.golden
- based on bracket.
-
- It is probably best to use target_df instead, as it is
- sometimes difficult to find a bracketing interval.
-
- INPUTS:
- y -- response variable
- x -- if None, uses self.x
- df -- target degrees of freedom
- weights -- optional array of weights
- tol -- (relative) tolerance for convergence
- brack -- an initial guess at the bracketing interval
-
- OUTPUTS: None
- The smoothing spline is determined by self.coef,
- subsequent calls of __call__ will be the smoothing spline.
-
- """
-
- def _gcv(pen, y, x):
- self.fit(y, x=x, pen=N.exp(pen))
- a = self.gcv()
- return a
-
- a = golden(_gcv, args=(y,x), brack=bracket, tol=tol)
Modified: trunk/scipy/stats/models/bspline.py
===================================================================
--- trunk/scipy/stats/models/bspline.py 2008-08-08 21:10:46 UTC (rev 4629)
+++ trunk/scipy/stats/models/bspline.py 2008-08-08 21:35:24 UTC (rev 4630)
@@ -20,8 +20,16 @@
from scipy.linalg import solveh_banded
from scipy.optimize import golden
-from scipy.stats.models import _bspline
+from scipy.stats.models import _hbspline
+
+# Issue warning regarding heavy development status of this module
+import warnings
+_msg = "The bspline code is technology preview and requires significant work\
+on the public API and documentation. The API will likely change in the future"
+warnings.warn(_msg, UserWarning)
+
+
def _band2array(a, lower=0, symmetric=False, hermitian=False):
"""
Take an upper or lower triangular banded matrix and return a
@@ -190,6 +198,10 @@
Bspline to avoid extra evaluation in the __call__ method
'''
+ # FIXME: update parameter names, replace single character names
+ # FIXME: `order` should be actual spline order (implemented as order+1)
+ ## FIXME: update the use of spline order in extension code (evaluate is recursively called)
+ # FIXME: eliminate duplicate M and m attributes (m is order, M is related to tau size)
def __init__(self, knots, order=4, M=None, coef=None, x=None):
@@ -277,7 +289,7 @@
x.shape = (N.product(_shape,axis=0),)
if i < self.tau.shape[0] - 1:
## TODO: OWNDATA flags...
- v = _bspline.evaluate(x, self.tau, self.m, d, i, i+1)
+ v = _hbspline.evaluate(x, self.tau, self.m, d, i, i+1)
else:
return N.zeros(x.shape, N.float64)
@@ -321,7 +333,7 @@
d = N.asarray(d)
if d.shape == ():
- v = _bspline.evaluate(x, self.tau, self.m, int(d), lower, upper)
+ v = _hbspline.evaluate(x, self.tau, self.m, int(d), lower, upper)
else:
if d.shape[0] != 2:
raise ValueError, "if d is not an integer, expecting a jx2 \
@@ -330,7 +342,7 @@
v = 0
for i in range(d.shape[1]):
- v += d[1,i] * _bspline.evaluate(x, self.tau, self.m, d[0,i], lower, upper)
+ v += d[1,i] * _hbspline.evaluate(x, self.tau, self.m, d[0,i], lower, upper)
v.shape = (upper-lower,) + _shape
if upper == self.tau.shape[0] - self.m:
@@ -373,7 +385,7 @@
d = N.squeeze(d)
if N.asarray(d).shape == ():
- self.g = _bspline.gram(self.tau, self.m, int(d), int(d))
+ self.g = _hbspline.gram(self.tau, self.m, int(d), int(d))
else:
d = N.asarray(d)
if d.shape[0] != 2:
@@ -385,7 +397,7 @@
self.g = 0
for i in range(d.shape[1]):
for j in range(d.shape[1]):
- self.g += d[1,i]* d[1,j] * _bspline.gram(self.tau, self.m, int(d[0,i]), int(d[0,j]))
+ self.g += d[1,i]* d[1,j] * _hbspline.gram(self.tau, self.m, int(d[0,i]), int(d[0,j]))
self.g = self.g.T
self.d = d
return N.nan_to_num(self.g)
@@ -425,6 +437,8 @@
fhat = ARGMIN_f SUM_i=1^n (y_i-f(x_i))^2 + pen * int f^(2)^2
+ int is integral. pen is lambda (from Hastie)
+
See Chapter 5 of
Hastie, Tibshirani and Friedman (2001). "The Elements of Statistical
@@ -553,7 +567,7 @@
"""
if self.pen > 0:
- _invband = _bspline.invband(self.chol.copy())
+ _invband = _hbspline.invband(self.chol.copy())
tr = _trace_symbanded(_invband, self.btb, lower=1)
return tr
else:
Deleted: trunk/scipy/stats/models/bspline_module.py
===================================================================
--- trunk/scipy/stats/models/bspline_module.py 2008-08-08 21:10:46 UTC (rev 4629)
+++ trunk/scipy/stats/models/bspline_module.py 2008-08-08 21:35:24 UTC (rev 4630)
@@ -1,381 +0,0 @@
-import numpy as N
-from scipy.weave import ext_tools
-import scipy.special.orthogonal
-
-def setup_bspline_module():
- """
- Builds an extension module with Bspline basis calculators using
- weave.
- """
-
- mod = ext_tools.ext_module('_bspline', compiler='gcc')
- knots = N.linspace(0,1,11).astype(N.float64)
- nknots = knots.shape[0]
- x = N.array([0.4,0.5], N.float64)
- nx = x.shape[0]
- m = 4
- d = 0
- lower = 0
- upper = 13
-
- # Bspline code in C
- eval_code = '''
- double *bspline(double **output, double *x, int nx,
- double *knots, int nknots,
- int m, int d, int lower, int upper)
- {
- int nbasis;
- int index, i, j, k;
- double *result, *b, *b0, *b1;
- double *f0, *f1;
- double denom;
-
- nbasis = upper - lower;
-
- result = *((double **) output);
- f0 = (double *) malloc(sizeof(*f0) * nx);
- f1 = (double *) malloc(sizeof(*f1) * nx);
-
- if (m == 1) {
- for(i=0; i<nbasis; i++) {
- index = i + lower;
-
- if(index < nknots - 1) {
- if ((knots[index] != knots[index+1]) && (d <= 0)) {
- for (k=0; k<nx; k++) {
-
- *result = (double) (x[k] >= knots[index]) * (x[k] < knots[index+1]);
- result++;
- }
- }
- else {
- for (k=0; k<nx; k++) {
- *result = 0.;
- result++;
- }
- }
- }
- else {
- for (k=0; k<nx; k++) {
- *result = 0.;
- result++;
- }
- }
- }
- }
- else {
- b = (double *) malloc(sizeof(*b) * (nbasis+1) * nx);
- bspline(&b, x, nx, knots, nknots, m-1, d-1, lower, upper+1);
-
- for(i=0; i<nbasis; i++) {
- b0 = b + nx*i;
- b1 = b + nx*(i+1);
-
- index = i+lower;
-
- if ((knots[index] != knots[index+m-1]) && (index+m-1 < nknots)) {
- denom = knots[index+m-1] - knots[index];
- if (d <= 0) {
- for (k=0; k<nx; k++) {
- f0[k] = (x[k] - knots[index]) / denom;
- }
- }
- else {
- for (k=0; k<nx; k++) {
- f0[k] = (m-1) / (knots[index+m-1] - knots[index]);
- }
- }
- }
- else {
- for (k=0; k<nx; k++) {
- f0[k] = 0.;
- }
- }
-
- index = i+lower+1;
- if ((knots[index] != knots[index+m-1]) && (index+m-1 < nknots)) {
- denom = knots[index+m-1] - knots[index];
- if (d <= 0) {
- for (k=0; k<nx; k++) {
- f1[k] = (knots[index+m-1] - x[k]) / denom;
- }
- }
- else {
- for (k=0; k<nx; k++) {
- f1[k] = -(m-1) / (knots[index+m-1] - knots[index]);
- }
- }
- }
- else {
- for (k=0; k<nx; k++) {
- f1[k] = 0.;
- }
- }
-
- for (k=0; k<nx; k++) {
- *result = f0[k]*(*b0) + f1[k]*(*b1);
- b0++; b1++; result++;
- }
- }
- free(b);
- }
- free(f0); free(f1);
- result = result - nx * nbasis;
-
- return(result);
- }
- '''
-
- eval_ext_code = '''
-
- npy_intp dim[2] = {upper-lower, Nx[0]};
- PyArrayObject *basis;
- double *data;
-
- basis = (PyArrayObject *) PyArray_SimpleNew(2, dim, PyArray_DOUBLE);
- data = (double *) basis->data;
- bspline(&data, x, Nx[0], knots, Nknots[0], m, d, lower, upper);
- return_val = (PyObject *) basis;
- Py_DECREF((PyObject *) basis);
-
- '''
-
- bspline_eval = ext_tools.ext_function('evaluate',
- eval_ext_code,
- ['x', 'knots',
- 'm', 'd', 'lower', 'upper'])
- mod.add_function(bspline_eval)
- bspline_eval.customize.add_support_code(eval_code)
-
- nq = 18
- qx, qw = scipy.special.orthogonal.p_roots(nq)
- dl = dr = 2
-
- gram_code = '''
-
- double *bspline_prod(double *x, int nx, double *knots, int nknots,
- int m, int l, int r, int dl, int dr)
- {
- double *result, *bl, *br;
- int k;
-
- if (fabs(r - l) <= m) {
- result = (double *) malloc(sizeof(*result) * nx);
- bl = (double *) malloc(sizeof(*bl) * nx);
- br = (double *) malloc(sizeof(*br) * nx);
-
- bl = bspline(&bl, x, nx, knots, nknots, m, dl, l, l+1);
- br = bspline(&br, x, nx, knots, nknots, m, dr, r, r+1);
-
- for (k=0; k<nx; k++) {
- result[k] = bl[k] * br[k];
- }
- free(bl); free(br);
- }
- else {
- for (k=0; k<nx; k++) {
- result[k] = 0.;
- }
- }
-
- return(result);
- }
-
-
- double bspline_quad(double *knots, int nknots,
- int m, int l, int r, int dl, int dr)
-
- /* This is based on scipy.integrate.fixed_quad */
-
- {
- double *y;
- double qx[%(nq)d]={%(qx)s};
- double qw[%(nq)d]={%(qw)s};
- double x[%(nq)d];
- int nq=%(nq)d;
- int k, kk;
- int lower, upper;
- double result, a, b, partial;
-
- result = 0;
-
- /* TO DO: figure out knot span more efficiently */
-
- lower = l - m - 1;
- if (lower < 0) { lower = 0;}
- upper = lower + 2 * m + 4;
- if (upper > nknots - 1) { upper = nknots-1; }
-
- for (k=lower; k<upper; k++) {
- partial = 0.;
- a = knots[k]; b=knots[k+1];
- for (kk=0; kk<nq; kk++) {
- x[kk] = (b - a) * (qx[kk] + 1) / 2. + a;
- }
-
- y = bspline_prod(x, nq, knots, nknots, m, l, r, dl, dr);
-
- for (kk=0; kk<nq; kk++) {
- partial += y[kk] * qw[kk];
- }
- free(y); /* bspline_prod malloc's memory, but does not free it */
-
- result += (b - a) * partial / 2.;
-
- }
-
- return(result);
- }
-
- void bspline_gram(double **output, double *knots, int nknots,
- int m, int dl, int dr)
-
- /* Presumes that the first m and last m knots are to be ignored, i.e.
- the interior knots are knots[(m+1):-(m+1)] and the boundary knots are
- knots[m] and knots[-m]. In this setting the first basis element of interest
- is the 1st not the 0th. Should maybe be fixed? */
-
- {
- double *result;
- int l, r, i, j;
- int nbasis;
-
- nbasis = nknots - m;
-
- result = *((double **) output);
- for (i=0; i<nbasis; i++) {
- for (j=0; j<m; j++) {
- l = i;
- r = l+j;
- *result = bspline_quad(knots, nknots, m, l, r, dl, dr);
- result++;
- }
- }
- }
-
- ''' % {'qx':`[q for q in N.real(qx)]`[1:-1], 'qw':`[q for q in qw]`[1:-1], 'nq':nq}
-
- gram_ext_code = '''
-
- npy_intp dim[2] = {Nknots[0]-m, m};
- double *data;
- PyArrayObject *gram;
-
- gram = (PyArrayObject *) PyArray_SimpleNew(2, dim, PyArray_DOUBLE);
- data = (double *) gram->data;
- bspline_gram(&data, knots, Nknots[0], m, dl, dr);
- return_val = (PyObject *) gram;
- Py_DECREF((PyObject *) gram);
-
- '''
-
- bspline_gram = ext_tools.ext_function('gram',
- gram_ext_code,
- ['knots',
- 'm', 'dl', 'dr'])
-
- bspline_gram.customize.add_support_code(gram_code)
- mod.add_function(bspline_gram)
-
- L = N.zeros((3,10), N.float64)
-
- invband_support_code = '''
-
- void invband_compute(double **dataptr, double *L, int n, int m) {
-
- /* Note: m is number of bands not including the diagonal so L is of size (m+1)xn */
-
- int i,j,k;
- int idx, idy;
- double *data, *odata;
- double diag;
-
- data = *((double **) dataptr);
-
- for (i=0; i<n; i++) {
- diag = L[i];
- data[i] = 1.0 / (diag*diag) ;
-
- for (j=0; j<=m; j++) {
- L[j*n+i] /= diag;
- if (j > 0) { data[j*n+i] = 0;}
- }
- }
-
- for (i=n-1; i>=0; i--) {
- for (j=1; j <= (m<n-1-i ? m:n-1-i); j++) {
- for (k=1; k<=(n-1-i<m ? n-1-i:m); k++) {
- idx = (j<k ? k-j:j-k); idy = (j<k ? i+j:i+k);
- data[j*n+i] -= L[k*n+i] * data[idx*n+idy];
- }
- }
-
- for (k=1; k<=(n-1-i<m ? n-1-i:m); k++) {
- data[i] -= L[k*n+i] * data[k*n+i];
- }
- }
-
- return;
- }
- '''
-
- invband_ext_code = '''
-
- npy_intp dim[2] = {NL[0], NL[1]};
- int i, j;
- double *data;
- PyArrayObject *invband;
-
- invband = (PyArrayObject *) PyArray_SimpleNew(2, dim, PyArray_DOUBLE);
- data = (double *) invband->data;
- invband_compute(&data, L, NL[1], NL[0]-1);
-
- return_val = (PyObject *) invband;
- Py_DECREF((PyObject *) invband);
-
- '''
-
- invband = ext_tools.ext_function('invband',
- invband_ext_code,
- ['L'])
- invband.customize.add_support_code(invband_support_code)
- mod.add_function(invband)
-
- return mod
-
-mod = setup_bspline_module()
-
-def build_bspline_module():
- mod.compile()
-
-# try:
-# import _bspline
-# except ImportError:
-# build_bspline_module()
-# import _bspline
-
-## if __name__ == '__main__':
-## knots = N.hstack([[0]*3, N.linspace(0,1,11).astype(N.float64), [1]*3])
-## x = N.array([0.4,0.5])
-## print bspline_ext.bspline_eval(x, knots, 4, 2, 0, 13)
-
-## knots = N.hstack([[0]*3, N.linspace(0,1,501).astype(N.float64), [1]*3])
-## nknots = knots.shape[0]
-## x = N.linspace(0,1,1000)
-## m = 4
-## d = 0
-
-
-## import time, gc
-## t = 0
-## for i in range(100):
-## lower = i
-## toc = time.time()
-## gc.collect()
-## y = bspline_ext.bspline_eval(x, knots, m, 2, 0, 503)
-## z = bspline_ext.bspline_prod(x, knots, m, 2, 1, 0, 0)
-## tic = time.time()
-## t += tic-toc
-## del(y); del(z)
-
-## print t / 100
Modified: trunk/scipy/stats/models/setup.py
===================================================================
--- trunk/scipy/stats/models/setup.py 2008-08-08 21:10:46 UTC (rev 4629)
+++ trunk/scipy/stats/models/setup.py 2008-08-08 21:35:24 UTC (rev 4630)
@@ -1,27 +1,21 @@
-def configuration(parent_package='',top_path=None, package_name='models'):
+def configuration(parent_package='',top_path=None):
from numpy.distutils.misc_util import Configuration
- config = Configuration(package_name,parent_package,top_path)
+ config = Configuration('models',parent_package,top_path)
- config.add_subpackage('*')
+ config.add_subpackage('family')
+ config.add_subpackage('robust')
config.add_data_dir('tests')
- try:
- from scipy.stats.models.bspline_module import mod
- n, s, d = weave_ext(mod)
- config.add_extension(n, s, **d)
- except ImportError: pass
-
+ config.add_extension('_hbspline',
+ sources=['src/bspline_ext.c',
+ 'src/bspline_impl.c'],
+ )
return config
-def weave_ext(mod):
- d = mod.setup_extension().__dict__
- n = d['name']; del(d['name'])
- s = d['sources']; del(d['sources'])
- return n, s, d
-
if __name__ == '__main__':
from numpy.distutils.core import setup
- setup(**configuration(top_path='', package_name='scipy.stats.models').todict())
+ setup(**configuration(top_path='').todict())
+
Modified: trunk/scipy/stats/models/smoothers.py
===================================================================
--- trunk/scipy/stats/models/smoothers.py 2008-08-08 21:10:46 UTC (rev 4629)
+++ trunk/scipy/stats/models/smoothers.py 2008-08-08 21:35:24 UTC (rev 4630)
@@ -9,10 +9,9 @@
from scipy.linalg import solveh_banded
from scipy.optimize import golden
-from scipy.stats.models import _bspline
-from scipy.stats.models.bspline import bspline, _band2array
+from scipy.stats.models import _hbspline
+from scipy.stats.models.bspline import BSpline, _band2array
-
class PolySmoother:
"""
Polynomial smoother up to a given order.
@@ -61,7 +60,7 @@
_y = y * _w
self.coef = N.dot(L.pinv(X).T, _y)
-class SmoothingSpline(bspline):
+class SmoothingSpline(BSpline):
penmax = 30.
@@ -153,7 +152,7 @@
"""
if self.pen > 0:
- _invband = _bspline.invband(self.chol.copy())
+ _invband = _hbspline.invband(self.chol.copy())
tr = _trace_symbanded(_invband, self.btb, lower=1)
return tr
else:
@@ -174,7 +173,7 @@
def __init__(self, knots, order=4, coef=None, M=None, target_df=None):
if target_df is not None:
self.target_df = target_df
- bspline.__init__(self, knots, order=order, coef=coef, M=M)
+ BSpline.__init__(self, knots, order=order, coef=coef, M=M)
self.target_reached = False
def fit(self, y, x=None, df=None, weights=None, tol=1.0e-03):
Copied: trunk/scipy/stats/models/src (from rev 4629, branches/stats_models/src)
Deleted: trunk/scipy/stats/models/src/bspline_ext.c
===================================================================
--- branches/stats_models/src/bspline_ext.c 2008-08-08 21:10:46 UTC (rev 4629)
+++ trunk/scipy/stats/models/src/bspline_ext.c 2008-08-08 21:35:24 UTC (rev 4630)
@@ -1,135 +0,0 @@
-#include "Python.h"
-#include "numpy/arrayobject.h"
-
-/* function prototypes */
-
-double *bspline(double**, double*, int, double *, int, int, int, int, int);
-void bspline_gram(double **, double *, int, int, int, int);
-void invband_compute(double **, double *, int, int);
-
-
-static PyObject *BSpline_Invband(PyObject *self, PyObject *args)
-{
-
- double *data;
- double *L_data;
- npy_intp *dims_invband;
- npy_intp *dims_L;
- PyArrayObject *L = NULL;
- PyArrayObject *invband = NULL;
-
- if(!PyArg_ParseTuple(args, "O", &L))
- goto exit;
-
- dims_L = PyArray_DIMS(L);
- L_data = (double *)PyArray_DATA(L);
-
- dims_invband = calloc(2, sizeof(npy_intp));
- dims_invband[0] = dims_L[0];
- dims_invband[1] = dims_L[1];
-
- invband = (PyArrayObject*)PyArray_SimpleNew(2, dims_invband, PyArray_DOUBLE);
- data = (double *)PyArray_DATA(invband);
- free(dims_invband);
-
- invband_compute(&data, L_data, (int)dims_L[0], (int)dims_L[1]);
-
-exit:
-
- return PyErr_Occurred() ? NULL : (PyObject*)Py_BuildValue("O", invband);
-
-}
-
-
-
-static PyObject *BSpline_Gram(PyObject *self, PyObject *args)
-{
-
- int m;
- int dl;
- int dr;
- double *knots;
- double *data;
- npy_intp *nknots;
- npy_intp *dims_gram;
- PyArrayObject *knots_array = NULL;
- PyArrayObject *gram_array = NULL;
-
- if(!PyArg_ParseTuple(args, "Oiii", &knots_array, &m, &dl, &dr))
- goto exit;
-
- nknots = PyArray_DIMS(knots_array);
- knots = (double *)PyArray_DATA(knots_array);
-
- dims_gram = calloc(2, sizeof(npy_intp));
- dims_gram[0] = (int)nknots[0] - m;
- dims_gram[1] = m;
-
- gram_array = (PyArrayObject*)PyArray_SimpleNew(2, dims_gram, PyArray_DOUBLE);
- data = (double *)PyArray_DATA(gram_array);
- free(dims_gram);
-
- bspline_gram(&data, knots, (int)nknots[0], m, dl, dr);
-
-exit:
-
- return PyErr_Occurred() ? NULL : (PyObject*)Py_BuildValue("O", gram_array);
-
-}
-
-
-static PyObject *BSpline_Evaluate(PyObject *self, PyObject *args)
-{
-
- int i;
- int upper;
- int lower;
- int m;
- int d;
- double *knots;
- double *x;
- double *data;
- npy_intp *nknots;
- npy_intp *nx;
- npy_intp dims_basis[2];
- PyArrayObject *knots_array = NULL;
- PyArrayObject *x_array = NULL;
- PyArrayObject *basis_array = NULL;
-
- if(!PyArg_ParseTuple(args, "OOiiii", &x_array, &knots_array, &m, &d, &lower, &upper))
- goto exit;
-
- nknots = PyArray_DIMS(knots_array);
- nx = PyArray_DIMS(x_array);
-
- knots = (double *)PyArray_DATA(knots_array);
- x = (double *)PyArray_DATA(x_array);
-
- dims_basis[0] = upper-lower;
- dims_basis[1] = (int)nx[0];
- basis_array = (PyArrayObject*)PyArray_SimpleNew(2, dims_basis, PyArray_DOUBLE);
- data = (double *)PyArray_DATA(basis_array);
-
- bspline(&data, x, (int)nx[0], knots, (int)nknots[0], m, d, lower, upper);
-
-exit:
-
- return PyErr_Occurred() ? NULL : (PyObject*)Py_BuildValue("O", basis_array);
-
-}
-
-
-static PyMethodDef BSplineMethods[] =
-{
- { "evaluate", BSpline_Evaluate, METH_VARARGS, NULL },
- { "gram", BSpline_Gram, METH_VARARGS, NULL },
- { "invband", BSpline_Invband, METH_VARARGS, NULL },
- { NULL, NULL, 0, NULL},
-};
-
-PyMODINIT_FUNC init_hbspline(void)
-{
- Py_InitModule("_hbspline", BSplineMethods);
- import_array();
-}
-
Copied: trunk/scipy/stats/models/src/bspline_ext.c (from rev 4629, branches/stats_models/src/bspline_ext.c)
Deleted: trunk/scipy/stats/models/src/bspline_impl.c
===================================================================
--- branches/stats_models/src/bspline_impl.c 2008-08-08 21:10:46 UTC (rev 4629)
+++ trunk/scipy/stats/models/src/bspline_impl.c 2008-08-08 21:35:24 UTC (rev 4630)
@@ -1,274 +0,0 @@
-
-#include <stdlib.h>
-
-/* function prototypes */
-
-double *bspline(double **, double *, int, double *, int, int, int, int, int);
-double bspline_quad(double *, int, int, int, int, int, int);
-double *bspline_prod(double *, int, double *, int, int, int, int, int, int);
-void bspline_gram(double **, double *, int, int, int, int);
-void invband_compute(double **, double *, int, int);
-
-
-double *bspline(double **output, double *x, int nx,
- double *knots, int nknots,
- int m, int d, int lower, int upper){
-
- int nbasis;
- int index, i, j, k;
- double *result, *b, *b0, *b1;
- double *f0, *f1;
- double denom;
-
- nbasis = upper - lower;
-
- result = *((double **) output);
- f0 = (double *) malloc(sizeof(*f0) * nx);
- f1 = (double *) malloc(sizeof(*f1) * nx);
-
- if (m == 1) {
- for(i=0; i<nbasis; i++) {
- index = i + lower;
-
- if(index < nknots - 1) {
- if ((knots[index] != knots[index+1]) && (d <= 0)) {
- for (k=0; k<nx; k++) {
-
- *result = (double) (x[k] >= knots[index]) * (x[k] < knots[index+1]);
- result++;
- }
- }
- else {
- for (k=0; k<nx; k++) {
- *result = 0.;
- result++;
- }
- }
- }
- else {
- for (k=0; k<nx; k++) {
- *result = 0.;
- result++;
- }
- }
- }
- }
- else {
- b = (double *) malloc(sizeof(*b) * (nbasis+1) * nx);
- bspline(&b, x, nx, knots, nknots, m-1, d-1, lower, upper+1);
-
- for(i=0; i<nbasis; i++) {
- b0 = b + nx*i;
- b1 = b + nx*(i+1);
-
- index = i+lower;
-
- if ((knots[index] != knots[index+m-1]) && (index+m-1 < nknots)) {
- denom = knots[index+m-1] - knots[index];
- if (d <= 0) {
- for (k=0; k<nx; k++) {
- f0[k] = (x[k] - knots[index]) / denom;
- }
- }
- else {
- for (k=0; k<nx; k++) {
- f0[k] = (m-1) / (knots[index+m-1] - knots[index]);
- }
- }
- }
- else {
- for (k=0; k<nx; k++) {
- f0[k] = 0.;
- }
- }
-
- index = i+lower+1;
- if ((knots[index] != knots[index+m-1]) && (index+m-1 < nknots)) {
- denom = knots[index+m-1] - knots[index];
- if (d <= 0) {
- for (k=0; k<nx; k++) {
- f1[k] = (knots[index+m-1] - x[k]) / denom;
- }
- }
- else {
- for (k=0; k<nx; k++) {
- f1[k] = -(m-1) / (knots[index+m-1] - knots[index]);
- }
- }
- }
- else {
- for (k=0; k<nx; k++) {
- f1[k] = 0.;
- }
- }
-
- for (k=0; k<nx; k++) {
- *result = f0[k]*(*b0) + f1[k]*(*b1);
- b0++; b1++; result++;
- }
- }
- free(b);
- }
- free(f0); free(f1);
- result = result - nx * nbasis;
-
- return(result);
-}
-
-
-double bspline_quad(double *knots, int nknots,
- int m, int l, int r, int dl, int dr)
-
- /* This is based on scipy.integrate.fixed_quad */
-
- {
- double *y;
- double qx[18]={-0.9915651684209309, -0.95582394957139838,
- -0.89260246649755604, -0.80370495897252303, -0.69168704306035333,
- -0.55977083107394743, -0.41175116146284346, -0.25188622569150576,
- -0.084775013041735417, 0.084775013041735306, 0.25188622569150554,
- 0.41175116146284246, 0.55977083107394743, 0.69168704306035189,
- 0.80370495897252314, 0.89260246649755637, 0.95582394957139616,
- 0.9915651684209319};
- double qw[18]={0.021616013526480963, 0.049714548894972385,
- 0.076425730254889301, 0.10094204410628659, 0.12255520671147889,
- 0.14064291467065104, 0.15468467512626605, 0.16427648374583206,
- 0.16914238296314324, 0.16914238296314299, 0.16427648374583295,
- 0.1546846751262658, 0.14064291467065093, 0.12255520671147752,
- 0.10094204410628753, 0.076425730254888483, 0.049714548894967854,
- 0.021616013526484387};
- double x[18];
- int nq=18;
- int k, kk;
- int lower, upper;
- double result, a, b, partial;
-
- result = 0;
-
- /* TO DO: figure out knot span more efficiently */
-
- lower = l - m - 1;
- if (lower < 0) { lower = 0;}
- upper = lower + 2 * m + 4;
- if (upper > nknots - 1) { upper = nknots-1; }
-
- for (k=lower; k<upper; k++) {
- partial = 0.;
- a = knots[k]; b=knots[k+1];
- for (kk=0; kk<nq; kk++) {
- x[kk] = (b - a) * (qx[kk] + 1) / 2. + a;
- }
-
- y = bspline_prod(x, nq, knots, nknots, m, l, r, dl, dr);
-
- for (kk=0; kk<nq; kk++) {
- partial += y[kk] * qw[kk];
- }
- free(y); /* bspline_prod malloc's memory, but does not free it */
-
- result += (b - a) * partial / 2.;
-
- }
-
- return(result);
-
-}
-
-
-
-double *bspline_prod(double *x, int nx, double *knots, int nknots,
- int m, int l, int r, int dl, int dr){
-
- double *result, *bl, *br;
- int k;
-
- if (abs(r - l) <= m) {
- result = (double *) malloc(sizeof(*result) * nx);
- bl = (double *) malloc(sizeof(*bl) * nx);
- br = (double *) malloc(sizeof(*br) * nx);
-
- bl = bspline(&bl, x, nx, knots, nknots, m, dl, l, l+1);
- br = bspline(&br, x, nx, knots, nknots, m, dr, r, r+1);
-
- for (k=0; k<nx; k++) {
- result[k] = bl[k] * br[k];
- }
- free(bl); free(br);
- }
- else {
- for (k=0; k<nx; k++) {
- result[k] = 0.;
- }
- }
-
- return(result);
-}
-
-void bspline_gram(double **output, double *knots, int nknots,
- int m, int dl, int dr){
-
- /* Presumes that the first m and last m knots are to be ignored, i.e.
- the interior knots are knots[(m+1):-(m+1)] and the boundary knots are
- knots[m] and knots[-m]. In this setting the first basis element of interest
- is the 1st not the 0th. Should maybe be fixed? */
-
-
- double *result;
- int l, r, i, j;
- int nbasis;
-
- nbasis = nknots - m;
-
- result = *((double **) output);
- for (i=0; i<nbasis; i++) {
- for (j=0; j<m; j++) {
- l = i;
- r = l+j;
- *result = bspline_quad(knots, nknots, m, l, r, dl, dr);
- result++;
- }
- }
-
-}
-
-
-void invband_compute(double **dataptr, double *L, int n, int m) {
-
- /* Note: m is number of bands not including the diagonal so L is of size (m+1)xn */
-
- int i,j,k;
- int idx, idy;
- double *data, *odata;
- double diag;
-
- data = *((double **) dataptr);
-
- for (i=0; i<n; i++) {
- diag = L[i];
- data[i] = 1.0 / (diag*diag) ;
-
- for (j=0; j<=m; j++) {
- L[j*n+i] /= diag;
- if (j > 0) { data[j*n+i] = 0;}
- }
- }
-
- for (i=n-1; i>=0; i--) {
- for (j=1; j <= (m<n-1-i ? m:n-1-i); j++) {
- for (k=1; k<=(n-1-i<m ? n-1-i:m); k++) {
- idx = (j<k ? k-j:j-k); idy = (j<k ? i+j:i+k);
- data[j*n+i] -= L[k*n+i] * data[idx*n+idy];
- }
- }
-
- for (k=1; k<=(n-1-i<m ? n-1-i:m); k++) {
- data[i] -= L[k*n+i] * data[k*n+i];
- }
- }
-
- return;
-
-}
-
-
-
Copied: trunk/scipy/stats/models/src/bspline_impl.c (from rev 4629, branches/stats_models/src/bspline_impl.c)
Modified: trunk/scipy/stats/models/tests/test_bspline.py
===================================================================
--- trunk/scipy/stats/models/tests/test_bspline.py 2008-08-08 21:10:46 UTC (rev 4629)
+++ trunk/scipy/stats/models/tests/test_bspline.py 2008-08-08 21:35:24 UTC (rev 4630)
@@ -2,26 +2,47 @@
Test functions for models.bspline
"""
-import numpy as N
+import numpy as np
from numpy.testing import *
+import scipy.stats.models.bspline as bsp
-import scipy.stats.models as S
-try:
- import scipy.stats.models.bspline as B
-except ImportError:
- B = None
-
-
class TestBSpline(TestCase):
def test1(self):
- if B:
- b = B.BSpline(N.linspace(0,10,11), x=N.linspace(0,10,101))
- old = b._basisx.shape
- b.x = N.linspace(0,10,51)
- new = b._basisx.shape
- self.assertEqual((old[0], 51), new)
+ b = bsp.BSpline(np.linspace(0,10,11), x=np.linspace(0,10,101))
+ old = b._basisx.shape
+ b.x = np.linspace(0,10,51)
+ new = b._basisx.shape
+ self.assertEqual((old[0], 51), new)
+ # FIXME: Have no idea what this test does. It's here to simply verify the
+ # C extension is working (in a technical sense, not functional).
+ def test_basis(self):
+ b = bsp.BSpline(np.linspace(0,1,11))
+ x = np.array([0.4, 0.5])
+ v = b.basis(x, lower=0, upper=13)
+ t = np.array([[ 0. , 0. ],
+ [ 0. , 0. ],
+ [ 0. , 0. ],
+ [ 0. , 0. ],
+ [ 0.16666667, 0. ],
+ [ 0.66666667, 0.16666667],
+ [ 0.16666667, 0.66666667],
+ [ 0. , 0.16666667],
+ [ 0. , 0. ],
+ [ 0. , 0. ],
+ [ 0. , 0. ],
+ [ 0. , 0. ],
+ [ 0. , 0. ]])
+ assert_array_almost_equal(v, t, decimal=6)
+ # FIXME: Have no idea what this test does. It's here to simply verify the
+ # C extension is working (in a technical sense, not functional).
+ def test_gram(self):
+ b = bsp.BSpline(np.linspace(0,1,11))
+ grm = b.gram()
+ assert grm.shape == (4, 13)
+
+
if __name__ == "__main__":
run_module_suite()
Modified: trunk/scipy/stats/models/tests/test_formula.py
===================================================================
--- trunk/scipy/stats/models/tests/test_formula.py 2008-08-08 21:10:46 UTC (rev 4629)
+++ trunk/scipy/stats/models/tests/test_formula.py 2008-08-08 21:35:24 UTC (rev 4630)
@@ -63,6 +63,7 @@
self.formula += self.terms[i]
self.formula.namespace = self.namespace
+ @dec.skipknownfailure
def test_namespace(self):
space1 = {'X':N.arange(50), 'Y':N.arange(50)*2}
space2 = {'X':N.arange(20), 'Y':N.arange(20)*2}
Modified: trunk/scipy/stats/models/tests/test_scale.py
===================================================================
--- trunk/scipy/stats/models/tests/test_scale.py 2008-08-08 21:10:46 UTC (rev 4629)
+++ trunk/scipy/stats/models/tests/test_scale.py 2008-08-08 21:35:24 UTC (rev 4630)
@@ -30,11 +30,16 @@
m = scale.MAD(X, axis=-1)
self.assertEquals(m.shape, (40,10))
+ # FIXME: Fix the axis length bug in stats.models.robust.scale.huber
+ # Then resolve ticket #587
+ @dec.skipknownfailure
def test_huber(self):
X = W((40,10))
m = scale.huber(X)
self.assertEquals(m.shape, (10,))
+ # FIXME: Fix the axis length bug in stats.models.robust.scale.huber
+ @dec.skipknownfailure
def test_huberaxes(self):
X = W((40,10,30))
m = scale.huber(X, axis=0)
More information about the Scipy-svn
mailing list