Skip to content

Commit adac56e

Browse files
eddelbuetteleddelbuettel
committed
Add files used in issue #474 
1 parent 41ea652 commit adac56e

6 files changed

+698
-0
lines changed

‎local/issue474_ebsc_rank.cpp

Lines changed: 115 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,115 @@
1+
2+
// The eBsc package on CRAN (at version 4.17) reports an error under the Intel compiler.
3+
// This error now also appears with RcppArmadillo 15.0.0. In the example for plot.eBsc,
4+
// the function drbasis() is called with n=250 and a second parameter p that varies from
5+
// 1 to 6. In the case of 5, NA values return leading an rank failure. The drbasis()
6+
// function calls an underlying C function drbasis from file drbasisC.cpp which has case
7+
// statements for the differenct values of p. We extrace the one for 5 here.
8+
9+
// [[Rcpp::depends(RcppArmadillo)]]
10+
11+
#include <RcppArmadillo/RcppArmadillo>
12+
13+
// [[Rcpp::export]]
14+
arma::mat compute(int nn) {
15+
using namespace arma;
16+
using namespace std;
17+
//Rcpp::List out;
18+
const double E = exp(1);
19+
arma::vec t(nn);
20+
arma::vec k(nn);
21+
typedef complex<double> dcomp;
22+
double Pi = M_PI;
23+
complex<double> comp0 (0,1);
24+
complex<double> comp1 (-1,0);
25+
complex<double> comp2 (0,-0.5);
26+
complex<double> comp3 (1,-1);
27+
complex<double> comp4 (-1,1);
28+
complex<double> comp5 (0,0.5);
29+
complex<double> comp6 (0,0.25);
30+
int const0 = 1;
31+
double const1 = -sqrt(3);
32+
double const2 = 2*sqrt(3);
33+
double const3 = -sqrt(5);
34+
double const4 = 6*sqrt(5);
35+
double const5 = -6*sqrt(5);
36+
double const6 = -sqrt(7);
37+
double const7 = 12*sqrt(7);
38+
double const8 = -30*sqrt(7);
39+
double const9 = 20*sqrt(7);
40+
double const10 = 3;
41+
double const11 = -60;
42+
double const12 = 270;
43+
double const13 = -420;
44+
double const14 = 210;
45+
double const15 = -sqrt(11);
46+
double const16 = 30*sqrt(11);
47+
double const17 = -210*sqrt(11);
48+
double const18 = 560*sqrt(11);
49+
double const19 = 630*sqrt(11);
50+
double const20 = 252*sqrt(11);
51+
52+
for (int j=0;j<nn;j++) {
53+
t(j)= j/double((nn-1));
54+
}
55+
for (int j=0;j<nn;j++) {
56+
k(j)= j+1;
57+
}
58+
59+
60+
arma::mat M5(nn,nn);
61+
for (int i=0;i<nn; i++) {
62+
for (int j=0;j<nn;j++) {
63+
if (j==0) {
64+
M5(i,j)=const0;
65+
} else if (j==1) {
66+
M5(i,j)= const1 + const2*t(i);
67+
} else if (j==2) {
68+
M5(i,j)= const3 + const4*t(i) + const5* pow(t(i),2);
69+
} else if (j==3) {
70+
M5(i,j)= const6 + const7*t(i) + const8* pow(t(i),2) + const9* pow(t(i),3);
71+
} else if (j==4) {
72+
M5(i,j)= const10 + const11*t(i) + const12* pow(t(i),2) + const13* pow(t(i),3) + const14* pow(t(i),4);
73+
} else {
74+
dcomp ans = (sqrt(2)*(1 + sqrt(5)/2) -
75+
(comp5*(sqrt(10 - 2*sqrt(5)) +
76+
sqrt(2*(5 +
77+
sqrt(5)))))/sqrt(2))*(pow(comp1,1 + k(j))/pow(E,pow(comp1,0.1)* (dcomp(-3 + k(j)))*Pi*(1 - t(i))) +
78+
pow(E,-(pow(comp1,0.1)*(dcomp(-3 + k(j)))* Pi*t(i)))) +
79+
(-(1/sqrt(2)) -
80+
(comp5*(sqrt(10 - 2*sqrt(5)) + sqrt(2*(5 + sqrt(5)))))/
81+
sqrt(2))*(pow(comp1,1 + k(j))/pow(E,pow(comp1,0.3)*(dcomp (-3 + k(j)))*Pi*(1 - t(i))) + pow(E,-(pow(comp1,0.3)*(dcomp (-3 + k(j)))*Pi*t(i)))) +
82+
(sqrt(2)*(1 + sqrt(5)/2) + (comp5*(sqrt(10 - 2*sqrt(5)) + sqrt(2*(5 + sqrt(5)))))/sqrt(2))
83+
*(pow(comp1,1 + k(j))/ pow(E,(sqrt(0.625 + sqrt(5)/8) - comp6*(-1 + sqrt(5)))*(dcomp (-3 + k(j)))*Pi*
84+
(1 - t(i))) +pow(E,-((sqrt(0.625 + sqrt(5)/8.) - comp6*(-1 + sqrt(5)))*(dcomp (-3 + k(j)))*
85+
Pi*t(i)))) +(-(1/sqrt(2)) + (comp5*(sqrt(10 - 2*sqrt(5)) + sqrt(2*(5 + sqrt(5)))))/sqrt(2))
86+
*(pow(comp1,1 + k(j))/pow(E,(dcomp(sqrt(0.625 - sqrt(5)/8)) - comp6*(1 + sqrt(5)))*
87+
(dcomp (-3 + k(j)))*Pi*(1 - t(i))) +pow(E,-((sqrt(0.625 - sqrt(5)/8) - comp6*(1 + sqrt(5)))* (dcomp (-3 + k(j)))*Pi*t(i))))- sqrt(2)*cos((-3 + k(j))*Pi*t(i));
88+
M5(i,j) =ans.real();
89+
}
90+
M5(i,j) = M5(i,j)/sqrt(nn);
91+
}
92+
}
93+
/*
94+
arma::mat vectors=M5;
95+
arma::mat vectorsQR=Turn(M5); arma::vec values=eigenvalues(nn,5);
96+
out = Rcpp::List::create(Rcpp::Named("eigenvectors") = vectors,
97+
Rcpp::Named("eigenvectorsQR") = vectorsQR,
98+
Rcpp::Named("eigenvalues") = values,
99+
Rcpp::Named("x") = t) ;
100+
return out;
101+
*/
102+
return M5;
103+
}
104+
105+
/*** R
106+
res <- compute(250)
107+
if (!any(is.na(res))) {
108+
cat("Kappa: ")
109+
kappa(res)
110+
} else {
111+
cat("Nb NAs: ")
112+
sum(is.na(res))
113+
print(res[1:10,241:250])
114+
}
115+
*/

