update, reenable old versions etc

DEploid-dev · Nov 16, 2018 · c30b6e6 · c30b6e6
1 parent 9e9de7f
commit c30b6e6
Show file tree

Hide file tree

Showing 5 changed files with 30 additions and 83 deletions.
diff --git a/.travis.yml b/.travis.yml
@@ -22,9 +22,9 @@ script:
   - make check
   - if [ $TRAVIS_OS_NAME == linux ]; then ./tests/test_binary.sh; fi
   - ./tests/testPOS.sh
-  #- ./tests/test_binaryVcfVsTxt.sh
-  #- ./tests/test-against-previous-version.sh
-  #- ./tests/test_binaryReproducible.sh
+  - ./tests/test_binaryVcfVsTxt.sh
+  - ./tests/test-against-previous-version.sh
+  - ./tests/test_binaryReproducible.sh
   - if [ $TRAVIS_OS_NAME == linux ]; then cd docs/doxygen; doxygen Doxyfile; cd ../..; fi
 
 after_success:

diff --git a/Makefile.am b/Makefile.am
@@ -22,6 +22,7 @@ common_LDADD = -lz
 common_src = src/random/fastfunc.cpp \
 			 src/random/random_generator.cpp \
 			 src/random/mersenne_twister.cpp \
+			 src/chooseK.cpp \
 			 src/ibd.cpp \
 			 src/mcmc.cpp \
 			 src/panel.cpp \

diff --git a/src/dEploid.cpp b/src/dEploid.cpp
@@ -26,7 +26,7 @@
 #include <iostream>  // std::cout
 #include "mcmc.hpp"
 #include "dEploidIO.hpp"
