[Lme4-commits] r1461 - in pkg: . Gqr Gqr/R Gqr/man Gqr/src Gqr/tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Nov 29 21:10:06 CET 2011
Author: dmbates
Date: 2011-11-29 21:10:06 +0100 (Tue, 29 Nov 2011)
New Revision: 1461
Added:
pkg/Gqr/
pkg/Gqr/DESCRIPTION
pkg/Gqr/NAMESPACE
pkg/Gqr/R/
pkg/Gqr/R/GaussQuad.R
pkg/Gqr/man/
pkg/Gqr/man/GaussQuad.Rd
pkg/Gqr/src/
pkg/Gqr/src/GQ.cpp
pkg/Gqr/src/GQ.h
pkg/Gqr/src/Makevars
pkg/Gqr/src/Makevars.win
pkg/Gqr/src/external.cpp
pkg/Gqr/tests/
pkg/Gqr/tests/GH_test.R
Log:
Added: pkg/Gqr/DESCRIPTION
===================================================================
--- pkg/Gqr/DESCRIPTION (rev 0)
+++ pkg/Gqr/DESCRIPTION 2011-11-29 20:10:06 UTC (rev 1461)
@@ -0,0 +1,18 @@
+Package: Gqr
+Version: 0.1.0
+Date: 2011-11-30
+Title: Determine knots and weights for Gaussian quadrature rules
+Author: Douglas Bates <bates at stat.wisc.edu>,
+ Martin Maechler <maechler at R-project.org> and
+ Ben Bolker <bbolker at gmail.com>
+Maintainer: <lme4-authors at R-forge.wu-wien.ac.at>
+Description: Determine the knots and the weights for Gaussian quadrature
+ rules based on Legendra, Chebyshev, Gegenbauer, Jacobit, Laguerre,
+ Hermite, Expontential or Rational kernels. Generalized forms of the
+ kernels can be used.
+Depends: R(>= 2.14.0)
+LinkingTo: Rcpp
+LazyLoad: yes
+License: GPL (>=2)
+URL: http://lme4.r-forge.r-project.org/
+Packaged: 2011-11-29 20:07:15 UTC; bates
Added: pkg/Gqr/NAMESPACE
===================================================================
--- pkg/Gqr/NAMESPACE (rev 0)
+++ pkg/Gqr/NAMESPACE 2011-11-29 20:10:06 UTC (rev 1461)
@@ -0,0 +1,2 @@
+useDynLib(Gqr, .registration=TRUE)
+export(GaussQuad)
Added: pkg/Gqr/R/GaussQuad.R
===================================================================
--- pkg/Gqr/R/GaussQuad.R (rev 0)
+++ pkg/Gqr/R/GaussQuad.R 2011-11-29 20:10:06 UTC (rev 1461)
@@ -0,0 +1,15 @@
+GaussQuad <-
+ function(n, rule, a=0, b=1, alpha=0, beta=0) {
+ rules <- c("Legendre", "Chebyshev", "Gegenbauer", "Jacobi",
+ "Laguerre", "Hermite", "Exponential", "Rational", "Type2")
+ if (is.na(rr <- pmatch(rule, rules)))
+ stop("Unknown quadrature rule, ", rule)
+ stopifnot(length(n <- as.integer(n)) == 1L,
+ length(a <- as.numeric(a)) == 1L,
+ length(b <- as.numeric(b)) == 1L,
+ length(alpha <- as.numeric(alpha)) == 1L,
+ length(beta <- as.numeric(beta)) == 1L,
+ n > 0L)
+ .Call(Gauss_Quad, n, rr, a, b, alpha, beta)
+ }
+
Added: pkg/Gqr/man/GaussQuad.Rd
===================================================================
--- pkg/Gqr/man/GaussQuad.Rd (rev 0)
+++ pkg/Gqr/man/GaussQuad.Rd 2011-11-29 20:10:06 UTC (rev 1461)
@@ -0,0 +1,66 @@
+\name{GaussQuad}
+\alias{GaussQuad}
+\title{
+ Create knot and weight vectors for Gaussian quadrature
+}
+\description{
+ Evaluate the positions (the \dQuote{knots}) and the weights for
+ various Gaussian quadrature rules. The quadrature rules are with
+ respect to a kernel as described in the details section. Optional
+ arguments \code{a}, \code{b}, \code{alpha} and \code{beta} are used to
+ generalize the rules.
+}
+\usage{
+GaussQuad(n, rule, a = 0, b = 1, alpha = 0, beta = 0)
+}
+\arguments{
+ \item{n}{integer - number of quadrature points}
+ \item{rule}{character - a partial match to one of the quadrature rule
+ names. See the details section.}
+ \item{a}{numeric scalar - an optional shift parameter or interval endpoint}
+ \item{b}{numeric scalar - an optional scale parameter or interval endpoint}
+ \item{alpha}{numeric scalar - an optional power}
+ \item{beta}{numeric scalar - another optional power}
+}
+\details{
+ The possible values for the \code{rule} character string and the
+ corresponding integrals are:
+ \describe{
+ \item{\dQuote{Legendre}}{\eqn{\int_a^b f(x)\,dx}{int_a^b f(x) dx}}
+ \item{\dQuote{Chebyshev}}{\eqn{\int_a^b f(x)((b-x)*(x-a))^{-0.5}\,dx}{int_a^b f(x) ((b-x)*(x-a))^(-0.5) dx}}
+ \item{\dQuote{Gegenbauer}}{\eqn{\int_a^b f(x)((b-x)*(x-a))^\alpha\,dx}{int_a^b f(x) ((b-x)*(x-a))^alpha dx}}
+ \item{\dQuote{Jacobi}}{\eqn{\int_a^b f(x)(b-x)^\alpha (x-a)^\beta\,dx}{int_a^b f(x) (b-x)^alpha*(x-a)^beta}}
+ \item{\dQuote{Laguerre}}{\eqn{\int_a^\infty f(x)(x-a)^\alpha*exp(-b*(x-a))\,dx}{int_a^inf f(x)(x-a)^alpha*exp(-b*(x-a))}}
+ \item{\dQuote{Hermite}}{\eqn{\int_{-\infty}^\infty f(x)|x-a|^\alpha \exp(-b*(x-a)^2)\,dx}{int_-inf^inf f(x)|x-a|^alpha*exp(-b*(x-a)^2) dx}}
+ \item{\dQuote{Exponential}}{\eqn{\int_a^b f(x)|x-(a+b)/2.0|^\alpha\,dx}{int_-inf^inf f(x)|x-(a+b)/2.0|^alpha dx}}
+ \item{\dQuote{Rational}}{\eqn{\int_a^\infty f(x)(x-a)^\alpha*(x+b)^\beta\,dx}{int_a^inf f(x)(x-a)^alpha*(x+b)^beta}}
+ }
+}
+\value{
+ A list with components
+ \item{knots}{The positions at which to evaluate the function to be integrated}
+ \item{weights}{The numeric weights to be applied.}
+}
+\references{
+ The original FORTRAN implementation is from
+ Sylvan Elhay, Jaroslav Kautsky (1987),
+ \dQuote{Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
+ Interpolatory Quadrature}, \bold{ACM Transactions on Mathematical Software},
+ Volume \bold{13}, Number 4, December 1987, pages 399--415.
+
+ One of the Eispack routines, an implicit QL algorithms from
+ Roger Martin and James Wilkinson (1968), \dQuote{The Implicit QL Algorithm},
+ \bold{Numerische Mathematik}, Volume \bold{12}, Number 5, December
+ 1968, pages 377--383, is used in a modified form.
+
+ The C++ implementation by John Burkardt was modified to a C++ class by
+ Douglas Bates.
+}
+%\author{}
+%\note{}
+%\seealso{}
+\examples{
+do.call(data.frame, GaussQuad(5L, "Hermite")) ## 5-point "physicist" Gauss-Hermite rule
+do.call(data.frame, GaussQuad(5L, "H", b=0.5)) ## 5-point "probabilist" Gauss-Hermite rule
+}
+\keyword{math}
Added: pkg/Gqr/src/GQ.cpp
===================================================================
--- pkg/Gqr/src/GQ.cpp (rev 0)
+++ pkg/Gqr/src/GQ.cpp 2011-11-29 20:10:06 UTC (rev 1461)
@@ -0,0 +1,388 @@
+#include "GQ.h"
+#include <limits>
+#include <cmath>
+#include <stdexcept>
+#include <algorithm>
+#include <iostream>
+
+namespace GaussQuad {
+ inline double r8_sign(const double& x) {return (x < 0) ? -1. : 1.;}
+
+ GQ::GQ(int n, Rule r, double a, double b, double alpha, double beta)
+ : d_n(n), d_rule(r), d_a(a), d_b(b), d_alpha(alpha), d_beta(beta),
+ d_t(n), d_wts(n), d_diag(n), d_sub(n) {
+ std::vector<int> mlt(d_n), ndx(d_n);
+ // Compute the Gauss quadrature formula for default a and b.
+ // Get the Jacobi matrix and zero-th moment.
+ class_matrix();
+ // Compute the knots and weights.
+ sgqf();
+ // Prepare to scale the quadrature formula to other weight
+ // function with valid A and B.
+ std::fill(mlt.begin(), mlt.end(), 1);
+ for (int i = 0; i < d_n; i++) ndx[i] = i + 1;
+ scqf(mlt, ndx);
+ }
+
+/**
+ * parchk checks parameters alpha and beta for classical weight functions.
+ *
+ * @param m the order of the highest moment to be calculated
+ */
+ void GQ::parchk(int m) {
+ // Check alpha for Gegenbauer, Jacobi, Laguerre, Hermite, Exponential.
+ if (Gegenbauer <= d_rule && d_alpha <= -1.0)
+ throw std::runtime_error("parchk: 3 <= kind and alpha <= -1.");
+ // check beta for Jacobi.
+ if (d_rule == Jacobi && d_beta <= -1.0)
+ throw std::runtime_error("parchk: Jacobi rule and beta <= -1.0");
+ // check alpha and beta for rational.
+ if (d_rule == Rational) {
+ double tmp = d_alpha + d_beta + double(m) + 1.0;
+ if (0.0 <= tmp || tmp <= d_beta)
+ throw std::runtime_error("parchk: Rational rule condition on alpha and beta fails");
+ }
+ }
+
+/**
+ * Compute the Jacobi matrix for a quadrature rule.
+ *
+ * This routine computes the diagonal and sub-diagonal
+ * elements of the order d_n tridiagonal symmetric Jacobi matrix
+ * associated with the polynomials orthogonal with respect to
+ * the weight function specified by d_rule.
+ *
+ * For weight functions other than Rational, d_n elements are defined in d_sub even
+ * though only d_n-1 are needed. For the Rational weight function, d_sub[d_n-1] is
+ * set to zero.
+ *
+ * Sets the zero-th moment of the weight function, d_zemu.
+ */
+ void GQ::class_matrix() {
+ const double pi = 3.14159265358979323846264338327950;
+ double a2b2 = 0., ab = 0., abi = 0, abp2 = 0., apone = d_alpha + 1.;
+
+ parchk(2 * d_n - 1);
+ d_zemu = 0.;
+
+ if (500.0 * std::numeric_limits<double>::epsilon() < std::abs(std::pow(tgamma(0.5), 2) - pi))
+ throw std::runtime_error("Gamma function does not match machine parameters.");
+ switch(d_rule) {
+ case Legendre:
+// ab = 0.0;
+// d_zemu = 2.0 / (ab + 1.0);
+ d_zemu = 2.;
+ std::fill(d_diag.begin(), d_diag.end(), double());
+ for (int i = 1; i <= d_n; i++) {
+ double abi = i;// + ab * (i % 2);
+ double abj = 2 * i;// + ab;
+ d_sub[i-1] = std::sqrt(abi * abi / (abj * abj - 1.0));
+ }
+ break;
+ case Chebyshev:
+ d_zemu = pi;
+ std::fill(d_diag.begin(), d_diag.end(), double());
+ std::fill(d_sub.begin(), d_sub.end(), 0.5);
+ d_sub[0] = std::sqrt(0.5);
+ break;
+ case Gegenbauer:
+ ab = d_alpha * 2.0;
+ d_zemu = std::pow(2.0, ab + 1.0) * std::pow(tgamma(d_alpha + 1.0), 2)
+ / tgamma (ab + 2.0);
+ std::fill(d_diag.begin(), d_diag.end(), double());
+ d_sub[0] = std::sqrt(1.0 / (2.0 * d_alpha + 3.0));
+ for(int i = 2; i <= d_n; i++) {
+ d_sub[i-1] = std::sqrt(i * (i + ab) / (4.0 * std::pow(i + d_alpha, 2) - 1.0));
+ }
+ break;
+ case Jacobi:
+ ab = d_alpha + d_beta;
+ abp2 = 2.0 + ab;
+ d_zemu = std::pow(2.0, ab + 1.0) * tgamma (d_alpha + 1.0)
+ * tgamma (d_beta + 1.0) / tgamma (abp2);
+ d_diag[0] = (d_beta - d_alpha) / abp2;
+ d_sub[0] = std::sqrt(4.0 * (1.0 + d_alpha) * (1.0 + d_beta)
+ / ((abi + 1.0) * abp2 * abp2));
+ a2b2 = d_beta * d_beta - d_alpha * d_alpha;
+ for(int i = 2; i <= d_n; i++) {
+ double abi = 2.0 * i + ab;
+ d_diag[i-1] = a2b2 / ((abi - 2.0) * abi);
+ abi = abi * abi;
+ d_sub[i-1] = std::sqrt(4.0 * i * (i + d_alpha) * (i + d_beta) * (i + ab)
+ / ((abi - 1.0) * abi));
+ }
+ break;
+ case Laguerre:
+ d_zemu = tgamma(d_alpha + 1.0);
+ for(int i = 1; i <= d_n; i++) {
+ d_diag[i-1] = 2.0 * i - 1.0 + d_alpha;
+ d_sub[i-1] = std::sqrt(i * (i + d_alpha));
+ }
+ break;
+ case Hermite:
+ d_zemu = tgamma((d_alpha + 1.0) / 2.0);
+ std::fill(d_diag.begin(), d_diag.end(), double());
+ for(int i = 1; i <= d_n; i++) {
+ d_sub[i-1] = std::sqrt((i + d_alpha * (i % 2)) / 2.);
+ }
+ break;
+ case Exponential:
+ ab = d_alpha;
+ d_zemu = 2.0 / (ab + 1.0);
+ std::fill(d_diag.begin(), d_diag.end(), double());
+ for(int i = 1; i <= d_n; i++) {
+ abi = i + ab * (i % 2);
+ double abj = 2 * i + ab;
+ d_sub[i-1] = std::sqrt(abi * abi / (abj * abj - 1.0));
+ }
+ break;
+ case Rational:
+ ab = d_alpha + d_beta;
+ d_zemu = tgamma(d_alpha + 1.0) * tgamma(-(ab + 1.0)) / tgamma(-d_beta);
+ d_diag[0] = -apone / (ab + 2.0);
+ d_sub[0] = -d_diag[0] * (d_beta + 1.0) / (ab + 2.0) / (ab + 3.0);
+ for(int i = 2; i <= d_n; i++) {
+ double abti = ab + 2.0 * i;
+ d_diag[i-1] = ab * apone + 2.0 * (ab + i) * (i - 1);
+ d_diag[i-1] = - d_diag[i-1] / abti / (abti - 2.0);
+ }
+
+ for(int i = 2; i <= d_n - 1; i++) {
+ double abti = ab + 2.0 * i;
+ d_sub[i-1] = i * (d_alpha + i) / (abti - 1.0) * (d_beta + i)
+ / (abti * abti) * (ab + i) / (abti + 1.0);
+ }
+ d_sub[d_n - 1] = 0.0;
+ std::transform(d_sub.begin(), d_sub.end(), d_sub.begin(), ::sqrt);
+ break;
+ case Type2:
+ // This case is not handled in John Burhardt's C++ code
+ break;
+ }
+ }
+
+/**
+ * Diagonalize a symmetric tridiagonal matrix.
+ *
+ * This routine is a slightly modified version of the EISPACK routine to
+ * perform the implicit QL algorithm on a symmetric tridiagonal matrix.
+ *
+ * The authors thank the authors of EISPACK for permission to use this
+ * routine.
+ *
+ * Reference:
+ *
+ * Sylvan Elhay, Jaroslav Kautsky,
+ * Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
+ * Interpolatory Quadrature,
+ * ACM Transactions on Mathematical Software,
+ * Volume 13, Number 4, December 1987, pages 399-415.
+ *
+ * Roger Martin, James Wilkinson,
+ * The Implicit QL Algorithm,
+ * Numerische Mathematik,
+ * Volume 12, Number 5, December 1968, pages 377-383.
+ *
+ * It has been modified to produce the product Q' * Z, where Z is an input
+ * vector and Q is the orthogonal matrix diagonalizing the input matrix.
+ * The changes consist (essentially) of applying the orthogonal transformations
+ * directly to Z as they are generated.
+ *
+ * @param z the value of Q' * Z where Q is the matrix that
+ * diagonalizes the input symmetric tridiagonal matrix
+ */
+ void GQ::imtqlx(std::vector<double>& z) {
+ double b;
+ double c;
+ double f;
+ double g;
+ int i;
+ int itn = 30;
+ int m = 1;
+ int mml;
+ double p;
+ double prec = std::numeric_limits<double>::epsilon();
+ double r;
+ double s;
+ std::vector<double> bi(d_sub);
+
+ if (d_n == 1) return;
+
+ bi[d_n - 1] = 0.0;
+
+ for (int l = 1; l <= d_n; l++) {
+ int j = 0;
+ for ( ; ;) {
+ for (m = l; m <= d_n; m++) {
+ if (m == d_n) break;
+
+ if (std::abs(bi[m-1]) <= prec * (std::abs(d_t[m-1]) + std::abs(d_t[m]))) break;
+ }
+ p = d_t[l-1];
+ if (m == l) break;
+ if (itn <= j) throw std::runtime_error("imtqlx: Iteration limit exceeded");
+ j++;
+ g = (d_t[l] - p) / (2.0 * bi[l-1]);
+ r = std::sqrt(g * g + 1.0);
+ g = d_t[m-1] - p + bi[l-1] / (g + std::abs(r) * r8_sign(g));
+ s = 1.0;
+ c = 1.0;
+ p = 0.0;
+ mml = m - l;
+
+ for(int ii = 1; ii <= mml; ii++) {
+ i = m - ii;
+ f = s * bi[i-1];
+ b = c * bi[i-1];
+
+ if (std::abs(g) <= std::abs(f)) {
+ c = g / f;
+ r = std::sqrt(c * c + 1.0);
+ bi[i] = f * r;
+ s = 1.0 / r;
+ c = c * s;
+ } else {
+ s = f / g;
+ r = std::sqrt(s * s + 1.0);
+ bi[i] = g * r;
+ c = 1.0 / r;
+ s = s * c;
+ }
+ g = d_t[i] - p;
+ r = (d_t[i-1] - g) * s + 2.0 * c * b;
+ p = s * r;
+ d_t[i] = g + p;
+ g = c * r - b;
+ f = z[i];
+ z[i] = s * z[i-1] + c * f;
+ z[i-1] = c * z[i-1] - s * f;
+ }
+ d_t[l-1] = d_t[l-1] - p;
+ bi[l-1] = g;
+ bi[m-1] = 0.0;
+ }
+ }
+ // Sorting.
+ for(int ii = 2; ii <= m; ii++) {
+ int i = ii - 1;
+ int k = i;
+ double p = d_t[i-1];
+
+ for (int j = ii; j <= d_n; j++) {
+ if (d_t[j-1] < p) {
+ k = j;
+ p = d_t[j-1];
+ }
+ }
+
+ if (k != i) {
+ d_t[k-1] = d_t[i-1];
+ d_t[i-1] = p;
+ p = z[i-1];
+ z[i-1] = z[k-1];
+ z[k-1] = p;
+ }
+ }
+ }
+
+
+/**
+ * Scale a quadrature formula to a nonstandard interval.
+ *
+ * @param mlt the multiplicity of the knots
+ * @param ndx
+ */
+ void GQ::scqf (const std::vector<int>& mlt, const std::vector<int>& ndx) {
+ double al = d_alpha;
+ double be = d_beta;
+ double shft = (d_a + d_b)/2.;
+ double slp = (d_b - d_a)/2.;
+
+ switch (d_rule) {
+ case Legendre:
+ case Chebyshev:
+ case Gegenbauer:
+ case Jacobi:
+ case Exponential:
+ case Type2:
+ if (std::abs(d_b - d_a) <= std::numeric_limits<double>::epsilon())
+ throw std::runtime_error("scqf: |B - A| too small");
+ break; // unnecessary but it can't hurt
+ case Laguerre:
+ case Hermite:
+ if (d_b <= 0.0) throw std::runtime_error("scqf: b <= 0");
+ break;
+ case Rational:
+ break;
+ }
+
+ switch (d_rule) {
+ case Legendre:
+ al = 0.0;
+ be = 0.0;
+ break;
+ case Chebyshev:
+ al = -0.5;
+ be = -0.5;
+ break;
+ case Gegenbauer:
+ be = d_alpha;
+ break;
+ case Jacobi:
+ break;
+ case Laguerre:
+ shft = d_a;
+ slp = 1.0 / d_b;
+ be = 0.0;
+ break;
+ case Hermite:
+ shft = d_a;
+ slp = 1.0 / std::sqrt(d_b);
+ be = 0.0;
+ break;
+ case Exponential:
+ be = 0.0;
+ break;
+ case Rational:
+ if (d_a + d_b <= 0.0) throw std::runtime_error("scqf: A + B <= 0.");
+ shft = d_a;
+ slp = d_a + d_b;
+ break;
+ case Type2:
+ al = 0.5;
+ be = 0.5;
+ break;
+ }
+
+ double p = std::pow(slp, al + be + 1.0);
+
+ for (int k = 0; k < d_n; k++) {
+ d_t[k] = shft + slp * d_t[k];
+ int l = std::abs(ndx[k]);
+ if (l != 0) {
+ double tmp = p;
+ for(int i = l - 1; i <= l - 1 + mlt[k] - 1; i++) {
+ d_wts[i] *= tmp;
+ tmp = tmp * slp;
+ }
+ }
+ }
+ }
+
+/**
+ * Compute knots and weights of a Gauss Quadrature formula.
+ *
+ * @param zemu the zero-th moment of the weight function
+ */
+ void GQ::sgqf() {
+ if (d_zemu <= 0.0) throw std::runtime_error("sgqf: zemu <= 0.");
+ // Set up vectors for imtqlx.
+ std::copy(d_diag.begin(), d_diag.end(), d_t.begin());
+ std::fill(d_wts.begin(), d_wts.end(), 0.);
+ d_wts[0] = std::sqrt(d_zemu);
+ // Diagonalize the Jacobi matrix.
+ imtqlx(d_wts);
+ for(int i = 0; i < d_n; i++) d_wts[i] *= d_wts[i];
+ }
+}
+
Added: pkg/Gqr/src/GQ.h
===================================================================
--- pkg/Gqr/src/GQ.h (rev 0)
+++ pkg/Gqr/src/GQ.h 2011-11-29 20:10:06 UTC (rev 1461)
@@ -0,0 +1,40 @@
+// -*- mode: C++; c-indent-level: 4; c-basic-offset: 4; tab-width: 8 -*-
+#ifndef LME4_GHQ_H
+#define LME4_GHQ_H
+
+#include <vector>
+
+namespace GaussQuad {
+ enum Rule{Legendre=1, Chebyshev, Gegenbauer, Jacobi, Laguerre,
+ Hermite, Exponential, Rational, Type2};
+ class GQ {
+ int d_n;
+ Rule d_rule;
+ double d_a, d_b, d_alpha, d_beta, d_zemu;
+ std::vector<double> d_t, d_wts, d_diag, d_sub;
+ public:
+ GQ(int, Rule, double, double, double, double);
+
+ const std::vector<double>& diag() const {return d_diag;}
+ const std::vector<double>& sub() const {return d_sub;}
+ const std::vector<double>& t() const {return d_t;}
+ const std::vector<double>& wts() const {return d_wts;}
+
+ double a() const {return d_a;}
+ double b() const {return d_b;}
+ double alpha() const {return d_alpha;}
+ double beta() const {return d_beta;}
+ double zemu() const {return d_zemu;}
+
+ int n() const {return d_n;}
+ int rule() const {return int(d_rule);}
+
+ void class_matrix();
+ void imtqlx(std::vector<double>&);
+ void scqf(const std::vector<int>&, const std::vector<int>&);
+ void sgqf();
+ void parchk(int);
+ };
+}
+
+#endif
Added: pkg/Gqr/src/Makevars
===================================================================
--- pkg/Gqr/src/Makevars (rev 0)
+++ pkg/Gqr/src/Makevars 2011-11-29 20:10:06 UTC (rev 1461)
@@ -0,0 +1,4 @@
+## -*- mode: makefile; -*-
+
+PKG_LIBS = `$(R_HOME)/bin/Rscript -e "Rcpp:::LdFlags()"`
+
Added: pkg/Gqr/src/Makevars.win
===================================================================
--- pkg/Gqr/src/Makevars.win (rev 0)
+++ pkg/Gqr/src/Makevars.win 2011-11-29 20:10:06 UTC (rev 1461)
@@ -0,0 +1,11 @@
+## -*- mode: makefile; -*-
+
+## This assumes that we can call Rscript to ask Rcpp about its locations
+## Use the R_HOME indirection to support installations of multiple R version
+
+PKG_LIBS = $(shell "${R_HOME}/bin${R_ARCH_BIN}/Rscript.exe" -e "Rcpp:::LdFlags()")
+
+PKG_CPPFLAGS = -I. # do we need this?
+
+
+
Added: pkg/Gqr/src/external.cpp
===================================================================
--- pkg/Gqr/src/external.cpp (rev 0)
+++ pkg/Gqr/src/external.cpp 2011-11-29 20:10:06 UTC (rev 1461)
@@ -0,0 +1,59 @@
+// external.cpp: externally .Call'able functions in GaussQuad package
+//
+// Copyright (C) 2011 Douglas Bates, Martin Maechler and Ben Bolker
+//
+// This file is part of the GaussQuad package
+
+#include <Rcpp.h>
+#include "GQ.h"
+
+extern "C" {
+ using GaussQuad::GQ;
+
+ SEXP Gauss_Quad(SEXP n, SEXP rule, SEXP a, SEXP b, SEXP alpha, SEXP beta) {
+ BEGIN_RCPP;
+ GaussQuad::Rule rr=GaussQuad::Legendre; // -Wall
+ switch(::Rf_asInteger(rule)) {
+ case 1: rr = GaussQuad::Legendre; break;
+ case 2: rr = GaussQuad::Chebyshev; break;
+ case 3: rr = GaussQuad::Gegenbauer; break;
+ case 4: rr = GaussQuad::Jacobi; break;
+ case 5: rr = GaussQuad::Laguerre; break;
+ case 6: rr = GaussQuad::Hermite; break;
+ case 7: rr = GaussQuad::Exponential; break;
+ case 8: rr = GaussQuad::Rational; break;
+ case 9: rr = GaussQuad::Type2; break;
+ default: throw std::invalid_argument("Unknown rule");
+ }
+ GQ gq(::Rf_asInteger(n), rr, ::Rf_asReal(a),
+ ::Rf_asReal(b), ::Rf_asReal(alpha), ::Rf_asReal(beta));
+ return Rcpp::List::create(Rcpp::Named("knots") = gq.t(),
+ Rcpp::Named("weights") = gq.wts()
+// , Rcpp::Named("diag") = gq.diag()
+// , Rcpp::Named("sub") = gq.sub()
+// , Rcpp::Named("zemu") = ::Rf_ScalarReal(gq.zemu())
+ );
+ END_RCPP;
+ }
+}
+
+#include <R_ext/Rdynload.h>
+
+#define CALLDEF(name, n) {#name, (DL_FUNC) &name, n}
+
+static R_CallMethodDef CallEntries[] = {
+ CALLDEF(Gauss_Quad, 6),
+ {NULL, NULL, 0}
+};
+
+/** Initializer for GaussQuad, called upon loading the package.
+ *
+ * Register routines that can be called directly from R.
+ */
+extern "C"
+void R_init_Gqr(DllInfo *dll)
+{
+ R_registerRoutines(dll, NULL, CallEntries, NULL, NULL);
+ R_useDynamicSymbols(dll, (Rboolean)FALSE);
+}
+
Added: pkg/Gqr/tests/GH_test.R
===================================================================
--- pkg/Gqr/tests/GH_test.R (rev 0)
+++ pkg/Gqr/tests/GH_test.R 2011-11-29 20:10:06 UTC (rev 1461)
@@ -0,0 +1,10 @@
+library(Gqr)
+wfmt <- "http://people.sc.fsu.edu/~jburkardt/datasets/quadrature_rules_hermite_probabilist/hermite_probabilist_%03d_w.txt "
+xfmt <- "http://people.sc.fsu.edu/~jburkardt/datasets/quadrature_rules_hermite_probabilist/hermite_probabilist_%03d_x.txt "
+nn <- 2^(1:7) - 1L
+for (n in nn) {
+ loc <- GaussQuad(n, "H", b=0.5)
+ webw <- scan(url(sprintf(wfmt, n)))
+ webx <- scan(url(sprintf(xfmt, n)))
+ stopifnot(all.equal(webw, loc$weights), all.equal(webx, loc$knots))
+}
More information about the Lme4-commits
mailing list