struct Params {
int no;
FILE *fi;
- bufsize_type size;
- int sp;
engine_type *engine;
double *mw;
};
Params *params = (Params*)arg;
FILE *fi = params->fi;
double *mw = params->mw;
- Buffer buffer(params->size, params->sp, nSamples, cvlen, tmpF);
cvec = new int[cvlen];
theta = new double[cvlen];
gammas = new gamma_dist*[cvlen];
rgs = new gamma_generator*[cvlen];
- float *vec = new float[cvlen];
+ float **vecs = new float*[nSpC];
+ for (int i = 0; i < nSpC; i++) vecs[i] = new float[cvlen];
int cnt = 0;
while (fscanf(fi, "%d", &cvec[0]) == 1) {
sum = 0.0;
- vec[0] = theta[0];
+ vecs[i][0] = theta[0];
for (int j = 1; j < cvlen; j++)
if (eel[j] >= EPSILON) {
- vec[j] = theta[j] / eel[j];
- sum += vec[j];
+ vecs[i][j] = theta[j] / eel[j];
+ sum += vecs[i][j];
}
else assert(theta[j] < EPSILON);
- assert(sum >= EPSILON);
- for (int j = 1; j < cvlen; j++) vec[j] /= sum;
- buffer.write(vec);
+ assert(sum >= EPSILON);
+ for (int j = 1; j < cvlen; j++) vecs[i][j] /= sum;
}
+ buffer->write(nSpC, vecs);
+
for (int j = 0; j < cvlen; j++) {
delete gammas[j];
delete rgs[j];
delete[] gammas;
delete[] rgs;
- delete[] vec;
+ for (int i = 0; i < nSpC; i++) delete[] vecs[i];
+ delete[] vecs;
return NULL;
}
model.read(modelF);
calcExpectedEffectiveLengths<ModelType>(model);
- char splitF[STRLEN];
- bufsize_type buf_maxcv = bufsize_type(nMB) * 1024 * 1024 / FLOATSIZE / cvlen;
- bufsize_type quotient, left;
- int sum = 0;
- sprintf(splitF, "%s.split", imdName);
- FILE *fi = fopen(splitF, "r");
- int num_threads;
- assert(fscanf(fi, "%d", &num_threads) == 1);
- assert(num_threads <= nThreads);
+ int num_threads = min(nThreads, nCV);
- quotient = buf_maxcv / num_threads;
- assert(quotient > 0);
- left = buf_maxcv % num_threads;
+ buffer = new Buffer(nMB, nSamples, cvlen, tmpF);
paramsArray = new Params[num_threads];
threads = new pthread_t[num_threads];
paramsArray[i].no = i;
sprintf(inpF, "%s%d", cvsF, i);
paramsArray[i].fi = fopen(inpF, "r");
-
- int num_samples;
- assert(fscanf(fi, "%d", &num_samples) == 1);
- num_samples *= nSpC;
- if (bufsize_type(num_samples) <= quotient) paramsArray[i].size = num_samples;
- else {
- paramsArray[i].size = quotient;
- if (left > 0) { ++paramsArray[i].size; --left; }
- }
- paramsArray[i].size *= cvlen * FLOATSIZE;
- paramsArray[i].sp = sum;
- sum += num_samples;
-
paramsArray[i].engine = engineFactory::new_engine();
paramsArray[i].mw = model.getMW();
}
- fclose(fi);
-
/* set thread attribute to be joinable */
pthread_attr_init(&attr);
pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE);
}
delete[] paramsArray;
+ delete buffer; // Must delete here, force the content left in the buffer be written into the disk
+
if (verbose) { printf("Sampling is finished!\n"); }
}
void calculate_credibility_intervals(char* imdName) {
FILE *fo;
char outF[STRLEN];
+ int num_threads = nThreads;
iso_tau = new CIType[M + 1];
gene_tau = new CIType[m];
- // nThreads must be intact here.
assert(M > 0);
- int quotient = M / nThreads;
- if (quotient < 1) { nThreads = M; quotient = 1; }
+ int quotient = M / num_threads;
+ if (quotient < 1) { num_threads = M; quotient = 1; }
int cur_gene_id = 0;
int num_isoforms = 0;
// A just so so strategy for paralleling
- ciParamsArray = new CIParams[nThreads];
- for (int i = 0; i < nThreads; i++) {
+ ciParamsArray = new CIParams[num_threads];
+ for (int i = 0; i < num_threads; i++) {
ciParamsArray[i].no = i;
ciParamsArray[i].start_gene_id = cur_gene_id;
num_isoforms = 0;
- while ((m - cur_gene_id > nThreads - i - 1) && (i == nThreads - 1 || num_isoforms < quotient)) {
+ while ((m - cur_gene_id > num_threads - i - 1) && (i == num_threads - 1 || num_isoforms < quotient)) {
num_isoforms += gi.spAt(cur_gene_id + 1) - gi.spAt(cur_gene_id);
++cur_gene_id;
}
ciParamsArray[i].end_gene_id = cur_gene_id;
}
- threads = new pthread_t[nThreads];
+ threads = new pthread_t[num_threads];
/* set thread attribute to be joinable */
pthread_attr_init(&attr);
pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE);
// paralleling
- for (int i = 0; i < nThreads; i++) {
+ for (int i = 0; i < num_threads; i++) {
rc = pthread_create(&threads[i], &attr, &calcCI_batch, (void*)(&ciParamsArray[i]));
pthread_assert(rc, "pthread_create", "Cannot create thread " + itos(i) + " (numbered from 0) in calculate_credibility_intervals!");
}
- for (int i = 0; i < nThreads; i++) {
+ for (int i = 0; i < num_threads; i++) {
rc = pthread_join(threads[i], &status);
pthread_assert(rc, "pthread_join", "Cannot join thread " + itos(i) + " (numbered from 0) in calculate_credibility_intervals!");
}
int main(int argc, char* argv[]) {
if (argc < 8) {
- printf("Usage: rsem-calculate-credibility-intervals reference_name sample_name sampleToken confidence nCV nSpC nMB [-p #Threads] [-q]\n");
+ printf("Usage: rsem-calculate-credibility-intervals reference_name imdName statName confidence nCV nSpC nMB [-p #Threads] [-q]\n");
exit(-1);
}
+ strcpy(imdName, argv[2]);
+ strcpy(statName, argv[3]);
+
confidence = atof(argv[4]);
nCV = atoi(argv[5]);
nSpC = atoi(argv[6]);
cvlen = M + 1;
assert(nSamples > 0 && cvlen > 1); // for Buffter.h: (bufsize_type)nSamples
- sprintf(imdName, "%s.temp/%s", argv[2], argv[3]);
- sprintf(statName, "%s.stat/%s", argv[2], argv[3]);
sprintf(tmpF, "%s.tmp", imdName);
sprintf(cvsF, "%s.countvectors", imdName);
projects / rsem.git / blobdiff