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]]
SEXP eta = PROTECT(allocVector(REALSXP,H_c)); n_prot++;
double *eta_c; eta_c = REAL(eta);
for (i=0;i<N_c;i++)
{
start_c = i * H_c;
for (j=0; j<H_c; j++)
{
eta_c[start_c + j] = ...
It looks you expect to be able to write into the N_c*H_c element
of eta, but you only allocated H_c elements for it.
Bill Dunlap
TIBCO Software
wdunlap tibco.com
On Fri, Feb 2, 2018 at 6:22 AM, Jesper Hybel Pedersen <jeshyb at dtu.dk>
wrote:
> 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]]
>
> ______________________________________________
> R-devel at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
>
[[alternative HTML version deleted]]