[Rcpp-commits] r3443 - in pkg/RcppSMC: . inst/include src src/compareGSLtoR
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Jan 22 00:50:23 CET 2012
Author: edd
Date: 2012-01-22 00:50:23 +0100 (Sun, 22 Jan 2012)
New Revision: 3443
Added:
pkg/RcppSMC/.Rbuildignore
pkg/RcppSMC/inst/include/rngR.h
pkg/RcppSMC/src/compareGSLtoR/
pkg/RcppSMC/src/compareGSLtoR/multinom.r
pkg/RcppSMC/src/compareGSLtoR/ugauss.r
Modified:
pkg/RcppSMC/ChangeLog
pkg/RcppSMC/DESCRIPTION
pkg/RcppSMC/inst/include/rng.hh
pkg/RcppSMC/inst/include/sampler.hh
pkg/RcppSMC/inst/include/smctc.hh
pkg/RcppSMC/src/Makevars
pkg/RcppSMC/src/pffuncs.cc
pkg/RcppSMC/src/rng.cc
pkg/RcppSMC/src/simfunctions.cc
Log:
release 0.0.2 which no longer needs the GSL
see ChangeLog for more aggregate details
Added: pkg/RcppSMC/.Rbuildignore
===================================================================
--- pkg/RcppSMC/.Rbuildignore (rev 0)
+++ pkg/RcppSMC/.Rbuildignore 2012-01-21 23:50:23 UTC (rev 3443)
@@ -0,0 +1 @@
+src/compareGSLtoR
Modified: pkg/RcppSMC/ChangeLog
===================================================================
--- pkg/RcppSMC/ChangeLog 2012-01-18 18:34:09 UTC (rev 3442)
+++ pkg/RcppSMC/ChangeLog 2012-01-21 23:50:23 UTC (rev 3443)
@@ -1,3 +1,20 @@
+2012-01-21 Dirk Eddelbuettel <edd at debian.org>
+
+ * DESCRIPTION: Version 0.0.2
+
+ * inst/include/rngR.h: New implementation of the RNG class, relying
+ solely on R thus removing the need to build against the GSL
+
+ * inst/include/rng.hh: #ifdef'ed out as GNU GSL RNGs no longer used
+ * src/rng.cc: Idem
+ * src/simfunctions.cc: Idem; and two calls changed from GSL to R
+
+ * inst/include/sampler.hh: Comment out ctor with GSL RNG type
+ * inst/include/smctc.hh: Include new rngR.h instead of rng.hh
+ * src/pffuncs.cc: Idem
+
+ * src/Makevars: Remove GSL vars for compilation
+
2012-01-15 Dirk Eddelbuettel <edd at debian.org>
* src/rareEventsEx.cpp: Adapted main() from rare-events/main.cc in
Modified: pkg/RcppSMC/DESCRIPTION
===================================================================
--- pkg/RcppSMC/DESCRIPTION 2012-01-18 18:34:09 UTC (rev 3442)
+++ pkg/RcppSMC/DESCRIPTION 2012-01-21 23:50:23 UTC (rev 3443)
@@ -1,7 +1,7 @@
Package: RcppSMC
Type: Package
Title: Rcpp bindings for Sequential Monte Carlo
-Version: 0.0.1
+Version: 0.0.2
Date: $Date$
Author: Dirk Eddelbuettel
Maintainer: Dirk Eddelbuettel <edd at debian.org>
Property changes on: pkg/RcppSMC/DESCRIPTION
___________________________________________________________________
Added: svn:keywords
+ Date
Modified: pkg/RcppSMC/inst/include/rng.hh
===================================================================
--- pkg/RcppSMC/inst/include/rng.hh 2012-01-18 18:34:09 UTC (rev 3442)
+++ pkg/RcppSMC/inst/include/rng.hh 2012-01-21 23:50:23 UTC (rev 3443)
@@ -1,3 +1,4 @@
+#if 0
// SMCTC: rng.hh
//
// Copyright Adam Johansen, 2008.
@@ -120,3 +121,4 @@
#endif
+#endif
Added: pkg/RcppSMC/inst/include/rngR.h
===================================================================
--- pkg/RcppSMC/inst/include/rngR.h (rev 0)
+++ pkg/RcppSMC/inst/include/rngR.h 2012-01-21 23:50:23 UTC (rev 3443)
@@ -0,0 +1,214 @@
+// -*- mode: C++; c-indent-level: 4; c-basic-offset: 4; indent-tabs-mode: nil; -*-
+//
+// rngR.h: Rcpp wrapper for SMC library -- R-only RNG class
+//
+// Copyright (C) 2012 Dirk Eddelbuettel
+//
+// This file is part of RcppSMC.
+//
+// RcppSMC is free software: you can redistribute it and/or modify it
+// under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 2 of the License, or
+// (at your option) any later version.
+//
+// RcppSMC is distributed in the hope that it will be useful, but
+// WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with RInside. If not, see <http://www.gnu.org/licenses/>.
+// from examples/pf/pfexample.cc; pffuncs.cc and pffuncs.hh also used
+
+// RcppSMC builds on top of, and wraps, SMCTC which is
+//
+// Copyright Adam Johansen, 2008.
+//
+// and released under GPL-3, see the copyright headers in inst/include/
+//
+// The main reason rngR.h was written was to provide a drop-in replacement,
+// maintaining the same API, but not requiring the GNU GSL -- as GNU R
+// already provides what we need. The text below is very lightly edited, and
+// if it refers to the GSL was probably written by Adam Johansen for SMCTC.
+// I added a few comments with handle 'DEdd' -- Dirk Eddelbuettel, Jan 2012
+
+//! \file
+//! \brief Random number generation.
+//!
+//! This file contains the definitions for the smc::rng and smc::rnginfo class.
+//! It wraps the random number generation facilities provided by the GSL and provides a convenient interfaces to access several of its more commonly-used features.
+
+#ifndef __SMC_RNGR_H
+#define __SMC_RNGR_H 1.0
+
+#include <Rcpp.h>
+
+namespace smc {
+ ///A gsl-rng information handling class (not templated)
+ ///DEdd: Now essentially empty
+ class gslrnginfo {
+ private:
+ ///This is a null terminated array of the available random number generators.
+ //const gsl_rng_type** typePtArray;
+ ///The number of available random number generators.
+ //int nNumber;
+
+ protected:
+ gslrnginfo() {
+ // DEdd: nothing to do for us
+ }
+
+ public:
+ ///Returns a reference to the sole static instance of this class.
+ static gslrnginfo* GetInstance();
+
+ ///Returns the number of available random number generators.
+ //int GetNumber();
+ ///Returns the name of random number generator number nIndex.
+ //const char* GetNameByIndex(int nIndex);
+ ///Returns a pointer to random number generator nIndex.
+ //const gsl_rng_type* GetPointerByIndex(int nIndex);
+ ///Returns a pointer to the random number generator with name szName (or null if it doesn't exist).
+ //const gsl_rng_type* GetPointerByName(const char* szName);
+ };
+
+ ///The global application instance of the gslrnginfo class:
+ extern gslrnginfo rngset;
+
+
+ ///A random number generator class.
+
+ /// At present this serves as a wrapper for the gsl random number generation code.
+ /// DEdd: Rewritten for R package using R's RNGs
+ class rng {
+ private:
+ ///This is the type of random number generator underlying the class.
+ //const gsl_rng_type* type;
+ ///This is a pointer to the internal workspace of the rng including its current state.
+ //gsl_rng* pWorkspace;
+
+ Rcpp::RNGScope *scopeptr; // DEdd: RNGScope saves and later restores R's RNG state
+
+ public:
+ ///Initialise the random number generator using default settings
+ rng(void) {
+ scopeptr = new Rcpp::RNGScope; // needed for state saving/restoring of R RNGs
+ }
+
+ ///Initialise the random number generator using the default seed for the type
+ //rng(const gsl_rng_type* Type);
+ ///Initialise the random number generator using specified type and seed
+ //rng(const gsl_rng_type* Type,unsigned long int lSeed);
+
+ ///Free the workspace allocated for random number generation
+ ~rng() {
+ delete scopeptr;
+ }
+
+ ///Provide access to the raw random number generator
+ //gsl_rng* GetRaw(void);
+
+ ///Generate a multinomial random vector with parameters (n,w[1:k]) and store it in X
+ void Multinomial(unsigned n, unsigned k, double* w, unsigned* X) {
+ Rcpp::IntegerVector v(k);
+ double sum = 0.0;
+ unsigned int i;
+ for (i=0; i<k; i++) sum += w[i];
+ for (i=0; i<k; i++) w[i] = w[i] / sum;
+
+ // R sources: rmultinom(int n, double* prob, int K, int* rN);
+ rmultinom(static_cast<int>(n), const_cast<double*>(w), static_cast<int>(k), &(v[0]));
+
+ for (i=0; i<k; i++) {
+ X[i] = static_cast<unsigned int>(v[i]);
+ }
+ }
+
+ ///Returns a random integer generated uniformly between the minimum and maximum values specified
+ long UniformDiscrete(long lMin, long lMax) {
+ return ::Rf_runif(static_cast<double>(lMin), static_cast<double>(lMax));
+ }
+
+ ///Returns a random number generated from a Beta distribution with the specified parameters.
+ double Beta(double da, double db) {
+ return ::Rf_rbeta(da, db);
+ }
+
+ ///Returns a random number generated from a Cauchy distribution with the specified scale parameter.
+ double Cauchy(double dScale) {
+ return ::Rf_rcauchy(0.0, dScale);
+ }
+
+ ///Returns a random number generated from an exponential distribution with the specified mean.
+ double Exponential(double dMean) {
+ return ::Rf_rexp(dMean);
+ }
+
+ ///Return a random number generated from a gamma distribution with shape alpha and scale beta.
+ double Gamma(double dAlpha, double dBeta) {
+ return ::Rf_rgamma(dAlpha, dBeta);
+ }
+
+ ///Returns a random number generated from a Laplacian distribution with the specified scale.
+ // "ported" from the gsl_ran_laplace functon
+ double Laplacian(double dScale) {
+ double u;
+ do {
+ u = 2 * ::Rf_runif(0.0,1.0) - 1.0;
+ } while (u == 0.0);
+
+ if (u < 0) {
+ return dScale * log (-u);
+ } else {
+ return -dScale * log (u);
+ }
+ }
+
+ ///Returns a random number generated from a Lognormal distribution of location mu and scale sigma
+ double Lognormal(double dMu, double dSigma) {
+ return ::Rf_rlnorm(dMu, dSigma);
+ }
+
+ ///Return a random number generated from a normal distribution with a specified mean and standard deviation
+ double Normal(double dMean, double dStd) {
+ return ::Rf_rnorm(dMean, dStd);
+ }
+
+ ///Return a random number generated from a standard normal distribution
+ double NormalS(void) {
+ return ::Rf_rnorm(0.0, 1.0);
+ }
+
+ ///Returns a random number from a normal distribution, conditional upon it exceeding the specified threshold.
+ // no Rcpp yet -- double NormalTruncated(double dMean, double dStd, double dThreshold);
+#if 0
+ ///This function simply calls gsl_ran_gaussian_tail with the specified parameters and performs appropriate shifting.
+ /// \param dMean The mean of the distribution.
+ /// \param dStd The standard deviation of the distribution
+ /// \param dThreshold The lower truncation threshold.
+ double rng::NormalTruncated(double dMean, double dStd, double dThreshold)
+ {
+ return dMean + gsl_ran_gaussian_tail(pWorkspace, dThreshold - dMean, dStd);
+ }
+#endif
+
+ ///Return a student-t random number generated with a specified number of degrees of freedom
+ double StudentT(double dDF) {
+ return ::Rf_rt(dDF);
+ }
+
+ ///Return a random number generated uniformly between dMin and dMax
+ double Uniform(double dMin, double dMax) {
+ return ::Rf_runif(dMin, dMax);
+ }
+
+ ///Returns a random number generated from the standard uniform[0,1) distribution
+ double UniformS(void) {
+ return ::Rf_runif(0.0, 1.0);
+ }
+
+ };
+}
+
+
+#endif
Modified: pkg/RcppSMC/inst/include/sampler.hh
===================================================================
--- pkg/RcppSMC/inst/include/sampler.hh 2012-01-18 18:34:09 UTC (rev 3442)
+++ pkg/RcppSMC/inst/include/sampler.hh 2012-01-21 23:50:23 UTC (rev 3443)
@@ -91,7 +91,7 @@
///Create an particle system containing lSize uninitialised particles with the specified mode.
sampler(long lSize, HistoryType htHistoryMode);
///Create an particle system constaining lSize uninitialised particles with the specified mode and random number generator.
- sampler(long lSize, HistoryType htHistoryMode, const gsl_rng_type* rngType, unsigned long nSeed);
+ // -- no GSL sampler(long lSize, HistoryType htHistoryMode, const gsl_rng_type* rngType, unsigned long nSeed);
///Dispose of a sampler.
~sampler();
///Calculates and Returns the Effective Sample Size.
@@ -172,6 +172,7 @@
dResampleThreshold = 0.5 * N;
}
+#if 0
/// The constructor prepares a sampler for use but does not assign any moves to the moveset, initialise the particles
/// or otherwise perform any sampling related tasks. Its main function is to allocate a region of memory in which to
/// store the particle set and to initialise a random number generator.
@@ -200,8 +201,8 @@
rtResampleMode = SMC_RESAMPLE_STRATIFIED;
dResampleThreshold = 0.5 * N;
}
+#endif
-
template <class Space>
sampler<Space>::~sampler()
{
Modified: pkg/RcppSMC/inst/include/smctc.hh
===================================================================
--- pkg/RcppSMC/inst/include/smctc.hh 2012-01-18 18:34:09 UTC (rev 3442)
+++ pkg/RcppSMC/inst/include/smctc.hh 2012-01-21 23:50:23 UTC (rev 3443)
@@ -98,7 +98,8 @@
#include <cstdlib>
#include <iostream>
-#include <gsl/gsl_rng.h>
+//#include <gsl/gsl_rng.h>
+#include "rngR.h"
#include "smc-exception.hh"
#include "sampler.hh"
Modified: pkg/RcppSMC/src/Makevars
===================================================================
--- pkg/RcppSMC/src/Makevars 2012-01-18 18:34:09 UTC (rev 3442)
+++ pkg/RcppSMC/src/Makevars 2012-01-21 23:50:23 UTC (rev 3443)
@@ -1,10 +1,8 @@
-# set by configure
-GSL_CFLAGS = -I/usr/include
-GSL_LIBS = -L/usr/lib -lgsl -lgslcblas -lm
-RCPP_LDFLAGS = -L/usr/local/lib/R/site-library/Rcpp/lib -lRcpp -Wl,-rpath,/usr/local/lib/R/site-library/Rcpp/lib
+## Use the R_HOME indirection to support installations of multiple R version
+PKG_LIBS = `$(R_HOME)/bin/Rscript -e "Rcpp:::LdFlags()"`
-# combine with standard arguments for R
-PKG_CPPFLAGS = $(GSL_CFLAGS) -I../inst/include
-PKG_LIBS += $(GSL_LIBS) $(RCPP_LDFLAGS)
+## combine with standard arguments for R
+PKG_CPPFLAGS = -I../inst/include
+PKG_LIBS += $(RCPP_LDFLAGS)
Added: pkg/RcppSMC/src/compareGSLtoR/multinom.r
===================================================================
--- pkg/RcppSMC/src/compareGSLtoR/multinom.r (rev 0)
+++ pkg/RcppSMC/src/compareGSLtoR/multinom.r 2012-01-21 23:50:23 UTC (rev 3443)
@@ -0,0 +1,53 @@
+
+library(inline)
+
+fgsl <- cxxfunction(signature(ns="integer", ws="numeric"),
+ plugin="RcppGSL",
+ include="#include <gsl/gsl_randist.h>", body='
+
+ unsigned int n = Rcpp::as<int>(ns);
+ Rcpp::NumericVector w = Rcpp::NumericVector(ws);
+ unsigned int k = w.size();
+
+ //void rng::Multinomial(unsigned n, unsigned k, const double* w, unsigned* X)
+ //{
+ // gsl_ran_multinomial(pWorkspace, k, n, w, X);
+ //}
+
+ const gsl_rng_type* type;
+ gsl_rng* pWorkspace;
+
+ gsl_rng_env_setup();
+ type = gsl_rng_default;
+ pWorkspace = gsl_rng_alloc(gsl_rng_default);
+
+ unsigned int *z = new unsigned int[k];
+
+ gsl_ran_multinomial(pWorkspace, k, n, w.begin(), z);
+
+ Rcpp::IntegerVector X(k);
+ for (unsigned int i=0; i<k; i++) X[i] = z[i];
+
+ gsl_rng_free(pWorkspace);
+ delete[] z;
+
+ return X;
+')
+
+fr <- cxxfunction(signature(ns="integer", ws="numeric"), plugin="Rcpp", body='
+
+ RNGScope tmp;
+ unsigned int n = Rcpp::as<int>(ns);
+ Rcpp::NumericVector w(ws);
+ unsigned int k = w.size();
+ Rcpp::IntegerVector x(k);
+
+ rmultinom(n, w.begin(), k, x.begin());
+
+ return x;
+')
+
+
+
+fgsl(30000, c(0.25, 0.25, 0.25, 0.25))
+fr(30000, c(0.25, 0.25, 0.25, 0.25))
Added: pkg/RcppSMC/src/compareGSLtoR/ugauss.r
===================================================================
--- pkg/RcppSMC/src/compareGSLtoR/ugauss.r (rev 0)
+++ pkg/RcppSMC/src/compareGSLtoR/ugauss.r 2012-01-21 23:50:23 UTC (rev 3443)
@@ -0,0 +1,16 @@
+
+library(inline)
+
+f1 <- cxxfunction(signature(xs="numeric"), plugin="Rcpp", body='
+ Rcpp::NumericVector x(xs);
+ return Rcpp::wrap(Rcpp::dnorm(x));
+')
+
+f2 <- cxxfunction(signature(xs="numeric"), plugin="RcppGSL",
+ include="#include <gsl/gsl_randist.h>", body='
+ double x = Rcpp::as<double>(xs);
+ double y = gsl_ran_ugaussian_pdf(x);
+ return Rcpp::wrap(y);
+')
+
+for (i in seq(-3,3,by=0.25)) print(c(f1(i), f2(i)))
Modified: pkg/RcppSMC/src/pffuncs.cc
===================================================================
--- pkg/RcppSMC/src/pffuncs.cc 2012-01-18 18:34:09 UTC (rev 3442)
+++ pkg/RcppSMC/src/pffuncs.cc 2012-01-21 23:50:23 UTC (rev 3443)
@@ -1,9 +1,10 @@
#include <iostream>
#include <cmath>
-#include <gsl/gsl_randist.h>
+//#include <gsl/gsl_randist.h>
#include "smctc.hh"
#include "pffuncs.hh"
+#include "rngR.h"
using namespace std;
Modified: pkg/RcppSMC/src/rng.cc
===================================================================
--- pkg/RcppSMC/src/rng.cc 2012-01-18 18:34:09 UTC (rev 3442)
+++ pkg/RcppSMC/src/rng.cc 2012-01-21 23:50:23 UTC (rev 3443)
@@ -1,3 +1,4 @@
+#if 0
#include <iostream>
#include <cstring>
@@ -244,3 +245,4 @@
return gsl_rng_uniform(pWorkspace);
}
}
+#endif
Modified: pkg/RcppSMC/src/simfunctions.cc
===================================================================
--- pkg/RcppSMC/src/simfunctions.cc 2012-01-18 18:34:09 UTC (rev 3442)
+++ pkg/RcppSMC/src/simfunctions.cc 2012-01-21 23:50:23 UTC (rev 3443)
@@ -1,6 +1,7 @@
#include <iostream>
#include <cmath>
-#include <gsl/gsl_randist.h>
+//#include <gsl/gsl_randist.h>
+#include <Rcpp.h>
#include "smctc.hh"
#include "simfunctions.hh"
@@ -18,10 +19,14 @@
mElement<double> *x = X.GetElement(0);
mElement<double> *y = x->pNext;
//Begin with the density exluding the effect of the potential
- lp = log(gsl_ran_ugaussian_pdf(x->value));
+ //lp = log(gsl_ran_ugaussian_pdf(x->value));
+ Rcpp::NumericVector z = Rcpp::dnorm( Rcpp::NumericVector(x->value), 0.0, 1.0);
+ lp = log( Rcpp::as<double>( z ));
while(y) {
- lp += log(gsl_ran_ugaussian_pdf(y->value - x->value));
+ //lp += log(gsl_ran_ugaussian_pdf(y->value - x->value));
+ z = Rcpp::dnorm( Rcpp::NumericVector(y->value - x->value), 0.0, 1.0);
+ lp += log( Rcpp::as<double>(z));
x = y;
y = x->pNext;
}
More information about the Rcpp-commits
mailing list