[Scipy-svn] r4626 - in branches/stats_models: . tests
scipy-svn at scipy.org
scipy-svn at scipy.org
Fri Aug 8 15:55:34 EDT 2008
Author: chris.burns
Date: 2008-08-08 14:55:31 -0500 (Fri, 08 Aug 2008)
New Revision: 4626
Removed:
branches/stats_models/_bspline.py
branches/stats_models/bspline_module.py
Modified:
branches/stats_models/tests/test_bspline.py
Log:
Delete old weave code. Fix bspline import in test.
Deleted: branches/stats_models/_bspline.py
===================================================================
--- branches/stats_models/_bspline.py 2008-08-08 19:52:18 UTC (rev 4625)
+++ branches/stats_models/_bspline.py 2008-08-08 19:55:31 UTC (rev 4626)
@@ -1,668 +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 _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
- 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
-
- '''
- # 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):
-
- 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 = _hbspline.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 = _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 \
- 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] * _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:
- 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 = _hbspline.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] * _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)
-
-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
-
- int is integral. pen is lambda (from Hastie)
-
- 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 = _hbspline.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)
-
-
-
-
-
Deleted: branches/stats_models/bspline_module.py
===================================================================
--- branches/stats_models/bspline_module.py 2008-08-08 19:52:18 UTC (rev 4625)
+++ branches/stats_models/bspline_module.py 2008-08-08 19:55:31 UTC (rev 4626)
@@ -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: branches/stats_models/tests/test_bspline.py
===================================================================
--- branches/stats_models/tests/test_bspline.py 2008-08-08 19:52:18 UTC (rev 4625)
+++ branches/stats_models/tests/test_bspline.py 2008-08-08 19:55:31 UTC (rev 4626)
@@ -4,7 +4,7 @@
import numpy as np
from numpy.testing import *
-import scipy.stats.models._bspline as bsp
+import scipy.stats.models.bspline as bsp
class TestBSpline(TestCase):
More information about the Scipy-svn
mailing list