Skip to content

Commit 3c46521

Browse files
prataspratas
committed
fixed bug having N
1 parent 0fe24a8 commit 3c46521

9 files changed

Lines changed: 96 additions & 76 deletions

File tree

‎src/buffer.c‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -49,6 +49,7 @@ void UpdateBuffer(BUF *B){
4949
//
5050
void ResetCBuffer(CBUF *B){
5151
Free(B->buf-B->guard);
52+
// B->buf = (uint8_t *) Calloc(B->size, sizeof(uint8_t));
5253
B->buf = (uint8_t *) Calloc(B->size+B->guard, sizeof(uint8_t));
5354
B->buf += B->guard;
5455

‎src/common.c‎

Lines changed: 6 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -8,6 +8,7 @@
88
#include "mem.h"
99
#include "common.h"
1010
#include "msg.h"
11+
#include "param.h"
1112

1213
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1314

@@ -757,20 +758,20 @@ void PrintArgs(Parameters *P, Threads T, char *ref, char *tar){
757758
uint32_t n;
758759

759760
fprintf(stderr, "==[ CONFIGURATION ]=================\n");
760-
fprintf(stderr, "Verbose mode ....................... %s\n", P->verbose == 0
761+
fprintf(stderr, "Verbose mode ....................... %s\n", P->verbose == 0
761762
? "no" : "yes");
762-
fprintf(stderr, "Force mode ......................... %s\n", P->force == 0 ?
763+
fprintf(stderr, "Force mode ......................... %s\n", P->force == 0 ?
763764
"no" : "yes");
764765
fprintf(stderr, "Compression level .................. %u\n", P->level);
765766
fprintf(stderr, "Sub-sampling ....................... %u\n", P->sample);
766767
fprintf(stderr, "Number of threads .................. %u\n", P->nThreads);
767768
for(n = 0 ; n < P->nModels ; ++n){
768769
fprintf(stderr, "Reference model %u:\n", n+1);
769-
fprintf(stderr, " [+] Context order ................ %u\n",
770+
fprintf(stderr, " [+] Context order ................ %u\n",
770771
T.model[n].ctx);
771-
fprintf(stderr, " [+] Alpha denominator ............ %u\n",
772+
fprintf(stderr, " [+] Alpha denominator ............ %u\n",
772773
T.model[n].den);
773-
fprintf(stderr, " [+] Inverted repeats ............. %s\n",
774+
fprintf(stderr, " [+] Inverted repeats ............. %s\n",
774775
T.model[n].ir == 0 ? "no" : "yes");
775776
fprintf(stderr, " [+] Allowable substitutions ...... %u\n",
776777
T.model[n].edits);

‎src/defs.h‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -9,6 +9,7 @@
99
#else
1010
#include <unistd.h>
1111
#endif
12+
#include <assert.h>
1213

1314
typedef uint64_t ULL;
1415
typedef uint64_t U64;
@@ -27,7 +28,7 @@ typedef int8_t I8;
2728

2829
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2930

30-
#define MAX_NAME_OUT 10000
31+
#define MAX_NAME_OUT 100000
3132
#define SBASE 65
3233
#define LOCAL_SIMILARITY 1
3334
#define MAX_LABEL 1024

‎src/magnet.c‎

Lines changed: 46 additions & 36 deletions
Original file line numberDiff line numberDiff line change
@@ -66,6 +66,7 @@ void CompressTarget(Threads T){
6666
Shadow = (CModel **) Calloc(P->nModels, sizeof(CModel *));
6767
for(n = 0 ; n < P->nModels ; ++n)
6868
Shadow[n] = CreateShadowModel(Models[n]);
69+
6970
pModel = (PModel **) Calloc(totModels, sizeof(PModel *));
7071
for(n = 0 ; n < totModels ; ++n)
7172
pModel[n] = CreatePModel(ALPHABET_SIZE);
@@ -79,28 +80,28 @@ void CompressTarget(Threads T){
7980
exit(1);
8081
}
8182

82-
char name_o[4096];
83-
sprintf(name_o, "%s.%u", P->output, T.id);
84-
FILE *Writer = Fopen(name_o, "w");
85-
83+
char name_o [10000];
84+
sprintf(name_o, "%s.%u", P->output, T.id);
85+
FILE *Writer = Fopen(name_o, "w");
8686
srand(T.id);
87-
Read *Read = CreateRead(10000, 40000);
88-
while((Read = GetRead(Reader, Read)) != NULL){
87+
Read *Read = CreateRead(50000, 100000);
8988

90-
if(PA->nRead % P->nThreads == T.id){
89+
while((Read = GetRead(Reader, Read)) != NULL)
90+
{
91+
if(PA->nRead % P->nThreads == T.id)
92+
{
9193
nBase = strlen(Read->bases) - 1; // IT ALSO LOADS '\n' AT THE END
9294
bits = 0;
9395

94-
for(idxPos = 0 ; idxPos < nBase ; ++idxPos){
95-
96+
for(idxPos = 0 ; idxPos < nBase ; ++idxPos)
97+
{
9698
sym = Read->bases[idxPos];
97-
98-
if(sym == 'N') sym = 0;// rand() % 4; // ASSUME A 'A' -> 0, CAUSE:
99-
// RANDOM BASE MODIFY THE RESULTS USING DIFFERENT THREADS
100-
else sym = DNASymToNum(sym);
99+
100+
if((sym = DNASymToNum(sym)) == 4)
101+
sym = rand() % 4;
101102

102103
symBuf->buf[symBuf->idx] = sym;
103-
memset((void *)PT->freqs, 0, ALPHABET_SIZE * sizeof(double));
104+
memset((void *)PT->freqs, 0, (ALPHABET_SIZE)* sizeof(double));
104105
n = 0;
105106
pos = &symBuf->buf[symBuf->idx-1];
106107
for(cModel = 0 ; cModel < P->nModels ; ++cModel){
@@ -126,22 +127,14 @@ void CompressTarget(Threads T){
126127
CorrectXModels(Shadow, pModel, sym, P->nModels);
127128
UpdateCBuffer(symBuf);
128129
}
129-
130+
130131
if(BPBB(bits, nBase) < P->threshold){
131-
if(P->invert){
132-
fputc('0', Writer); // IGNORE READ
133-
}
134-
else{
135-
fputc('1', Writer); // WRITE READ
136-
}
132+
if(P->invert) fputc('0', Writer); // IGNORE READ
133+
else fputc('1', Writer); // WRITE READ
137134
}
138135
else{
139-
if(P->invert){
140-
fputc('1', Writer); // WRITE READ
141-
}
142-
else{
143-
fputc('0', Writer); // IGNORE READ
144-
}
136+
if(P->invert) fputc('1', Writer); // WRITE READ
137+
else fputc('0', Writer); // IGNORE READ
145138
}
146139
ResetModelsAndParam(symBuf, Shadow, CMW);
147140
}
@@ -197,10 +190,12 @@ void LoadReference(char *refName){
197190
idx = 0;
198191
continue;
199192
}
193+
200194
if(sym == 'N') // WE CAN RAND HERE CAUSE IS ALWAYS IN ONE THREAD
201-
symBuf->buf[symBuf->idx] = sym = (rand() % 4);
195+
symBuf->buf[symBuf->idx] = sym = 0; //(rand() % 4);
202196
else
203197
symBuf->buf[symBuf->idx] = sym = DNASymToNum(sym);
198+
204199
for(n = 0 ; n < P->nModels ; ++n){
205200
CModel *CM = Models[n];
206201
GetPModelIdx(symBuf->buf+symBuf->idx-1, CM);
@@ -248,20 +243,25 @@ void CompressAction(Threads *T, char *refName, char *baseName){
248243
fprintf(stderr, "Done!\n");
249244

250245
fprintf(stderr, " [+] Joinning streams ............. ");
251-
FILE *OUT = Fopen(P->output, "w");
252-
FILE *IN = Fopen(P->base, "r");
253-
FILE **TMP = (FILE **) Calloc(P->nThreads, sizeof(FILE *));
246+
FILE *OUT = Fopen(P->output, "w");
247+
FILE *OUT2 = Fopen(P->output2, "w");
248+
FILE *IN = Fopen(P->base, "r");
249+
FILE **TMP = (FILE **) Calloc(P->nThreads, sizeof(FILE *));
254250
for(n = 0 ; n < P->nThreads ; ++n){
255251
char name_o[MAX_NAME_OUT];
256252
sprintf(name_o, "%s.%u", P->output, n);
257253
TMP[n] = Fopen(name_o, "r");
258254
}
255+
259256
Read *Read = CreateRead(10000, 40000);
260257
n = 0;
261258
while((Read = GetRead(IN, Read)) != NULL){
262259
if(fgetc(TMP[n++ % P->nThreads]) == '1')
263260
PutRead(Read, OUT);
261+
else
262+
PutRead(Read, OUT2);
264263
}
264+
265265
for(n = 0 ; n < P->nThreads ; ++n){
266266
fclose(TMP[n]);
267267
char name_o[MAX_NAME_OUT];
@@ -270,6 +270,7 @@ void CompressAction(Threads *T, char *refName, char *baseName){
270270
}
271271
fclose(IN);
272272
fclose(OUT);
273+
fclose(OUT2);
273274
fprintf(stderr, "Done!\n");
274275
}
275276

@@ -282,7 +283,7 @@ void CompressAction(Threads *T, char *refName, char *baseName){
282283

283284
int32_t main(int argc, char *argv[]){
284285

285-
char **p = *&argv, **xargv, *xpl = NULL;
286+
char **p = *&argv, **xargv, *xpl = NULL;
286287
int32_t xargc = 0;
287288
uint32_t n, k, col, ref;
288289
double gamma;
@@ -344,12 +345,16 @@ int32_t main(int argc, char *argv[]){
344345
P->threshold = fabs(ArgsDouble (0.9, p, argc, "-t", "--threshold"));
345346
P->gamma = ((int)(P->gamma * 65536)) / 65536.0;
346347
P->output = ArgsFileGen(p, argc, "-o", "filtered", ".fq");
348+
P->output2 = ArgsFileGen(p, argc, "-2", "not-filtered", ".fq");
347349

348350
FILE *OUTPUT = NULL;
349351
if(!P->force)
350352
FAccessWPerm(P->output);
351353
OUTPUT = Fopen(P->output, "w");
352354

355+
FILE *OUTPUT2 = NULL;
356+
OUTPUT2 = Fopen(P->output2, "w");
357+
353358
if(P->nModels == 0){
354359
fprintf(stderr, "Error: at least you need to use a context model!\n");
355360
return EXIT_FAILURE;
@@ -389,11 +394,16 @@ int32_t main(int argc, char *argv[]){
389394
fprintf(stderr, "\n");
390395

391396
fprintf(stderr, "==[ RESULTS ]=======================\n");
392-
fprintf(stderr, "Results have been sent to file:%s\n", P->output);
397+
fprintf(stderr, "Filtered reads in : %s\n", P->output);
398+
fprintf(stderr, "Unfiltered reads in : %s\n", P->output2);
399+
fprintf(stderr, "\n");
393400

394-
//fprintf(stderr, "==[ STATISTICS ]====================\n");
395-
//StopCalcAll(Time, clock());
396-
//fprintf(stderr, "\n");
401+
if(P->verbose)
402+
{
403+
fprintf(stderr, "==[ STATISTICS ]====================\n");
404+
StopCalcAll(Time, clock());
405+
fprintf(stderr, "\n");
406+
}
397407

398408
RemoveClock(Time);
399409
for(ref = 0 ; ref < P->nThreads ; ++ref)

‎src/models.c‎

Lines changed: 23 additions & 23 deletions
Original file line numberDiff line numberDiff line change
@@ -83,41 +83,41 @@ static void InsertKey(HashTable *H, U32 hi, U64 idx, U8 s){
8383

8484
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
8585

86-
void GetFreqsFromHCC(HCC c, uint32_t a, PModel *P){
87-
P->sum = (P->freqs[0] = 1 + a * ( c & 0x0f)); // NO +1: N ARITHCODER
88-
P->sum += (P->freqs[1] = 1 + a * ((c & (0x0f<<4 ))>>4)); // NO +1: N ARITHCODER
89-
P->sum += (P->freqs[2] = 1 + a * ((c & (0x0f<<8 ))>>8)); // NO +1: N ARITHCODER
90-
P->sum += (P->freqs[3] = 1 + a * ((c & (0x0f<<12))>>12)); // NO +1: N ARITHCODER
86+
void GetFreqsFromHCC(HCC c, uint32_t a, PModel *PM){
87+
PM->sum = (PM->freqs[0] = 1 + a * ( c & 0x0f));
88+
PM->sum += (PM->freqs[1] = 1 + a * ((c & (0x0f<<4 ))>>4));
89+
PM->sum += (PM->freqs[2] = 1 + a * ((c & (0x0f<<8 ))>>8));
90+
PM->sum += (PM->freqs[3] = 1 + a * ((c & (0x0f<<12))>>12));
9191
}
9292

9393
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9494

95-
void GetHCCounters(HashTable *H, U64 key, PModel *P, uint32_t a){
95+
void GetHCCounters(HashTable *H, U64 key, PModel *PM, uint32_t a){
9696
U32 n, hIndex = key % HASH_SIZE;
9797
U16 b = key & 0xffff;
9898

9999
U32 pos = H->index[hIndex];
100100
// FROM INDEX-1 TO 0
101101
for(n = pos+1 ; n-- ; ){
102102
if(H->entries[hIndex][n].key == b){
103-
GetFreqsFromHCC(H->entries[hIndex][n].counters, a, P);
103+
GetFreqsFromHCC(H->entries[hIndex][n].counters, a, PM);
104104
return;
105105
}
106106
}
107107
// FROM MAX_COLISIONS TO INDEX
108108
for(n = (H->maxC-1) ; n > pos ; --n){
109109
if(H->entries[hIndex][n].key == b){
110-
GetFreqsFromHCC(H->entries[hIndex][n].counters, a, P);
110+
GetFreqsFromHCC(H->entries[hIndex][n].counters, a, PM);
111111
return;
112112
}
113113
}
114114

115115
// TODO: MAKE THIS ALREADY DONE!
116-
P->freqs[0] = 1;
117-
P->freqs[1] = 1;
118-
P->freqs[2] = 1;
119-
P->freqs[3] = 1;
120-
P->sum = 4;
116+
PM->freqs[0] = 1;
117+
PM->freqs[1] = 1;
118+
PM->freqs[2] = 1;
119+
PM->freqs[3] = 1;
120+
PM->sum = 4;
121121
return;
122122
}
123123

@@ -156,7 +156,7 @@ void UpdateCModelCounter(CModel *M, U32 sym, U64 im){
156156
}
157157
else{
158158
AC = &M->array.counters[idx << 2];
159-
if(++AC[sym] == M->maxCount){
159+
if(++AC[sym] >= M->maxCount){
160160
AC[0] >>= 1;
161161
AC[1] >>= 1;
162162
AC[2] >>= 1;
@@ -328,8 +328,8 @@ void HitSUBS(CModel *M){
328328

329329
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
330330

331-
void CorrectCModelSUBS(CModel *M, PModel *P, uint8_t sym){
332-
int32_t best = BestId(P->freqs, P->sum);
331+
void CorrectCModelSUBS(CModel *M, PModel *PM, uint8_t sym){
332+
int32_t best = BestId(PM->freqs, PM->sum);
333333
switch(best){
334334
case -2: // IT IS A ZERO COUNTER [NOT SEEN BEFORE]
335335
if(M->SUBS.in != 0)
@@ -369,19 +369,19 @@ void CorrectXModels(CModel **Shadow, PModel **PM, uint8_t sym, uint32_t nModels)
369369

370370
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
371371

372-
void ComputePModel(CModel *M, PModel *P, uint64_t idx, uint32_t aDen){
372+
void ComputePModel(CModel *M, PModel *PM, uint64_t idx, uint32_t aDen){
373373
ACC *ac;
374374
switch(M->mode){
375375
case HASH_TABLE_MODE:
376-
GetHCCounters(&M->hTable, ZHASH(idx), P, aDen);
376+
GetHCCounters(&M->hTable, ZHASH(idx), PM, aDen);
377377
break;
378378
case ARRAY_MODE:
379379
ac = &M->array.counters[idx<<2];
380-
P->freqs[0] = 1 + aDen * ac[0]; // +1 IS NOT NEEDED BECAUSE THERE IS NO AC
381-
P->freqs[1] = 1 + aDen * ac[1]; // +1 IS NOT NEEDED BECAUSE THERE IS NO AC
382-
P->freqs[2] = 1 + aDen * ac[2]; // +1 IS NOT NEEDED BECAUSE THERE IS NO AC
383-
P->freqs[3] = 1 + aDen * ac[3]; // +1 IS NOT NEEDED BECAUSE THERE IS NO AC
384-
P->sum = P->freqs[0] + P->freqs[1] + P->freqs[2] + P->freqs[3];
380+
PM->freqs[0] = 1 + aDen * ac[0];
381+
PM->freqs[1] = 1 + aDen * ac[1];
382+
PM->freqs[2] = 1 + aDen * ac[2];
383+
PM->freqs[3] = 1 + aDen * ac[3];
384+
PM->sum = PM->freqs[0] + PM->freqs[1] + PM->freqs[2] + PM->freqs[3];
385385
break;
386386
default:
387387
fprintf(stderr, "Error: not implemented!\n");

‎src/msg.c‎

Lines changed: 6 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -71,7 +71,12 @@ void PrintMenu(void){
7171
" \n"
7272
" -o [OUT_FILE] \n"
7373
" Output file name. \n"
74-
" The OUT_FILE will contain the filtered reads. \n"
74+
" The OUT_FILE will contain the filtered FASTQ reads. \n"
75+
" \n"
76+
" -2 [OUT_FILE] \n"
77+
" Output file name. \n"
78+
" The OUT_FILE will contain the unfiltered FASTQ reads. \n"
79+
" \n"
7580
" \n"
7681
" -n [NUMBER], --threads [NUMBER] \n"
7782
" Number of threads. \n"

‎src/param.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -33,6 +33,7 @@ typedef struct{
3333
U8 nFiles;
3434
U64 *size;
3535
char *output;
36+
char *output2;
3637
char *base;
3738
}
3839
Parameters;

‎src/pmodels.c‎

Lines changed: 9 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -43,6 +43,7 @@ void RenormalizeWeights(CMWeight *CMW){
4343

4444
void CalcDecayment(CMWeight *CMW, PModel **PM, uint8_t sym, double gamma){
4545
uint32_t n;
46+
assert(sym < 4);
4647
CMW->totalWeight = 0;
4748
for(n = 0 ; n < CMW->totModels ; ++n){
4849
CMW->weight[n] = Power(CMW->weight[n], gamma) * (double) PM[n]->freqs[sym]
@@ -99,18 +100,18 @@ void RemoveFPModel(FloatPModel *FM){
99100

100101
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
101102

102-
void ComputeWeightedFreqs(double w, PModel *P, FloatPModel *PT){
103-
double f = w / P->sum;
104-
PT->freqs[0] += (double) P->freqs[0] * f;
105-
PT->freqs[1] += (double) P->freqs[1] * f;
106-
PT->freqs[2] += (double) P->freqs[2] * f;
107-
PT->freqs[3] += (double) P->freqs[3] * f;
103+
void ComputeWeightedFreqs(double w, PModel *PM, FloatPModel *PT){
104+
double f = w / PM->sum;
105+
PT->freqs[0] += (double) PM->freqs[0] * f;
106+
PT->freqs[1] += (double) PM->freqs[1] * f;
107+
PT->freqs[2] += (double) PM->freqs[2] * f;
108+
PT->freqs[3] += (double) PM->freqs[3] * f;
108109
}
109110

110111
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
111112

112-
double PModelSymbolLog(PModel *P, U32 s){
113-
return log((double) P->sum / P->freqs[s]) / M_LN2; // TODO: FLOG2 ?
113+
double PModelSymbolLog(PModel *PM, U32 s){
114+
return log((double) PM->sum / PM->freqs[s]) / M_LN2; // TODO: FLOG2 ?
114115
}
115116

116117
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

0 commit comments

Comments
 (0)