add --min-P

Zilong-Li · Feb 12, 2024 · f08e02e · f08e02e
1 parent b520950
commit f08e02e
Show file tree

Hide file tree

Showing 3 changed files with 14 additions and 4 deletions.
diff --git a/src/admixture.cpp b/src/admixture.cpp
@@ -38,7 +38,7 @@ double Admixture::runOptimalWithBigAss(int ind, const std::unique_ptr<BigAss> &
                 {
                     c12 = c1 * C + c2;
                     double xz = cl(c12, s);
-                    // if(cf(c1, s) < tol || cf(c2, s) < tol) xz = 0.0;
+                    if(AE(c1, s) < magicTol || AE(c2, s) < magicTol) xz = 0.0;
                     double zy = Hz(c1) * Hz(c2);
                     tmp += xz * zy;
                 }
@@ -92,6 +92,7 @@ double Admixture::runNativeWithBigAss(int ind, const std::unique_ptr<BigAss> & g
                 {
                     c12 = c1 * C + c2;
                     double xz = cl(c12, s);
+                    if(AE(c1, s) < magicTol || AE(c2, s) < magicTol) xz = 0.0;
                     for(k1 = 0; k1 < K; k1++)
                     {
                         for(k2 = 0; k2 < K; k2++)
@@ -208,8 +209,11 @@ void Admixture::setStartPoint(const std::unique_ptr<BigAss> & genome, std::strin
     if(!qfile.empty()) load_csv(Q, qfile);
 }
 
-void Admixture::setFlags(bool debug_, bool nonewQ_, bool cF_)
+void Admixture::setFlags(double cftol, double Ftol, double Qtol, bool debug_, bool nonewQ_, bool cF_)
 {
+    magicTol = cftol;
+    clusterFreqThreshold = Ftol;
+    admixtureThreshold = Qtol;
     debug = debug_;
     nonewQ = nonewQ_;
     cF = cF_;
@@ -239,7 +243,7 @@ int run_admix_main(Options & opts)
 
     Admixture admixer(genome->nsamples, genome->G, genome->C, opts.K, opts.seed);
     cao.warn(tim.date(), "-> running admixture with seed =", opts.seed);
-    admixer.setFlags(opts.debug, opts.nQ, opts.cF);
+    admixer.setFlags(opts.ptol, opts.ftol, opts.qtol, opts.debug, opts.nQ, opts.cF);
     admixer.setStartPoint(genome, opts.in_qfile);
     vector<future<double>> llike;
     if(!opts.noaccel)

diff --git a/src/admixture.hpp b/src/admixture.hpp
@@ -18,6 +18,7 @@ class Admixture
     // BOUNDING
     double clusterFreqThreshold{1e-6}; // threshold for F
     double admixtureThreshold{1e-6}; // threshold for Q
+    double magicTol{1e-2}; // threshold for cluster frequency estimated from fastphase
     bool debug = false;
     bool nonewQ = false;
     bool cF = false;
@@ -49,7 +50,7 @@ class Admixture
     void updateIteration();
     void protectPars();
     void constrainF();
-    void setFlags(bool, bool, bool);
+    void setFlags(double, double, double, bool, bool, bool);
     void setStartPoint(const std::unique_ptr<BigAss> & genome, std::string qfile);
     double runNativeWithBigAss(int ind, const std::unique_ptr<BigAss> & genome);
     double runOptimalWithBigAss(int ind, const std::unique_ptr<BigAss> & genome);

diff --git a/src/main.cpp b/src/main.cpp
@@ -217,6 +217,10 @@ int main(int argc, char * argv[])
     cmd_admix.add_argument("-F", "--constrain-F")
         .help("apply constraint on F so that it is not smaller than cluster frequency in fastphase model")
         .flag();
+    cmd_admix.add_argument("-P", "--min-P")
+        .help("set cluster likelihood to zeros if P (in fastphase) < min-P")
+        .default_value(0.0)
+        .scan<'g', double>();
 
     argparse::ArgumentParser cmd_convert("convert", VERSION, default_arguments::help);
     cmd_convert.add_description("different file format converter");
@@ -308,6 +312,7 @@ int main(int argc, char * argv[])
         }
         else if(program.is_subcommand_used(cmd_admix))
         {
+            opts.ptol = cmd_admix.get<double>("--min-P");
             opts.in_bin.assign(cmd_admix.get("--bin"));
             opts.out.assign(cmd_admix.get("--out"));
             opts.seed = cmd_admix.get<int>("--seed");