[Rcpp-devel] R.e. Fwd: problem with rmultinom function

[Christian Gunning]

On Tue, May 31, 2011 at 10:55 AM,
<rcpp-devel-request at r-forge.wu-wien.ac.at> wrote:
> Hi,
> I sent this message a couple of days ago, but I haven't heard anything yet.
> I think you didnt receive it, that's why I send it again.

I didn't see anything.  If you're not subscribed to the list,
submissions silently fail.  Plus, you must send from the subscribing
address.  I've made this mistake a few times...

> I have a strange problem that I spent several days trying to fix it, but
> unfortunately I wasn't able to catch it. I'm running the rcpp code for 3000
> iterations(running the same code for 3000 times) on different datasets. but
> for some datasets, in the final iterations (like 2667 th iteration) the code
> stops running and it jumps out of the R environment. and if I run the same
> code (without any changes) on the same dataset(still without any changes)
> the problem might never happen and it goes to the last iteration without any
> problem(or the problem may occur in another iteration). and the interesting
> part is that error code is different each time, I got error codes like:
>
> what():  invalid partial string match
> what():  'prob' is missing
> what():  unused argument(s) (<environment>)
> and etc...
> I'm pretty sure that problem occurs in the line of code that I am using
> rmultinom function, but I included the necessary function
> Function rmultinom = stats["rmultinom"];
>
> to be more clear I wrote a simple function and I attached it in this email.
> in this code I didn't read any dataset, I just used some variables to define
> the size of structures that I need to get from dataset. if you run the code
> with these values for P, K, unum,qnum, the error occurs in the first 10
> iterations, but If you change the values of those parameters to smaller
> values, the error may never happen.
> I tried to used an cpp implementation for rmultinom, it works but the
> program would be as slow as R version.
>
>
> do you have any idea why this problem occur? or did anybody ever have the
> same problem or can give me some tips to fix it?
>
> ps: all datasets are generated through the sampe process(with a single
> simulation code)

I'm attaching a cleanup of the code using inline -- you can now just
source this in R.  I'm not sure, but I think you might be working
under some mistaken assumptions.  Hopefully the attached code is more
conceptually clear.

There were a few offhand things that I modified, e.g. maxItr wasn't
used in the attached code.  There are a number of things that I
*didn't* change for which there seem to be faster/better methods, but
I'm not sure of your purpose.  For example, can any of the in-loop
assignments (qnum, qlength) be moved to the top?

Lastly, are you *just* using rmultinom and Cwhich to roll a K or P
sided die?  If so, what about something like as<int>(runif(1, 0, K))?
Even if you're rolling a weighted die, it seems like there's a better
way...

Hopefully this sheds at least a bit of light on your question.
Apologies if I've misunderstood.
-Christian

## Compile/run from R using the following (assuming inline package is
installed):
source('gibbs.R.txt')
fun1(10)

-- 
A man, a plan, a cat, a ham, a yak, a yam, a hat, a canal – Panama!
-------------- next part --------------
#include "SamplingCode.h"


inc = '
//using namespace Rcpp;
int Cwhich(IntegerVector binVec){
  int vecsize = binVec.size();
  for ( int i =0; i < vecsize; i++)
    if (binVec(i) == 1)
      return i;
  std::cout << "soemthing went wrong, no value equal to one was found !" << std::endl;
  return -1;
}
'

#RcppExport SEXP doGibbs(SEXP RmaxItr){//doGibbs <- function(Questions, K, P, gibbs.max.iter, Z.Q=NULL, Y.Q = NULL){

src = '
	Environment stats("package:stats"); 
	Function rmultinom = stats["rmultinom"];
	Rcpp::RNGScope scope;
//	Rcpp::List questions(Rquestions);
	
	unsigned int K =50;
	unsigned int P = 50;
    NumericVector probs1(P, 1.0/P); //added -- instead of rep
    NumericVector probs2(K, 1.0/K); //added
	unsigned int maxItr = Rcpp::as<int>(RmaxItr);
	Rcpp::List ZQ;
	Rcpp::List YQ;
	unsigned int unum=300;

	std::cout<<unum<<std::endl;
	std::cout<<"point1"<<std::endl;
	for(unsigned itr=0; itr<maxItr; itr++){ //changed
			std::cout<<"itr"<<itr<<std::endl;
		for ( unsigned int u=0; u < unum; u++) {
			Rcpp::List QZQ;
			unsigned int qnum= 100;
            IntegerVector tmpu(qnum); //added
			//YQ.insert(u, rep(0, qnum ));	

			for( unsigned int q=0; q<qnum; q++){
				//as<Rcpp::IntegerVector>(YQ(u))(q) = Cwhich(rmultinom(1,1, probs1 ));
				tmpu(q) = Cwhich(rmultinom(1,1, probs1 ));
				unsigned int qlength=200;
				IntegerVector tmpq(qlength); //added
				//QZQ.insert(q, rep(0, qlength ));

				for( unsigned int n=0; n<qlength; n++){   
						tmpq(n)= Cwhich(rmultinom(1,1, probs2 ));	
						//as<Rcpp::IntegerVector>(QZQ(q))(n)= Cwhich(rmultinom(1,1, probs2 ));	
				}
				QZQ.insert(q, tmpq);
			}
			YQ.insert(u, tmpu);
			ZQ.insert(u, QZQ);
		}
    }
'

require(inline)
fun1 = cxxfunction(signature(RmaxItr='integer'), body=src, includes=inc, plugin='Rcpp')
#}