-
- delete[] pdf;
- delete[] cdf;
- delete[] clen;
-}
-
-template<class ModelType>
-void writeEstimatedParameters(char* modelF, char* imdName) {
- ModelType model;
- double denom;
- char outF[STRLEN];
- FILE *fo;
-
- model.read(modelF);
-
- calcExpectedEffectiveLengths<ModelType>(model);
-
- denom = pme_theta[0];
- for (int i = 1; i <= M; i++)
- if (eel[i] < EPSILON) pme_theta[i] = 0.0;
- else denom += pme_theta[i];
-
- general_assert(denom >= EPSILON, "No Expected Effective Length is no less than " + ftos(MINEEL, 6) + "?!");
-
- for (int i = 0; i <= M; i++) pme_theta[i] /= denom;
-
- denom = 0.0;
- double *mw = model.getMW();
- for (int i = 0; i <= M; i++) {
- pme_theta[i] = (mw[i] < EPSILON ? 0.0 : pme_theta[i] / mw[i]);
- denom += pme_theta[i];
- }
- assert(denom >= EPSILON);
- for (int i = 0; i <= M; i++) pme_theta[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] = pme_theta[i] / eel[i];
- denom += tau[i];
- }
-
- general_assert(denom >= EPSILON, "No alignable reads?!");
-
- for (int i = 1; i <= M; i++) {
- tau[i] /= denom;
- }
-
- //isoform level results
- sprintf(outF, "%s.iso_res", imdName);
- fo = fopen(outF, "a");
- general_assert(fo != NULL, "Cannot open " + cstrtos(outF) + "!");
-
- for (int i = 1; i <= M; i++)
- fprintf(fo, "%.2f%c", pme_c[i], (i < M ? '\t' : '\n'));
- for (int i = 1; i <= M; i++)
- fprintf(fo, "%.15g%c", tau[i], (i < M ? '\t' : '\n'));
-
- fclose(fo);
-
- //gene level results
- sprintf(outF, "%s.gene_res", imdName);
- fo = fopen(outF, "a");
- general_assert(fo != NULL, "Cannot open " + cstrtos(outF) + "!");
-
- for (int i = 0; i < m; i++) {
- double sumC = 0.0; // sum of pme counts
- int b = gi.spAt(i), e = gi.spAt(i + 1);
- for (int j = b; j < e; j++) {
- sumC += pme_c[j];
- }
- fprintf(fo, "%.15g%c", sumC, (i < m - 1 ? '\t' : '\n'));
- }
- for (int i = 0; i < m; i++) {
- double sumT = 0.0; // sum of tau values
- int b = gi.spAt(i), e = gi.spAt(i + 1);
- for (int j = b; j < e; j++) {
- sumT += tau[j];
- }
- fprintf(fo, "%.15g%c", sumT, (i < m - 1 ? '\t' : '\n'));
- }
- fclose(fo);
-
- delete[] tau;
-
- if (verbose) { printf("Gibbs based expression values are written!\n"); }