diff --git a/.github/workflows/make.yml b/.github/workflows/make.yml
index 071b491..dfd4ed1 100644
--- a/.github/workflows/make.yml
+++ b/.github/workflows/make.yml
@@ -25,4 +25,7 @@ jobs:
       uses: actions/upload-artifact@v4
       with:
         name: test-output
-        path: test/output/
+        path: |
+          test/output
+          test/testAll.sh
+          test/testAll.sh.log
diff --git a/.gitmodules b/.gitmodules
deleted file mode 100644
index b6e1f14..0000000
--- a/.gitmodules
+++ /dev/null
@@ -1,3 +0,0 @@
-[submodule "htslib"]
-	path = htslib
-	url = https://github.com/samtools/htslib
diff --git a/Makefile b/Makefile
index fa754e7..e181c43 100644
--- a/Makefile
+++ b/Makefile
@@ -1,144 +1,138 @@
-# make FLAGS=-D__REGRESSION__
-#c++ -I /opt/homebrew/Cellar/eigen/3.4.0_1/include/ -c -I/Users/fvr124/metaDMG-cpp/htslib   regression.cpp -std=c++14  -D__REGRESSION__
+# --- Flags og kompileringsvalg ---
+FLAGS     := -O2
+CFLAGS    := $(FLAGS)
+CXXFLAGS  := $(FLAGS)
+CPPFLAGS  := $(CPPFLAGS) -Wall -Wextra
+LDFLAGS   := -lgsl
 
 CC  ?= gcc
 CXX ?= g++
 
-LIBS = -lz -lm -lbz2 -llzma -lpthread -lcurl -lgsl -lgslcblas
-
-CRYPTO_TRY=$(shell echo 'int main(){}'|$(CXX) -x c++ - -lcrypto 2>/dev/null -o /dev/null; echo $$?)
-ifeq "$(CRYPTO_TRY)" "0"
-$(info Crypto library is available to link; adding -lcrypto to LIBS)
-LIBS += -lcrypto
-else
-$(info Crypto library is not available to link; will not use -lcrypto)
-endif
+CXXSRC := $(wildcard *.cpp)
+CSRC   := $(wildcard *.c)
+OBJ    := $(CSRC:.c=.o) $(CXXSRC:.cpp=.o)
 
 
-#if htslib source is defined
-ifdef HTSSRC
-
-#if hts source is set to systemwide
-ifeq ($(HTSSRC),systemwide)
-$(info HTSSRC set to systemwide; assuming systemwide installation)
-LIBS += -lhts
-
-else
-
-#if hts source path is given
-# Adjust $(HTSSRC) to point to your top-level htslib directory
-$(info HTSSRC defined: $(HTSSRC))
-CPPFLAGS += -I"$(realpath $(HTSSRC))"
-LIBHTS := $(HTSSRC)/libhts.a
-LIBS := $(LIBHTS) $(LIBS)
 
+# --- Brugervalgte paths og biblioteker ---
+ifdef PREFIX
+CPPFLAGS += -I$(PREFIX)/include
+LDFLAGS  := -L$(PREFIX)/lib $(LDFLAGS)
 endif
 
-#if htssrc not defined
-else
-
-$(info HTSSRC not defined; using htslib submodule)
-$(info Use `make HTSSRC=/path/to/htslib` to build metadamage using a local htslib installation)
-$(info Use `make HTSSRC=systemwide` to build metadamage using the systemwide htslib installation)
-
-
-HTSSRC := $(CURDIR)/htslib
-CPPFLAGS += -I$(HTSSRC)
-LIBHTS := $(HTSSRC)/libhts.a
-LIBS := $(LIBHTS) $(LIBS)
-
-all: .activate_module
-
+ifdef EXTRA_INC
+CPPFLAGS += $(addprefix -I,$(EXTRA_INC))
 endif
 
-.PHONY: .activate_module
-
-.activate_module:
-	git submodule update --init --recursive
-	$(MAKE) -C $(HTSSRC)
-
-
-
-#modied from htslib makefile
-FLAGS = -O3
-CPPFLAGS := $(filter-out -DNDEBUG,$(CPPFLAGS))
-FLAGS2 = $(CPPFLAGS) $(FLAGS) $(LDFLAGS)
-
-CFLAGS := $(FLAGS2) $(CFLAGS)
-CXXFLAGS := $(FLAGS2) $(CXXFLAGS)
+ifdef EXTRA_LIB
+LDFLAGS  := $(addprefix -L,$(EXTRA_LIB)) $(LDFLAGS)
+endif
 
-OPT=-Wl,-O2
-# filter -O2 optimization flag CXXFLAGS
-CXXFLAGS := $(filter-out -O2,$(CXXFLAGS))
-CXXFLAGS := $(subst $(OPT),,$(CXXFLAGS))
+ifdef EXTRA_LIBS
+LIBS += $(EXTRA_LIBS)
+endif
 
-CSRC = $(wildcard *.c)
-CXXSRC = $(wildcard *.cpp)
-OBJ = $(CSRC:.c=.o) $(CXXSRC:.cpp=.o)
+# --- Crypto library detektion ---
+HAVE_CRYPTO := $(shell echo 'int main(){}'|$(CXX) -x c++ - -lcrypto -o /dev/null 2>/dev/null && echo 0 || echo 1)
+LIBS := -lz -lm -lbz2 -llzma -lpthread -lcurl
 
-prefix      = /usr/local
-exec_prefix = $(prefix)
-bindir      = $(exec_prefix)/bin
+# --- Mål og bygning ---
+PROGRAM = metaDMG-cpp
+all: version.h $(PROGRAM) misc
 
-INSTALL = install
-INSTALL_DIR = $(INSTALL) -dm0755
-INSTALL_PROGRAM = $(INSTALL) -m0755
+ifeq ($(HAVE_CRYPTO),0)
+  $(info Crypto library is available to link; adding -lcrypto)
+  LIBS += -lcrypto
+else
+  $(info Crypto library is not available to link; skipping -lcrypto)
+endif
 
-PROGRAMS = metaDMG-cpp
+# --- HTSLIB håndtering ---
+# HTSSRC is an absolute path (e.g., $(CURDIR)/htslib or user-specified)
 
-all: $(PROGRAMS) misc
+ifndef HTSSRC
+  $(info HTSSRC not defined; cloning htslib from GitHub)
+  HTSSRC := $(CURDIR)/htslib
+  ABSPATH=$(HTSSRC) #donkykong
+endif
 
-PACKAGE_VERSION  = 0.4
+ABSPATH=$(HTSSRC) #donkykong
+ifeq ($(HTSSRC),systemwide)
+  $(info HTSSRC set to systemwide; using systemwide installation)
+  LIBS += -lhts
+  LIBHTS :=
+else
+  # Use HTSSRC directly for include path
+  CPPFLAGS += -I$(HTSSRC)
+  LIBHTS := $(HTSSRC)/libhts.a
+  LIBS := $(LIBHTS) $(LIBS)
+  $(PROGRAM): $(LIBHTS)
+
+  ifneq ($(filter /%,$(HTSSRC)),$(HTSSRC))
+    ABSPATH=../$(HTSSRC)
+  endif
+endif
 
