+// assign threads
+void init() {
+ int quotient, left;
+ seedType seed;
+ set<seedType> seedSet;
+ set<seedType>::iterator iter;
+ char outF[STRLEN];
+
+ quotient = NSAMPLES / nThreads;
+ left = NSAMPLES % nThreads;
+
+ sprintf(cvsF, "%s.countvectors", imdName);
+ seedSet.clear();
+ params = new Params[nThreads];
+ threads = new pthread_t[nThreads];
+ for (int i = 0; i < nThreads; i++) {
+ params[i].no = i;
+
+ params[i].nsamples = quotient;
+ if (i < left) params[i].nsamples++;
+
+ if (i == 0) {
+ params[i].fo = fopen(cvsF, "w");
+ fprintf(params[i].fo, "%d %d\n", NSAMPLES, M + 1);
+ }
+ else {
+ sprintf(outF, "%s%d", cvsF, i);
+ params[i].fo = fopen(outF, "w");
+ }
+
+ do {
+ seed = engine();
+ iter = seedSet.find(seed);
+ } while (iter != seedSet.end());
+ params[i].rng = new engine_type(seed);
+ seedSet.insert(seed);
+
+ params[i].pme_c = new double[M + 1];
+ memset(params[i].pme_c, 0, sizeof(double) * (M + 1));
+ params[i].pme_theta = new double[M + 1];
+ memset(params[i].pme_theta, 0, sizeof(double) * (M + 1));
+ }
+
+ /* set thread attribute to be joinable */
+ pthread_attr_init(&attr);
+ pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE);
+
+}
+
+void sampleTheta(engine_type& engine, vector<double>& theta) {
+ distribution_type gm(1);
+ generator_type gmg(engine, gm);
+ double denom;
+
+ theta.clear();
+ theta.resize(M + 1);
+ denom = 0.0;
+ for (int i = 0; i <= M; i++) {
+ theta[i] = gmg();
+ denom += theta[i];
+ }
+ assert(denom > EPSILON);
+ for (int i = 0; i <= M; i++) theta[i] /= denom;
+}
+