Rcpp Version 1.0.14
Loading...
Searching...
No Matches
lmGSL.R
Go to the documentation of this file.
1#
2# lm() via GNU GSL -- building on something from the Intro to HPC tutorials
3#
4# Copyright (C) 2010 Dirk Eddelbuettel and Romain Francois
5#
6# This file is part of Rcpp.
7#
8# Rcpp is free software: you can redistribute it and/or modify it
9# under the terms of the GNU General Public License as published by
10# the Free Software Foundation, either version 2 of the License, or
11# (at your option) any later version.
12#
13# Rcpp is distributed in the hope that it will be useful, but
14# WITHOUT ANY WARRANTY; without even the implied warranty of
15# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16# GNU General Public License for more details.
17#
18# You should have received a copy of the GNU General Public License
19# along with Rcpp. If not, see <http://www.gnu.org/licenses/>.
20
21suppressMessages(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
27suppressMessages(require(inline))
28
29lmGSL_old <- function() {
30
31 src <- '
32
33 Rcpp::NumericVector Yr(Ysexp);
34 Rcpp::NumericMatrix Xr(Xsexp);
35
36 int i,j,n = Xr.nrow(), k = Xr.ncol();
37 double chisq;
38
39 gsl_matrix *X = gsl_matrix_alloc (n, k);
40 gsl_vector *y = gsl_vector_alloc (n);
41 gsl_vector *c = gsl_vector_alloc (k);
42 gsl_matrix *cov = gsl_matrix_alloc (k, k);
43 for (i = 0; i < n; i++) {
44 for (j = 0; j < k; j++)
45 gsl_matrix_set (X, i, j, Xr(i,j));
46 gsl_vector_set (y, i, Yr(i));
47 }
48
49 gsl_multifit_linear_workspace *work = gsl_multifit_linear_alloc (n, k);
50 gsl_multifit_linear (X, y, c, cov, &chisq, work);
51 gsl_multifit_linear_free (work);
52
53 Rcpp::NumericVector coefr(k), stderrestr(k);
54 for (i = 0; i < k; i++) {
55 coefr(i) = gsl_vector_get(c,i);
56 stderrestr(i) = sqrt(gsl_matrix_get(cov,i,i));
57 }
58 gsl_matrix_free (X);
59 gsl_vector_free (y);
60 gsl_vector_free (c);
61 gsl_matrix_free (cov);
62
63
64 return Rcpp::List::create( Rcpp::Named( "coef", coefr),
65 Rcpp::Named( "stderr", stderrestr));
66 '
67
68 ## turn into a function that R can call
69 ## compileargs redundant on Debian/Ubuntu as gsl headers are found anyway
70 fun_old <- cxxfunction(signature(Ysexp="numeric", Xsexp="numeric"),
71 src,
72 includes="#include <gsl/gsl_multifit.h>",
73 plugin="RcppGSL")
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
79lmGSL <- function() {
80
81sourceCpp(code='
82#include <RcppGSL.h>
83#include <gsl/gsl_multifit.h>
84// [[Rcpp::depends(RcppGSL)]]
85
86// [[Rcpp::export]]
87Rcpp::List fun(Rcpp::NumericVector Yr, Rcpp::NumericMatrix Xr){
88
89 int i, j, n = Xr.nrow(), k = Xr.ncol();
90 double chisq;
91
92 RcppGSL::Matrix X(n, k); // allocate a gsl_matrix<double> of dim n, k
93 RcppGSL::Vector y(n); // allocate a gsl_vector<double> of length n
94 RcppGSL::Vector c(k); // allocate a gsl_vector<double> of length k
95 RcppGSL::Matrix cov(k, k); // allocate a gsl_matrix<double> of dim k, k
96
97 for (i = 0; i < n; i++) {
98 for (j = 0; j < k; j++)
99 X(i, j) = Xr(i, j);
100 y[i] = Yr(i); // Note vector requires [] not ()
101 }
102
103 gsl_multifit_linear_workspace *work = gsl_multifit_linear_alloc (n, k);
104 gsl_multifit_linear (X, y, c, cov, &chisq, work);
105 gsl_multifit_linear_free (work);
106
107 Rcpp::NumericVector coefr(k), stderrestr(k);
108 for (i = 0; i < k; i++) {
109 coefr(i) = c[i];
110 stderrestr(i) = sqrt(cov(i,i));
111 }
112
113
114 return Rcpp::List::create( Rcpp::Named("coef") = coefr,
115 Rcpp::Named("stderr") = stderrestr);
116}')
117fun
118}