Update: I did eventually discover an error in the C code which is probably
ultimately responsible for the bug (the variable eff_size is allowed to get
too large, and overrun the end of lev2[]).  However the behaviour of R (or
Rstudio) in response to this is still somewhat mystifying to me!

Robin

On 21 May 2013 14:10, Robin Evans <rje42@cam.ac.uk> wrote:

> 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.weylandt@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 <rje42@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 <rje42@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
>

	[[alternative HTML version deleted]]

