From e3678c264edd1f12b8e27a7dfce85d78de7b1821 Mon Sep 17 00:00:00 2001
From: Adrian Tan <atks@umich.edu>
Date: Tue, 21 Mar 2017 13:21:50 -0400
Subject: [PATCH] fixed bug in bcf_ip2g

---
 hts_utils.cpp |  40 +++++++++--
 hts_utils.h   |   5 ++
 test.cpp      | 193 +++++++++++++++++++++++++++++---------------------
 3 files changed, 153 insertions(+), 85 deletions(-)

diff --git a/hts_utils.cpp b/hts_utils.cpp
index e8c67638..40fc7cfc 100644
--- a/hts_utils.cpp
+++ b/hts_utils.cpp
@@ -925,6 +925,11 @@ uint32_t bcf_ag2p(uint32_t no_alleles, uint32_t no_genotypes)
 
 /**
  * Gets genotype from genotype index and ploidy.
+ *
+ * The genotype index is computed by a summation of a series which is 
+ * monotonically decreasing.  This allows you to compute the inverse function
+ * from index to the ordered genotypes by using a "water rapids algorithm" with
+ * decreasing height of each mini water fall.
  */
 std::vector<int32_t> bcf_ip2g(int32_t genotype_index, uint32_t no_ploidy)
 {
@@ -933,22 +938,23 @@ std::vector<int32_t> bcf_ip2g(int32_t genotype_index, uint32_t no_ploidy)
     int32_t pth = no_ploidy;
     int32_t max_allele_index = genotype_index;
     int32_t leftover_genotype_index = genotype_index;
-
+    
     while (pth>0)
     {
+        int32_t old_pth = pth;
+        
         for (int32_t allele_index=0; allele_index <= max_allele_index; ++allele_index)
         {
             int32_t i = choose(pth+allele_index-1, pth);
 
-            if (i>=leftover_genotype_index)
+            if (i>=leftover_genotype_index || allele_index==max_allele_index)
             {
                 if (i>leftover_genotype_index) --allele_index;
                 leftover_genotype_index -= choose(pth+allele_index-1, pth);
                 --pth;
                 max_allele_index = allele_index;
-                genotype[pth] = allele_index;
+                genotype[pth] = allele_index;                
                 break;
-
             }
         }
     }
@@ -980,6 +986,32 @@ uint32_t bcf_g2i(int32_t* g, uint32_t n)
     }
 }
 
+/**
+ * Gets index of a genotype of n ploidy.
+ */
+uint32_t bcf_g2i(std::vector<int32_t>& g)
+{
+    int32_t n = g.size();
+    
+    if (n==1)
+    {
+        return g[0];
+    }
+    if (n==2)
+    {
+        return g[0] + (((g[1]+1)*(g[1]))>>1);
+    }
+    else
+    {
+        uint32_t index = 0;
+        for (uint32_t i=0; i<n; ++i)
+        {
+            index += bcf_ap2g(g[i], i+1);
+        }
+        return index;
+    }
+}
+
 /**
  * Gets a string representation of a variant.
  */
diff --git a/hts_utils.h b/hts_utils.h
index 9fd134c2..2dcb03e3 100644
--- a/hts_utils.h
+++ b/hts_utils.h
@@ -348,6 +348,11 @@ std::vector<int32_t> bcf_ip2g(int32_t genotype_index, uint32_t no_ploidy);
  */
 uint32_t bcf_g2i(int32_t* g, uint32_t p);
 
+/**
+ * Gets index of a genotype of n ploidy.
+ */
+uint32_t bcf_g2i(std::vector<int32_t>& g);
+
 /**
  * Gets index of a diploid genotype.
  */
diff --git a/test.cpp b/test.cpp
index 93ca8768..ff1fde25 100644
--- a/test.cpp
+++ b/test.cpp
@@ -75,83 +75,83 @@ class Igor : Program
         
     };
 
