[Scipy-svn] r4626 - in branches/stats_models: . tests

scipy-svn@scip... scipy-svn@scip...
Fri Aug 8 14:55:34 CDT 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