Re: [Rcpp-devel] Detecting Singular Matrices in RcppArmadillo
There's an R object that has the machine precision, which could be a
reasonable threshold.
.Machine$double.eps
I believe there is a similar constant in the C++ standard library.
You might also try checking the condition of the matrix instead of the
determinant, but you might take a performance hit. EG:
https://www.quora.com/Why-is-a-condition-number-more-useful-than-a-determinant-for-spotting-problems-with-the-solution-to-ax-b
- Neal
On Fri, Jan 28, 2022 at 8:54 AM Alex Ilich wrote:
> Hi, I have some code I've written to calculate regression parameters using
> ordinarily least squares using RcppArmadillo. This is meant to be performed
> iteratively over many small subsets (sliding window over spatial raster
> data) so I'd like to automatically detect when the regression can't be
> performed and just return a vector of NA's. Currently I check to see if the
> determinant is 0 in order to see if the XtX matrix can be inverted and it
> seems to work in the vast majority of cases but I did run into a case where
> a small amount of precision error was introduced and for some reason this
> causes the determinant to be non-zero in C++ however, the matrix still is
> considered singular (uninvertible) leading to an error. I was wondering if
> anyone had any suggestions? Maybe there's a minimum value non-zero value
> where it's considered essentially singular that I could use instead of 0,
> or maybe there's a way to return a vector of NAs if an error is detected
> rather than stopping execution and printing an error message. I've also
> seen that pinv can be used, but from some tests I ran it's substantially
> slower, so I'd like to avoid using that. I've included some example code
> below. Any suggestions would be greatly appreciated.
>
> Thanks,
> Alex
>
> library(Rcpp)
> cppFunction(depends="RcppArmadillo", "NumericVector C_OLS(arma::mat X,
> arma::mat Y){
> arma::mat Xt = trans(X);
> arma::mat XtX = Xt * X;
> double d = det(XtX);
> Rcout << d; //print determinant
> if(d==0){
> NumericVector B2(X.n_cols, NA_REAL);
> return B2;
> } else{
> arma::mat XtX_inv= inv(XtX);
> NumericVector B = Rcpp::as(wrap(XtX_inv * (Xt *
> Y)));
> return B;
> }}")
>
> X<- matrix(data=c(-250, -250, -250, -250, -250, -250, 250, 200, 150,
> 100, 50, 0, 1, 1, 1, 1, 1, 1), ncol=3, nrow=6)
>
>
> X2<- X
>
> X2[1,2]<- X[1,2] + 5.684342e-14 #Add a slight numeric precision issue
>
> Y<- matrix(1:6, ncol=1)
>
> #Determinant of XtX in R is 0 for both
> det(t(X) %*% X) ==0 #TRUE
> det(t(X2) %*% X2) ==0 #TRUE
>
> C_OLS(X,Y) #Determinant is 0 in C++ and correctly returns a vector of NA's
> #0[1] NA NA NA
>
> C_OLS(X2,Y) #Determinant is not 0 in C++ so the if statement is not
> triggered but the matrix is still uninvertible
> # 4.57764e-05
> # error: inv(): matrix is singular
> # Error in C_OLS(X2, Y) : inv(): matrix is singular
> ___
> Rcpp-devel mailing list
> Rcpp-devel@lists.r-forge.r-project.org
> https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
___
Rcpp-devel mailing list
Rcpp-devel@lists.r-forge.r-project.org
https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
Re: [Rcpp-devel] Rcpp - libtorrent
Rcpp should make this very straightforward; assuming you are putting this in a package (because why wouldn't you?) you would just need to set some flags in src/Makefile so that the external library gets picked up correctly by R, and write a few Rcpp wrapper functions for interacting with it. See also Writing R Extensions for "offical" docs and Dirk's book for plain english explanations. Feel free to share a link eg if you get something up on github - I've wanted to share data via torrent before and didn't find a nice solution. On Sun, Feb 10, 2019 at 11:21 AM Morgan Morgan wrote: > Hi All, > > I hope you are well. > > I wanted to have your view regarding the feasibility of a project. > > There is c++ library called libtorrent. I was wondering if it would be > possible to use Rcpp in order to call some of libtorrent's functionality > from R. > > It seems that there are some python bindings for this library however > nothing for R. > > Do you think it would be something possible? Is it difficult from your > point of view? > > Thank you > Best regards > Morgan > ___ > Rcpp-devel mailing list > Rcpp-devel@lists.r-forge.r-project.org > https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel ___ Rcpp-devel mailing list Rcpp-devel@lists.r-forge.r-project.org https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
Re: [Rcpp-devel] Manipulating json with Rcpp
(Resending to list, sorry) | I was pleasantly surprised how easy it was to make a header-only package > | (by cannabalizing BH and copying the hpp file to inst/include). > > Why? Is LinkingTo: BH not good enough? > I'm not using Boost for json, instead using this https://nlohmann.github.io/json/ for the syntax sugar. But I wanted to start with something that was working instead of starting from scratch, in case there were package configuration that had to be set to make // [[Rcpp::depends(RcppJson)]] compile correctly. But Rcpp::depends seems to find everything automatically, it "just works". > | Maybe someone else will find this useful some day. > > There are also R (at least) cpp-using package ndjson and validatejsonr on > CRAN. > Those are good finds. Not sure if you can use them from inside C++ though. It looks like Bob Rudis is also using the nlohmann/json.h, maybe I can ask him to move json.h from src/ to inst/include for his next release. I think that's all that needed to make Rcpp::depends / LinkingTo work. On Thu, Oct 11, 2018 at 12:18 PM Dirk Eddelbuettel wrote: > > On 11 October 2018 at 11:50, Neal Fultz wrote: > | I recently needed to deal with json on the C++ side of things, and found > | the nlohmann json library. > | > | I made a small wrapper package: https://github.com/nfultz/RcppJson > > Cool. > > | I was pleasantly surprised how easy it was to make a header-only package > | (by cannabalizing BH and copying the hpp file to inst/include). > > Why? Is LinkingTo: BH not good enough? > | > | Maybe someone else will find this useful some day. > > There are also R (at least) cpp-using package ndjson and validatejsonr on > CRAN. > > Dirk > > -- > http://dirk.eddelbuettel.com | @eddelbuettel | e...@debian.org > ___ Rcpp-devel mailing list Rcpp-devel@lists.r-forge.r-project.org https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
[Rcpp-devel] Manipulating json with Rcpp
Hi all, I recently needed to deal with json on the C++ side of things, and found the nlohmann json library. I made a small wrapper package: https://github.com/nfultz/RcppJson I was pleasantly surprised how easy it was to make a header-only package (by cannabalizing BH and copying the hpp file to inst/include). Maybe someone else will find this useful some day. Best, Neal Fultz ___ Rcpp-devel mailing list Rcpp-devel@lists.r-forge.r-project.org https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
Re: [Rcpp-devel] rwishart Rcpp implementation
I had ported the bayesm version a while a back, here it is:
List rwishart(int nu, NumericMatrix const& V){// function to draw
from Wishart (nu,V) and IW// // W ~ W(nu,V) E[W]=nuV// // WI=W^-1
E[WI]=V^-1/(nu-m-1)
RNGScope rngscope; int m = V.nrow();mat Vm(V.begin(), m, m,
false);// Can't vectorise because Rcpp-sugar rchisq doesnt
vectorise df argument // T has sqrt chisqs on diagonal and normals
below diagonal and garbage above diagonal mat T(m,m);for(int i =
0; i < m; i++) {T(i,i) = sqrt(R::rchisq(nu-i)); } for(int j = 0;
j < m; j++) {for(int i = j+1; i < m; i++) { T(i,j) =
R::norm_rand(); }}// explicitly declare T as triangular // top
triangular is NaN ### mat C = trans(trimatl(T)) * chol(Vm); mat CI =
C.i();// C is the upper triangular root of Wishart // therefore, W=C'C
this is the LU decomposition // Inv(W) = CICI'this is the UL
decomp *not LU*!// W is Wishart draw, IW is W^-1 return List::create(
_["W"] = trans(C) * C,_["IW"] = CI * trans(CI),_["C"] = C,
_["CI"] = CI);}
On Thu, Jun 11, 2015 at 12:59 PM, Yue Li wrote:
> Dear List,
>
> I wonder if anyone could share their Rcpp code for an equivalence of
> rWishart. I’ve come across a similar question in the forum:
>
>
> http://stackoverflow.com/questions/23463852/generate-wishart-distribution-using-rwishart-within-rcpp
>
> Based on Dirk’s suggestion, I have my following attempt:
>
> // [[Rcpp::export]]
> mat mvrnormArma(int n, vec mu, mat sigma) {
>
> int ncols = sigma.n_cols;
>
> mat Y = randn(n, ncols);
>
> return repmat(mu, 1, n).t() + Y * chol(sigma);
> }
>
>
> // [[Rcpp::export]]
> mat rWishartCpp(mat sigma, int n) {
>
> mat X = mvrnormArma(1, zeros(sigma.n_rows), sigma);
>
> return X.t() * X;
> }
>
> But this implementation does not take into account the degree of freedom.
> In fact, I’m very fuzzy on how to incorporate df into the sampling process:
>
> I also checked the R code from baysm, namely rwishart, which samples from
> rchisq instead (code pasted below). But I’m not sure how to relate the
> simple naive code I have above with the baysm code below …
>
> function baysm_rwishart(nu, V)
> {
> m = nrow(V)
>
> df = (nu + nu - m + 1) - (nu - m + 1):nu
>
> if (m > 1) {
> T = diag(sqrt(rchisq(c(rep(1, m)), df)))
> T[lower.tri(T)] = rnorm((m * (m + 1)/2 - m))
> }
> else {
> T = sqrt(rchisq(1, df))
> }
>
> U = chol(V)
>
> C = t(T) %*% U
>
> CI = backsolve(C, diag(m))
>
> return(list(W = crossprod(C),
> IW = crossprod(t(CI)), C = C,
> CI = CI))
> }
>
> Yue
>
>
>
> ___
> Rcpp-devel mailing list
> Rcpp-devel@lists.r-forge.r-project.org
> https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
>
___
Rcpp-devel mailing list
Rcpp-devel@lists.r-forge.r-project.org
https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
Re: [Rcpp-devel] rwishart Rcpp implementation
I had ported the bayesm version a while a back, here it is:
List rwishart(int nu, NumericMatrix const& V){// function to draw
from Wishart (nu,V) and IW// // W ~ W(nu,V) E[W]=nuV// // WI=W^-1
E[WI]=V^-1/(nu-m-1)
RNGScope rngscope; int m = V.nrow();mat Vm(V.begin(), m, m,
false);// Can't vectorise because Rcpp-sugar rchisq doesnt
vectorise df argument // T has sqrt chisqs on diagonal and normals
below diagonal and garbage above diagonal mat T(m,m);for(int i =
0; i < m; i++) {T(i,i) = sqrt(R::rchisq(nu-i)); } for(int j = 0;
j < m; j++) {for(int i = j+1; i < m; i++) { T(i,j) =
R::norm_rand(); }}// explicitly declare T as triangular // top
triangular is NaN ### mat C = trans(trimatl(T)) * chol(Vm); mat CI =
C.i();// C is the upper triangular root of Wishart // therefore, W=C'C
this is the LU decomposition // Inv(W) = CICI'this is the UL
decomp *not LU*!// W is Wishart draw, IW is W^-1 return List::create(
_["W"] = trans(C) * C,_["IW"] = CI * trans(CI),_["C"] = C,
_["CI"] = CI);}
On Thu, Jun 11, 2015 at 12:59 PM, Yue Li wrote:
> Dear List,
>
> I wonder if anyone could share their Rcpp code for an equivalence of
> rWishart. I’ve come across a similar question in the forum:
>
>
> http://stackoverflow.com/questions/23463852/generate-wishart-distribution-using-rwishart-within-rcpp
>
> Based on Dirk’s suggestion, I have my following attempt:
>
> // [[Rcpp::export]]
> mat mvrnormArma(int n, vec mu, mat sigma) {
>
> int ncols = sigma.n_cols;
>
> mat Y = randn(n, ncols);
>
> return repmat(mu, 1, n).t() + Y * chol(sigma);
> }
>
>
> // [[Rcpp::export]]
> mat rWishartCpp(mat sigma, int n) {
>
> mat X = mvrnormArma(1, zeros(sigma.n_rows), sigma);
>
> return X.t() * X;
> }
>
> But this implementation does not take into account the degree of freedom.
> In fact, I’m very fuzzy on how to incorporate df into the sampling process:
>
> I also checked the R code from baysm, namely rwishart, which samples from
> rchisq instead (code pasted below). But I’m not sure how to relate the
> simple naive code I have above with the baysm code below …
>
> function baysm_rwishart(nu, V)
> {
> m = nrow(V)
>
> df = (nu + nu - m + 1) - (nu - m + 1):nu
>
> if (m > 1) {
> T = diag(sqrt(rchisq(c(rep(1, m)), df)))
> T[lower.tri(T)] = rnorm((m * (m + 1)/2 - m))
> }
> else {
> T = sqrt(rchisq(1, df))
> }
>
> U = chol(V)
>
> C = t(T) %*% U
>
> CI = backsolve(C, diag(m))
>
> return(list(W = crossprod(C),
> IW = crossprod(t(CI)), C = C,
> CI = CI))
> }
>
> Yue
>
>
>
> ___
> Rcpp-devel mailing list
> Rcpp-devel@lists.r-forge.r-project.org
> https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
>
___
Rcpp-devel mailing list
Rcpp-devel@lists.r-forge.r-project.org
https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
Re: [Rcpp-devel] Difference between runif and unif_rand?
I use the debian packages. Upgrading r-cran-rcpp from testing to unstable seems to have fixed this for me, so good call. On Fri, Jul 12, 2013 at 10:25:17PM -0400, Krzysztof Sakrejda wrote: > On Fri, Jul 12, 2013 at 8:50 PM, Dirk Eddelbuettel wrote: > > > > On 12 July 2013 at 12:28, Neal Fultz wrote: > > | I've been updating the C in an r package to use Rcpp and started seeing > > | so odd results from a function that samples from a discrete > > | distribution. > > For what it's worth, I can't reproduce this so it may be specific to > your versions of R/Rcpp. Krzysztof ___ Rcpp-devel mailing list Rcpp-devel@lists.r-forge.r-project.org https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
[Rcpp-devel] Difference between runif and unif_rand?
I've been updating the C in an r package to use Rcpp and started seeing
so odd results from a function that samples from a discrete
distribution.
I think I've narrowed it down to the sugar runif. The two programs below
are identical except f_works uses unif_rand and f_broke goes through runif.
---
library(Rcpp)
cppFunction('
int f_works(NumericVector p) {
int i = 0;
double u = unif_rand();
do { u -= p[i++]; } while(u > 0);
return i;
}
')
cppFunction('
int f_broke(NumericVector p) {
int i = 0;
double u = Rcpp::runif(1)[0];
do { u -= p[i++]; } while(u > 0);
return i;
}
')
---
When I run f_broke, sometimes there is a number that doesn't belong:
R>p <- 1:4 / 10;
R>table(replicate(1, f_works(p)))
1234
961 2045 2984 4010
R>table(replicate(1, f_broke(p)))
0.427225098479539 0.458629200235009 0.687304416438565 0.735608365153894
0.978602966060862
1 1 1 1
1
1 2 3 4
1053 2010 3021 3911
I'm pretty stumped about why I consistently see odd numbers every couple
thousand draws, becuase runif just wraps unif_rand anyway. My
functions are declared as int, so it should be impossible to
return a double.
What is happening here that I don't get?
Thanks much,
Neal
___
Rcpp-devel mailing list
Rcpp-devel@lists.r-forge.r-project.org
https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
