|
19 | 19 | # along with Rcpp. If not, see <http://www.gnu.org/licenses/>. |
20 | 20 |
|
21 | 21 | suppressMessages(require(Rcpp)) |
| 22 | + |
| 23 | +## NOTE: This is the old way to compile Rcpp code inline. |
| 24 | +## The code here has left as a historical artifact and tribute to the old way. |
| 25 | +## Please use the code under the "new" inline compilation section. |
| 26 | + |
22 | 27 | suppressMessages(require(inline)) |
23 | 28 |
|
24 | | -lmGSL <- function() { |
| 29 | +lmGSL_old <- function() { |
25 | 30 |
|
26 | 31 | src <- ' |
27 | 32 |
|
@@ -62,8 +67,55 @@ lmGSL <- function() { |
62 | 67 |
|
63 | 68 | ## turn into a function that R can call |
64 | 69 | ## compileargs redundant on Debian/Ubuntu as gsl headers are found anyway |
65 | | - fun <- cxxfunction(signature(Ysexp="numeric", Xsexp="numeric"), |
| 70 | + fun_old <- cxxfunction(signature(Ysexp="numeric", Xsexp="numeric"), |
66 | 71 | src, |
67 | 72 | includes="#include <gsl/gsl_multifit.h>", |
68 | 73 | plugin="RcppGSL") |
69 | 74 | } |
| 75 | + |
| 76 | +## NOTE: Within this section, the new way to compile Rcpp code inline has been |
| 77 | +## written. Please use the code next as a template for your own project. |
| 78 | + |
| 79 | +lmGSL <- function() { |
| 80 | + |
| 81 | +sourceCpp(code=' |
| 82 | +#include <RcppGSL.h> |
| 83 | +#include <gsl/gsl_multifit.h> |
| 84 | +// [[Rcpp::depends(RcppGSL)]] |
| 85 | +
|
| 86 | +// [[Rcpp::export]] |
| 87 | +Rcpp::List fun(Rcpp::NumericVector Yr, Rcpp::NumericMatrix Xr){ |
| 88 | + |
| 89 | + int i,j,n = Xr.nrow(), k = Xr.ncol(); |
| 90 | + double chisq; |
| 91 | + |
| 92 | + gsl_matrix *X = gsl_matrix_alloc (n, k); |
| 93 | + gsl_vector *y = gsl_vector_alloc (n); |
| 94 | + gsl_vector *c = gsl_vector_alloc (k); |
| 95 | + gsl_matrix *cov = gsl_matrix_alloc (k, k); |
| 96 | + for (i = 0; i < n; i++) { |
| 97 | + for (j = 0; j < k; j++) |
| 98 | + gsl_matrix_set (X, i, j, Xr(i,j)); |
| 99 | + gsl_vector_set (y, i, Yr(i)); |
| 100 | + } |
| 101 | + |
| 102 | + gsl_multifit_linear_workspace *work = gsl_multifit_linear_alloc (n, k); |
| 103 | + gsl_multifit_linear (X, y, c, cov, &chisq, work); |
| 104 | + gsl_multifit_linear_free (work); |
| 105 | + |
| 106 | + Rcpp::NumericVector coefr(k), stderrestr(k); |
| 107 | + for (i = 0; i < k; i++) { |
| 108 | + coefr(i) = gsl_vector_get(c,i); |
| 109 | + stderrestr(i) = sqrt(gsl_matrix_get(cov,i,i)); |
| 110 | + } |
| 111 | + gsl_matrix_free (X); |
| 112 | + gsl_vector_free (y); |
| 113 | + gsl_vector_free (c); |
| 114 | + gsl_matrix_free (cov); |
| 115 | + |
| 116 | + |
| 117 | + return Rcpp::List::create( Rcpp::Named( "coef", coefr), |
| 118 | + Rcpp::Named( "stderr", stderrestr)); |
| 119 | +}') |
| 120 | +fun |
| 121 | +} |
0 commit comments