-ifneq "$(wildcard .git)" ""
-PACKAGE_VERSION := $(shell git describe --always --dirty)
-version.h: $(if $(wildcard version.h),$(if $(findstring "$(PACKAGE_VERSION)",$(shell cat version.h)),,force))
+# Ensure htslib is cloned and built only if libhts.a is missing
+$(LIBHTS): .clone_htslib
+
+.clone_htslib:
+	@if [ ! -d "$(HTSSRC)" ]; then \
+		echo "Cloning htslib into $(HTSSRC) with submodules..."; \
+		git clone --recursive https://github.com/samtools/htslib.git $(HTSSRC) || { echo "Clone failed"; exit 1; }; \
+	else \
+		echo "$(HTSSRC) already exists, skipping clone."; \
+	fi
+	@if [ ! -f "$(LIBHTS)" ]; then \
+		$(MAKE) -C $(HTSSRC) libhts.a || { echo "Failed to build libhts.a"; exit 1; }; \
+	else \
+		echo "$(LIBHTS) already exists, skipping build."; \
+	fi
+	$(info CPPFLAGS=$(CPPFLAGS) HTSSRC=$(HTSSRC) (absolute path))
+
+# --- Versionsnummer og version.h ---
+PACKAGE_VERSION := 0.4
+ifneq ("$(wildcard .git)","")
+  PACKAGE_VERSION := $(shell git describe --always --dirty)
 endif
 
 version.h:
-	echo '#define METADAMAGE_VERSION "$(PACKAGE_VERSION)"' > $@
+	@echo '#define METADAMAGE_VERSION "$(PACKAGE_VERSION)"' > version.h.tmp
+	@cmp -s version.h.tmp version.h || mv version.h.tmp version.h
+	@rm -f version.h.tmp
 
-.PHONY: all clean install install-all install-misc misc test
+$(PROGRAM): version.h $(OBJ) $(LIBHTS)
+	$(CXX) $(CXXFLAGS) $(CPPFLAGS) -o $@ $(OBJ) $(LIBS) $(LDFLAGS)
 
-misc:
-	$(MAKE) -C misc HTSSRC="$(realpath $(HTSSRC))"
+.PHONY: misc
+misc: $(LIBHTS) $(OBJ)
+	$(MAKE) -C misc HTSSRC=$(ABSPATH)
 
+# --- Automatisk afhængighedshåndtering ---
 -include $(OBJ:.o=.d)
 
-%.o: %.c
-	$(CC) -c  $(CFLAGS) $*.c
-	$(CC) -MM $(CFLAGS) $*.c >$*.d
-
-%.o: %.cpp
-	$(CXX) -c  $(CXXFLAGS) $*.cpp
-	$(CXX) -MM $(CXXFLAGS) $*.cpp >$*.d
+%.o: %.c $(LIBHTS)
+	$(CC) -c $(CPPFLAGS) $(CFLAGS) $< -o $@
+	$(CC) -MM $(CPPFLAGS) $(CFLAGS) $< > $*.d
 
+%.o: %.cpp $(LIBHTS)
+	$(CXX) -c $(CPPFLAGS) $(CXXFLAGS) $< -o $@
+	$(CXX) -MM $(CPPFLAGS) $(CXXFLAGS) $< > $*.d
 
-metaDMG-cpp: version.h $(OBJ)
-	$(CXX) $(FLAGS) -o metaDMG-cpp *.o $(LIBS)
+# --- Rens og tests ---
+.PHONY: clean test testclean force
 
+clean: 
+	rm -f *.o *.d $(PROGRAM) version.h *~
+	rm -rf $(HTSSRC)
+	$(MAKE) -C misc clean
 
 testclean:
