[Pomp-commits] r134 - in pkg: . R man src

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

Author: kingaa
Date: 2009-05-28 14:05:02 +0200 (Thu, 28 May 2009)
New Revision: 134

Modified:
   pkg/DESCRIPTION
   pkg/R/bsplines.R
   pkg/man/bsplines.Rd
   pkg/src/bspline.c
Log:
fix up the bspline.basis function

Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION	2009-05-26 09:15:40 UTC (rev 133)
+++ pkg/DESCRIPTION	2009-05-28 12:05:02 UTC (rev 134)
@@ -1,8 +1,8 @@
 Package: pomp
 Type: Package
 Title: Statistical inference for partially observed Markov processes
-Version: 0.24-2
-Date: 2009-05-26
+Version: 0.24-3
+Date: 2009-05-28
 Author: Aaron A. King, Edward L. Ionides, Carles Martinez Breto, Steve Ellner, Bruce Kendall
 Maintainer: Aaron A. King <kingaa at umich.edu>
 Description: Inference methods for partially-observed Markov processes

Modified: pkg/R/bsplines.R
===================================================================
--- pkg/R/bsplines.R	2009-05-26 09:15:40 UTC (rev 133)
+++ pkg/R/bsplines.R	2009-05-28 12:05:02 UTC (rev 134)
@@ -1,15 +1,23 @@
 bspline <- function (x, i, degree, knots)
   .Call(bspline_basis_function,x,as.integer(i),as.integer(degree),knots)
 
-bspline.basis <- function (x, degree = 3, knots)
-  .Call(bspline_basis,x,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)
+}
 
 periodic.bspline.basis <- function (x, nbasis, degree = 3, period = 1) {
   if (nbasis<degree)
-    stop("periodic.bspline.basis error: must have nbasis >= degree",call.=FALSE)
+    stop("periodic.bspline.basis error: must have ",sQuote("nbasis")," >= ",sQuote("degree"),call.=FALSE)
   dx <- period/nbasis
   knots <- seq(-degree*dx,period+degree*dx,by=dx)
-  y <- bspline.basis(x%%period,degree,knots)
+  y <- .Call(bspline_basis,x%%period,degree,knots)
   if (degree>0)
     y[,1:degree] <- y[,1:degree]+y[,-(1:nbasis)]
   shift <- floor((degree-1)/2)

Modified: pkg/man/bsplines.Rd
===================================================================
--- pkg/man/bsplines.Rd	2009-05-26 09:15:40 UTC (rev 133)
+++ pkg/man/bsplines.Rd	2009-05-28 12:05:02 UTC (rev 134)
@@ -3,38 +3,39 @@
 \alias{periodic.bspline.basis}
 \title{B-spline bases}
 \description{
-  These functions generate B-spline basis functions as lookup tables.
+  These functions generate B-spline basis functions.
   \code{bspline.basis} gives a set of basis functions.
   \code{periodic.bspline.basis} gives a basis of periodic spline functions.
 }
 \usage{
-bspline.basis(x, degree = 3, knots)
+bspline.basis(x, nbasis, degree = 3)
 periodic.bspline.basis(x, nbasis, degree = 3, period = 1)
 }
 \arguments{
   \item{x}{Vector at which the spline functions are to be evaluated.}
+  \item{nbasis}{The number of basis functions to return.}
   \item{degree}{Degree of requested B-splines.}
-  \item{knots}{Vector of positions of the knots.}
-  \item{nbasis}{The number of basis functions to return.}
   \item{period}{The period of the requested periodic B-splines.}
 }
 \value{
   \item{bspline.basis}{
-    Returns a matrix with \code{length(x)} rows and \code{nbasis=length(knots)-degree-1} columns.
-    Each column contains the values one of the \code{nbasis} spline functions.
+    Returns a matrix with \code{length(x)} rows and \code{nbasis} columns.
+    Each column contains the values one of the spline basis functions.
   }
   \item{periodic.bspline.basis}{
     Returns a matrix with \code{length(x)} rows and \code{nbasis} columns.
+    The basis functions returned are periodic with period \code{period}.
   }
 }
 \author{Aaron A. King (kingaa at umich dot edu)}
 \examples{
-x <- seq(-0.2,1.2,by=0.01)
-y <- bspline.basis(x,degree=7,seq(0,1,length=14))
-matplot(x,y)
+x <- seq(0,2,by=0.01)
+y <- bspline.basis(x,degree=3,nbasis=9)
+matplot(x,y,type='l',ylim=c(0,1.1))
+lines(x,apply(y,1,sum),lwd=2)
 
 x <- seq(-1,2,by=0.01)
 y <- periodic.bspline.basis(x,nbasis=5)
-matplot(x,y)
+matplot(x,y,type='l')
 }
 \keyword{smooth}

Modified: pkg/src/bspline.c
===================================================================
--- pkg/src/bspline.c	2009-05-26 09:15:40 UTC (rev 133)
+++ pkg/src/bspline.c	2009-05-28 12:05:02 UTC (rev 134)
@@ -8,18 +8,23 @@
 
 SEXP bspline_basis (SEXP x, SEXP degree, SEXP knots) {
   int nprotect = 0;
-  SEXP y;
+  SEXP y, xr, kr, di;
   int nx = length(x);
   int nknots = length(knots);
-  int deg = INTEGER_VALUE(degree);
-  int nbasis = nknots-deg-1;
+  int deg;
+  int nbasis;
   double *ydata;
   int i;
+  PROTECT(di = AS_INTEGER(degree)); nprotect++;
+  deg = INTEGER_VALUE(di);
   if (deg < 0) error("must have degree > 0 in 'bspline.basis'");
+  nbasis = nknots-deg-1;
+  PROTECT(xr = AS_NUMERIC(x)); nprotect++;
+  PROTECT(kr = AS_NUMERIC(knots)); nprotect++;
   PROTECT(y = allocMatrix(REALSXP,nx,nbasis)); nprotect++;
   ydata = REAL(y);
   for (i = 0; i < nbasis; i++) {
-    bspline_internal(ydata,REAL(x),nx,i,deg,REAL(knots),nknots);
+    bspline_internal(ydata,REAL(xr),nx,i,deg,REAL(kr),nknots);
     ydata += nx;
   }
   UNPROTECT(nprotect);