[Scipy-svn] r4605 - in branches/stats_models: . src
scipy-svn@scip...
scipy-svn@scip...
Wed Aug 6 17:14:34 CDT 2008
Author: tom.waite
Date: 2008-08-06 17:14:32 -0500 (Wed, 06 Aug 2008)
New Revision: 4605
Added:
branches/stats_models/src/
branches/stats_models/src/bspline_ext.c
branches/stats_models/src/bspline_impl.c
Log:
adding c-ext files for bspline
Added: branches/stats_models/src/bspline_ext.c
===================================================================
--- branches/stats_models/src/bspline_ext.c 2008-08-06 21:46:50 UTC (rev 4604)
+++ branches/stats_models/src/bspline_ext.c 2008-08-06 22:14:32 UTC (rev 4605)
@@ -0,0 +1,142 @@
+#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;
+ PyObject *L = NULL;
+ PyObject *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 = (PyObject*)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]);
+
+ Py_DECREF(invband);
+
+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;
+ PyObject *knots_array = NULL;
+ PyObject *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 = (PyObject*)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);
+
+ Py_DECREF(gram_array);
+
+exit:
+
+ return PyErr_Occurred() ? NULL : (PyObject*)Py_BuildValue("O", gram_array);
+
+}
+
+
+static PyObject *BSpline_Evaluate(PyObject *self, PyObject *args)
+{
+
+ int upper;
+ int lower;
+ int m;
+ int d;
+ double *knots;
+ double *x;
+ double *data;
+ npy_intp *nknots;
+ npy_intp *nx;
+ npy_intp *dims_basis;
+ PyObject *knots_array = NULL;
+ PyObject *x_array = NULL;
+ PyObject *basis_array = NULL;
+
+ if(!PyArg_ParseTuple(args, "OOiiii", &knots_array, &x_array, &lower, &upper, &m, &d))
+ 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 = calloc(2, sizeof(npy_intp));
+ dims_basis[0] = upper-lower;
+ dims_basis[1] = (int)nx[0];
+ basis_array = (PyObject*)PyArray_SimpleNew(2, dims_basis, PyArray_DOUBLE);
+ data = (double *)PyArray_DATA(basis_array);
+ free(dims_basis);
+
+ bspline(data, x, (int)nx[0], knots, (int)nknots[0], m, d, lower, upper);
+
+ Py_DECREF(basis_array);
+
+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_segment(void)
+{
+ Py_InitModule("_hbspline", BSplineMethods);
+ import_array();
+}
+
Added: branches/stats_models/src/bspline_impl.c
===================================================================
--- branches/stats_models/src/bspline_impl.c 2008-08-06 21:46:50 UTC (rev 4604)
+++ branches/stats_models/src/bspline_impl.c 2008-08-06 22:14:32 UTC (rev 4605)
@@ -0,0 +1,272 @@
+
+/* 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 (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);
+}
+
+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;
+
+}
+
+
+
More information about the Scipy-svn
mailing list