-
+#include "chooseK.hpp"
 
 int main(int argc, char *argv[]) {
     try {
@@ -86,7 +86,7 @@ int main(int argc, char *argv[]) {
                 delete lassoMcmcSample;
             }
             dEploidIO.writeHap(hap, "lasso");
-        } else if (dEploidIO.useIBD()) {  // ibd
+        } else if (dEploidIO.useBestPractice()) {  // best practice
             MersenneTwister rg(dEploidIO.randomSeed());
             if (dEploidIO.usePanel()) {
                 // #############################################################
@@ -229,15 +229,30 @@ int main(int argc, char *argv[]) {
 
             dEploidIO.paintIBD();
             dEploidIO.writeHap(hap, "final");
-        } else {
+        } else { // classic version, with IBD flag
+            if (dEploidIO.useIBD()) {  // ibd
+                McmcSample * ibdMcmcSample = new McmcSample();
+                MersenneTwister ibdRg(dEploidIO.randomSeed());
+                McmcMachinery ibdMcmcMachinery(&dEploidIO.plaf_,
+                                               &dEploidIO.refCount_,
+                                               &dEploidIO.altCount_,
+                                               dEploidIO.panel,
+                                               &dEploidIO,
+                                               "DEploid-IBD",
+                                               "DEploid-IBD",
+                                               ibdMcmcSample,
+                                               &ibdRg,
+                                               true);
+            }
+
             McmcSample * mcmcSample = new McmcSample();
             MersenneTwister rg(dEploidIO.randomSeed());
             McmcMachinery mcmcMachinery(&dEploidIO.plaf_,
                                 &dEploidIO.refCount_,
                                 &dEploidIO.altCount_,
                                 dEploidIO.panel,
                                 &dEploidIO,
-                                "DEploid first version",
+                                "DEploid classic version",
                                 "",
                                 mcmcSample,
                                 &rg,

diff --git a/src/dEploidIO.cpp b/src/dEploidIO.cpp
@@ -123,6 +123,7 @@ void DEploidIO::init() {
     this->setUseIBD( false );
     this->setUseIbdOnly(false);
     this->setUseLasso( false );
+    this->setUseBestPractice(false);
     this->setDoExportSwitchMissCopy ( true );
     this->setDoAllowInbreeding( false );
     this->mcmcBurn_ = 0.5;
@@ -464,6 +465,8 @@ void DEploidIO::parse () {
             this->readInitialHaps();
         } else if ( *argv_i == "-ibd" ) {
             this->setUseIBD(true);
+        } else if (*argv_i == "-best") {
+            this->setUseBestPractice(true);
         } else if ( *argv_i == "-ibdonly" ) {
             this->setUseIBD(true);
             this->setUseIbdOnly(true);
@@ -1154,59 +1157,3 @@ void DEploidIO::ibdTrimming() {
 }
 
 
-void ChooseK::appendProportions(const vector <double> &p) {
-    this->proportions_.push_back(vector<double>(p.begin(), p.end()));
-}
-
-
-void ChooseK::gatherKs() {
-    assert(this->ks.size() == 0);
-    for (auto const& vec : this->proportions_) {
-        size_t k = 0;
-        for (auto const& v : vec) {
-            if (v > 0.01) {
-                k++;
-            }
-        }
-        ks.push_back(k);
-    }
-    assert(this->ks.size() == this->proportions_.size());
-}
-
-
-void ChooseK::findKmode() {
-    this->gatherKs();
-
-    vector <size_t> uniqueK;
-    for (size_t i = 0; i < this->proportions_[0].size(); i++) {
-        uniqueK.push_back(i+1);
-    }
-    size_t maxCount = 0;
-    this->max_at_ = 0;
-    vector <size_t> kCount(uniqueK.size(), 0);
-    for (auto const& k : ks) {
-        kCount[k-1]++;
-        cout <<k<<" " <<kCount[k-1]<<endl;
-        if (kCount[k-1] > maxCount) {
-            max_at_ = k-1;
-            maxCount = kCount[max_at_];
-        }
-    }
-
-    cout << "k most frequent: " << max_at_+1
-         << " with " << maxCount << "count" << endl;
-}
-
-
-vector <double> ChooseK::chosenP() {
-    this->findKmode();
-    size_t idx = max_at_;
-    for (size_t i = idx; i < this->ks.size(); i++) {
-        if (ks[i] == ks[idx]) {
-            idx = i;
-        }
-    }
-    this->haveChosenP_ = true;
-    return this->proportions_[idx];
-}
-
diff --git a/src/dEploidIO.hpp b/src/dEploidIO.hpp
@@ -73,6 +73,7 @@ class DEploidIO{
     bool useIBD() const { return this->useIBD_;}
     bool useIbdOnly() const { return this->useIbdOnly_;}
     bool useLasso() const { return this->useLasso_;}
+    bool useBestPractice() const { return this->useBestPractice_;}
     void paintIBD();
     double ibdLLK_;
     void getIBDprobsIntegrated(vector < vector <double> > &prob);
@@ -164,9 +165,11 @@ class DEploidIO{
     bool doAllowInbreeding_;
     bool doLsPainting_;
     bool doIbdPainting_;
+
     bool useIBD_;
     bool useIbdOnly_;
     bool useLasso_;
+    bool useBestPractice_;
 
     double vqslod_;
     vector <double> obsWsaf_;
@@ -329,6 +332,7 @@ class DEploidIO{
     void setUseIBD( const bool setTo) { this->useIBD_ = setTo; }
     void setUseIbdOnly(const bool setTo) { this->useIbdOnly_ = setTo;}
     void setUseLasso( const bool setTo) { this->useLasso_ = setTo; }
+    void setUseBestPractice ( const bool setTo) {this->useBestPractice_ = setTo;}
 
     bool initialPropWasGiven() const { return initialPropWasGiven_; }
 
@@ -434,24 +438,4 @@ class DEploidIO{
     void computeObsWsaf();
 };
 
-
-class ChooseK {
- public:
-    ChooseK() {this->haveChosenP_ = false;}
-    ~ChooseK() {}
-    void appendProportions(const vector <double> &p);
-    size_t chosenK() const {
-        assert(this->haveChosenP_ == true);
-        return this->max_at_+1; }
-    vector <double> chosenP();
-
- private:
-    void findKmode();
-    vector < vector<double> > proportions_;
-    vector < size_t > ks;
-    void gatherKs();
-    size_t max_at_;
-    bool haveChosenP_;
-};
-
 #endif