Hi
I'm trying to develop some C code to find the fixpoint of a contraction mapping, the code compiles and gives the right results when executed in R.
However R-gui session is frequently terminated. I'm suspecting some access violation error due to the exception code 0xc0000005
In the error report windows 10 gives me.
It is the first time I'm writing any C-code so I'm guessing I have done something really stupid. I have been trying to debug the code for a couple of days now,
But I simply can't figure out what generates the problem. Could it be something particular to my windows set up and security stuff?
I'm in the process of reading Writing R Extensions and Hadley Wickham's Advanced R but might have missed something.
The windows error report:
Faulting application name: Rgui.exe, version: 3.33.6774.0, time stamp: 0x58bd6d26
Faulting module name: R.dll, version: 3.33.6774.0, time stamp: 0x58bd6d0b
Exception code: 0xc0000005
Fault offset: 0x000000000010b273
Faulting process id: 0x1d14
Faulting application start time: 0x01d39aede45c96e9
Faulting application path: C:\Program Files\R\R-3.3.3\bin\x64\Rgui.exe
Faulting module path: C:\Program Files\R\R-3.3.3\bin\x64\R.dll
Report Id: c78d7c52-72c5-40f3-a3cc-927323d2af07
Faulting package full name:
Faulting package-relative application ID:
####### How I call the C-function in R
dyn.load("C://users//jeshyb//desktop//myC//lce_fixpoint_cc.dll")
N = 10
H = 3
v <- rnorm(N*H)
mu <- 0.1
psi <- matrix(c(1,0,0.5,0.5,0,1),nrow=2)
K <- dim(psi)[1]
p <- rep(1/H,N*H)
error <- 1e-10
f<-function(p,v,mu,psi,N,H,K)
{
.Call("lce_fixpoint_cc",p, v, mu, psi, as.integer(N), as.integer(H), as.integer(K),error)
}
for (i in 1:100)
{
v <- rnorm(N*H)
p <- rep(1/H,N*H)
a<-f(p,v,mu,psi,N,H,K)
}
a<-f(p,v,mu,psi,N,H,K)
plot(a)
######## The C- function
#include <R.h>
#include <Rinternals.h>
SEXP lce_fixpoint_cc(SEXP q, SEXP v, SEXP mu, SEXP psi, SEXP N,SEXP H, SEXP K, SEXP err)
{
int n_prot = 0;
/* Make ready integer and double constants */
PROTECT(N); n_prot++;
PROTECT(H); n_prot++;
PROTECT(K); n_prot++;
int N_c = asInteger(N);
int H_c = asInteger(H);
int K_c = asInteger(K);
double mu_c = asReal(mu);
double mu2_c = 1.0 - mu_c;
double error_c = asReal(err);
double lowest_double = 1e-15;
double tmp_c;
double denom;
double error_temp;
double error_i_c;
/* Make ready vector froms input */
PROTECT(q); n_prot++;
PROTECT(v); n_prot++;
PROTECT(psi); n_prot++;
double *v_c; v_c = REAL(v);
double *psi_c; psi_c = REAL(psi);
/* Initialize new vectors not given as input */
SEXP q_copy = PROTECT(duplicate(q)); n_prot++;
double *q_c; q_c = REAL(q_copy);
SEXP q_new = PROTECT(allocVector(REALSXP,length(q))); n_prot++;
double *q_new_c; q_new_c = REAL(q_new);
SEXP eta = PROTECT(allocVector(REALSXP,H_c)); n_prot++;
double *eta_c; eta_c = REAL(eta);
SEXP exp_eta = PROTECT(allocVector(REALSXP,H_c)); n_prot++;
double *exp_eta_c; exp_eta_c = REAL(exp_eta);
SEXP psi_ln_psiq = PROTECT(allocVector(REALSXP,H_c)); n_prot++;
double *psi_ln_psiq_c; psi_ln_psiq_c = REAL(psi_ln_psiq);
int not_converged;
int maxIter = 10000;
int iter;
int start_c;
/* loop indeces */
int i;
int j;
int k;
/* loop over observational units to find choice probabilities for i=1,...,N */
for (i=0;i<N_c;i++)
{
start_c = i * H_c;
not_converged = 1;
iter = 0;
while(iter < maxIter && not_converged)
{
/* make v_ij + (1-mu)*ln(q_ij) */
for (j=0; j<H_c; j++)
{
eta_c[start_c + j] = v_c[start_c + j] + mu2_c * log(q_c[start_c + j]);
psi_ln_psiq_c[j] = 0.0;
}
/* Make psi_ln_psiq_c vector for individual i */
for (k=0;k<K_c;k++)
{
tmp_c = 0.0;
/* Calculate row k of psi %*% q */
for (j=0;j<H_c;j++)
{
tmp_c += psi_c[k + j*K_c] * q_c[start_c +j];
}
tmp_c = mu2_c * log(tmp_c);
for (j=0;j<H_c;j++)
{
psi_ln_psiq_c[j] += psi_c[k + j*K_c] * tmp_c;
}
}
denom = 0.0;
for (j=0;j<H_c;j++)
{
exp_eta_c[start_c + j] = exp( eta_c[start_c + j] - psi_ln_psiq_c[j] ) + lowest_double;
denom += exp_eta_c[start_c + j];
}
error_i_c = 0.0;
error_temp = 0.0;
/* calculate error and update choice prob */
for (j=0;j<H_c;j++)
{
q_new_c[start_c + j] = exp_eta_c[start_c + j]/denom;
error_temp = fabs(q_new_c[start_c + j] - q_c[start_c + j]);
if (error_temp>error_i_c)
{
error_i_c = error_temp;
}
q_c[start_c + j] = q_new_c[start_c + j];
}
not_converged = error_i_c > error_c;
iter++;
} /* End while loop for individual i to solve for q_i */
} /* End loop over individuals */
UNPROTECT(n_prot);
return(q_new);
}
Best regards
Jesper Hybel
[[alternative HTML version deleted]]