[Rcpp-devel] tracking volatile bug
Hi everybody,
I am facing a volatile bug which appears and disappears in
identical calls on identical data.
My RcppArmadilla code is expm_cpp() (it is obtained by rex2arma tool
which I have presented before, cf. attached file). It is compared and
benchmarked vs its R counterpart expm.higham() (also attached). Both
are computing matrix exponential by a same particular method.
The results are identical for R and Rcpp on small and medium sized
matrices nxn n=10:100 but on large matrices (n 800) various error messages
can interrupt (or not) a normal run of expm_cpp().
Sometimes message says (in French) type 'char' indisponible dans 'eval',
I suppose in English it must be unimplemented type 'char' in 'eval'
I have seen (thanks to Google) a similar error message in the following snippet:
/* call this R command: source(FileName) */
int errorOccurred;
SEXP e = lang2(install(source), mkString(FileName));
/* mkChar instead of mkString would lead to this runtime error:
* Error in source(FileName) : unimplemented type 'char' in 'eval' */
R_tryEval(e, R_GlobalEnv, errorOccurred);
which suggests that somewhere in Rcpp or RcppArmadillo there is
a mkChar() call instead of mkString().
Other times, error message can say something like
argument type[1]='x' must be one of 'M','1','O','I','F' or 'E'
or argument type[1]='character' must be a one-letter character string
This latter message is somewhat volatile per se. The part of message
just after type[1]= can be 'character' (as above) or 'method' or 'ANY' etc.
I have found these messages in the Matrix package
https://r-forge.r-project.org/scm/viewvc.php/pkg/Matrix/src/Mutils.c?view=markuproot=matrixpathrev=2614
function char La_norm_type(const char *typstr)
Seemingly, the argument *typstr is getting corrupted somewhere
on the road.
It is useless to say that debugging is of no help here.
If run in a debugger, the program stops normally with
or without error messages cited above.
I have also tried to low the level of gcc optimization
both in compiling R and the Rcpp code but it didn't help.
Anybody has an experience in tracking down similar cases?
This example can be run as:
library(expm)
library(Rcpp)
source(expm_cpp.R)
source(expm.higham.R)
n=1000
As=matrix(rnorm(n*n), n)
stopifnot(diff(range(expm.higham(As, balancing=TRUE)-expm_cpp(As,
balancing=TRUE))) 1.e-14)
The last command may be run several times before an error shows up.
Cheers,
Serguei.
cppFunction(depends='RcppArmadillo', rebuild=TRUE, includes='
template typename T
inline unsigned which_max(T v) {
unsigned i;
v.max(i);
return i+1;
}
templatetypename T
inline unsigned which_min(T v) {
unsigned i;
v.min(i);
return i+1;
}
',
using namespace arma;
using namespace Rcpp;
SEXP expm_cpp(
NumericMatrix A_in_,
bool balancing) {
// auxiliary functions
Environment base_env_r_=Environment::base_env();
Function rep_r_=base_env_r_[\rep\];
Function c_r_=base_env_r_[\c\];
// External R function declarations
Function expm_balance_r_=Environment(\package:expm\)[\balance\];
Function Matrix_norm_r_=Environment(\package:Matrix\)[\norm\];
// Input variable declarations and conversion
mat A(A_in_.begin(), A_in_.nrow(), A_in_.ncol(), false);
// Output and intermediate variable declarations
mat A2;
mat B;
mat B2;
mat B4;
mat B6;
List baP;
List baS;
mat C;
vec c_;
ivec d;
vec dd;
int i;
mat I;
int k;
int l;
int n;
double nA;
mat P;
ivec pp;
double s;
vec t;
mat tt;
mat U;
mat V;
mat X;
// Translated code starts here
d=ivec(IntegerVector::create(A.n_rows, A.n_cols));
if (d.size() != 2 || d.at(0) != d.at(1)) stop(\'A' must be a square
matrix\);
n=d.at(0);
if (n = 1) return wrap(exp(A));
if (balancing) {
baP=asList(expm_balance_r_(A, \P\));
baS=asList(expm_balance_r_(asmat(baP[\z\]), \S\));
A=asmat(baS[\z\]);
}
nA=asdouble(Matrix_norm_r_(A, \1\));
I=eyemat(n, n);
if (nA = 2.1) {
t=vec({0.015, 0.25, 0.95, 2.1});
l=which_max(nA = t);
C=join_vert(join_vert(join_vert(vec({120, 60, 12, 1, 0, 0, 0, 0, 0,
0}).st(), vec({30240, 15120, 3360, 420, 30, 1, 0, 0, 0, 0}).st()),
vec({17297280, 8648640, 1995840, 277200, 25200, 1512, 56, 1, 0, 0}).st()),
vec({17643225600, 8821612800, 2075673600, 302702400, 30270240, 2162160, 110880,
3960, 90, 1}).st());
A2=A*A;
P=I;
U=C.at((l)-1, 1) * I;
V=C.at((l)-1, 0) * I;
for (int incr_arma_=(1 = l ? 1 : -1), k=1; k != l+incr_arma_;
k+=incr_arma_) {
P=P*A2;
U=U + C((l)-1, ((2 * k) + 2)-1) * P;
V=V + C((l)-1, ((2 * k) + 1)-1) * P;
}
U=A*U;
X=solve(V - U, V + U);
} else {
s=log2(nA / 5.4);
B=A;
if (s 0) {
s=ceil(s);
B=B / (pow(2, s));
}
c_=vec({6476475253248, 3238237626624, 7771770303897600,
1187353796428800, 129060195264000, 10559470521600,
Re: [Rcpp-devel] tracking volatile bug
Serguei, First, I think we mentioned to you after your first post on rex2arma that all serious Rcpp development is generally done in packages. You should probably embrace that mode too. Many of us prefer it. Second, 'expm' does not need to be reimplemented. See my RcppKalman package (on GitHub only) to see how to access the expm code (from the package of the same name by Goulet, Dutang,Maechler et all -- all serious people) directly from RcppArmadillo. Third, your code calls back into R via Rcpp::Environment. As a first test I'd try to do without. The recommendation usually is minimally self-contained reproducible examples for bug report. I, like many other people, have found that a good number of times the bug becomes apparent while preparing such a report. Maybe someone will have time to go through your somewhat-non-standard in detail. Cheers, 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
Re: [Rcpp-devel] tracking volatile bug
Le 18 nov. 2014 à 14:42, Serguei Sokol serguei.so...@gmail.com a écrit : Hi everybody, I am facing a volatile bug which appears and disappears in identical calls on identical data. My RcppArmadilla code is expm_cpp() (it is obtained by rex2arma tool which I have presented before, cf. attached file). It is compared and benchmarked vs its R counterpart expm.higham() (also attached). Both are computing matrix exponential by a same particular method. The results are identical for R and Rcpp on small and medium sized matrices nxn n=10:100 but on large matrices (n 800) various error messages can interrupt (or not) a normal run of expm_cpp(). Sometimes message says (in French) type 'char' indisponible dans 'eval', I suppose in English it must be unimplemented type 'char' in 'eval' I have seen (thanks to Google) a similar error message in the following snippet: /* call this R command: source(FileName) */ int errorOccurred; SEXP e = lang2(install(source), mkString(FileName)); /* mkChar instead of mkString would lead to this runtime error: * Error in source(FileName) : unimplemented type 'char' in 'eval' */ R_tryEval(e, R_GlobalEnv, errorOccurred); e is not protected here. Does the problem go away if you protect it: SEXP e = PROTECT( lang2(install(source), mkString(FileName))) ; Or more R/C++ idiomatically, which is both nicer and safer. Language e( source, FileName ) ; Romain which suggests that somewhere in Rcpp or RcppArmadillo there is a mkChar() call instead of mkString(). Other times, error message can say something like argument type[1]='x' must be one of 'M','1','O','I','F' or 'E' or argument type[1]='character' must be a one-letter character string This latter message is somewhat volatile per se. The part of message just after type[1]= can be 'character' (as above) or 'method' or 'ANY' etc. That kind of problem usually is related to premature GC. I have found these messages in the Matrix package https://r-forge.r-project.org/scm/viewvc.php/pkg/Matrix/src/Mutils.c?view=markuproot=matrixpathrev=2614 function char La_norm_type(const char *typstr) Seemingly, the argument *typstr is getting corrupted somewhere on the road. It is useless to say that debugging is of no help here. If run in a debugger, the program stops normally with or without error messages cited above. I have also tried to low the level of gcc optimization both in compiling R and the Rcpp code but it didn't help. Anybody has an experience in tracking down similar cases? This example can be run as: library(expm) library(Rcpp) source(expm_cpp.R) source(expm.higham.R) n=1000 As=matrix(rnorm(n*n), n) stopifnot(diff(range(expm.higham(As, balancing=TRUE)-expm_cpp(As, balancing=TRUE))) 1.e-14) The last command may be run several times before an error shows up. Cheers, Serguei. expm_cpp.Rexpm.higham.R___ 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] tracking volatile bug
Le 18/11/2014 16:39, Romain François a écrit : Le 18 nov. 2014 à 14:42, Serguei Sokol serguei.so...@gmail.com a écrit : /* call this R command: source(FileName) */ int errorOccurred; SEXP e = lang2(install(source), mkString(FileName)); /* mkChar instead of mkString would lead to this runtime error: * Error in source(FileName) : unimplemented type 'char' in 'eval' */ R_tryEval(e, R_GlobalEnv, errorOccurred); e is not protected here. Does the problem go away if you protect it: SEXP e = PROTECT( lang2(install(source), mkString(FileName))) ; you mean SEXP e = PROTECT(lang2(install(source), mkChar(expm.higham.R))); (because mkString() had no problem)? No, it does not, the problem is still there. Or more R/C++ idiomatically, which is both nicer and safer. Language e( source, FileName ) ; This works as expected. We can suppose that behind the stage FileName is wrap()-ed here with mkString() and not with mkChar(). That's why it's working. which suggests that somewhere in Rcpp or RcppArmadillo there is a mkChar() call instead of mkString(). Other times, error message can say something like argument type[1]='x' must be one of 'M','1','O','I','F' or 'E' or argument type[1]='character' must be a one-letter character string This latter message is somewhat volatile per se. The part of message just after type[1]= can be 'character' (as above) or 'method' or 'ANY' etc. That kind of problem usually is related to premature GC. I will try gctorture(TRUE) suggested by Martin. Thanks, Serguei. ___ 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] tracking volatile bug
I will try gctorture(TRUE) suggested by Martin. Try running it under valgrind as well. E.g., % R --debugger=valgrind ==15338== Memcheck, a memory error detector ==15338== Copyright (C) 2002-2011, and GNU GPL'd, by Julian Seward et al. ==15338== Using Valgrind-3.7.0 and LibVEX; rerun with -h for copyright info ==15338== Command: /opt/sw/R/R-3.1.1.atlas1/lib/R/bin/exec/R ==15338== ==15338== Conditional jump or move depends on uninitialised value(s) ==15338==at 0x40188B6: index (strchr.S:56) ==15338==by 0x4007C92: expand_dynamic_string_token (dl-load.c:431) ==15338==by 0x40085EE: _dl_map_object (dl-load.c:2501) ==15338==by 0x400188D: map_doit (rtld.c:633) ==15338==by 0x400F175: _dl_catch_error (dl-error.c:178) ... gctorture(TRUE) # call your test code here Bill Dunlap TIBCO Software wdunlap tibco.com On Tue, Nov 18, 2014 at 10:10 AM, Serguei Sokol serguei.so...@gmail.com wrote: Le 18/11/2014 16:39, Romain François a écrit : Le 18 nov. 2014 à 14:42, Serguei Sokol serguei.so...@gmail.com a écrit : /* call this R command: source(FileName) */ int errorOccurred; SEXP e = lang2(install(source), mkString(FileName)); /* mkChar instead of mkString would lead to this runtime error: * Error in source(FileName) : unimplemented type 'char' in 'eval' */ R_tryEval(e, R_GlobalEnv, errorOccurred); e is not protected here. Does the problem go away if you protect it: SEXP e = PROTECT( lang2(install(source), mkString(FileName))) ; you mean SEXP e = PROTECT(lang2(install(source), mkChar(expm.higham.R))); (because mkString() had no problem)? No, it does not, the problem is still there. Or more R/C++ idiomatically, which is both nicer and safer. Language e( source, FileName ) ; This works as expected. We can suppose that behind the stage FileName is wrap()-ed here with mkString() and not with mkChar(). That's why it's working. which suggests that somewhere in Rcpp or RcppArmadillo there is a mkChar() call instead of mkString(). Other times, error message can say something like argument type[1]='x' must be one of 'M','1','O','I','F' or 'E' or argument type[1]='character' must be a one-letter character string This latter message is somewhat volatile per se. The part of message just after type[1]= can be 'character' (as above) or 'method' or 'ANY' etc. That kind of problem usually is related to premature GC. I will try gctorture(TRUE) suggested by Martin. Thanks, Serguei. ___ 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
[Rcpp-devel] Bug in loadRcppClass/loadModule?
Dear list,
In the process of writing a comprehensive unit testing application for Rcpp
I may have come across a bug in the code. It seems to me that the following
block should execute just fine, exporting the C++ class norm to the
global environment:
require(inline)
require(Rcpp)
inc - '
using namespace Rcpp;
double norm( double x, double y ) {
return sqrt( x*x + y*y );
}
RCPP_MODULE(mod) {
function( norm, norm );
}'
fx - cxxfunction(signature(), plugin=Rcpp, include=inc)
mod - Module(mod, getDynLib(fx))
loadRcppClass('norm', 'norm', mod)
What happens, though, is that the following error is returned:
Error in as.environment(pos) :
no item called moduleName on the search list
D
igging in to
loadRcppClass
, I find that the function fails at the line:
mod - loadModule(module, NULL, env = where, loadNow = TRUE)
Entering
loadModule
, the function fails here, at the get statement:
loadM - as.environment(module)
module - get(loadM, moduleName)
Isn't this backward? get syntax is:
get(x, pos = -1, envir = as.environment(pos), mode = any,
inherits = TRUE)
Where x is the object sought in the specified environment. In this case,
the function is failing because it can't find moduleName in the
environment's search list, but the reason for this seems to be that the
current statement is search for an environment within a character string,
rather than a character string representing a named object within an
environment.
Is this in need of a patch, or am I missing something obvious?
Cheers,
Aaron
___
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
