Code cleanup

Makefile

@@ -1,9 +1,9 @@
 # --- Flags og kompileringsvalg ---
-FLAGS     := -O2
+FLAGS     := -O3
 CFLAGS    := $(FLAGS)
 CXXFLAGS  := $(FLAGS)
 CPPFLAGS  := $(CPPFLAGS) -Wall -Wextra
-LDFLAGS   := -lgsl
+LDFLAGS   := -lgsl -lgslcblas
 
 CC  ?= gcc
 CXX ?= g++

README.md

@@ -49,7 +49,7 @@ To install metaDMG-cpp do:
 ```
 git clone https://github.com/metaDMG-dev/metaDMG-cpp.git
 cd metaDMG-cpp
-make HTTSRC=../htslib
+make HTSSRC=../htslib
 ```
 
 ### Updating to latest version

dfit_optim.cpp

@@ -60,7 +60,7 @@ double **read1_ugly_matrix(const char *fname){
   ret[0][1] = 0;//<- function counter
   ret[1] = new double[xcol.size()];
   ret[2] = new double[xcol.size()];
-  for(int i=0;i<xcol.size();i++){
+  for(size_t i=0;i<xcol.size();i++){
     ret[0][i+2] = xcol[i];
     ret[1][i] = kcol[i];
     ret[2][i] = ncol[i];
@@ -161,10 +161,10 @@ void getstat(double **dat,double *pars,double *statpars){
     double q = pars[1];
     double c = pars[2];
     double phi = pars[3];
-    double llh = pars[4];
+    //double llh = pars[4];
     int NUMROWS = dat[0][0];
     double *XCOL = dat[0]+2;
-    double *KCOL = dat[1];
+    //double *KCOL = dat[1];
     double *NCOL = dat[2];
 
     double N_pos = dat[2][0];

main_dfit.cpp

@@ -339,7 +339,7 @@ void slave_block(std::map<int, mydataD> &retmap,int howmany,sam_hdr_t *hdr,int2c
   dat[3] = new double [2*howmany];
 
   for (std::map<int, mydataD>::iterator it = retmap.begin(); it != retmap.end(); it++) {
-    int taxid = it->first;
+    //int taxid = it->first;
     mydataD md = it->second;
     if (it->second.nal == 0)
       continue;
@@ -378,7 +378,6 @@ void slave_block(std::map<int, mydataD> &retmap,int howmany,sam_hdr_t *hdr,int2c
       double** bootcidata = (double**)malloc(nbootstrap * sizeof(double*));
 
       double acumtmp = 0; double qcumtmp = 0; double ccumtmp = 0; double phicumtmp = 0;
-      double astdtmp = 0; double qstdtmp = 0; double cstdtmp = 0; double phistdtmp = 0;
 
       //store the confidence interval  values
       for (int j = 0; j < npars; j++){
@@ -670,7 +669,7 @@ int main_dfit(int argc, char **argv) {
     int nbootstrap = 1;
     int doboot = 0;
     int seed = time(NULL);
-    int nthreads = 1;
+    size_t nthreads = 1;
     double CI = 0.95;
     int doCI = 2;
     int rng_type = -1;
@@ -853,7 +852,7 @@ int main_dfit(int argc, char **argv) {
         }
 
 	      ding *dings = new ding[nthreads];
-        for(int i=0;i<nthreads;i++){
+        for(size_t i=0;i<nthreads;i++){
           dings[i].retmap = &(ary[i]);
           dings[i].howmany = howmany;
           dings[i].hdr = hdr;
@@ -867,8 +866,8 @@ int main_dfit(int argc, char **argv) {
           dings[i].seed = (seed+i)*100;
           dings[i].rng_type = rng_type;
           dings[i].doboot = doboot;
-          fprintf(stderr, "\t-> Initiating thread %d with thread specific seed %d inferred from global seed %d with pseudo-random number generator type  %d\n",i,dings[i].seed/100,seed,rng_type);
-          //fprintf(stderr,"INITIATED THREAD WHAT %d WITH SEED VALUE WHAT %d WITH SPECIFIC SEED %d and seedtype %d \n",i,seed,dings[i].seed,rng_type);
+          fprintf(stderr, "\t-> Initiating thread %lu with thread specific seed %d inferred from global seed %d with pseudo-random number generator type  %d\n",i,dings[i].seed/100,seed,rng_type);
+          //fprintf(stderr,"INITIATED THREAD WHAT %lu WITH SEED VALUE WHAT %d WITH SPECIFIC SEED %d and seedtype %d \n",i,seed,dings[i].seed,rng_type);
 
           kstring_t *kstr = new kstring_t;
           kstr->s = NULL; kstr->l = kstr->m = 0;
@@ -879,13 +878,13 @@ int main_dfit(int argc, char **argv) {
           dings[i].showfits = showfits;
         }
         pthread_t mythreads[nthreads];
-        for(int i=0;i<nthreads;i++)
+        for(size_t i=0;i<nthreads;i++)
           pthread_create(&mythreads[i],NULL,slaveslave, &dings[i]);
 
-        for(int i=0;i<nthreads;i++)
+        for(size_t i=0;i<nthreads;i++)
           pthread_join(mythreads[i],NULL);
 
-        for(int i=0;i<nthreads;i++){
+        for(size_t i=0;i<nthreads;i++){
           ksprintf(kstr,"%s",dings[i].kstr->s);
           ksprintf(bootkstr,"%s",dings[i].bootkstr->s);
         }

main_pmd.cpp

@@ -244,8 +244,9 @@ void parse_sequencingdata(char *refName, char *fname, int mapped_only, int se_on
 }
 
 int usage(FILE *fp, int val) {
-    fprintf(stderr, "./metaDMG-cpp pmd [-T ref.fa -@ threads -a se_only -q minmapQ -v VERBOSE] file.bam\n");
-    fprintf(stderr, "-a is an integer zero or one, indicating if paired end reads should be discarded\n");
+    fprintf(fp, "ERROR (%d)\nUsage:", val);
+    fprintf(fp, "\t./metaDMG-cpp pmd [-T ref.fa -@ threads -a se_only -q minmapQ -v VERBOSE] file.bam\n");
+    fprintf(fp, "\t-a is an integer zero or one, indicating if paired end reads should be discarded\n");
     return 0;
 }
 
@@ -264,7 +265,6 @@ int main_pmd(int argc, char **argv) {
     }
 
     int VERBOSE = 0;
-    unsigned long int seed = 0;
 
     clock_t t = clock();
     time_t t2 = time(NULL);
@@ -343,7 +343,8 @@ int main_pmd(int argc, char **argv) {
         fname = strdup(argv[optind]);
 
     if (!fname) {
-        fprintf(stderr, "\t-> No input file specified\n");
+        if(VERBOSE)
+          fprintf(stderr, "\t-> No input file specified\n");
         usage(stdout, 0);
         return 0;
     }

merge_bdamage.cpp

@@ -34,8 +34,8 @@ void merge_bdamage(const std::vector<std::string> &bdamage_files, const char* ou
 
   // now merge all but the first one, into the first one. Either adding or making new ids
   fprintf(stderr,"\t-> Doing merging old size: %lu\n",myvec[0].size());
-  for(int i=1;i<myvec.size();i++){
-    fprintf(stderr,"\t-> Merging %d into 0\n",i);
+  for(size_t i=1;i<myvec.size();i++){
+    fprintf(stderr,"\t-> Merging %lu into 0\n",i);
     std::map<int, mydataD> &small= myvec[i];
 
     for(auto slave =small.begin();slave!=small.end();slave++){

metaDMG.cpp

@@ -64,7 +64,7 @@ double *getval(std::map<int, double *> &retmap, int2intvec &child, int taxid, in
         int2intvec::iterator it2 = child.find(taxid);
         if (it2 != child.end()) {
             std::vector<int> &avec = it2->second;
-            for (int i = 0; i < avec.size(); i++) {
+            for (size_t i = 0; i < avec.size(); i++) {
                 //	fprintf(stderr,"%d/%d %d\n",i,avec.size(),avec[i]);
                 double *tmp = getval(retmap, child, avec[i], howmany);
                 for (int i = 0; i < 3 * howmany + 1; i++)
@@ -150,12 +150,12 @@ mydata2 getval_stats(std::map<int, mydata2> &retmap, int2intvec &child, int taxi
         int2intvec::iterator it2 = child.find(taxid);
         if (it2 != child.end()) {
             std::vector<int> &avec = it2->second;
-            for (int i = 0; i < avec.size(); i++) {
+            for (size_t i = 0; i < avec.size(); i++) {
                 mydata2 tmp = getval_stats(retmap, child, avec[i]);
                 ret.nreads += tmp.nreads;
             }
 
-            for (int i = 0; i < avec.size(); i++) {
+            for (size_t i = 0; i < avec.size(); i++) {
                 //	fprintf(stderr,"%d/%d %d\n",i,avec.size(),avec[i]);
                 mydata2 tmp = getval_stats(retmap, child, avec[i]);
                 if (tmp.nreads == 0)

misc/MAP.cpp

@@ -154,8 +154,6 @@ N MAX IS  28921
 
 void M3_stat_file(double M3[MAX_ROWS][MAX_COLS],gzFile file){
 
-    int N_max = 0;
-
     // N_x=1_forward  N_x=1_reverse
     gzprintf(file,"%d \t %d \t",(int) M3[0][NCOL],(int)M3[NUMROWS-1][NCOL]);
 
@@ -220,7 +218,7 @@ void M3_stat_file(double M3[MAX_ROWS][MAX_COLS],gzFile file){
 
 const char** getColumnNames(int* colnumber) {
     const int maxColumnNames = 200;
-    const char** columnNames = (const char**)malloc(200 * sizeof(const char*));
+    const char** columnNames = (const char**)malloc(maxColumnNames * sizeof(const char*));
     if (columnNames == NULL) {
         fprintf(stderr,"Failed to allocate memory for column names");
         return NULL;
@@ -384,11 +382,6 @@ void M3Print_to_OutStat(int argc,char **argv){
             free(columnNames);
         }
 
-        char* Id_copy = strdup(M3file);
-        char* sample_id = strtok(Id_copy,".");
-        free(Id_copy);
-
-
         gzprintf(gz,"%s \t %s \t","tmp",tax_id[0]);
         M3_stat_file(M3,gz);
         MAP_stat_file(M3,LlhRes,gz);

misc/dfit_optim.cpp

@@ -58,7 +58,7 @@ double **read1_ugly_matrix(const char *fname){
   ret[0][1] = 0;//<- function counter
   ret[1] = new double[xcol.size()];
   ret[2] = new double[xcol.size()];
-  for(int i=0;i<xcol.size();i++){
+  for(size_t i=0;i<xcol.size();i++){
     ret[0][i+2] = xcol[i];
     ret[1][i] = kcol[i];
     ret[2][i] = ncol[i];
@@ -159,10 +159,10 @@ void getstat(double **dat,double *pars,double *statpars){
     double q = pars[1];
     double c = pars[2];
     double phi = pars[3];
-    double llh = pars[4];
+    //double llh = pars[4];
     int NUMROWS = dat[0][0];
     double *XCOL = dat[0]+2;
-    double *KCOL = dat[1];
+    //double *KCOL = dat[1];
     double *NCOL = dat[2];
     double N_pos = dat[2][0];
 

misc/extract_reads.cpp

@@ -122,7 +122,7 @@ int2int getkeysint(const char *key,int value){
 void doflush(queue *myq,int2int &keeplist,bam_hdr_t *hdr,samFile *outhts,int strict){
   // fprintf(stderr,"flush: %lu strictk:%d\n",myq->l,strict);
   if(strict==1){//will only print specific match
-    for(int i=0;i<myq->l;i++){
+    for(size_t i=0;i<myq->l;i++){
       int2int::iterator it=keeplist.find(myq->ary[i]->core.tid);
       if(it!=keeplist.end())
 	assert(sam_write1(outhts, hdr,myq->ary[i])>=0);
@@ -131,15 +131,15 @@ void doflush(queue *myq,int2int &keeplist,bam_hdr_t *hdr,samFile *outhts,int str
   if(strict==0){
     //if writedata=0 then no aligments will be printed, otherwise all
     int writedata = 0;
-    for(int i=0;i<myq->l;i++){
+    for(size_t i=0;i<myq->l;i++){
       int2int::iterator it=keeplist.find(myq->ary[i]->core.tid);
       if(it!=keeplist.end()){
 	writedata=1;
 	break;
       }
     }
     if(writedata>0){
-      for(int i=0;i<myq->l;i++){
+      for(size_t i=0;i<myq->l;i++){
 	assert(sam_write1(outhts, hdr,myq->ary[i])>=0);
       }
     }
@@ -248,7 +248,7 @@ void gettaxids_to_use(int taxid,int2intvec &child,int2int &i2i){
     fprintf(stderr,"\t-> Problem finding taxid: %d from nodesfile\n",taxid);
   }else{ 
     std::vector<int> &avec = it->second;
-    for(int i=0;i<avec.size();i++){
+    for(size_t i=0;i<avec.size();i++){
       //	fprintf(stderr,"%d/%d %d\n",i,avec.size(),avec[i]);
       gettaxids_to_use(avec[i],child,i2i);