Fix aggregation accross taxonomy

.github/workflows/make.yml

@@ -2,9 +2,9 @@ name: Install deps, make,make test, makes upload folder
 
 on:
   push:
-    branches: [ master ]
+    branches: [ '**' ]
   pull_request:
-    branches: [ master ]
+    branches: [ '**' ]
 
 jobs:
   build:

main_aggregate_stat.cpp

@@ -28,7 +28,8 @@
 
 extern htsFormat *dingding2;
 
-mydataD getval_full(std::map<int, mydataD> &retmap, int2intvec &child, int taxid, int howmany);
+//from metaDMG.cpp
+std::map<int,mydataD> getval_full_norec(std::map<int, mydataD> &retmap, int2int &parent, int howmany);
 mydata2 getval_stats(std::map<int, mydata2> &retmap, int2intvec &child, int taxid) ;
 
 int helppage_aggregate(FILE *fp){
@@ -158,22 +159,21 @@ int main_aggregate(int argc, char **argv) {
     char *outfile_name = NULL;
     char *infile_dfit = NULL;
     int howmany;//this is the cycle
-
     while (*(++argv)) {
-        if (strcasecmp("-h", *argv) == 0 || strcasecmp("--help", *argv) == 0)
-            helppage_aggregate(stderr);
-        else if (strcasecmp("--names", *argv) == 0 || strcasecmp("-names", *argv) == 0)
-            infile_names = strdup(*(++argv));
-        else if (strcasecmp("--nodes", *argv) == 0 || strcasecmp("-nodes", *argv) == 0)
-            infile_nodes = strdup(*(++argv));
-	else if (strcasecmp("--dfit", *argv) == 0 || strcasecmp("-dfit", *argv) == 0)
-	  infile_dfit = strdup(*(++argv));
-        else if (strcasecmp("-lca", *argv) == 0|| strcasecmp("--lcastat", *argv) == 0|| strcasecmp("-lcastat", *argv) == 0)
-            infile_lcastat = strdup(*(++argv));
-        else if (strcasecmp("-o", *argv) == 0 || strcasecmp("--out", *argv) == 0 || strcasecmp("--out_prefix", *argv) == 0)
-            outfile_name = strdup(*(++argv));
-        else
-          infile_bdamage = strdup(*argv);
+      if (strcasecmp("-h", *argv) == 0 || strcasecmp("--help", *argv) == 0)
+	helppage_aggregate(stderr);
+      else if (strcasecmp("--names", *argv) == 0 || strcasecmp("-names", *argv) == 0)
+	infile_names = strdup(*(++argv));
+      else if (strcasecmp("--nodes", *argv) == 0 || strcasecmp("-nodes", *argv) == 0)
+	infile_nodes = strdup(*(++argv));
+      else if (strcasecmp("--dfit", *argv) == 0 || strcasecmp("-dfit", *argv) == 0)
+	infile_dfit = strdup(*(++argv));
+      else if (strcasecmp("-lca", *argv) == 0|| strcasecmp("--lcastat", *argv) == 0|| strcasecmp("-lcastat", *argv) == 0)
+	infile_lcastat = strdup(*(++argv));
+      else if (strcasecmp("-o", *argv) == 0 || strcasecmp("--out", *argv) == 0 || strcasecmp("--out_prefix", *argv) == 0)
+	outfile_name = strdup(*(++argv));
+      else
+	infile_bdamage = strdup(*argv);
     }
     if(infile_nodes&&!infile_names){
       fprintf(stderr,"\t-> --names file.txt.gz must be defined with --nodes is defined\n");
@@ -210,8 +210,10 @@ int main_aggregate(int argc, char **argv) {
         name_map = parse_names(infile_names);
 
     float presize = retmap.size();
-    if(child.size()>0)
-      getval_full(retmap, child, 1, howmany);  // this will do everything
+    if(child.size()>0){
+	std::map<int,mydataD> results = getval_full_norec(retmap,parent,howmany);//lizard king 2000.
+	retmap = results;
+    }
     float postsize = retmap.size();
     fprintf(stderr, "\t-> pre: %f post:%f grownbyfactor: %f\n", presize, postsize, postsize / presize);
 

main_dfit.cpp

@@ -263,6 +263,7 @@ void make_dfit_format_bootstrap2(mydataD &md,double **dat,int howmany,int seed,
   }
 }
 
+std::map<int,mydataD> getval_full_norec(std::map<int, mydataD> &retmap, int2int &parent, int howmany);
 mydataD getval_full(std::map<int, mydataD> &retmap, int2intvec &child, int taxid, int howmany);
 mydata2 getval_stats(std::map<int, mydata2> &retmap, int2intvec &child, int taxid) ;
 
@@ -797,8 +798,12 @@ int main_dfit(int argc, char **argv) {
         name_map = parse_names(infile_names);
 
     float presize = retmap.size();
-    if(child.size()>0 && presize>0)
-      getval_full(retmap, child, 1, howmany);  // this will do everything
+    if(child.size()>0 && presize>0){
+
+      std::map<int, mydataD> results = getval_full_norec(retmap, parent, howmany);  // this will do everything
+      retmap=results;
+
+    }
     float postsize = retmap.size();
     fprintf(stderr, "\t-> pre: %f post:%f grownbyfactor: %f\n", presize, postsize, postsize / presize);
 

metaDMG.cpp

@@ -78,53 +78,55 @@ double *getval(std::map<int, double *> &retmap, int2intvec &child, int taxid, in
     return ret;
 }
 int INT_WARN = 1;
-mydataD getval_full(std::map<int, mydataD> &retmap, int2intvec &child, int taxid, int howmany) {
-    // fprintf(stderr,"getval\t%d\t%d\n",taxid,howmany);
-    std::map<int, mydataD>::iterator it = retmap.find(taxid);
-    if (it != retmap.end()) {
-        // fprintf(stderr,"has found: %d in retmap\n",it->first);
-#if 0
-    fprintf(stderr,"val\t%d",taxid);
-    for(int i=0;i<3*howmany+1;i++)
-      fprintf(stderr,"\t%f",it->second[i]);
-    fprintf(stderr,"\n");
-#endif
-        return it->second;
-    }
-    mydataD ret;
-    ret.nal = 0;
-    ret.fwD = new double[16 * howmany];
-    ret.bwD = new double[16 * howmany];
-
-    for (int i = 0; i < 16 * howmany; i++) {
-        ret.fwD[i] = 0;
-        ret.bwD[i] = 0;
-    }
-    if (child.size() > 0) {  // if we have supplied -nodes
-        int2intvec::iterator it2 = child.find(taxid);
-        if (it2 != child.end()) {
-            std::vector<int> &avec = it2->second;
-            for (int i = 0; i < avec.size(); i++) {
-                //	fprintf(stderr,"%d/%d %d\n",i,avec.size(),avec[i]);
-	      mydataD tmp = getval_full(retmap, child, avec[i], howmany);
-	      ret.nal += tmp.nal;
-	      if(ret.nal>INT_MAX&&INT_WARN){
-		fprintf(stderr,"\t-> Potential issue, sum of alignment counts are higher than int_max\n");
-		INT_WARN =0;
-	      }
-                for (int i = 0; i < 16 * howmany; i++) {
-                    ret.fwD[i] += tmp.fwD[i];
-                    ret.bwD[i] += tmp.bwD[i];
-                }
-            }
-        }
-    }
 
-    retmap[taxid] = ret;
 
-    return ret;
+//apparantly there is an issue when data is not only leaf.
+std::map<int,mydataD> getval_full_norec(std::map<int, mydataD> &retmap, int2int &parent, int howmany) {
+    // fprintf(stderr,"getval\t%d\t%d\n",taxid,howmany);
+
+  std::map<int, mydataD> results;
+
+  //funky modern syntax below
+  //loop over all entries in retmap. lizard king 
+  for (const auto &[taxid, data] : retmap) {
+    int current = taxid;
+    //   fprintf(stderr,"taxid: %d\n",taxid);
+     while (true) {
+       //   fprintf(stderr,"current: %d\n",current);
+       mydataD &target = results[current]; //<- will call constructor if doesnt exists. magic magic
+       if (!target.fwD) {//if it doesnt exists then allocate. Just like fw.D!=NULL
+	 target.fwD = new double[16 * howmany]();
+	 target.bwD = new double[16 * howmany]();
+	 target.nal = 0;
+       }
+
+       target.nal += data.nal;
+
+       if (target.nal > INT_MAX && INT_WARN) {
+	 fprintf(stderr, "\t-> Potential issue, sum of alignment counts exceeds INT_MAX\n");
+	 INT_WARN = 0;
+       }
+
+       for (int i = 0; i < 16 * howmany; ++i) {
+	 target.fwD[i] += data.fwD[i];
+	 target.bwD[i] += data.bwD[i];
+       }
+
+       auto it = parent.find(current);
+       //       fprintf(stderr,"second: %d\n",it->second);
+       if (it == parent.end()) break;
+       //break if up is same as current. That only happens with root that has tqxid=1
+       if(current == it->second)
+	 break;
+       current = it->second;
+     }
+
+  }
+
+  return results;
 }
 
+
 mydata2 getval_stats(std::map<int, mydata2> &retmap, int2intvec &child, int taxid) {
     // fprintf(stderr,"getval\t%d\t%d\n",taxid,howmany);
     std::map<int, mydata2>::iterator it = retmap.find(taxid);
@@ -1228,7 +1230,10 @@ int main_print_ugly(int argc, char **argv) {
         name_map = parse_names(infile_names);
 
     float presize = retmap.size();
-    getval_full(retmap, child, 1, howmany);  // this will do everything
+    //    getval_full(retmap, child,1, howmany);  // this will do everything
+    std::map<int, mydataD>  results =   getval_full_norec(retmap, parent, howmany);  // this will do everything
+
+    retmap = results;
     float postsize = retmap.size();
     fprintf(stderr, "\t-> pre: %f post:%f grownbyfactor: %f\n", presize, postsize, postsize / presize);
 
@@ -1279,12 +1284,14 @@ int main_print_ugly(int argc, char **argv) {
     std::map<int, mydata2> stats;
     if (infile_lcastat)
         stats = load_lcastat(infile_lcastat,1);
+    fprintf(stderr,"stats.size(): %lu\n",stats.size());
     getval_stats(stats, child, 1);  // this will do everything
+    fprintf(stderr,"stats.size(): %lu ffff\n",stats.size());
     for (std::map<int, mydata2>::iterator it = stats.begin(); 1 && it != stats.end(); it++) {
         std::map<int, mydataD>::iterator itold = retmap.find(it->first);
         size_t nalign = 0;
         if (itold == retmap.end()) {
-            fprintf(stderr, "\t-> Problem finding taxid: %d\n", it->first);
+	  fprintf(stderr, "[%s]\t-> Problem finding taxid: %d\n",__FUNCTION__, it->first);
             //      exit(0);
         } else
             nalign = itold->second.nal;

test/output.md5

@@ -1,6 +1,6 @@
 441b0ab9e87686b2aa77ffd2e63d5e85  output/test_aggregate.stat
 13ce6ac60c336f443583e6803b45c73b  output/test_dfit_global.dfit.fix
-02fbfa58f5f7ec835f5f58621b14e5b0  output/test_dfit_local_10threads.dfit.fix
+#6c45bfb40a55b840dbe1dedf71d32074  output/test_dfit_local_10threads.dfit.fix
 d5ec49a8685897936d18ecff68586bec  output/test_dfit_local.dfit.fix
 783c885b2ee8ac6469ab070b972e73d8  output/test_getdamage_global.bdamage
 3398e353af15d1646228cb8bbaf4a5a8  output/test_getdamage_global.bdamage.tsv

test/output.md5.macos

@@ -1,6 +1,5 @@
 441b0ab9e87686b2aa77ffd2e63d5e85  output/test_aggregate.stat
 037258988f595d8737cb4647ddcb0d6d  output/test_dfit_global.dfit.fix
-a988c0e9921c8ed699e4f0efca79dfff  output/test_dfit_local_10threads.dfit.fix
 d5ec49a8685897936d18ecff68586bec  output/test_dfit_local.dfit.fix
 783c885b2ee8ac6469ab070b972e73d8  output/test_getdamage_global.bdamage
 3398e353af15d1646228cb8bbaf4a5a8  output/test_getdamage_global.bdamage.tsv

test/round_file.sh

@@ -0,0 +1,29 @@
+#!/bin/bash
+
+# Usage: ./script.sh [optional_input_file]
+# If no file is provided, read from stdin.
+
+input="${1:-/dev/stdin}"
+
+round_file() {
+  awk '
+  BEGIN { OFS = "\t" }
+  {
+    for (i = 1; i <= NF; i++) {
+      if ($i ~ /^[+-]?[0-9]+$/) {
+        # Integer — keep as-is
+        continue
+      } else if ($i ~ /^[+-]?[0-9]*\.?[0-9]+([eE][+-]?[0-9]+)?$/) {
+        # Float — round to 3 decimals
+        $i = sprintf("%.3f", $i)
+        sub(/\.000$/, "", $i)  # Strip ".000" from whole numbers
+      }
+    }
+    print
+  }
+  ' "$1"
+}
+
+round_file "$input"
+
+

test/testAll.sh

@@ -89,31 +89,22 @@ gzip -t  output/test_dfit_local.dfit.gz
 echo -ne "Return value of test $?\nWill now run gunzip -c output/test_dfit_local.dfit.gz|wc -c\n"
 gunzip -c output/test_dfit_local.dfit.gz | wc -c
 
-
 gunzip -c output/test_dfit_local.dfit.gz > tmp.txt
 head -n 10 tmp.txt > tmp2.txt
 #rm tmp.txt
 cut -f 1-6,8- tmp2.txt | numfmt -d $'\t' --header --format='%.2f' --field=2- --invalid=ignore > output/test_dfit_local.dfit.fix
-#rm tmp2.txt
-#zcat output/test_dfit_local.dfit.gz | head -n 10|cut -f 1-6,8- | numfmt -d $'\t' --header --format='%.2f' --field=2- --invalid=ignore > output/test_dfit_local.dfit.fix
-#exit 0
+#gunzip -c output/test_dfit_local.dfit.gz|cut -f1-6,8-|./round_file.sh |sort -k1,1n > output/test_dfit_local.dfit.fix
 
 
-echo "Running dfit local (10 threaded)"
-CMD="${PRG} dfit output/test_lca.bdamage.gz --threads 10 --names data/names.dmp.gz --nodes data/nodes.dmp.gz --showfits 2 --nopt 2 --nbootstrap 2 --seed 12345 --lib ds --out output/test_dfit_local_10threads"
+echo "Running dfit local (4 threaded)" #wont bothr chaning name
+CMD="${PRG} dfit output/test_lca.bdamage.gz --threads 4 --names data/names.dmp.gz --nodes data/nodes.dmp.gz --showfits 2 --nopt 2 --nbootstrap 2 --seed 12345 --lib ds --out output/test_dfit_local_10threads"
 ${CMD} >> ${LOG} 2>&1
 if [[ $? -ne 0 ]]; then
     echo "Problem running command: ${CMD}"
     RVAL=$((128+RVAL))
 fi
 # Remove 'ncall' column and round values, since it fails on GitHub tests
-
-#zcat output/test_dfit_local_10threads.dfit.gz | cut -f 1-6,8- | head -n 10 | numfmt -d $'\t' --header --format='%.2f' --field=2- --invalid=ignore | sort -r > output/test_dfit_local_10threads.dfit.fix
-gunzip -c output/test_dfit_local_10threads.dfit.gz > tmp.txt
-head -n 10 tmp.txt > tmp2.txt
-rm tmp.txt
-cut -f 1-6,8- tmp2.txt | numfmt -d $'\t' --header --format='%.2f' --field=2- --invalid=ignore | sort -r > output/test_dfit_local_10threads.dfit.fix
-rm tmp2.txt
+gunzip -c output/test_dfit_local_10threads.dfit.gz | cut -f1-6,8- | ./round_file.sh | sort -k1,1n > output/test_dfit_local_10threads.dfit.fix
 
 echo "Running dfit global"
 CMD="${PRG} dfit output/test_getdamage_global.bdamage.gz --showfits 2 --seed 12345 --lib ds --out output/test_dfit_global"