-	rm -f test/acc2taxid.map.gzf570b1db7c.dedup.filtered.rname.bam.bin
-	rm -f test/data/f570b1db7c.dedup.filtered.rname.bam
+	rm -f test/*.bam.bin test/data/*.bam
 	rm -rf test/output test/logfile version.h
 
-clean: testclean
-	rm  -f *.o *.d $(PROGRAMS) version.h *~
-	$(MAKE) -C misc clean
-
 test:
-	echo "Running unit tests for metaDMG"
-	cd test;./testAll.sh
+	@echo "Running unit tests for metaDMG"
+	cd test ; ./testAll.sh
 
 force:
-
-install: all
-	$(INSTALL_DIR) $(DESTDIR)$(bindir)
-	$(INSTALL_PROGRAM) $(PROGRAMS) $(DESTDIR)$(bindir)
-	$(MAKE) -C misc HTSSRC="$(realpath $(HTSSRC))" install
-
-install-misc: misc
-	$(MAKE) -C misc HTSSRC="$(realpath $(HTSSRC))" install-misc
-
-install-all: install install-misc
diff --git a/htslib b/htslib
deleted file mode 160000
index e13611a..0000000
--- a/htslib
+++ /dev/null
@@ -1 +0,0 @@
-Subproject commit e13611a942095dfea0944fe28836934c17c2ef6c
diff --git a/main_aggregate_stat.cpp b/main_aggregate_stat.cpp
index 9590b0a..9d135ed 100644
--- a/main_aggregate_stat.cpp
+++ b/main_aggregate_stat.cpp
@@ -66,11 +66,11 @@ void to_root(int from,int to,std::map<int,mydata2> &stats,int2int &parent,int nr
     mydata2 &md2 = stats.find(to)->second;
     md2.data[0] = ((double) md1.data[0]*md1.nreads+md2.data[0]*md2.nreads)/((double) md1.nreads+md2.nreads); //weighted mean of read length
     md2.data[2] = ((double) md1.data[2]*md1.nreads+md2.data[2]*md2.nreads)/((double) md1.nreads+md2.nreads); //weighted mean of gc
-    
-    double variance1;
-    int nreads1;
-    double variance2;
-    int nreads2;
+
+#if 0
+    double variance1,variance2;
+    int nreads1, nreads2;
+#endif
     if(((double) md1.nreads+md2.nreads)>2){//pooled variance of length and GC
       /*variance1 = ((double) md1.nreads-1)*md1.data[1];
       nreads1 += md1.nreads;
@@ -130,7 +130,7 @@ std::map<int,char *> read_dfit(char *fname){
 }
 
 
-void aggr_stat2000(std::map<int, mydata2> &stats,int2int &parent){
+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;
@@ -217,14 +217,6 @@ int main_aggregate(int argc, char **argv) {
     float postsize = retmap.size();
     fprintf(stderr, "\t-> pre: %f post:%f grownbyfactor: %f\n", presize, postsize, postsize / presize);
 
-    for (std::map<int, mydataD>::iterator it = retmap.begin(); it != retmap.end(); it++) {
-        int taxid = it->first;
-        mydataD md = it->second;
-        if (it->second.nal == 0)
-	  continue;
-	//	fprintf(stderr,"retmap taxid:%d nreads: %d\n",it->first,it->second.nreads);
-    }
-
     std::map<int, mydata2> stats;
     if (infile_lcastat){
       stats = load_lcastat(infile_lcastat,1);
@@ -248,7 +240,7 @@ int main_aggregate(int argc, char **argv) {
     if(child.size()>0)//this will not work if we have data at internal nodes 
       getval_stats(stats, child, 1);  // this will do everything
 #endif
-    aggr_stat2000(stats,parent);
+    aggr_stat3000(stats,parent);
     //    exit(0);
     if(stats.size()>0){
       postsize = stats.size();
diff --git a/merge_bdamage.cpp b/merge_bdamage.cpp
index 9e33ef9..3b592d0 100644
--- a/merge_bdamage.cpp
+++ b/merge_bdamage.cpp
@@ -89,7 +89,7 @@ void merge_bdamage(const std::vector<std::string> &bdamage_files, const char* ou
 void merge_rlens(const std::vector<std::string> &rlens_files, const char* outname) {
   
   // results 
-  std::map<int, std::array<size_t, 200>> rlens_merged;
+  std::map<int, std::array<size_t, 500>> rlens_merged;
 
   // each file 
   for (const std::string &fname : rlens_files) {
@@ -109,7 +109,7 @@ void merge_rlens(const std::vector<std::string> &rlens_files, const char* outnam
 
     while ((ret = bgzf_getline(fp, '\n', &line)) >= 0) {
       int taxid;
-      size_t tmp[200];
+      size_t tmp[500];
       char *ptr = line.s;
       char *endptr;
 
@@ -117,15 +117,15 @@ void merge_rlens(const std::vector<std::string> &rlens_files, const char* outnam
       taxid = strtol(ptr, &endptr, 10); 
       ptr = endptr;
 
-      // each new cell (200)
-      for (int j = 0; j < 200; j++) {
+      // each new cell (500)
+      for (int j = 0; j < 500; j++) {
         tmp[j] = strtoull(ptr, &endptr, 10);
         ptr = endptr;
       }
 
       // sum and merge
       auto &entry = rlens_merged[taxid];
-      for (int j = 0; j < 200; j++) {
+      for (int j = 0; j < 500; j++) {
         entry[j] += tmp[j];
       }
     }
@@ -143,15 +143,15 @@ void merge_rlens(const std::vector<std::string> &rlens_files, const char* outnam
 
   kstring_t kstr = {0, 0, NULL};
   ksprintf(&kstr, "id");
-  for (int i = 0; i < 200; i++)
+  for (int i = 0; i < 500; i++)
     ksprintf(&kstr, "\trlen%d", i);
   ksprintf(&kstr, "\n");
 
   for (const auto &it : rlens_merged) {
     ksprintf(&kstr, "%d", it.first);
-    for (int i = 0; i < 199; i++)
+    for (int i = 0; i < 499; i++)
       ksprintf(&kstr, "\t%lu", it.second[i]);
-    ksprintf(&kstr, "\t%lu\n", it.second[199]);
+    ksprintf(&kstr, "\t%lu\n", it.second[499]);
 
     if (kstr.l > 1000000) {
       assert(bgzf_write(fp, kstr.s, kstr.l) == (ssize_t)kstr.l);
diff --git a/metaDMG.cpp b/metaDMG.cpp
index 2eb0c6d..877dcba 100644
--- a/metaDMG.cpp
+++ b/metaDMG.cpp
@@ -128,6 +128,7 @@ std::map<int,mydataD> getval_full_norec(std::map<int, mydataD> &retmap, int2int
 
 
 mydata2 getval_stats(std::map<int, mydata2> &retmap, int2intvec &child, int taxid) {
+  fprintf(stderr,"PAS PAA SATAN\n");
     // fprintf(stderr,"getval\t%d\t%d\n",taxid,howmany);
     std::map<int, mydata2>::iterator it = retmap.find(taxid);
     if (it != retmap.end()) {
@@ -318,7 +319,7 @@ int main_getdamage(int argc, char **argv) {
 
     gzFile fpstat = NULL;
     assert((fpstat = gzopen(buf, "wb")) != NULL);
-    gzprintf(fpstat,"taxid\tnreads\tmea_len\tvar_len\tmean_gc\tvar_gc\tlca\trank\n");
+    gzprintf(fpstat,"taxid\tnreads\tmean_len\tvar_len\tmean_gc\tvar_gc\tlca\trank\n");
     for (std::map<int, std::vector<float> >::iterator it = gcconts.begin(); it != gcconts.end(); it++) {
         std::map<int, triple>::iterator it2 = dmg->assoc.find(it->first);
         assert(it2 != dmg->assoc.end());
@@ -1284,9 +1285,13 @@ 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());
+#if 1
+    if(child.size()>0)
+       getval_stats(stats, child, 1);  // this will do everything
+#endif
+    void aggr_stat3000(std::map<int, mydata2> &stats,int2int &parent);
+    if(0&&parent.size()>0)
+      aggr_stat3000(stats,parent);
     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;
diff --git a/misc/Makefile b/misc/Makefile
index db71412..5a6569d 100644
--- a/misc/Makefile
+++ b/misc/Makefile
@@ -1,71 +1,94 @@
-#assumed HTSSRC is a absolute path
+# Makefile for the misc directory
+# Builds utilities for metaDMG-cpp, using htslib and shared.o from parent
+
 CC  ?= gcc
 CXX ?= g++
 
-LIBS = -lz -lm -lbz2 -llzma -lpthread -lcurl
-
-CRYPTO_TRY=$(shell echo 'int main(){}'|g++ -x c++ - -lcrypto 2>/dev/null; echo $$?)
-ifeq "$(CRYPTO_TRY)" "0"
-$(info Crypto library is available to link; adding -lcrypto to LIBS)
-LIBS += -lcrypto
+FLAGS     := -O2
+CFLAGS    := $(FLAGS)
+CXXFLAGS  := $(FLAGS)
+CPPFLAGS  := $(CPPFLAGS) -Wall -Wextra
+LDFLAGS   :=
+LIBS      := -lz -lm -lbz2 -llzma -lpthread -lcurl
+
+# --- Crypto library detection ---
+HAVE_CRYPTO := $(shell echo 'int main(){}'|$(CXX) -x c++ - -lcrypto -o /dev/null 2>/dev/null && echo 0 || echo 1)
+ifeq ($(HAVE_CRYPTO),0)
+  $(info Crypto library is available to link; adding -lcrypto)
+  LIBS += -lcrypto
 else
-$(info Crypto library is not available to link; will not use -lcrypto)
+  $(info Crypto library is not available to link; skipping -lcrypto)
 endif
 
-
-ifdef HTSSRC
-$(info misc-HTSSRC defined $(HTSSRC))
-CPPFLAGS += -I"$(HTSSRC)"
-LIBS := "$(HTSSRC)/libhts.a" $(LIBS)
+ifeq ($(HTSSRC),systemwide)
+  $(info Using systemwide htslib installation)
+  LIBS += -lhts
+  LIBHTS :=
 else
-$(info HTSSRC not defined, assuming systemwide installation)
-LIBS += -lhts
+  # Convert relative HTSSRC to absolute if needed
+  ifneq ($(filter /%,$(HTSSRC)),$(HTSSRC))
+    override HTSSRC := $(shell [ -d "$(HTSSRC)" ] && cd "$(HTSSRC)" && pwd || echo "__INVALID__")
+    ifeq ($(HTSSRC),__INVALID__)
+      $(error HTSSRC is not a valid directory in misc folder)
+    endif
+    $(info Converted HTSSRC to absolute path: $(HTSSRC))
+  endif
+  CPPFLAGS += -I$(HTSSRC)
+  LIBHTS := $(HTSSRC)/libhts.a
+  LIBS := $(LIBHTS) $(LIBS)
 endif
 
-FLAGS = $(CPPFLAGS) -O3 $(LDFLAGS)
-
-.PHONY: all clean install
+PROGRAMS := compressbam extract_reads
+PROGRAMS_MISC := MAP dfit_optim
+OBJ := kmath.o Likelihood.o M3Load.o
 
-prefix      = /usr/local
-exec_prefix = $(prefix)
-bindir      = $(exec_prefix)/bin
+all: $(PROGRAMS) $(PROGRAMS_MISC)
 
-INSTALL = install
-INSTALL_DIR = $(INSTALL) -dm0755
-INSTALL_PROGRAM = $(INSTALL) -Dm0755
+# Ensure htslib and shared.o are available
+$(PROGRAMS) $(PROGRAMS_MISC): $(LIBHTS)
+$(OBJ): $(LIBHTS)
 
-PROGRAMS = compressbam extract_reads
-PROGRAMS_MISC = MAP dfit_optim
+# Check if libhts.a exists
+$(LIBHTS):
+	@if [ ! -f "$(LIBHTS)" ]; then \
+		echo "Error: $(LIBHTS) not found; ensure parent Makefile has built htslib"; \
+		exit 1; \
+	fi
 
-all: $(PROGRAMS) $(PROGRAMS_MISC)
+# Check if shared.o exists for extract_reads
+../shared.o:
+	@if [ ! -f "../shared.o" ]; then \
+		echo "Error: ../shared.o not found; ensure parent Makefile has built it"; \
+		exit 1; \
+	fi
 
 kmath.o: bfgs.cpp bfgs.h
-	$(CC) $(FLAGS) -c bfgs.cpp
+	$(CXX) -c $(CPPFLAGS) $(CXXFLAGS) bfgs.cpp -o kmath.o
+	$(CXX) -MM $(CPPFLAGS) $(CXXFLAGS) bfgs.cpp > kmath.d
 
 Likelihood.o: Likelihood.cpp
-	$(CC) $(FLAGS) -c Likelihood.cpp
+	$(CXX) -c $(CPPFLAGS) $(CXXFLAGS) Likelihood.cpp -o Likelihood.o
+	$(CXX) -MM $(CPPFLAGS) $(CXXFLAGS) Likelihood.cpp > Likelihood.d
 
 M3Load.o: M3Load.cpp
-	$(CC) $(FLAGS) -c M3Load.cpp
+	$(CXX) -c $(CPPFLAGS) $(CXXFLAGS) M3Load.cpp -o M3Load.o
+	$(CXX) -MM $(CPPFLAGS) $(CXXFLAGS) M3Load.cpp > M3Load.d
 
 compressbam: compressbam.cpp
-	$(CXX) $(FLAGS) compressbam.cpp -o compressbam $(LIBS)
+	$(CXX) $(CXXFLAGS) $(CPPFLAGS) compressbam.cpp -o compressbam $(LIBS) $(LDFLAGS)
 
 extract_reads: extract_reads.cpp ../shared.o
-	$(CXX) $(FLAGS) extract_reads.cpp -o extract_reads ../shared.o $(LIBS)
+	$(CXX) $(CXXFLAGS) $(CPPFLAGS) extract_reads.cpp ../shared.o -o extract_reads $(LIBS) $(LDFLAGS)
 
-MAP: MAP.cpp M3Load.o bfgs.o Likelihood.o
-	$(CXX) $(FLAGS) -std=c++11 MAP.cpp M3Load.o bfgs.o Likelihood.o -o MAP $(LIBS) -D__WITH_MAIN__
+MAP: MAP.cpp M3Load.o kmath.o Likelihood.o
+	$(CXX) $(CXXFLAGS) $(CPPFLAGS) -std=c++11 MAP.cpp M3Load.o kmath.o Likelihood.o -o MAP $(LIBS) $(LDFLAGS) -D__WITH_MAIN__
 
-dfit_optim: dfit_optim.cpp bfgs.o
-	$(CXX) $(FLAGS) dfit_optim.cpp bfgs.o -o dfit_optim $(LIBS) -D__WITH_MAIN__
-clean:
-	rm -f *.o $(PROGRAMS) $(PROGRAMS_MISC) *~
+dfit_optim: dfit_optim.cpp kmath.o
+	$(CXX) $(CXXFLAGS) $(CPPFLAGS) dfit_optim.cpp kmath.o -o dfit_optim $(LIBS) $(LDFLAGS) -D__WITH_MAIN__
 
-install: $(PROGRAMS)
-	$(INSTALL_DIR) $(DESTDIR)$(bindir)
-	$(INSTALL_PROGRAM) $(PROGRAMS) $(DESTDIR)$(bindir)
+# --- Automatisk afhængighedshåndtering ---
+-include $(OBJ:.o=.d)
 
-install-misc: $(PROGRAMS_MISC)
-	$(INSTALL_DIR) $(DESTDIR)$(bindir)
-	$(INSTALL_PROGRAM) $(PROGRAMS_MISC) $(DESTDIR)$(bindir)
+.PHONY: clean
+clean:
+	rm -f *.o *.d $(PROGRAMS) $(PROGRAMS_MISC) *~
diff --git a/mrand.h b/mrand.h
index ec02f5a..4030997 100644
--- a/mrand.h
+++ b/mrand.h
@@ -25,7 +25,7 @@ struct Ran {
 	inline Uint int32() { return (Uint)int64(); }
 };
 */
