[Pomp-commits] r200 - in pkg: R man src

[noreply at r-forge.r-project.org]

Author: kingaa
Date: 2010-01-27 23:20:11 +0100 (Wed, 27 Jan 2010)
New Revision: 200

Modified:
   pkg/R/bsplines.R
   pkg/man/bsplines.Rd
   pkg/src/bspline.c
Log:
- some improvements in the B-spline codes


Modified: pkg/R/bsplines.R
===================================================================
--- pkg/R/bsplines.R	2010-01-26 15:11:44 UTC (rev 199)
+++ pkg/R/bsplines.R	2010-01-27 22:20:11 UTC (rev 200)
@@ -1,15 +1,5 @@
-bspline <- function (x, i, degree, knots)
-  .Call(bspline_basis_function,x,as.integer(i),as.integer(degree),knots)
-
 bspline.basis <- function (x, nbasis, degree = 3) {
-  if (nbasis<=degree)
-    stop("bspline.basis error: must have ",sQuote("nbasis")," > ",sQuote("degree"),call.=FALSE)
-  min.x <- min(x,na.rm=TRUE)
-  max.x <- max(x,na.rm=TRUE)
-  dx <- (max.x-min.x)/(nbasis-degree)
-  tails <- degree*dx
-  knots <- seq(from=min.x-tails,to=max.x+tails,by=dx)
-  .Call(bspline_basis,x,degree,knots)
+  .Call(bspline_basis,x,nbasis,degree)
 }
 
 periodic.bspline.basis <- function (x, nbasis, degree = 3, period = 1) {

Modified: pkg/man/bsplines.Rd
===================================================================
--- pkg/man/bsplines.Rd	2010-01-26 15:11:44 UTC (rev 199)
+++ pkg/man/bsplines.Rd	2010-01-27 22:20:11 UTC (rev 200)
@@ -27,6 +27,10 @@
     The basis functions returned are periodic with period \code{period}.
   }
 }
+\details{
+  Direct access to the underlying C routines is available.
+  See the header file \dQuote{pomp.h} for details.
+}
 \author{Aaron A. King \email{kingaa at umich dot edu}}
 \examples{
 x <- seq(0,2,by=0.01)

Modified: pkg/src/bspline.c
===================================================================
--- pkg/src/bspline.c	2010-01-26 15:11:44 UTC (rev 199)
+++ pkg/src/bspline.c	2010-01-27 22:20:11 UTC (rev 200)
@@ -7,24 +7,32 @@
 
 // B-spline basis
 
-SEXP bspline_basis (SEXP x, SEXP degree, SEXP knots) {
+SEXP bspline_basis (SEXP x, SEXP nbasis, SEXP degree) {
   int nprotect = 0;
-  SEXP y, xr, kr;
+  SEXP y, xr;
   int nx = length(x);
-  int nknots = length(knots);
+  int nb = INTEGER_VALUE(nbasis);
   int deg = INTEGER_VALUE(degree);
-  int nb;
-  double *ydata;
+  int nk = nb+deg+1;
+  double dx, minx, maxx;
+  double knots[nk];
+  double *xdata, *ydata;
   int i;
-  nb = nknots-deg-1;
   if (deg < 0) error("bspline.basis error: must have degree > 0");
-  if (nb<=deg) error("bspline.basis error: must have nbasis > degree");
+  if (nb <= deg) error("bspline.basis error: must have nbasis > degree");
   PROTECT(xr = AS_NUMERIC(x)); nprotect++;
-  PROTECT(kr = AS_NUMERIC(knots)); nprotect++;
   PROTECT(y = allocMatrix(REALSXP,nx,nb)); nprotect++;
+  xdata = REAL(xr);
   ydata = REAL(y);
+  for (i = 1, minx = maxx = xdata[0]; i < nx; i++) {
+    minx = (minx > xdata[i]) ? xdata[i] : minx;
+    maxx = (maxx < xdata[i]) ? xdata[i] : maxx;
+  }
+  dx = (maxx-minx)/((double) (nb-deg));
+  knots[0] = minx-deg*dx;
+  for (i = 1; i < nk; i++) knots[i] = knots[i-1]+dx;
   for (i = 0; i < nb; i++) {
-    bspline_internal(ydata,REAL(xr),nx,i,deg,REAL(kr),nknots);
+    bspline_internal(ydata,xdata,nx,i,deg,&knots[0],nk);
     ydata += nx;
   }
   UNPROTECT(nprotect);