Remove assert statements and general cleanup

Makefile

@@ -128,6 +128,9 @@ misc: $(LIBHTS) $(OBJ)
 # --- Automatisk afhængighedshåndtering ---
 -include $(OBJ:.o=.d)
 
+bfgs.o: bfgs.cpp $(LIBHTS)
+	$(CXX) -c $(CPPFLAGS) $(filter-out -Wall -Wextra,$(CXXFLAGS)) -w $< -o $@
+
 %.o: %.c $(LIBHTS)
 	$(CC) -c $(CPPFLAGS) $(CFLAGS) $< -o $@
 

dfit_helppage.cpp

@@ -2,7 +2,6 @@
 #include <cstring>
 #include <cstdlib>
 #include <ctime>
-#include <cassert>
 #include <cstdint>
 #include <iostream>
 #include <math.h>

dfit_optim.cpp

@@ -1,5 +1,4 @@
 #include <cmath>
-#include <cassert>
 #include <cstring>
 #include <vector>
 #include <cstdio>
@@ -12,7 +11,10 @@
 double **read1_ugly_matrix(const char *fname){
   fprintf(stderr,"\t-> Reading file: \'%s\'\n",fname);
   gzFile gz = Z_NULL;
-  assert((gz=gzopen(fname,"rb"))!=Z_NULL);
+  if((gz=gzopen(fname,"rb"))==Z_NULL){
+    fprintf(stderr,"\t-> Problem opening file: %s will exit\n",fname);
+    exit(1);
+  }
   char buf[4096];
   char *taxid = NULL;
   gzgets(gz,buf,4096);
@@ -24,9 +26,12 @@ double **read1_ugly_matrix(const char *fname){
     char *tok = strtok(buf,"\t\n ");
     if(taxid==NULL)
       taxid = strdup(tok);
-    else
-      assert(strcmp(taxid,tok)==0);
-
+    else{
+      if (strcmp(taxid, tok) != 0) {
+	fprintf(stderr, "\t-> Error: taxid mismatch (%s vs %s) will exit\n", taxid, tok);
+	exit(1);
+      }
+    }
     tok = strtok(NULL,"\t\n ");
     int is5 = 1;
     if(strcmp(tok,"3'")==0)

main_aggregate_stat.cpp

@@ -13,7 +13,6 @@
 #include <sys/time.h>
 #include <math.h>
 #include <algorithm> // for std::sort
-#include <cassert>
 
 #include <iostream>
 #include "profile.h"
@@ -173,7 +172,11 @@ std::map<int,char *> read_dfit(char *fname){
 
   BGZF *fpfpfp = NULL;
   fpfpfp = bgzf_open(fname, "rb");
-  assert(fpfpfp!=NULL);
+  if(fpfpfp==NULL){
+    fprintf(stderr,"\t-> Problem opening file: %s will exit\n",fname);
+    exit(1);
+  }
+
   kstring_t *kstr = new kstring_t;
   kstr->s=NULL;kstr->l=kstr->m = 0;
   while(bgzf_getline(fpfpfp,'\n',kstr)){
@@ -200,25 +203,6 @@ std::map<int,char *> read_dfit(char *fname){
   return ret;
 }
 
-/*
-void aggr_stat3000(std::map<int, mydata2> &stats,int2int &parent){
-  std::map<int,int> dasmap;
-  for(std::map<int,mydata2>::iterator it = stats.begin();it!=stats.end();it++)
-    dasmap[it->first] = it->second.nreads;
-  //fprintf(stderr,"dasmap.size(): %lu stats.size():%lu\n",dasmap.size(), stats.size());
-
-  for(std::map<int,int>::iterator itt=dasmap.begin();itt!=dasmap.end();itt++){
-  //  for(int i=0;i<dasvector.size();i++){
-    int focal_taxid = itt->first;
-    int2int::iterator it = parent.find(focal_taxid);
-    assert(it!=parent.end());
-    int target = it->second;
-    to_root(focal_taxid,target,stats,parent,itt->second);
-  }
-
-}
-
-*/
 void aggr_stat3000(std::map<int, mydata2> &stats,int2int &parent){
   std::map<int,int> dasmap;
   std::map<int,double *> datamap;
@@ -234,7 +218,10 @@ void aggr_stat3000(std::map<int, mydata2> &stats,int2int &parent){
   for(std::map<int,int>::iterator itt=dasmap.begin();itt!=dasmap.end();itt++){
     int focal_taxid = itt->first;
     int2int::iterator it = parent.find(focal_taxid);
-    assert(it!=parent.end());
+    if (it == parent.end()) {
+      fprintf(stderr, "\t-> Error: iterator reached end (key:%d not found),will exit\n",focal_taxid);
+      exit(1);
+    }
     int target = it->second;
 
     if(target!=focal_taxid)
@@ -276,10 +263,26 @@ int main_aggregate(int argc, char **argv) {
       fprintf(stderr,"\t-> --names file.txt.gz must be defined with --nodes is defined\n");
       exit(1);
     }
-    fprintf(stderr,"aggregate infile_bdamage: %s infile_names: %s infile_nodes: %s infile_lcastat: %s infile_dfit: %s outfile_name: %s\n",infile_bdamage,infile_names,infile_nodes,infile_lcastat,infile_dfit,outfile_name);
+    fprintf(stderr,
+	    "aggregate infile_bdamage: %s infile_names: %s infile_nodes: %s infile_lcastat: %s infile_dfit: %s outfile_name: %s\n",
+	    infile_bdamage ? infile_bdamage : "NULL",
+	    infile_names   ? infile_names   : "NULL",
+	    infile_nodes   ? infile_nodes   : "NULL",
+	    infile_lcastat ? infile_lcastat : "NULL",
+	    infile_dfit    ? infile_dfit    : "NULL",
+	    outfile_name   ? outfile_name   : "NULL"
+	    );
     if(outfile_name==NULL)
       outfile_name = strdup(infile_bdamage);
-    fprintf(stderr,"aggregate infile_bdamage: %s infile_names: %s infile_nodes: %s infile_lcastat: %s infile_dfit: %s outfile_name: %s\n",infile_bdamage,infile_names,infile_nodes,infile_lcastat,infile_dfit,outfile_name);
+    fprintf(stderr,
+	    "aggregate infile_bdamage: %s infile_names: %s infile_nodes: %s infile_lcastat: %s infile_dfit: %s outfile_name: %s\n",
+	    infile_bdamage ? infile_bdamage : "NULL",
+	    infile_names   ? infile_names   : "NULL",
+	    infile_nodes   ? infile_nodes   : "NULL",
+	    infile_lcastat ? infile_lcastat : "NULL",
+	    infile_dfit    ? infile_dfit    : "NULL",
+	    outfile_name   ? outfile_name   : "NULL"
+	    );
     char buf[1024];
     snprintf(buf, 1024, "%s.stat.gz", outfile_name);
     fprintf(stderr, "\t-> Dumping file: \'%s\'\n", buf);
@@ -328,7 +331,10 @@ int main_aggregate(int argc, char **argv) {
     }
     if(dfit_int_char.size()>0){
       std::map<int, char *>::iterator it = dfit_int_char.find(-1);
-      assert(it!=dfit_int_char.end());
+      if (it == dfit_int_char.end()) {
+	fprintf(stderr, "\t-> Error: iterator reached end in dfit_int_char, will exit\n");
+	exit(1);
+      }
       ksprintf(kstr,"\t%s",it->second);
     }
     ksprintf(kstr,"\n");

main_dfit.cpp

@@ -79,6 +79,7 @@ void make_bootstrap_data(double **in,double **out,int howmany,mrand_t *rand_allo
 }
 
 void print_dat(double **dat, int howmany, int libprep) {
+  (void) libprep;
     int rows = 4;  // [metadata, kvec, nvec, freq]
     int cols = 2 * howmany; 
     printf("printing d:\n");
@@ -503,7 +504,13 @@ void slave_block(std::map<int, mydataD> &retmap,int howmany,sam_hdr_t *hdr,int2c
     }
 
     // Sigma and Z
-    double stats[2+2*(int)dat[0][0]];
+    int n = 2 + 2 * (int)dat[0][0];
+
+    double *stats = (double *)malloc(n * sizeof(double));
+    if (stats == NULL) {
+      fprintf(stderr, "\t-> Error: failed to allocate memory for stats, will exit\n");
+      exit(1);
+    }
     getstat(dat,pars,stats);
 
     // stats contains standard deviation, then significance, then the calculated Dx for each position then the normalized 
@@ -626,7 +633,7 @@ void slave_block(std::map<int, mydataD> &retmap,int howmany,sam_hdr_t *hdr,int2c
     }
     ksprintf(kstr,"\n");
     delete[] CI_val;  // Free the array of pointers
-
+    free(stats);
   }
 
   // Free the memory
@@ -772,9 +779,9 @@ int main_dfit(int argc, char **argv) {
       snprintf(bootbuf, 1024, "%s.boot.stat.gz", outfile_name);
       bootfp = bgzf_open(bootbuf, "wb");    
       ksprintf(bootkstr,"taxid\tA_b\tq_b\tc_b\tphi_b\n");
-      if(bgzf_write(bootfp,bootkstr->s,bootkstr->l)!=bootkstr->l){
+      if(bgzf_write(bootfp,bootkstr->s,bootkstr->l)!=(ssize_t)bootkstr->l){
 	fprintf(stderr,"\t-> Problem writing boot file\n");
-	exit(0);
+	exit(1);
       }
       bootkstr->l = 0;
     }
@@ -884,7 +891,12 @@ int main_dfit(int argc, char **argv) {
           dings[i].bootkstr = bootkstr;
           dings[i].showfits = showfits;
         }
-        pthread_t mythreads[nthreads];
+
+	pthread_t *mythreads = (pthread_t *)malloc(nthreads * sizeof(pthread_t));
+	if (mythreads == NULL) {
+	  fprintf(stderr, "\t-> Error: failed to allocate memory for threads, will exit\n");
+	  exit(1);
+	}
         for(size_t i=0;i<nthreads;i++)
           pthread_create(&mythreads[i],NULL,slaveslave, &dings[i]);
 
@@ -898,7 +910,7 @@ int main_dfit(int argc, char **argv) {
       }
     }
 
-    if(bgzf_write(fpfpfp,kstr->s,kstr->l) != kstr->l){
+    if(bgzf_write(fpfpfp,kstr->s,kstr->l) != (ssize_t) kstr->l){
       fprintf(stderr, "\t-> Problems write to output BGZ file\n");
       exit(1);
     }
@@ -907,7 +919,7 @@ int main_dfit(int argc, char **argv) {
 
 
     if(doboot>0){
-      if(bgzf_write(bootfp,bootkstr->s,bootkstr->l) != bootkstr->l){
+      if(bgzf_write(bootfp,bootkstr->s,bootkstr->l) != (ssize_t) bootkstr->l){
         fprintf(stderr, "\t-> Problemst write to output BGZ file\n");
         exit(1);
       }

main_pmd.cpp

@@ -15,10 +15,6 @@
 
 #include "profile.h"  // for reconstructRefWithPosHTS
 
-//extern htsFormat *dingding2;
-//htsFormat *dingding3 = dingding2;
-
-
 int nproc = 0;  // number of reads processed
 
 extern char refToChar[256];
@@ -34,7 +30,7 @@ double pmd_stat(char *seq, char *ref, int len, char *qs) {
     // make revcom
     static char compref1[512];
     for (int i = 0; i < len; i++)
-        compref1[i] = revTable[ref[len - i - 1]];
+      compref1[i] = revTable[(unsigned char)ref[len - i - 1]];
 
     // do flipcount
     int count[2] = {0, 0};
@@ -45,17 +41,11 @@ double pmd_stat(char *seq, char *ref, int len, char *qs) {
             count[1] = count[1] + 1;
     }
 
-    char *myref = NULL;
-    if (count[0] > count[1])
-        myref = compref1;
-    else
-        myref = ref;
-
     // do stat
     double llh1 = 0;
     double llh2 = 0;
     for (int i = 0; i < len; i++) {
-        double eps = qsToProp[qs[i]];
+      double eps = qsToProp[(unsigned char)qs[i]];
         if (ref[i] == 'C') {
             double Dz = C0 + p * pow(1 - p, i + 1);
             double pm1 = (1 - ppi) * (1 - eps) * (1 - Dz) + ppi * eps / 3 * (1 - Dz) + ppi / 3 * (1 - eps) * Dz + (1 - ppi) * eps / 3 * Dz;
@@ -83,7 +73,7 @@ double pmd_stat(char *seq, char *ref, int len, char *qs) {
     return llh1 - llh2;
 }
 
-void wrapper(const bam1_t *b, const char *reconstructedReference, const std::vector<int> &reconstructedReferencePos, const int &minQualBase, int MAXLENGTH, float **mm5p, float **mm3p, float incval, char myread[512], char myref[512]) {
+void wrapper(const bam1_t *b, const char *reconstructedReference, const std::vector<int> /*&reconstructedReferencePos*/, const int &minQualBase, int /*MAXLENGTH*/, float /***mm5p*/, float /***mm3p*/, float /*incval*/, char myread[512], char myref[512]) {
     const char *alphabetHTSLIB = "NACNGNNNTNNNNNNN";
     char refeBase;
     char readBase;
@@ -225,7 +215,7 @@ void parse_sequencingdata(char *refName, char *fname, int mapped_only, int se_on
         memset(myread, 'N', 512);
         memset(myrefe, 'N', 512);
         reconstructRefWithPosHTS(b, mypair, reconstructedRef);
-        wrapper(b, mypair.first->s, mypair.second, 0, 0, NULL, NULL, 1, myread, myrefe);
+        wrapper(b, mypair.first->s, mypair.second, 0, 0, 0, 0, 1, myread, myrefe);
 #if 0
     fprintf(stderr,"\nmyread:\n%.*s\nmyReference:\n%.*s\n",b->core.l_qseq,myread,b->core.l_qseq,myrefe);
     fprintf(stderr,"---read[%d]----\n",nproc);
@@ -237,7 +227,6 @@ void parse_sequencingdata(char *refName, char *fname, int mapped_only, int se_on
     }
 
     bam_destroy1(b);
-    free(dingding3);
     sam_hdr_destroy(hdr);
     sam_close(in);
     free(fname);
@@ -359,6 +348,7 @@ int main_pmd(int argc, char **argv) {
             "\t[ALL done] cpu-time used =  %.2f sec\n"
             "\t[ALL done] walltime used =  %.2f sec\n",
             (float)(clock() - t) / CLOCKS_PER_SEC, (float)(time(NULL) - t2));
+    (void)nthreads;//<- avoid compiler complain.
     return 0;
 }