‎local/issue474_ebsc_rank_full.cpp

Lines changed: 193 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,193 @@
1+
2+
// The eBsc package on CRAN (at version 4.17) reports an error under the Intel compiler.
3+
// This error now also appears with RcppArmadillo 15.0.0. In the example for plot.eBsc,
4+
// the function drbasis() is called with n=250 and a second parameter p that varies from
5+
// 1 to 6. In the case of 5, NA values return leading an rank failure. The drbasis()
6+
// function calls an underlying C function drbasis from file drbasisC.cpp which has case
7+
// statements for the differenct values of p. We extrace the one for 5 here.
8+
9+
// [[Rcpp::depends(RcppArmadillo)]]
10+
11+
#include <RcppArmadillo/RcppArmadillo>
12+
13+
// three helpers used below
14+
double Qfunc(double x, int qq=999) {
15+
using namespace arma;
16+
using namespace std;
17+
double Qf=1;
18+
double pi = M_PI;
19+
20+
switch (qq)
21+
{
22+
case 1:
23+
Qf = 1;
24+
break;
25+
case 2:
26+
Qf =1.0/3.0 + 2*pow(cos(pi*x),2.0)/3.0;
27+
break;
28+
case 3:
29+
Qf = 2.0/15.0 + 11*pow(cos(pi*x),2.0)/15.0 + 2*pow(cos(pi*x),4.0)/15.0;
30+
break;
31+
case 4:
32+
Qf =(17+180*pow(cos(pi*x),2)+114*pow(cos(pi*x),4)+4*pow(cos(pi*x),6))/315.0;
33+
break;
34+
case 5:
35+
Qf =(62+1072*pow(cos(pi*x),2)+1452*pow(cos(pi*x),4)+247*pow(cos(pi*x),6)+2*pow(cos(pi*x),8))/2835.0;
36+
break;
37+
case 6:
38+
Qf =(1382+35396*pow(cos(pi*x),2)+83021*pow(cos(pi*x),4)+34096*pow(cos(pi*x),6)+2026*pow(cos(pi*x),8)+4*pow(cos(pi*x),10))/155295.0;
39+
break;
40+
}
41+
42+
return Qf;
43+
}
44+
45+
// [[Rcpp::export]]
46+
arma::vec eigenvalues(int nn, int qq=999){
47+
using namespace arma;
48+
double pi = M_PI;
49+
arma::vec svec(nn-qq);
50+
arma::vec atten(nn-qq);
51+
arma::vec atten1(2*nn);
52+
arma::vec ev(nn);
53+
for (int i=0;i<nn-qq;i++){
54+
svec(i)=pow(pi*((i+1)+1.0/2*((qq+1)%2)+floor((qq-1)/2.0)),(2*qq))/nn;}
55+
for (int j=0;j<2*nn;j++){
56+
atten1(j)= pow(sin((j+1)*pi/(2*nn))/((j+1)*pi/(2*nn)),2*qq) / Qfunc((j+1)/(2.0*nn),qq);}
57+
for (int j=0;j<nn-qq;j++){
58+
atten(j)= atten1(j+qq);}
59+
svec = svec%atten;
60+
ev = zeros<vec>(nn);
61+
for (int i=qq;i<nn;i++) {
62+
ev(i)=svec(i-qq);
63+
}
64+
return ev;
65+
}
66+
67+
// [[Rcpp::export]]
68+
arma::mat Turn(arma::mat M){
69+
int nn =arma::rank(M);
70+
arma::mat Mo,R,CMo;
71+
arma::qr(Mo,R,M);
72+
arma::mat A=M+Mo;
73+
arma::mat B=M-Mo;
74+
arma::vec a(nn),b(nn);//contain the norm of the coloumns of A and B matrices
75+
arma::vec x(nn),y(nn); //auxillary vectors
76+
arma::mat D(nn,2);// same as the R code
77+
arma::vec d(nn);// same as the R code
78+
arma::mat dd; // diag(d)
79+
for(int i=0;i<nn;i++){
80+
for(int j=0;j<nn;j++){
81+
x(j)=A(j,i);
82+
y(j)=B(i,j);
83+
}
84+
a(i) = norm(x,2);
85+
b(i) = norm(y,2);
86+
}
87+
88+
for(int i=0;i<nn;i++){
89+
D(i,0)=a(i);
90+
D(i,1)=b(i);
91+
}
92+
for(int i=0;i<nn;i++){
93+
if(D(i,0) <= D(i,1)){d(i)=-1;}
94+
else d(i)=1;
95+
}
96+
dd=diagmat(d);
97+
CMo = Mo*dd;
98+
return CMo;
99+
}
100+
101+
102+
// [[Rcpp::export]]
103+
Rcpp::List compute(int nn) {
104+
using namespace arma;
105+
using namespace std;
106+
Rcpp::List out;
107+
const double E = exp(1);
108+
arma::vec t(nn);
109+
arma::vec k(nn);
110+
typedef complex<double> dcomp;
111+
double Pi = M_PI;
112+
complex<double> comp0 (0,1);
113+
complex<double> comp1 (-1,0);
114+
complex<double> comp2 (0,-0.5);
115+
complex<double> comp3 (1,-1);
116+
complex<double> comp4 (-1,1);
117+
complex<double> comp5 (0,0.5);
118+
complex<double> comp6 (0,0.25);
119+
int const0 = 1;
120+
double const1 = -sqrt(3);
121+
double const2 = 2*sqrt(3);
122+
double const3 = -sqrt(5);
123+
double const4 = 6*sqrt(5);
124+
double const5 = -6*sqrt(5);
125+
double const6 = -sqrt(7);
126+
double const7 = 12*sqrt(7);
127+
double const8 = -30*sqrt(7);
128+
double const9 = 20*sqrt(7);
129+
double const10 = 3;
130+
double const11 = -60;
131+
double const12 = 270;
132+
double const13 = -420;
133+
double const14 = 210;
134+
double const15 = -sqrt(11);
135+
double const16 = 30*sqrt(11);
136+
double const17 = -210*sqrt(11);
137+
double const18 = 560*sqrt(11);
138+
double const19 = 630*sqrt(11);
139+
double const20 = 252*sqrt(11);
140+
141+
for (int j=0;j<nn;j++) {
142+
t(j)= j/double((nn-1));
143+
}
144+
for (int j=0;j<nn;j++) {
145+
k(j)= j+1;
146+
}
147+
148+
149+
arma::mat M5(nn,nn);
150+
for (int i=0;i<nn; i++) {
151+
for (int j=0;j<nn;j++) {
152+
if (j==0) {
153+
M5(i,j)=const0;
154+
} else if (j==1) {
155+
M5(i,j)= const1 + const2*t(i);
156+
} else if (j==2) {
157+
M5(i,j)= const3 + const4*t(i) + const5* pow(t(i),2);
158+
} else if (j==3) {
159+
M5(i,j)= const6 + const7*t(i) + const8* pow(t(i),2) + const9* pow(t(i),3);
160+
} else if (j==4) {
161+
M5(i,j)= const10 + const11*t(i) + const12* pow(t(i),2) + const13* pow(t(i),3) + const14* pow(t(i),4);
162+
} else {
163+
dcomp ans = (sqrt(2)*(1 + sqrt(5)/2) -
164+
(comp5*(sqrt(10 - 2*sqrt(5)) +
165+
sqrt(2*(5 +
166+
sqrt(5)))))/sqrt(2))*(pow(comp1,1 + k(j))/pow(E,pow(comp1,0.1)* (dcomp(-3 + k(j)))*Pi*(1 - t(i))) +
167+
pow(E,-(pow(comp1,0.1)*(dcomp(-3 + k(j)))* Pi*t(i)))) +
168+
(-(1/sqrt(2)) -
169+
(comp5*(sqrt(10 - 2*sqrt(5)) + sqrt(2*(5 + sqrt(5)))))/
170+
sqrt(2))*(pow(comp1,1 + k(j))/pow(E,pow(comp1,0.3)*(dcomp (-3 + k(j)))*Pi*(1 - t(i))) + pow(E,-(pow(comp1,0.3)*(dcomp (-3 + k(j)))*Pi*t(i)))) +
171+
(sqrt(2)*(1 + sqrt(5)/2) + (comp5*(sqrt(10 - 2*sqrt(5)) + sqrt(2*(5 + sqrt(5)))))/sqrt(2))
172+
*(pow(comp1,1 + k(j))/ pow(E,(sqrt(0.625 + sqrt(5)/8) - comp6*(-1 + sqrt(5)))*(dcomp (-3 + k(j)))*Pi*
173+
(1 - t(i))) +pow(E,-((sqrt(0.625 + sqrt(5)/8.) - comp6*(-1 + sqrt(5)))*(dcomp (-3 + k(j)))*
174+
Pi*t(i)))) +(-(1/sqrt(2)) + (comp5*(sqrt(10 - 2*sqrt(5)) + sqrt(2*(5 + sqrt(5)))))/sqrt(2))
175+
*(pow(comp1,1 + k(j))/pow(E,(dcomp(sqrt(0.625 - sqrt(5)/8)) - comp6*(1 + sqrt(5)))*
176+
(dcomp (-3 + k(j)))*Pi*(1 - t(i))) +pow(E,-((sqrt(0.625 - sqrt(5)/8) - comp6*(1 + sqrt(5)))* (dcomp (-3 + k(j)))*Pi*t(i))))- sqrt(2)*cos((-3 + k(j))*Pi*t(i));
177+
M5(i,j) =ans.real();
178+
}
179+
M5(i,j) = M5(i,j)/sqrt(nn);
180+
}
181+
}
182+
arma::mat vectors=M5; arma::mat vectorsQR=Turn(M5); arma::vec values=eigenvalues(nn,5);
183+
out = Rcpp::List::create(Rcpp::Named("eigenvectors") = vectors,
184+
Rcpp::Named("eigenvectorsQR") = vectorsQR,
185+
Rcpp::Named("eigenvalues") = values,
186+
Rcpp::Named("x") = t) ;
187+
return out;
188+
}
189+
190+
/*** R
191+
res <- compute(250)
192+
summary(res)
193+
*/

0 commit comments

Comments
 (0)