-typedef struct{
+typedef struct mrand{
   int type;
   #if defined(__linux__) || defined(__unix__)
     struct drand48_data buf0;
diff --git a/ngsLCA.cpp b/ngsLCA.cpp
index 9be49c9..a644516 100644
--- a/ngsLCA.cpp
+++ b/ngsLCA.cpp
@@ -412,7 +412,7 @@ std::vector<int> purge(std::vector<int> &taxids, std::vector<int> &editdist) {
     return tmpnewvec;
 }
 
-void hts(gzFile fp, samFile *fp_in, int2int &i2i, int2int &parent, bam_hdr_t *hdr, int2char &rank, int2char &name_map, int minmapq, int discard, int editMin, int editMax, double scoreLow, double scoreHigh, int minlength, int lca_rank, char *prefix, int howmany, samFile *fp_usedreads, int skipnorank, int2int &rank2level, int nthreads, int weighttype,long maxreads) {
+void hts(gzFile fp, samFile *fp_in, int2int &i2i, int2int &parent, bam_hdr_t *hdr, int2char &rank, int2char &name_map, int minmapq, int discard, int editMin, int editMax, double scoreLow, double scoreHigh, int minlength, int lca_rank, char *prefix, int howmany, samFile *fp_usedreads, int skipnorank, int2int &rank2level, int nthreads, int weighttype,long maxreads,samFile *fp_famout) {
   fprintf(stderr, "[%s] \t-> editMin:%d editmMax:%d scoreLow:%f scoreHigh:%f minlength:%d discard: %d prefix: %s howmany: %d skipnorank: %d weighttype: %d maxreads: %ld\n", __FUNCTION__, editMin, editMax, scoreLow, scoreHigh, minlength, discard, prefix, howmany, skipnorank, weighttype,maxreads);
     assert(fp_in != NULL);
     damage *dmg = new damage(howmany, nthreads, 13);
@@ -496,7 +496,14 @@ void hts(gzFile fp, samFile *fp_in, int2int &i2i, int2int &parent, bam_hdr_t *hd
                     // lca_rank is integer and is different from minus one
                     int2int::iterator myit = rank2level.find(lca);
                     assert(myit != rank2level.end());
-		    
+		    //write to family out
+		    if (myit->second != -1 && (myit->second <= 16)) {
+		      //family is 16, see rank2level.txt
+		       for (int i = 0; i < myq->l; i++)
+			 if(fp_famout)
+			   assert(sam_write1(fp_famout, hdr, myq->ary[i]) >= 0);
+		    }
+		    //standard analyses
 		    if (myit->second != -1 && (myit->second <= lca_rank)) {
                         adder(lca, strlen(seq), gccontent(seq));
 			            //fprintf(stderr,"Looping through alignments we have :%d \n",myq->l);
@@ -614,6 +621,13 @@ void hts(gzFile fp, samFile *fp_in, int2int &i2i, int2int &parent, bam_hdr_t *hd
             // lca_rank is integer and is different from minus one
             int2int::iterator myit = rank2level.find(lca);
             assert(myit != rank2level.end());
+	    //write to family out
+	    if (myit->second != -1 && (myit->second <= 16)) {
+	      //family is 16, see rank2level.txt
+	      for (int i = 0; i < myq->l; i++)
+		if(fp_famout)
+		  assert(sam_write1(fp_famout, hdr, myq->ary[i]) >= 0);
+	    }
             if (myit->second != -1 && (myit->second <= lca_rank)) {
                 adder(lca, strlen(seq), gccontent(seq));
                 //      if(correct_rank(lca_rank,lca,rank,norank2species)){
@@ -718,7 +732,7 @@ int2node makeNodes(int2int &parent){
 int main_lca(int argc, char **argv) {
     htsFormat *dingding2 = (htsFormat *)calloc(1, sizeof(htsFormat));
     if (argc == 1) {
-        fprintf(stderr, "\t-> ./metaDMG-cpp lca --names --nodes --acc2tax [-edit_dist_[min/max] --sim_score_[low/high] --min_mapq --bam --lca_rank  --filtered_acc2tax --used_reads [0,1] --weight_type [0,1] #%s \n", METADAMAGE_VERSION);
+        fprintf(stderr, "\t-> ./metaDMG-cpp lca --names --nodes --acc2tax [-edit_dist_[min/max] --sim_score_[low/high] --min_mapq --bam --lca_rank  --filtered_acc2tax --used_reads [0,1] --weight_type [0,1] --famout [0,1] #%s \n", METADAMAGE_VERSION);
         return 0;
     }
     catchkill();
@@ -738,9 +752,7 @@ int main_lca(int argc, char **argv) {
     int2int *i2i = NULL;
     //  fprintf(stderr,"p->header: %p\n",p->header);
     if (p->header)
-      i2i = bamRefId2tax(p->header, p->acc2taxfile, p->htsfile, errmap, p->tempfolder, p->reallyDump, p->filteredAcc2taxfile,NULL);
-
-
+      i2i = bamRefId2tax(p->header, p->acc2taxfile, p->htsfile, errmap, p->tempfolder, p->useDump, p->filteredAcc2taxfile,NULL);
     // map of taxid -> taxid
     int2int parent;
     // map of taxid -> rank
@@ -771,7 +783,7 @@ int main_lca(int argc, char **argv) {
     }
     samFile *usedreads_sam = NULL;
     if (p->usedreads_sam != NULL) {  // p->usedreads sam is const *, sorry this is confusing
-        char out_mode[5] = "ws";
+        char out_mode[5] = "wb";
         if ((usedreads_sam = sam_open_format(p->usedreads_sam, out_mode, dingding2)) == 0) {
             fprintf(stderr, "Error opening file for writing\n");
             exit(0);
@@ -779,8 +791,18 @@ int main_lca(int argc, char **argv) {
         if (sam_hdr_write(usedreads_sam, p->header) < 0)
             fprintf(stderr, "writing headers to %s", p->usedreads_sam);
     }
+    samFile *famout_sam = NULL;
+    if (p->famout_sam != NULL) {  // p->usedreads sam is const *, sorry this is confusing
+      char out_mode[5] = "wb";
+      if ((famout_sam = sam_open_format(p->famout_sam, out_mode, dingding2)) == 0) {
+	fprintf(stderr, "Error opening file for writing\n");
+	exit(0);
+      }
+      if (sam_hdr_write(famout_sam, p->header) < 0)
+	fprintf(stderr, "writing headers to %s", p->famout_sam);
+    }
     gzprintf(p->fp1,"queryid\tseq\tlen\tnaln\tgc\tlca\ttaxa_path\n");
-    hts(p->fp1, p->hts, *i2i, parent, p->header, rank, name_map, p->minmapq, p->discard, p->editdistMin, p->editdistMax, p->simscoreLow, p->simscoreHigh, p->minlength, lca_rank, p->outnames, p->howmany, usedreads_sam, p->skipnorank, tax2level, p->nthreads, p->weighttype,p->maxreads);
+    hts(p->fp1, p->hts, *i2i, parent, p->header, rank, name_map, p->minmapq, p->discard, p->editdistMin, p->editdistMax, p->simscoreLow, p->simscoreHigh, p->minlength, lca_rank, p->outnames, p->howmany, usedreads_sam, p->skipnorank, tax2level, p->nthreads, p->weighttype,p->maxreads,famout_sam);
 
     fprintf(stderr, "\t-> Number of species with reads that map uniquely: %lu\n", specWeight.size());
 
@@ -798,8 +820,10 @@ int main_lca(int argc, char **argv) {
 
     if (usedreads_sam != NULL)
         sam_close(usedreads_sam);
+    if (famout_sam != NULL)
+      sam_close(famout_sam);
     if (p->fp_lcadist) {
-        gzprintf(p->fp_lcadist,"taxid\tnreads\tmea_len\tvar_len\tmean_gc\tvar_gc\tlca\trank\n");
+        gzprintf(p->fp_lcadist,"taxid\tnreads\tmean_len\tvar_len\tmean_gc\tvar_gc\tlca\trank\n");
         for (std::map<int, lcatriplet>::iterator it = lcastat.begin(); it != lcastat.end(); it++) {
             lcatriplet tmp = it->second;
             gzprintf(p->fp_lcadist, "%d\t%d\t%f\t%f\t%f\t%f", it->first, tmp.nalignments, mean(tmp.readlengths), var(tmp.readlengths), mean(tmp.gccontents), var(tmp.gccontents));
diff --git a/ngsLCA_cli.cpp b/ngsLCA_cli.cpp
index 2070f12..7b6c622 100644
--- a/ngsLCA_cli.cpp
+++ b/ngsLCA_cli.cpp
@@ -38,12 +38,13 @@ pars *pars_init() {
     p->skipnorank = 1;
     p->howmany = 5;
     p->usedreads_sam = NULL;
+    p->famout_sam = NULL;
     p->fixdb = 1;
     p->nthreads = 4;
     p->weighttype = 0;
     p->tempfolder = strdup("");
     p->ignore_errors = 0;
-    p->reallyDump = 0;
+    p->useDump = 0;
     p->maxreads = -1;
     return p;
 }
@@ -188,7 +189,8 @@ pars *get_pars(int argc, char **argv) {
         return NULL;
     }
 
-    int make_used_reads = 0;
+    int make_used_reads = 1;
+    int make_famout_reads = 1;
 
     while (*argv) {
         char *key = *argv;
@@ -236,6 +238,8 @@ pars *get_pars(int argc, char **argv) {
             p->maxreads = atol(val);
         else if (!strcasecmp("--used_reads", key))
             make_used_reads = atoi(val);
+	else if (!strcasecmp("--famout", key))
+            make_famout_reads = atoi(val);
         else if (!strcasecmp("--no_rank2species", key))
             p->norank2species = atoi(val);
         else if (!strcasecmp("--skip_no_rank", key))
@@ -244,8 +248,8 @@ pars *get_pars(int argc, char **argv) {
             p->nthreads = atoi(val);
         else if (!strcasecmp("--weight_type", key))
             p->weighttype = atoi(val);
-	else if (!strcasecmp("--reallyDump", key))
-	   p->reallyDump = atoi(val);
+	else if (!strcasecmp("--useDump", key))
+	   p->useDump = atoi(val);
 	else if (!strcasecmp("--ignore_errors", key)||!strcasecmp("-i", key))
 	  p->ignore_errors++;
         else if (!strcasecmp("--temp", key)) {
@@ -286,10 +290,15 @@ pars *get_pars(int argc, char **argv) {
     assert(p->fp3);
 #endif
     if (make_used_reads) {
-        snprintf(buf, 1024, "%s.usedreads.sam", p->outnames);
+        snprintf(buf, 1024, "%s.usedreads.bam", p->outnames);
         fprintf(stderr, "\t-> Will output the reads that are used for damage file:\t\'%s\'\n", buf);
         p->usedreads_sam = strdup(buf);
     }
+    if (make_used_reads) {
+      snprintf(buf, 1024, "%s.famoutreads.bam", p->outnames);
+      fprintf(stderr, "\t-> Will output the reads that has lca below family :\t\'%s\'\n", buf);
+      p->famout_sam = strdup(buf);
+    }
 
     pthread_mutex_lock(&mutex1);
     pthread_mutex_lock(&mutex2);
@@ -317,7 +326,6 @@ void print_pars(FILE *fp, pars *p) {
     fprintf(fp, "\t-> --ignore_errors\t%d\n", p->ignore_errors);
     fprintf(fp, "\t-> --temp\t%s\n", p->tempfolder);
     fprintf(fp, "\t-> --filtered_acc2tax\t%s\n", p->filteredAcc2taxfile);
-
 }
 
 #ifdef __WITH_MAIN__
diff --git a/ngsLCA_cli.h b/ngsLCA_cli.h
index e061f82..6801a63 100644
--- a/ngsLCA_cli.h
+++ b/ngsLCA_cli.h
@@ -35,12 +35,13 @@ typedef struct {
     int norank2species;
     int howmany;
     char *usedreads_sam;
+  char *famout_sam;
     int fixdb;  // used for disabling mod_db function
     int nthreads;
     int weighttype;
     char *tempfolder;
     int ignore_errors;
-  int reallyDump;
+  int useDump;
   long maxreads;
 } pars;
 
diff --git a/profile.cpp b/profile.cpp
index 4673fc4..4408ca5 100644
--- a/profile.cpp
+++ b/profile.cpp
@@ -287,8 +287,8 @@ inline void increaseCounters(const bam1_t *b, const char *reconstructedReference
 int damage::damage_analysis(bam1_t *b, int which, float incval) {
   //  fprintf(stderr,"\t-> incval: %f\n",incval);
   if (assoc.find(which) == assoc.end()) {
-    triple val = {0, getmatrix(MAXLENGTH, 16), getmatrix(MAXLENGTH, 16),new size_t[200]};
-    for(int i=0;i<200;i++)
+    triple val = {0, getmatrix(MAXLENGTH, 16), getmatrix(MAXLENGTH, 16),new size_t[500]};
+    for(int i=0;i<500;i++)
       val.rlens[i] = 0;
     assoc[which] = val;
     mm5pF = val.mm5pF;
@@ -298,7 +298,7 @@ int damage::damage_analysis(bam1_t *b, int which, float incval) {
     std::map<int, triple>::iterator it = assoc.find(which);
     it->second.nreads++;
     //    fprintf(stderr,"[%s] it->first:%d it->second.nreads:%d\n",__FUNCTION__,it->first,it->second.nreads);
-    assert(b->core.l_qseq<200);
+    assert(b->core.l_qseq<500);
     it->second.rlens[b->core.l_qseq] =  it->second.rlens[b->core.l_qseq] + 1; 
     if (b->core.l_qseq - 10 > temp_len) {
         temp_len = b->core.l_qseq;
@@ -374,30 +374,30 @@ void damage::bwrite(char *fname) {
     snprintf(onam, 1024, "%s.rlens.gz", fname);
     fprintf(stderr, "\t-> Will dump: \'%s\' this contains read length distributions for: %lu items\n", onam, assoc.size());
     fp = my_bgzf_open(onam, nthreads);
-    kstring_t kstr2000;
-    kstr2000.s = NULL;
-    kstr2000.l = kstr2000.m = 0;
-    ksprintf(&kstr2000,"id");
-    for(int i=0;i<200;i++)
-      ksprintf(&kstr2000,"\trlen%d",i);
-    ksprintf(&kstr2000,"\n");
+    kstring_t kstr3000;
+    kstr3000.s = NULL;
+    kstr3000.l = kstr3000.m = 0;
+    ksprintf(&kstr3000,"id");
+    for(int i=0;i<500;i++)
+      ksprintf(&kstr3000,"\trlen%d",i);
+    ksprintf(&kstr3000,"\n");
     for (std::map<int, triple>::iterator it = assoc.begin(); it != assoc.end(); it++) {
         if (it->second.nreads == 0)  // should never happen
             continue;
-	ksprintf(&kstr2000,"%d",it->first);
-	for(int i=0;i<200-1;i++)
-	  ksprintf(&kstr2000,"\t%lu",it->second.rlens[i]);
-	ksprintf(&kstr2000,"\t%lu\n",it->second.rlens[199]);
-	if(kstr2000.l>1000000){
-	  assert(bgzf_write(fp,kstr2000.s,kstr2000.l)==kstr2000.l);
-	  kstr2000.l  = 0;
+	ksprintf(&kstr3000,"%d",it->first);
+	for(int i=0;i<500-1;i++)
+	  ksprintf(&kstr3000,"\t%lu",it->second.rlens[i]);
+	ksprintf(&kstr3000,"\t%lu\n",it->second.rlens[499]);
+	if(kstr3000.l>1000000){
+	  assert(bgzf_write(fp,kstr3000.s,kstr3000.l)==kstr3000.l);
+	  kstr3000.l  = 0;
 	}
 	  
     }
-    assert(bgzf_write(fp,kstr2000.s,kstr2000.l)==kstr2000.l);
-    if(kstr2000.l>0)
-      free(kstr2000.s);
-    kstr2000.l  = 0;
+    assert(bgzf_write(fp,kstr3000.s,kstr3000.l)==kstr3000.l);
+    if(kstr3000.l>0)
+      free(kstr3000.s);
+    kstr3000.l  = 0;
     bgzf_close(fp);
 }
 
diff --git a/shared.cpp b/shared.cpp
index 98b4504..e8daf5b 100644
--- a/shared.cpp
+++ b/shared.cpp
@@ -15,6 +15,14 @@
 #include <map>      // for operator!=, __map_iterator, operator==
 #include <vector>   // for vector
 
+
+int file_is_older(const char *file1, const char *file2) {
+    struct stat st1, st2;
+    if (stat(file1, &st1) != 0) return -1;
+    if (stat(file2, &st2) != 0) return -1;
+    return st1.st_mtime < st2.st_mtime; // return 1 hvis file1 er ældre
+}
+
 BGZF *getbgzf(const char *str1, const char *mode, int nthreads) {
     BGZF *fp = NULL;
     fp = bgzf_open(str1, mode);
@@ -44,6 +52,15 @@ BGZF *getbgzf3(const char *str1, const char *str2, const char *str3, const char
     return getbgzf(tmp, mode, nthreads);
 }
 
+char *getfilename4(const char *str1, const char *str2, const char *str3, const char *str4) {
+    unsigned tmp_l = strlen(str1) + strlen(str2) + strlen(str3) + strlen(str4) + 5;
+    char *tmp=(char*)calloc(tmp_l,1);
+    snprintf(tmp, tmp_l, "%s%s%s%s", str1, str2, str3, str4);
+    return tmp;//getbgzf(tmp, mode, nthreads);
+}
+
+
+
 BGZF *getbgzf4(const char *str1, const char *str2, const char *str3, const char *str4, const char *mode, int nthreads) {
     unsigned tmp_l = strlen(str1) + strlen(str2) + strlen(str3) + strlen(str4) + 5;
     char tmp[tmp_l];
@@ -229,102 +246,127 @@ void parse_nodes2(int2int &parent, int2intvec &child) {
 }
 
 int SIG_COND = 1;
-int2int *bamRefId2tax(bam_hdr_t *hdr, char *acc2taxfile, char *bamfile, int2int &errmap, char *tempfolder, int forceDump, char *filteredAcc2taxfile,char2int *acc2taxidmap) {
-    fprintf(stderr, "\t-> Starting to extract (acc->taxid) from binary file: \'%s\'\n", acc2taxfile);
+int2int *bamRefId2tax(bam_hdr_t *hdr, char *acc2taxfile, char *bamfile, int2int &errmap,
+                      char *tempfolder, int usedump, char *filteredAcc2taxfile,
+                      char2int *acc2taxidmap) {
+    fprintf(stderr, "\t-> Starting to extract (acc->taxid) from binary file: '%s'\n", acc2taxfile);
     fflush(stderr);
-    int dodump = !fexists4(tempfolder, basename(acc2taxfile), basename(bamfile), ".bin");
-    if(acc2taxidmap!=NULL)
-      forceDump = 1;
-    dodump += forceDump;
-    
-    fprintf(stderr, "\t-> Checking if need to reload acc2tax and dump. dodump=%d forcedump=%d acc2taxidmap: %p\n", dodump,forceDump,acc2taxidmap);
+
+    int binExists = fexists4(tempfolder, basename(acc2taxfile), basename(bamfile), ".bin");
+    int useBinary = 0;
+
+    if (acc2taxidmap != NULL) {
+        useBinary = 0;
+    } else if (binExists && usedump) {
+        useBinary = 1;
+    } else if (!binExists && usedump) {
+        useBinary = 0;
+    } else if (binExists && !usedump) {
+        useBinary = 1;
+    } else {
+        useBinary = 0;
+    }
+
+    fprintf(stderr, "\t-> Binary exists: %d, usedump: %d, acc2taxidmap: %p => useBinary: %d\n",
+            binExists, usedump, acc2taxidmap, useBinary);
 
     time_t t = time(NULL);
     BGZF *fp = NULL;
-    if (dodump)
-      fp = getbgzf4(tempfolder, basename(acc2taxfile), basename(bamfile), ".bin", "wb", 4);
-    else
+    if (!useBinary && usedump)
+        fp = getbgzf4(tempfolder, basename(acc2taxfile), basename(bamfile), ".bin", "wb", 4);
+    else if (useBinary){
+      char *binfile = getfilename4(tempfolder, basename(acc2taxfile), basename(bamfile), ".bin"); //not cleaned up
+      
+      int binIsOlder = 0;
+      if (binExists && usedump) {
+	binIsOlder = file_is_older(binfile, bamfile);
+	if (binIsOlder == 1) {
+	  fprintf(stderr, "\n\t-> WARNING: Binary cache file '%s' is older than BAM file '%s'.\n", binfile, bamfile);
+	  fprintf(stderr, "\t-> This may indicate out-of-date accession-taxid mapping. Refusing to use it.\n");
+	  exit(1);
+	}
+      }
       fp = getbgzf4(tempfolder, basename(acc2taxfile), basename(bamfile), ".bin", "rb", 4);
-
-    // This contains refname(as int) -> taxid
+    }
     int2int *am = new int2int;
 
-    // Open the filtered output file if specified
     gzFile filteredFile = Z_NULL;
-    if (filteredAcc2taxfile != NULL && dodump) {
+    if (filteredAcc2taxfile != NULL && !useBinary) {
         filteredFile = gzopen(filteredAcc2taxfile, "wb");
         if (!filteredFile) {
             fprintf(stderr, "Error opening file '%s' for writing\n", filteredAcc2taxfile);
             exit(1);
         }
-        gzprintf(filteredFile, "BAM_Reference\tTaxID\tBam_hdr_index\n"); // Write header
+        gzprintf(filteredFile, "BAM_Reference\tTaxID\tBam_hdr_index\n");
     }
 
-    if (dodump) {
-        char buf[4096];
+    if (!useBinary) {
         int at = 0;
-
         kstring_t *kstr = (kstring_t *)malloc(sizeof(kstring_t));
         kstr->l = kstr->m = 0;
         kstr->s = NULL;
         BGZF *fp2 = getbgzf(acc2taxfile, "rb", 2);
-        bgzf_getline(fp2, '\n', kstr);  // Skip header
+        bgzf_getline(fp2, '\n', kstr);  // skip header
         kstr->l = 0;
-	int nprocs =0;
+        int nprocs = 0;
         while (SIG_COND && bgzf_getline(fp2, '\n', kstr)) {
-            if (kstr->l == 0)
-                break;
+            if (kstr->l == 0) break;
             if (!((at++ % 100000)) && isatty(fileno(stderr)))
-                fprintf(stderr, "\r\t-> At linenr: %d in \'%s\'      ", at, acc2taxfile);
+                fprintf(stderr, "\r\t-> At linenr: %d in '%s'      ", at, acc2taxfile);
+
             char *tok = strtok(kstr->s, "\t\n ");
             char *key = strtok(NULL, "\t\n ");
             tok = strtok(NULL, "\t\n ");
             int val = atoi(tok);
-	    if(acc2taxidmap!=NULL){
-	      if(!acc2taxidmap->insert(std::pair<char*,int>(strdup(key),val)).second){
-		fprintf(stderr,"\t-> Problem inserting key: %s with value: %d\n",key,val);
-		exit(0);
-	      } 
-	    }
-            int valinbam = bam_name2id(hdr, key);
+            if (!key) continue;
 
-            if (valinbam == -1)
-                continue;
-	    nprocs++;
-            if (fp != NULL) {
-                assert(bgzf_write(fp, &valinbam, sizeof(int)) == sizeof(int));
-                assert(bgzf_write(fp, &val, sizeof(int)) == sizeof(int));
+            if (acc2taxidmap != NULL) {
+                if (!acc2taxidmap->insert(std::pair<char*, int>(strdup(key), val)).second) {
+                    fprintf(stderr, "\t-> Problem inserting key: %s with value: %d\n", key, val);
+                    exit(1);
+                }
             }
 
+            int valinbam = bam_name2id(hdr, key);
+            if (valinbam == -1) continue;
+            nprocs++;
+
+            if (fp != NULL)
+                assert(bgzf_write(fp, &valinbam, sizeof(int)) == sizeof(int) &&
+                       bgzf_write(fp, &val, sizeof(int)) == sizeof(int));
+
             if (am->find(valinbam) != am->end())
-                fprintf(stderr, "\t-> Duplicate entries found \'%s\'\n", key);
+                fprintf(stderr, "\t-> Duplicate entries found '%s'\n", key);
             (*am)[valinbam] = val;
 
-            // Write to the filtered file if specified
             if (filteredFile != Z_NULL)
-	      gzprintf(filteredFile, "%s\t%d\t%d\n", key, val,valinbam);
+                gzprintf(filteredFile, "%s\t%d\t%d\n", key, val, valinbam);
             kstr->l = 0;
         }
         bgzf_close(fp2);
-	fprintf(stderr,"\t-> Number of items from acc2tax file that is relevant from the bamheader: %d\n",nprocs);
-	fflush(stderr);
+        free(kstr->s);
+        free(kstr);
+
+        fprintf(stderr, "\t-> Number of items from acc2tax file that is relevant from the bamheader: %d\n", nprocs);
+        fflush(stderr);
     } else {
         int valinbam, val;
         while (bgzf_read(fp, &valinbam, sizeof(int))) {
             assert(bgzf_read(fp, &val, sizeof(int)) == sizeof(int));
             (*am)[valinbam] = val;
-
         }
     }
 
-    // Close the filtered file if it was opened
     if (filteredFile != Z_NULL) {
         gzclose(filteredFile);
         fprintf(stderr, "\n\t-> Filtered acc2taxfile saved to '%s'\n", filteredAcc2taxfile);
     }
 
-    bgzf_close(fp);
-    fprintf(stderr, "\t-> Number of entries to use from accession to taxid: %lu, time taken: %.2f sec acc2taxidmap.size(): %lu\n", am->size(), (float)(time(NULL) - t),acc2taxidmap!=NULL?acc2taxidmap->size():0);
+    if (fp != NULL)
+        bgzf_close(fp);
+
+    fprintf(stderr, "\t-> Number of entries to use from accession to taxid: %lu, time taken: %.2f sec acc2taxidmap.size(): %lu\n",
+            am->size(), (float)(time(NULL) - t), acc2taxidmap != NULL ? acc2taxidmap->size() : 0);
     return am;
 }
 
diff --git a/test/output.md5 b/test/output.md5
index 0db493b..2ebc225 100644
--- a/test/output.md5
+++ b/test/output.md5
@@ -1,28 +1,28 @@
 441b0ab9e87686b2aa77ffd2e63d5e85  output/test_aggregate.stat
 13ce6ac60c336f443583e6803b45c73b  output/test_dfit_global.dfit.fix
 #6c45bfb40a55b840dbe1dedf71d32074  output/test_dfit_local_10threads.dfit.fix
-d5ec49a8685897936d18ecff68586bec  output/test_dfit_local.dfit.fix
+#d5ec49a8685897936d18ecff68586bec  output/test_dfit_local.dfit.fix
 783c885b2ee8ac6469ab070b972e73d8  output/test_getdamage_global.bdamage
 3398e353af15d1646228cb8bbaf4a5a8  output/test_getdamage_global.bdamage.tsv
 f0882ca11587cea63aa120d1a64c8d14  output/test_getdamage_global.res
-96caf31ec344dbe1343ea11bb1a27898  output/test_getdamage_global.rlens
-7e2d05365143cd0fde70256fac00a010  output/test_getdamage_global.stat
+3fda91168f2a07d82312e4589841df1e  output/test_getdamage_global.stat
 71308047bbc8f703fdd65c798de91cfa  output/test_getdamage_local.bdamage
 93ec5114d069f984ccaaa5079c15db1e  output/test_getdamage_local.bdamage.tsv
 d95e54a3886b9f411bde8d96de4af32e  output/test_getdamage_local.res
-c9d41aa644ba3d89075888b534c1b6c8  output/test_getdamage_local.rlens
-91a833ba74ca25355f9613233b9d51af  output/test_getdamage_local.stat
+3401f8b78adb9ab36f609d8e9e96b84a  output/test_getdamage_local.stat
 443b35425bc868259fde18dbe86011e3  output/test_lca.bdamage
 #d4532cede0b424946f92559b2e737899  output/test_lca.lca##this has been commented out because eithput is random between reqad and the reverse complement
-de0751cd16aad49a49b0063226d3fa2a  output/test_lca.rlens
-caf50bb495f772b5b7588cc41a5ce275  output/test_lca.stat
+e6a34153420df3f15a7c00eb6bf082fc  output/test_lca.stat
 eaa8bb5f32aa0ef05bcf923d0c2bd711  output/test_lca_taxa.uglyprint.mismatch
 823ac0bcf53f3a69329ff7c3eaf6a25f  output/test_lca_taxa.uglyprint.stat
 70363c8585d69e8d4457fc1d1ba7cd92  output/test_lca.uglyprint.mismatch
 1f9c4d0ab804d6a6e077573e6c8e736a  output/test_lca.uglyprint.stat
 a23c34eaca4a1d1742dd65f98d70efbf  output/test_mergedamage.bdamage
-1e70c3fbbceb07f1046bf735e39e89a4  output/test_mergedamage.rlens
 dc5bf1edd27b494e23f72d41d718600d  data/f570b1db7c.dedup.filtered.bam
 #c73c56f8a3e2dc9b13c949ec1404bd50  data/f570b1db7c.dedup.filtered.rname.bam
 #3568c671c0fa09797284746bb0528ddd  data/f570b1db7c.dedup.filtered.rname.bam ##apparantly samtools 1.20 has changed the within qname order when sorting by readname. This is of cause fine but some seqs would be reverse complemented. When switching to the new order then this entry should be used instead of the above
-2ae4149efce77b91786d4a13b5451c3b  output/test_lca.lca.sub
\ No newline at end of file
+2ae4149efce77b91786d4a13b5451c3b  output/test_lca.lca.sub
+53943f1ceab5c86c6d815fc55d0363c7  output/test_getdamage_global.rlens
+57013f84420d6460ae5e15aa0092dc1e  output/test_getdamage_local.rlens
+abc7a9ec248d48204d5f417f4c30c59d  output/test_lca.rlens
+9574c8051682dc50c040e5e8c16704a6  output/test_mergedamage.rlens
diff --git a/test/output.md5.macos b/test/output.md5.macos
index d2c81d2..dd3ab04 100644
--- a/test/output.md5.macos
+++ b/test/output.md5.macos
@@ -1,25 +1,25 @@
 441b0ab9e87686b2aa77ffd2e63d5e85  output/test_aggregate.stat
 037258988f595d8737cb4647ddcb0d6d  output/test_dfit_global.dfit.fix
-d5ec49a8685897936d18ecff68586bec  output/test_dfit_local.dfit.fix
+#d5ec49a8685897936d18ecff68586bec  output/test_dfit_local.dfit.fix
 783c885b2ee8ac6469ab070b972e73d8  output/test_getdamage_global.bdamage
 3398e353af15d1646228cb8bbaf4a5a8  output/test_getdamage_global.bdamage.tsv
 f0882ca11587cea63aa120d1a64c8d14  output/test_getdamage_global.res
-96caf31ec344dbe1343ea11bb1a27898  output/test_getdamage_global.rlens
-7e2d05365143cd0fde70256fac00a010  output/test_getdamage_global.stat
+3fda91168f2a07d82312e4589841df1e  output/test_getdamage_global.stat
 71308047bbc8f703fdd65c798de91cfa  output/test_getdamage_local.bdamage
 93ec5114d069f984ccaaa5079c15db1e  output/test_getdamage_local.bdamage.tsv
 d95e54a3886b9f411bde8d96de4af32e  output/test_getdamage_local.res
-c9d41aa644ba3d89075888b534c1b6c8  output/test_getdamage_local.rlens
-91a833ba74ca25355f9613233b9d51af  output/test_getdamage_local.stat
+3401f8b78adb9ab36f609d8e9e96b84a  output/test_getdamage_local.stat
 443b35425bc868259fde18dbe86011e3  output/test_lca.bdamage
 #d4532cede0b424946f92559b2e737899  output/test_lca.lca##this has been commented out because eithput is random between reqad and the reverse complement
-de0751cd16aad49a49b0063226d3fa2a  output/test_lca.rlens
-caf50bb495f772b5b7588cc41a5ce275  output/test_lca.stat
+e6a34153420df3f15a7c00eb6bf082fc  output/test_lca.stat
 eaa8bb5f32aa0ef05bcf923d0c2bd711  output/test_lca_taxa.uglyprint.mismatch
 823ac0bcf53f3a69329ff7c3eaf6a25f  output/test_lca_taxa.uglyprint.stat
 70363c8585d69e8d4457fc1d1ba7cd92  output/test_lca.uglyprint.mismatch
 1f9c4d0ab804d6a6e077573e6c8e736a  output/test_lca.uglyprint.stat
 a23c34eaca4a1d1742dd65f98d70efbf  output/test_mergedamage.bdamage
-1e70c3fbbceb07f1046bf735e39e89a4  output/test_mergedamage.rlens
 dc5bf1edd27b494e23f72d41d718600d  data/f570b1db7c.dedup.filtered.bam
-2ae4149efce77b91786d4a13b5451c3b  output/test_lca.lca.sub
\ No newline at end of file
+2ae4149efce77b91786d4a13b5451c3b  output/test_lca.lca.sub
+53943f1ceab5c86c6d815fc55d0363c7  output/test_getdamage_global.rlens
+57013f84420d6460ae5e15aa0092dc1e  output/test_getdamage_local.rlens
+abc7a9ec248d48204d5f417f4c30c59d  output/test_lca.rlens
+9574c8051682dc50c040e5e8c16704a6  output/test_mergedamage.rlens
\ No newline at end of file
diff --git a/test/testAll.sh b/test/testAll.sh
index a54e85f..9b050ba 100755
--- a/test/testAll.sh
+++ b/test/testAll.sh
@@ -16,6 +16,14 @@ echo "Testing Existence of ${PRG}"
 if [[ ! -f "${PRG}" ]]; then
     echo "Problem finding program: ${PRG}"
     RVAL=1
+    exit 1;
+fi
+
+echo "Testing Existence of samtools"
+if ! command -v samtools &> /dev/null ; then
+    echo "Problem finding program: samtools"
+    RVAL=1
+    exit 1;
 fi
 
 echo "Testing Existence of ${BAM1}"
@@ -89,12 +97,17 @@ 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
-#gunzip -c output/test_dfit_local.dfit.gz|cut -f1-6,8-|./round_file.sh |sort -k1,1n > output/test_dfit_local.dfit.fix
-
+gunzip -c output/test_dfit_local.dfit.gz | \
+cut -f1-6,8- | \
+awk 'BEGIN{OFS=FS="\t"} 
+     NR==1 {print; next} 
+     {
+       for(i=2;i<=NF;i++) {
+         if($i ~ /^[0-9.eE+-]+$/) $i = sprintf("%.2f", $i)
+       }
+       print
+     }' \
+> output/test_dfit_local.dfit.fix
 
 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"