X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=simulation.cpp;h=e516cc15c32848e143aaf4b3a56a1e738fbe3456;hb=fc69cf6af24c0550e55447fc82f01cb6f90c1c42;hp=ecc0a11a20eb6b91f3bdb4b906cc187bf0cac7aa;hpb=a95154919f950f86de9104b2b9dcf1f0c7e83387;p=rsem.git

diff --git a/simulation.cpp b/simulation.cpp
index ecc0a11..e516cc1 100644
--- a/simulation.cpp
+++ b/simulation.cpp
@@ -27,7 +27,7 @@
 #include "Transcript.h"
 #include "Transcripts.h"
 
-#include "simul_mersenne.h"
+#include "simul.h"
 
 using namespace std;
 
@@ -48,7 +48,7 @@ char outReadF[2][STRLEN];
 char refF[STRLEN], groupF[STRLEN], tiF[STRLEN];
 char geneResF[STRLEN], isoResF[STRLEN];
 
-simul_mersenne sampler;
+simul sampler;
 
 void genOutReadStreams(int type, char *outFN) {
 	switch(type) {
@@ -127,7 +127,7 @@ void simulate(char* modelF, char* resultsF) {
 	for (int i = 1; i <= M; i++) {
 	  getline(fin, line);
 	  size_t pos = 0;
-	  for (int j = 0; j < 3; j++) pos = line.find_first_of('\t', pos) + 1;
+	  for (int j = 0; j < 2; j++) pos = line.find_first_of('\t', pos) + 1;
 	  size_t pos2 = line.find_first_of('\t', pos);
 	  if (pos2 == string::npos) pos2 = line.length();
 	  tau = atof(line.substr(pos, pos2 - pos).c_str());
@@ -138,49 +138,37 @@ void simulate(char* modelF, char* resultsF) {
 	fin.close();
 	for (int i = 1; i <= M; i++) theta[i] = theta[i] / denom * (1.0 - theta[0]);
 	
-	int tmp2140 = 0;
+	int resimulation_count = 0;
 
 	//simulating...
 	model.startSimulation(&sampler, theta);
 	for (int i = 0; i < N; i++) {
-	  while (!model.simulate(i, read, sid)) { ++tmp2140; }//;
+		while (!model.simulate(i, read, sid)) { ++resimulation_count; }
 		read.write(n_os, os);
 		++counts[sid];
 		if ((i + 1) % 1000000 == 0 && verbose) printf("GEN %d\n", i + 1);
 	}
 	model.finishSimulation();
 
-	printf("Total number of resimulation is %d\n", tmp2140);
+	printf("Total number of resimulation is %d\n", resimulation_count);
 }
 
 void writeResFiles(char* outFN) {
 	FILE *fo;
 	double denom;
 
-	//calculate normalized read fraction
-	double *nrf = new double[M + 1];
-	memset(nrf, 0, sizeof(double) * (M + 1));
-	denom = 0.0;
-	for (int i = 1; i <= M; i++)
-	  if (eel[i] > EPSILON) {
-	    nrf[i] = counts[i];
-	    denom += nrf[i];
-	  }
-	  else {
-	    if (counts[i] > EPSILON) { printf("Warning: An isoform which EEL is less than %.6g gets sampled!\n", MINEEL); }
-	  }
-	assert(denom > 0.0);
-	for (int i = 1; i <= M; i++) nrf[i] /= denom;
-
 	//calculate tau values
 	double *tau = new double[M + 1];
 	memset(tau, 0, sizeof(double) * (M + 1));
 	denom = 0.0;
 	for (int i = 1; i <= M; i++) 
-	  if (eel[i] > EPSILON) {
-	    tau[i] = nrf[i] / eel[i];
-	    denom += tau[i];
-	  }   
+		if (eel[i] > EPSILON) {
+			tau[i] = counts[i] / eel[i];
+			denom += tau[i];
+		}
+		else {
+		    if (counts[i] > EPSILON) { printf("Warning: An isoform which EEL is less than %.6g gets sampled!\n", MINEEL); }
+		}
 	assert(denom > 0.0);
 	for (int i = 1; i <= M; i++) tau[i] /= denom;
 
@@ -189,7 +177,7 @@ void writeResFiles(char* outFN) {
 	fo = fopen(isoResF, "w");
 	for (int i = 1; i <= M; i++) {
 		const Transcript& transcript = transcripts.getTranscriptAt(i);
-		fprintf(fo, "%s\t%.2f\t%.15g\t%.15g", transcript.getTranscriptID().c_str(), counts[i], nrf[i], tau[i]);
+		fprintf(fo, "%s\t%.2f\t%.15g", transcript.getTranscriptID().c_str(), counts[i], tau[i]);
 		
 		if (transcript.getLeft() != "") { fprintf(fo, "\t%s", transcript.getLeft().c_str()); }
 		fprintf(fo, "\n");
@@ -200,22 +188,20 @@ void writeResFiles(char* outFN) {
 	sprintf(geneResF, "%s.sim.genes.results", outFN);
 	fo = fopen(geneResF, "w");
 	for (int i = 0; i < m; i++) {
-	  double sum_c = 0.0, sum_nrf = 0.0, sum_tau = 0.0;
+	  double sum_c = 0.0, sum_tau = 0.0;
 		int b = gi.spAt(i), e = gi.spAt(i + 1);
 		for (int j = b; j < e; j++) {
 			sum_c += counts[j];
-			sum_nrf += nrf[j];
 			sum_tau += tau[j];
 		}
 		const string& gene_id = transcripts.getTranscriptAt(b).getGeneID();
-		fprintf(fo, "%s\t%.2f\t%.15g\t%.15g\t", gene_id.c_str(), sum_c, sum_nrf, sum_tau);
+		fprintf(fo, "%s\t%.2f\t%.15g\t", gene_id.c_str(), sum_c, sum_tau);
 		for (int j = b; j < e; j++) {
 			fprintf(fo, "%s%c", transcripts.getTranscriptAt(j).getTranscriptID().c_str(), (j < e - 1 ? ',' : '\n'));
 		}
 	}
 	fclose(fo);
 
-	delete[] nrf;
 	delete[] tau;
 }