[Rcpp-commits] r4357 - in pkg/RcppGSL/inst/examples: . bSpline
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Jun 20 03:19:30 CEST 2013
Author: edd
Date: 2013-06-20 03:19:29 +0200 (Thu, 20 Jun 2013)
New Revision: 4357
Added:
pkg/RcppGSL/inst/examples/bSpline/
pkg/RcppGSL/inst/examples/bSpline/bSpline.R
pkg/RcppGSL/inst/examples/bSpline/bSpline.cpp
pkg/RcppGSL/inst/examples/bSpline/plotHelper.R
Log:
actually adding + committing bSpline which should have gone in in rev4198
Added: pkg/RcppGSL/inst/examples/bSpline/bSpline.R
===================================================================
--- pkg/RcppGSL/inst/examples/bSpline/bSpline.R (rev 0)
+++ pkg/RcppGSL/inst/examples/bSpline/bSpline.R 2013-06-20 01:19:29 UTC (rev 4357)
@@ -0,0 +1,27 @@
+
+## This example illustrated use of RcppGSL using the 'Rcpp attributes' feature
+##
+## The example comes from Section 39.7 of the GSL Reference manual, and constructs
+## a data set from the curve y(x) = \cos(x) \exp(-x/10) on the interval [0, 15] with
+## added Gaussian noise --- which is then fit via linear least squares using a cubic
+## B-spline basis functions with uniform breakpoints.
+##
+## Obviously all this could be done in R too as R can both generate data, and fit
+## models including (B-)splines. But the point to be made here is that we can very
+## easily translate a given GSL program (thanks to RcppGSL), and get it into R with
+## ease thanks to Rcpp and Rcpp attributes.
+
+require(Rcpp) # load Rcpp
+
+sourceCpp("bSpline.cpp") # compile two functions
+
+dat <- genData() # generate the data
+fit <- fitData(dat) # fit the model, returns matrix and gof measures
+
+X <- fit[["X"]] # extract vectors
+Y <- fit[["Y"]]
+
+op <- par(mar=c(3,3,1,1))
+plot(dat[,"x"], dat[,"y"], pch=19, col="#00000044")
+lines(X, Y, col="orange", lwd=2)
+par(op)
Added: pkg/RcppGSL/inst/examples/bSpline/bSpline.cpp
===================================================================
--- pkg/RcppGSL/inst/examples/bSpline/bSpline.cpp (rev 0)
+++ pkg/RcppGSL/inst/examples/bSpline/bSpline.cpp 2013-06-20 01:19:29 UTC (rev 4357)
@@ -0,0 +1,137 @@
+// -*- mode: C++; c-indent-level: 4; c-basic-offset: 4; indent-tabs-mode: nil; -*-
+
+// [[Rcpp::depends(RcppGSL)]]
+#include <RcppGSL.h>
+
+#include <gsl/gsl_bspline.h>
+#include <gsl/gsl_multifit.h>
+#include <gsl/gsl_rng.h>
+#include <gsl/gsl_randist.h>
+#include <gsl/gsl_statistics.h>
+
+const int N = 200; // number of data points to fit
+const int NCOEFFS = 12; // number of fit coefficients */
+const int NBREAK = (NCOEFFS - 2); // nbreak = ncoeffs + 2 - k = ncoeffs - 2 since k = 4 */
+
+// [[Rcpp::export]]
+Rcpp::List genData() {
+
+ const size_t n = N;
+ size_t i;
+ double dy;
+ gsl_rng *r;
+ RcppGSL::vector<double> w(n), x(n), y(n);
+
+ gsl_rng_env_setup();
+ r = gsl_rng_alloc(gsl_rng_default);
+
+ //printf("#m=0,S=0\n");
+ /* this is the data to be fitted */
+
+ for (i = 0; i < n; ++i) {
+ double sigma;
+ double xi = (15.0 / (N - 1)) * i;
+ double yi = cos(xi) * exp(-0.1 * xi);
+
+ sigma = 0.1 * yi;
+ dy = gsl_ran_gaussian(r, sigma);
+ yi += dy;
+
+ gsl_vector_set(x, i, xi);
+ gsl_vector_set(y, i, yi);
+ gsl_vector_set(w, i, 1.0 / (sigma * sigma));
+
+ //printf("%f %f\n", xi, yi);
+ }
+
+ Rcpp::DataFrame res = Rcpp::DataFrame::create(Rcpp::Named("x") = x,
+ Rcpp::Named("y") = y,
+ Rcpp::Named("w") = w);
+
+ x.free();
+ y.free();
+ w.free();
+ gsl_rng_free(r);
+
+ return(res);
+}
+
+
+
+// [[Rcpp::export]]
+Rcpp::List fitData(Rcpp::DataFrame ds) {
+
+ const size_t ncoeffs = NCOEFFS;
+ const size_t nbreak = NBREAK;
+
+ const size_t n = N;
+ size_t i, j;
+
+ Rcpp::DataFrame D(ds); // construct the data.frame object
+ RcppGSL::vector<double> y = D["y"]; // access columns by name,
+ RcppGSL::vector<double> x = D["x"]; // assigning to GSL vectors
+ RcppGSL::vector<double> w = D["w"];
+
+ gsl_bspline_workspace *bw;
+ gsl_vector *B;
+ gsl_vector *c;
+ gsl_matrix *X, *cov;
+ gsl_multifit_linear_workspace *mw;
+ double chisq, Rsq, dof, tss;
+
+ bw = gsl_bspline_alloc(4, nbreak); // allocate a cubic bspline workspace (k = 4)
+ B = gsl_vector_alloc(ncoeffs);
+
+ X = gsl_matrix_alloc(n, ncoeffs);
+ c = gsl_vector_alloc(ncoeffs);
+ cov = gsl_matrix_alloc(ncoeffs, ncoeffs);
+ mw = gsl_multifit_linear_alloc(n, ncoeffs);
+
+ gsl_bspline_knots_uniform(0.0, 15.0, bw); // use uniform breakpoints on [0, 15]
+
+ for (i = 0; i < n; ++i) { // construct the fit matrix X
+ double xi = gsl_vector_get(x, i);
+
+ gsl_bspline_eval(xi, B, bw); // compute B_j(xi) for all j
+
+ for (j = 0; j < ncoeffs; ++j) { // fill in row i of X
+ double Bj = gsl_vector_get(B, j);
+ gsl_matrix_set(X, i, j, Bj);
+ }
+ }
+
+ gsl_multifit_wlinear(X, w, y, c, cov, &chisq, mw); // do the fit
+
+ dof = n - ncoeffs;
+ tss = gsl_stats_wtss(w->data, 1, y->data, 1, y->size);
+ Rsq = 1.0 - chisq / tss;
+
+ Rcpp::NumericVector FX(151), FY(151); // output the smoothed curve
+ double xi, yi, yerr;
+ for (xi = 0.0, i=0; xi < 15.0; xi += 0.1, i++) {
+ gsl_bspline_eval(xi, B, bw);
+ gsl_multifit_linear_est(B, c, cov, &yi, &yerr);
+ FX[i] = xi;
+ FY[i] = yi;
+ }
+
+ Rcpp::List res =
+ Rcpp::List::create(Rcpp::Named("X")=FX,
+ Rcpp::Named("Y")=FY,
+ Rcpp::Named("chisqdof")=Rcpp::wrap(chisq/dof),
+ Rcpp::Named("rsq")=Rcpp::wrap(Rsq));
+
+ gsl_bspline_free(bw);
+ gsl_vector_free(B);
+ gsl_matrix_free(X);
+ gsl_vector_free(c);
+ gsl_matrix_free(cov);
+ gsl_multifit_linear_free(mw);
+
+ y.free();
+ x.free();
+ w.free();
+
+ return(res);
+}
+
Added: pkg/RcppGSL/inst/examples/bSpline/plotHelper.R
===================================================================
--- pkg/RcppGSL/inst/examples/bSpline/plotHelper.R (rev 0)
+++ pkg/RcppGSL/inst/examples/bSpline/plotHelper.R 2013-06-20 01:19:29 UTC (rev 4357)
@@ -0,0 +1,21 @@
+
+library(Rcpp)
+sourceCpp("bSpline.cpp") # compile two functions
+dat <- genData() # generate the data
+fit <- fitData(dat) # fit the model, returns matrix and gof measures
+
+op <- par(mar=c(3,3,1,1))
+plot(dat[,"x"], dat[,"y"], pch=19, col="#00000044")
+lines(fit[[1]], fit[[2]], col="orange", lwd=2)
+par(op)
+
+#spl <- read.table("/tmp/spline.txt")
+#splA <- spl[1:200,]
+#splB <- spl[-c(1:200),]
+#pdf("/tmp/plotHelper.pdf")
+#op <- par(mar=c(3,3,1,1))
+#plot(splA, pch=19, col="#00000044")
+#lines(splB, col="orange")
+#par(op)
+#dev.off()
+
More information about the Rcpp-commits
mailing list