add ModelTamer original subsample method

Author: HuaiyanRen

alignment/alignment.cpp

@@ -6407,9 +6407,32 @@ void createSUAlignment(Params &params) {
     std::mt19937 gen;
     gen.seed(params.ran_seed);
     int n_site = alignment->getNSite();
-    int n_target_site = static_cast<int>(ceil((params.model_tamer/100.0)*n_site));
+    int n_target_site = static_cast<int>(ceil(params.model_tamer * n_site / 100.0)); //method 1: directly subsample sites
 
     for (int i=0; i<n_sub; i++ ) {
+        if (params.model_tamer_method == 0) {
+            // original ModelTamer method
+            //1. estimated how many distinct pattern are needeed
+            int n_ptn = alignment->getNPattern();
+            int n_target_ptn = static_cast<int>(ceil(params.model_tamer * n_ptn / 100.0));
+
+            //2. initially subsample the estimated needed number of pattern
+            vector<int> init_sub_sites(n_site);
+            std::iota(init_sub_sites.begin(), init_sub_sites.end(), 0);
+            std::shuffle(init_sub_sites.begin(), init_sub_sites.end(), gen);
+            init_sub_sites.resize(n_target_ptn);
+            Alignment *init_sub_alignment = NULL;
+            init_sub_alignment = new Alignment();
+            init_sub_alignment ->extractSites(alignment, init_sub_sites);
+
+            //3. compute how many site are needed to subsample enough pattern
+            int n_init_ptn = init_sub_alignment->getNPattern();
+            n_target_site = (n_target_ptn * n_target_ptn + n_init_ptn - 1) / n_init_ptn; //ceil( (n_target_ptn/n_init_ptn) * n_target_ptn )
+
+            delete init_sub_alignment;
+        }
+        // else: method 1 has been done
+
         // subsample sites
         vector<int> sub_sites(n_site);
         std::iota(sub_sites.begin(), sub_sites.end(), 0);

utils/tools.cpp

@@ -1290,6 +1290,7 @@ void parseArg(int argc, char *argv[], Params &params) {
     params.model_tamer = 100;
     params.model_tamer_sub = 1;
     params.model_tamer_up = 1;
+    params.model_tamer_method = 0;
     params.gamma_shape = -1.0;
     params.min_gamma_shape = MIN_GAMMA_SHAPE;
     params.gamma_median = false;
@@ -3797,6 +3798,15 @@ void parseArg(int argc, char *argv[], Params &params) {
                 if (params.model_tamer_up < 1)
                     throw "Wrong number of ModelTamer upsampling time for -mt-up. Must be at least 1";
                 continue;
+            }
+            if (strcmp(argv[cnt], "-mt-method") == 0 || strcmp(argv[cnt], "--mt-method") == 0) {
+                cnt++;
+                if (cnt >= argc)
+                    throw "Use -mt-method <0|1>";
+                params.model_tamer_method = convert_int(argv[cnt]);
+                if (params.model_tamer_method < 0 || params.model_tamer_method > 1)
+                    throw "Wrong option for -mt-method. Only 0 or 1 is allowed.";
+                continue;
             }
 			if (strcmp(argv[cnt], "-a") == 0) {
 				cnt++;
@@ -7746,6 +7756,7 @@ void Params::setDefault() {
     model_tamer = 100;
     model_tamer_sub = 1;
     model_tamer_up = 1;
+    model_tamer_method = 0;
     gamma_shape = -1.0;
     min_gamma_shape = MIN_GAMMA_SHAPE;
     gamma_median = false;

utils/tools.h

@@ -1765,15 +1765,20 @@ class Params {
     double model_tamer;
 
     /**
-            Model tamer subsampling time.
+            ModelTamer subsampling time.
      */
     int model_tamer_sub;
 
     /**
-            Model tamer upsampling time.
+            ModelTamer upsampling time.
      */
     int model_tamer_up;
 
+    /**
+            ModelTamer subsampling method.
+     */
+    int model_tamer_method;
+
     /**
             shape parameter (alpha) of the Gamma distribution for site rates
      */