+ sprintf(thetaF, "%s.theta", statName);
+ fo = fopen(thetaF, "w");
+ fprintf(fo, "%d\n", M + 1);
+
+ // output theta'
+ for (int i = 0; i < M; i++) fprintf(fo, "%.15g ", theta[i]);
+ fprintf(fo, "%.15g\n", theta[M]);
+
+ //calculate expected effective lengths for each isoform
+ calcExpectedEffectiveLengths<ModelType>(model);
+
+ //correct theta vector
+ sum = theta[0];
+ for (int i = 1; i <= M; i++)
+ if (eel[i] < EPSILON) { theta[i] = 0.0; }
+ else sum += theta[i];
+ if (sum < EPSILON) { fprintf(stderr, "No Expected Effective Length is no less than %.6g?!\n", MINEEL); exit(-1); }
+ for (int i = 0; i <= M; i++) theta[i] /= sum;
+