-    /**
-     * n choose r.
-     */
-    uint32_t choose(uint32_t n, uint32_t r)
-    {
-        if (r>n)
-        {
-            return 0;
-        }
-        else if (r==n)
-        {
-            return 1;
-        }
-        else if (r==0)
-        {
-            return 1;
-        }
-        else
-        {
-            if (r>(n>>1))
-            {
-                r = n-r;
-            }
-
-            uint32_t num = n;
-            uint32_t denum = 1;
-
-            for (uint32_t i=1; i<r; ++i)
-            {
-                num *= n-i;
-                denum *= i+1;
-            }
-
-            return num/denum;
-        }
-    }
-
-    
-    /**
-     * Gets number of genotypes from number of alleles and ploidy.
-     */
-    uint32_t bcf_ap2g(uint32_t no_allele, uint32_t no_ploidy)
-    {
-        return choose(no_ploidy+no_allele-1, no_allele-1);
-    }
-
-    /**
-     * Gets number of genotypes from number of alleles and ploidy.
-     */
-    uint32_t bcf_g2i(std::string genotype)
-    {
-        uint32_t index = 0;
-        for (uint32_t i=0; i<genotype.size(); ++i)
-        {
-            uint32_t allele = genotype.at(i)-65;
-            index += bcf_ap2g(allele, i+1);
-        }
-        return index;
-    }
-
-    void print_genotypes(uint32_t A, uint32_t P, std::string genotype)
-    {
-        if (genotype.size()==P)
-        {
-            std::cerr << no << ") " << genotype << " " << bcf_g2i(genotype) << "\n";
-            ++no;
-        }
-        else
-        {
-            for (uint32_t a=0; a<A; ++a)
-            {
-                std::string s(1,(char)(a+65));
-                s.append(genotype);
-                print_genotypes(a+1, P, s);
-            }
-        }
-    }
+//    /**
+//     * n choose r.
+//     */
+//    uint32_t choose(uint32_t n, uint32_t r)
+//    {
+//        if (r>n)
+//        {
+//            return 0;
+//        }
+//        else if (r==n)
+//        {
+//            return 1;
+//        }
+//        else if (r==0)
+//        {
+//            return 1;
+//        }
+//        else
+//        {
+//            if (r>(n>>1))
+//            {
+//                r = n-r;
+//            }
+//
+//            uint32_t num = n;
+//            uint32_t denum = 1;
+//
+//            for (uint32_t i=1; i<r; ++i)
+//            {
+//                num *= n-i;
+//                denum *= i+1;
+//            }
+//
+//            return num/denum;
+//        }
+//    }
+//
+//    
+//    /**
+//     * Gets number of genotypes from number of alleles and ploidy.
+//     */
+//    uint32_t bcf_ap2g(uint32_t no_allele, uint32_t no_ploidy)
+//    {
+//        return choose(no_ploidy+no_allele-1, no_allele-1);
+//    }
+//
+//    /**
+//     * Gets number of genotypes from number of alleles and ploidy.
+//     */
+//    uint32_t bcf_g2i(std::string genotype)
+//    {
+//        uint32_t index = 0;
+//        for (uint32_t i=0; i<genotype.size(); ++i)
+//        {
+//            uint32_t allele = genotype.at(i)-65;
+//            index += bcf_ap2g(allele, i+1);
+//        }
+//        return index;
+//    }
+//
+//    void print_genotypes(uint32_t A, uint32_t P, std::string genotype)
+//    {
+//        if (genotype.size()==P)
+//        {
+//            std::cerr << no << ") " << genotype << " " << bcf_g2i(genotype) << "\n";
+//            ++no;
+//        }
+//        else
+//        {
+//            for (uint32_t a=0; a<A; ++a)
+//            {
+//                std::string s(1,(char)(a+65));
+//                s.append(genotype);
+//                print_genotypes(a+1, P, s);
+//            }
+//        }
+//    }
 
     /**
      * Computes composition and entropy ofrepeat tract.
@@ -274,9 +274,7 @@ class Igor : Program
         }
             
     }
-    
-    
-    
+        
     void test(int argc, char ** argv)
     {   
         version = "0.5";
@@ -601,6 +599,39 @@ class Igor : Program
     };
     
 
+    void test_ip2g(int argc, char ** argv)
+    {   
+        printf("allele\tploidy\tindex\tverified\n");
+        
+        for (int32_t no_allele=2; no_allele<=6; ++no_allele)
+        {
+            for (int32_t no_ploidy=2; no_ploidy<=6; ++no_ploidy)
+            {
+                int32_t max_index = bcf_ap2g(no_allele, no_ploidy);
+                
+                for (int32_t index=0; index<=max_index; ++index)
+                {
+                    fprintf(stderr, "%d\t%d\t%d\t", no_allele, no_ploidy, index);
+                    
+                    
+                    std::vector<int32_t> genotypes = bcf_ip2g(index, no_ploidy);
+                    
+                    for (int32_t j=0; j<no_ploidy; ++j)
+                    {
+                        if (j) printf("/");
+                        printf("%d", genotypes[j]);
+                    }
+                    
+                    uint32_t i = bcf_g2i(genotypes);  
+                
+                    std::string v = (i==index) ? "correct" : "wrong";
+                    
+                    printf("\t%s\n", v.c_str());
+                    
+                }
+            }
+        }
+    };
     ~Igor() {};
 
     private:
@@ -613,7 +644,7 @@ void test(int argc, char ** argv)
     Igor igor(argc, argv);
     
     
-    igor.test(argc, argv);
+    igor.test_ip2g(argc, argv);
 //    printf ("std::isnan(0.0)       : %d\n",std::isnan(0.0));
 //    printf ("std::isnan(1.0/0.0)   : %d\n",std::isnan(1.0/0.0));
 //    printf ("std::isnan(-1.0/0.0)  : %d\n",std::isnan(-1.0/0.0));