- 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;
-}
-
-// arr should be cumulative!
-// interval : [,)
-// random number should be in [0, arr[len - 1])
-// If by chance arr[len - 1] == 0.0, one possibility is to sample uniformly from 0...len-1
-int sample(uniform_generator& rg, vector<double>& arr, int len) {
- int l, r, mid;
- double prb = rg() * arr[len - 1];
-
- l = 0; r = len - 1;
- while (l <= r) {
- mid = (l + r) / 2;
- if (arr[mid] <= prb) l = mid + 1;
- else r = mid - 1;
- }
-
- if (l >= len) { printf("%d %lf %lf\n", len, arr[len - 1], prb); }
- assert(l < len);
-
- return l;
-}
-
-void writeForSee(FILE* fo, vector<int>& counts) {
- for (int i = 1; i < M; i++) fprintf(fo, "%d ", counts[i]);
- fprintf(fo, "%d\n", counts[M]);
-}
-
-void writeCountVector(FILE* fo, vector<int>& counts) {
- for (int i = 0; i < M; i++) {
- fprintf(fo, "%d ", counts[i]);
- }
- fprintf(fo, "%d\n", counts[M]);
-}
-
-void* Gibbs(void* arg) {