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=markup&root=matrix&pathrev=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;
}
template<typename 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=as<List>(expm_balance_r_(A, \"P\"));
baS=as<List>(expm_balance_r_(as<mat>(baP[\"z\"]), \"S\"));
A=as<mat>(baS[\"z\"]);
}
nA=as<double>(Matrix_norm_r_(A, \"1\"));
I=eye<mat>(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({64764752532480000, 32382376266240000, 7771770303897600,
1187353796428800, 129060195264000, 10559470521600, 670442572800, 33522128640,
1323241920, 40840800, 960960, 16380, 182, 1});
B2=B*B;
B4=B2*B2;
B6=B2*B4;
U=B*(B6*(c_.at(13) * B6 + c_.at(11) * B4 + c_.at(9) * B2) + c_.at(7) * B6
+ c_.at(5) * B4 + c_.at(3) * B2 + c_.at(1) * I);
V=B6*(c_.at(12) * B6 + c_.at(10) * B4 + c_.at(8) * B2) + c_.at(6) * B6 +
c_.at(4) * B4 + c_.at(2) * B2 + c_.at(0) * I;
X=solve(V - U, V + U);
if (s > 0) for (int incr_arma_=(1 <= s ? 1 : -1), t=1; t !=
s+incr_arma_; t+=incr_arma_) X=X*X;
}
if (balancing) {
dd=as<vec>(baS[\"scale\"]);
X=X % mat(vec((as<vec>(rep_r_(dd, n)) % as<vec>(rep_r_(1 / dd,
_[\"each\"]=n)))).begin(), X.n_rows, X.n_cols, false);
pp=conv_to<ivec>::from(as<vec>(baP[\"scale\"]));
if (as<int>(baP[\"i1\"]) > 1) {
for (int incr_arma_=((as<int>(baP[\"i1\"]) - 1) <= 1 ? 1 : -1),
i=(as<int>(baP[\"i1\"]) - 1); i != 1+incr_arma_; i+=incr_arma_) {
tt=X(span(), (i)-1);
X(span(), (i)-1)=X(span(), (pp.at((i)-1))-1);
X(span(), (pp.at((i)-1))-1)=tt;
tt=X((i)-1, span());
X((i)-1, span())=X((pp.at((i)-1))-1, span());
X((pp.at((i)-1))-1, span())=tt;
}
}
if (as<int>(baP[\"i2\"]) < n) {
for (int incr_arma_=((as<int>(baP[\"i2\"]) + 1) <= n ? 1 : -1),
i=(as<int>(baP[\"i2\"]) + 1); i != n+incr_arma_; i+=incr_arma_) {
tt=X(span(), (i)-1);
X(span(), (i)-1)=X(span(), (pp.at((i)-1))-1);
X(span(), (pp.at((i)-1))-1)=tt;
tt=X((i)-1, span());
X((i)-1, span())=X((pp.at((i)-1))-1, span());
X((pp.at((i)-1))-1, span())=tt;
}
}
}
return wrap(X);
}"
)
expm.higham=function (A, balancing = TRUE)
{
d <- dim(A)
if (length(d) != 2 || d[1] != d[2])
stop("'A' must be a square matrix")
n <- d[1]
if (n <= 1)
return(exp(A))
if (balancing) {
baP <- balance(A, "P")
baS <- balance(baP$z, "S")
A <- baS$z
}
nA <- Matrix::norm(A, "1")
#print(nA);
I <- diag(n)
if (nA <= 2.1) {
t <- c(0.015, 0.25, 0.95, 2.1)
l <- which.max(nA <= t)
C <- rbind(c(120, 60, 12, 1, 0, 0, 0, 0, 0, 0), c(30240,
15120, 3360, 420, 30, 1, 0, 0, 0, 0), c(17297280,
8648640, 1995840, 277200, 25200, 1512, 56, 1, 0,
0), c(17643225600, 8821612800, 2075673600, 302702400,
30270240, 2162160, 110880, 3960, 90, 1))
#print(C);
A2 <- A %*% A
P <- I
U <- C[l, 2] * I
V <- C[l, 1] * I
for (k in 1:l) {
P <- P %*% A2
U <- U + C[l, (2 * k) + 2] * P
V <- V + C[l, (2 * k) + 1] * P
}
U <- A %*% U
X <- solve(V - U, V + U)
}
else {
s <- log2(nA/5.4)
B <- A
if (s > 0) {
s <- ceiling(s)
B <- B/(2^s)
}
#print(B)
c. <- c(64764752532480000, 32382376266240000, 7771770303897600,
1187353796428800, 129060195264000, 10559470521600,
670442572800, 33522128640, 1323241920, 40840800,
960960, 16380, 182, 1)
B2 <- B %*% B
B4 <- B2 %*% B2
B6 <- B2 %*% B4
U <- B %*% (B6 %*% (c.[14] * B6 + c.[12] * B4 + c.[10] *
B2) + c.[8] * B6 + c.[6] * B4 + c.[4] * B2 + c.[2] *
I)
V <- B6 %*% (c.[13] * B6 + c.[11] * B4 + c.[9] * B2) +
c.[7] * B6 + c.[5] * B4 + c.[3] * B2 + c.[1] * I
X <- solve(V - U, V + U)
#print(U)
#print(V)
#print(X)
#print(s)
if (s > 0)
for (t in 1:s) X <- X %*% X
#print(X)
}
if (balancing) {
dd <- baS$scale
X <- X * (rep(dd, n) * rep(1/dd, each = n))
pp <- as.integer(baP$scale)
if (baP$i1 > 1) {
for (i in (baP$i1 - 1):1) {
tt <- X[, i]
X[, i] <- X[, pp[i]]
X[, pp[i]] <- tt
tt <- X[i, ]
X[i, ] <- X[pp[i], ]
X[pp[i], ] <- tt
}
}
if (baP$i2 < n) {
for (i in (baP$i2 + 1):n) {
tt <- X[, i]
X[, i] <- X[, pp[i]]
X[, pp[i]] <- tt
tt <- X[i, ]
X[i, ] <- X[pp[i], ]
X[pp[i], ] <- tt
}
}
}
X
}
_______________________________________________
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
