Sure! C code is below if it helps. The gist is that the function
oneMargin forms two matrices M and C, mostly by repeatedly taking
Kronecker products.
Robin
void kronecker (int *A, int *B, int *dima, int *dimb, int *C) {
int k = 0;
int i1,i2,j1,j2;
for (i2 = 0; i2 < dima[1]; i2++) {
for (j2 = 0; j2 < dimb[1]; j2++) {
for (i1 = 0; i1 < dima[0]; i1++) {
for (j1 = 0; j1 < dimb[0]; j1++) {
C[k] = A[i1 + dima[0]*i2]*B[j1 + dimb[0]*j2];
k++;
}}}}
}
void iterate (int *x, int *lev) {
bool ok = FALSE;
int i=0;
do {
if(x[i] < lev[i]-1) {
x[i]++;
ok = TRUE;
}
else {
x[i] = 0;
}
i++;
}
while (!ok);
}
void oneMargin(int *Mar, int *Eff, int *neff, int *lev, int *nvar, int
*M, int *C) {
int *M2, *mult, *Cj, *Cj2;
int i,j=1,k,l,m;
/* determine size of M to output */
for (i = 0; i < nvar[0]; i++) {
j *= (Mar[i] == 1 ? lev[i]*lev[i] : lev[i]);
}
M2 = malloc(sizeof(int)*j);
mult = malloc(sizeof(int)*j);
M[0] = 1;
int M_size[2] = {1,1}, mult_size[2], Cj_size[2], C_size[2] = {0,0};
for (i = nvar[0]-1; i >= 0; i--) {
if (Mar[i] == 1) {
for (j = 0; j < lev[i]*lev[i]; j++) {
mult[j] = ((j % (lev[i]+1)) == 0) ? 1 : 0;
}
mult_size[0] = mult_size[1] = lev[i];
}
else {
for (j = 0; j < lev[i]; j++) {
mult[j] = 1;
}
mult_size[0] = 1;
mult_size[1] = lev[i];
}
kronecker(M, mult, M_size, mult_size, M2);
M_size[0] *= mult_size[0];
M_size[1] *= mult_size[1];
/* copy M2 over to M */
for (j=0; j < M_size[0]*M_size[1]; j++) {
M[j] = M2[j];
}
}
free(M2);
/* Create C matrix */
int index[nvar[0]];
int prod;
int eff_size = 0;
/*swapped = 0; */
mult_size[0] = 1;
C_size[0] = 1;
for (j = 0; j < nvar[0]; j++) {
C_size[0] *= (Mar[j] == 1 ? lev[j] : 1);
}
Cj = malloc(sizeof(int)*C_size[0]);
Cj2 = malloc(sizeof(int)*C_size[0]);
int *lev2;
lev2 = malloc(sizeof(int)*nvar[0]);
/* loop over effects */
for (i = 0; i < neff[0]; i++) {
prod = 1;
for (j = 0; j < nvar[0]; j++) {
if (Eff[j + i*nvar[0]]) {
prod *= lev[j]-1;
lev2[eff_size] = lev[j]-1;
index[eff_size] = 0;
eff_size++;
}
}
/* loop over states of effect (excluding corner points) */
for (j = 0; j < prod; j++) {
Cj[0] = 1;
Cj_size[0] = Cj_size[1] = 1;
k = eff_size;
/* loop over variables */
for (l = nvar[0]-1; l >= 0; l--) {
/* skip variables not in margin */
if (Mar[l] == 0) continue;
mult_size[1] = lev[l];
/* kronecker factor depends on whether or not variable is in effect */
if (Eff[i*nvar[0]+l] == 0) {
/* for variables not in effect, just repeat matrix */
for (m = 0; m < lev[l]; m++) {
mult[m] = 1;
}
}
else {
/* otherwise multiply based on state */
k--;
for (m = 0; m < lev[l]; m++) {
mult[m] = (m == index[k]) ? lev[l] - 1 : -1;
}
}
kronecker(Cj, mult, Cj_size, mult_size, Cj2);
Cj_size[1] *= mult_size[1];
/* copy Cj2 over to Cj */
for (k=0; k < Cj_size[0]*Cj_size[1]; k++) {
Cj[k] = Cj2[k];
}
if (Cj_size[0]*Cj_size[1] > C_size[0]) Rprintf("pointer screwup");
}
/* copy row over to C */
for (m = 0; m < C_size[0]; m++) C[C_size[0]*C_size[1]+m] = Cj[m];
C_size[1] += 1;
iterate(index, lev2);
}
}
free(lev2);
free(Cj);
free(Cj2);
free(mult);
}
On 21 May 2013 14:04, R. Michael Weylandt <michael.weyla...@gmail.com> wrote:
>
> It might also help if you can point us to the C code to help debug.
>
> MW
>
> On Tue, May 21, 2013 at 10:53 AM, Robin Evans <rj...@cam.ac.uk> wrote:
> > I should add to this that I'm running on Scientific Linux 6. I later
> > noticed that the bug only seems to occur when I run the code from Rstudio,
> > and not if I use the terminal directly, so this may be the key to the
> > problem.
> >
> > Robin
> >
> > On 20 May 2013 16:12, Robin Evans <rj...@cam.ac.uk> wrote:
> >
> >> Hello,
> >>
> >> I've run into a problem which is both maddening and rather hard to
> >> replicate, therefore I wondered if someone might know of a plausible
> >> explanation for it. I couldn't find anything in the archives, though
> >> maybe I'm searching for the wrong thing.
> >>
> >> I'm calling some C code using .C, and get the vector I'm interested in
> >> back as the 7th location in the returned list. However I find that if
> >> I try to inspect the contents of this entry in the list in some ways,
> >> I get one answer, and if I look at it in others I get a different
> >> answer. It's quite possible that there's something wrong with the C
> >> code, but this doesn't seem to explain why this phenomenon would occur
> >> in R.
> >>
> >> The problem does not always occur - I have to run the code a few times
> >> and then call the console when it does, but the commands below show
> >> what can happen when it does. I apologise for not being able to get a
> >> cleaner example. Full code and output is below, but here is a
> >> stylised version:
> >>
> >> The following all give one answer (which is the wrong answer as far as
> >> I'm concerned) :
> >> * printing the whole list :
> >> .C(...) # and looking at the 7th entry
> >> * applying c() to the 7th element of the list
> >> c(.C(...)[[7]])
> >> * assigning the 7th element to a vector:
> >> x = .C(...)[[7]];
> >> x
> >>
> >> these give a different answer (which is the answer I would hope the C
> >> code returns):
> >> * using dput on the 7th entry:
> >> dput(.C(...)[[7]])
> >> * applying c() and then dput()
> >> dput(c(.C(...)[[7]]))
> >> * just printing the 7th entry of the list
> >> .C(...)[[7]]
> >>
> >> The answers are consistent in the sense that I always get the same
> >> answers from running the same command in the console. I have tried
> >> inspecting the returned objects to see if the objects are somehow of a
> >> different class than I expect, or are just being printed oddly, but
> >> have not found anything untoward.
> >>
> >> Any suggestions would be much appreciated!
> >>
> >> Regards,
> >>
> >> Robin
> >>
> >>
> >> # THESE COMMANDS GIVE ONE ANSWER
> >> # [the correct answer always begins with 1, the incorrect with -1]
> >>
> >> > .C("oneMargin", c(1L,1L,1L,1L,1L), c(1L,1L,1L,1L,1L), 1L,
> >> c(2L,2L,2L,2L,2L), 5L, ptr1, ptr2)[[7]]
> >> [1] 1 -1 -1 1 -1 1 1 -1 -1 1 1 -1 1 -1 -1 1 -1 1 1 -1 1 -1
> >> -1 1 1 -1 -1 1 -1 1 1 -1
> >>
> >> > dput(.C("oneMargin", c(1L,1L,1L,1L,1L), c(1L,1L,1L,1L,1L), 1L,
> >> c(2L,2L,2L,2L,2L), 5L, ptr1, ptr2)[[7]])
> >> c(1L, -1L, -1L, 1L, -1L, 1L, 1L, -1L, -1L, 1L, 1L, -1L, 1L, -1L,
> >> -1L, 1L, -1L, 1L, 1L, -1L, 1L, -1L, -1L, 1L, 1L, -1L, -1L, 1L,
> >> -1L, 1L, 1L, -1L)
> >>
> >> > x=dput(c(.C("oneMargin", c(1L,1L,1L,1L,1L), c(1L,1L,1L,1L,1L), 1L,
> >> c(2L,2L,2L,2L,2L), 5L, ptr1, ptr2)[[7]]))
> >> c(1L, -1L, -1L, 1L, -1L, 1L, 1L, -1L, -1L, 1L, 1L, -1L, 1L, -1L,
> >> -1L, 1L, -1L, 1L, 1L, -1L, 1L, -1L, -1L, 1L, 1L, -1L, -1L, 1L,
> >> -1L, 1L, 1L, -1L)
> >> > x
> >> [1] 1 -1 -1 1 -1 1 1 -1 -1 1 1 -1 1 -1 -1 1 -1 1 1 -1 1 -1
> >> -1 1 1 -1 -1 1 -1 1 1 -1
> >>
> >> # THESE ALL GIVE A DIFFERENT ONE!
> >>
> >> > .C("oneMargin", c(1L,1L,1L,1L,1L), c(1L,1L,1L,1L,1L), 1L,
> >> c(2L,2L,2L,2L,2L), 5L, ptr1, ptr2)
> >>
> >> # (OTHER ELEMENTS OF LIST REMOVED)
> >> [[7]]
> >> [1] -1 -1 1 1 1 1 -1 -1 1 1 -1 -1 -1 -1 1 1 1 1 -1 -1 -1 -1
> >> 1 1 -1 -1 1 1 1 1 -1 -1
> >>
> >> > c(.C("oneMargin", c(1L,1L,1L,1L,1L), c(1L,1L,1L,1L,1L), 1L,
> >> c(2L,2L,2L,2L,2L), 5L, ptr1, ptr2)[[7]])
> >> [1] -1 -1 1 1 1 1 -1 -1 1 1 -1 -1 -1 -1 1 1 1 1 -1 -1 -1 -1
> >> 1 1 -1 -1 1 1 1 1 -1 -1
> >> > x = .C("oneMargin", c(1L,1L,1L,1L,1L), c(1L,1L,1L,1L,1L), 1L,
> >> c(2L,2L,2L,2L,2L), 5L, ptr1, ptr2)[[7]]
> >> > x
> >> [1] -1 -1 1 1 1 1 -1 -1 1 1 -1 -1 -1 -1 1 1 1 1 -1 -1 -1 -1
> >> 1 1 -1 -1 1 1 1 1 -1 -1
> >>
> >>
> >> --
> >> Robin Evans
> >> Statistical Laboratory
> >> University of Cambridge
> >> blog: itsastatlife.blogspot.com
> >> web: www.statslab.cam.ac.uk/~rje42
> >>
> >> Causal Inference Workshop July 15th:
> >> http://www.statslab.cam.ac.uk/~rje42/uai13/main.htm
> >>
> >
> >
> >
> > --
> > Robin Evans
> > Statistical Laboratory
> > University of Cambridge
> > blog: itsastatlife.blogspot.com
> > web: www.statslab.cam.ac.uk/~rje42 <http://www.stat.washington.edu/~rje42>
> >
> > Causal Inference Workshop July 15th:
> > http://www.statslab.cam.ac.uk/~rje42/uai13/main.htm
> >
> > [[alternative HTML version deleted]]
> >
> > ______________________________________________
> > R-devel@r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-devel
______________________________________________
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
