- case 0 : writeEstimatedParameters<SingleModel>(modelF, imdName); break;
- case 1 : writeEstimatedParameters<SingleQModel>(modelF, imdName); break;
- case 2 : writeEstimatedParameters<PairedEndModel>(modelF, imdName); break;
- case 3 : writeEstimatedParameters<PairedEndQModel>(modelF, imdName); break;
+ case 0 : init_model_related<SingleModel>(modelF); break;
+ case 1 : init_model_related<SingleQModel>(modelF); break;
+ case 2 : init_model_related<PairedEndModel>(modelF); break;
+ case 3 : init_model_related<PairedEndQModel>(modelF); break;
+ }
+
+ if (verbose) printf("Gibbs started!\n");
+
+ init();
+ for (int i = 0; i < nThreads; i++) {
+ rc = pthread_create(&threads[i], &attr, Gibbs, (void*)(¶msArray[i]));
+ pthread_assert(rc, "pthread_create", "Cannot create thread " + itos(i) + " (numbered from 0)!");
+ }
+ for (int i = 0; i < nThreads; i++) {
+ rc = pthread_join(threads[i], NULL);
+ pthread_assert(rc, "pthread_join", "Cannot join thread " + itos(i) + " (numbered from 0)!");
+ }
+ release();
+
+ if (verbose) printf("Gibbs finished!\n");
+
+ writeResultsGibbs(M, refName, imdName, pme_c, pme_fpkm, pme_tpm);
+
+ if (var_opt) {
+ char varF[STRLEN];
+
+ // Load group info
+ int m;
+ GroupInfo gi;
+ char groupF[STRLEN];
+ sprintf(groupF, "%s.grp", refName);
+ gi.load(groupF);
+ m = gi.getm();
+
+ sprintf(varF, "%s.var", statName);
+ FILE *fo = fopen(varF, "w");
+ general_assert(fo != NULL, "Cannot open " + cstrtos(varF) + "!");
+ for (int i = 0; i < m; i++) {
+ int b = gi.spAt(i), e = gi.spAt(i + 1), number_of_isoforms = e - b;
+ for (int j = b; j < e; j++) {
+ fprintf(fo, "%s\t%d\t%.15g\t%.15g\n", refs.getRef(j).getName().c_str(), number_of_isoforms, pme_c[j], pve_c[j]);
+ }
+ }
+ fclose(fo);