add support for median category rates in discrete GAMMA model

examples/heterotachy/heterotachy.c

@@ -101,7 +101,7 @@ int main(int argc, char * argv[])
 
   /* fixed rate categories */
   double rate_cats[4];
-  pll_compute_gamma_cats(alpha, 4, rate_cats);
+  pll_compute_gamma_cats(alpha, 4, rate_cats, PLL_GAMMA_RATES_MEAN);
 
   /* set frequencies to the parameter sets */
   pll_set_frequencies(partition, 0, frequencies[0]);

examples/lg4/lg4.c

@@ -287,7 +287,7 @@ int main(int argc, char * argv[])
 
   /* compute the discretized category rates from a gamma distribution
      with alpha shape 1 and store them in rate_cats  */
-  pll_compute_gamma_cats(1.00000, 4, rate_cats);
+  pll_compute_gamma_cats(1.00000, 4, rate_cats, PLL_GAMMA_RATES_MEAN);
 
   /* set rate categories */
   pll_set_category_rates(partition, rate_cats);

examples/newick-fasta-rooted/newick-fasta-rooted.c

@@ -288,7 +288,7 @@ int main(int argc, char * argv[])
 
   /* compute the discretized category rates from a gamma distribution
      with alpha shape 1 and store them in rate_cats  */
-  pll_compute_gamma_cats(1, 4, rate_cats);
+  pll_compute_gamma_cats(1, 4, rate_cats, PLL_GAMMA_RATES_MEAN);
 
   /* set frequencies at model with index 0 (we currently have only one model) */
   pll_set_frequencies(partition, 0, frequencies);

examples/newick-fasta-unrooted/newick-fasta-unrooted.c

@@ -234,7 +234,7 @@ int main(int argc, char * argv[])
 
   /* compute the discretized category rates from a gamma distribution
      with alpha shape 1 and store them in rate_cats  */
-  pll_compute_gamma_cats(1, 4, rate_cats);
+  pll_compute_gamma_cats(1, 4, rate_cats, PLL_GAMMA_RATES_MEAN);
 
   /* set frequencies at model with index 0 (we currently have only one model) */
   pll_set_frequencies(partition, 0, frequencies);

examples/newick-phylip-unrooted/newick-phylip-unrooted.c

@@ -220,7 +220,7 @@ int main(int argc, char * argv[])
 
   /* compute the discretized category rates from a gamma distribution
      with alpha shape 1 and store them in rate_cats  */
-  pll_compute_gamma_cats(1, 4, rate_cats);
+  pll_compute_gamma_cats(1, 4, rate_cats, PLL_GAMMA_RATES_MEAN);
 
   /* set frequencies at model with index 0 (we currently have only one model) */
   pll_set_frequencies(partition, 0, frequencies);

examples/newton/newton.c

@@ -132,7 +132,7 @@ int main(int argc, char * argv[])
 
   /* discretized category rates from a gamma distribution with alpha shape 1 */
   double rate_cats[RATES];
-  pll_compute_gamma_cats(alpha, RATES, rate_cats);
+  pll_compute_gamma_cats(alpha, RATES, rate_cats, PLL_GAMMA_RATES_MEAN);
 
   /* set frequencies */
   pll_set_frequencies(partition, 0, frequencies);

examples/partial-traversal/partial.c

@@ -286,7 +286,7 @@ int main(int argc, char * argv[])
 
   /* compute the discretized category rates from a gamma distribution
      with alpha shape 1 and store them in rate_cats  */
-  pll_compute_gamma_cats(1, 4, rate_cats);
+  pll_compute_gamma_cats(1, 4, rate_cats, PLL_GAMMA_RATES_MEAN);
 
   /* set frequencies at model with index 0 (we currently have only one model) */
   pll_set_frequencies(partition, 0, frequencies);

examples/protein-list/protein-list.c

@@ -362,7 +362,7 @@ int main(int argc, char * argv[])
 
   /* compute the discretized category rates from a gamma distribution
      with alpha shape 1 and store them in rate_cats  */
-  pll_compute_gamma_cats(1, 4, rate_cats);
+  pll_compute_gamma_cats(1, 4, rate_cats, PLL_GAMMA_RATES_MEAN);
 
   /* set rate categories */
   pll_set_category_rates(partition, rate_cats);

examples/unrooted/unrooted.c

@@ -54,7 +54,7 @@ int main(int argc, char * argv[])
 
   /* discretized category rates from a gamma distribution with alpha shape 1 */
   double rate_cats[4];
-  pll_compute_gamma_cats(alpha, 4, rate_cats);
+  pll_compute_gamma_cats(alpha, 4, rate_cats, PLL_GAMMA_RATES_MEAN);
 
   /* set frequencies */
   pll_set_frequencies(partition, 0, frequencies);

src/gamma.c

@@ -219,7 +219,8 @@ static double PointChi2 (double prob, double v)
 
 PLL_EXPORT int pll_compute_gamma_cats(double alpha,
                                       unsigned int categories,
-                                      double * output_rates)
+                                      double * output_rates,
+                                      int rates_mode)
 {
   unsigned int i;
 
@@ -245,7 +246,21 @@ PLL_EXPORT int pll_compute_gamma_cats(double alpha,
   {
     output_rates[0] = 1.0;
   }
-  else
+  else if (rates_mode == PLL_GAMMA_RATES_MEDIAN)
+  {
+    double
+      middle = 1.0 / (2.0 * categories),
+      t = 0.0;
+
+    for(i = 0; i < categories; i++)
+      output_rates[i] = POINT_GAMMA((double)(i * 2 + 1) * middle, alfa, beta);
+
+    for (i = 0; i < categories; i++)
+      t += output_rates[i];
+     for( i = 0; i < categories; i++)
+       output_rates[i] *= factor / t;
+  }
+  else if (rates_mode == PLL_GAMMA_RATES_MEAN)
   {
     gammaProbs = (double *)malloc(categories * sizeof(double));
 
@@ -266,6 +281,12 @@ PLL_EXPORT int pll_compute_gamma_cats(double alpha,
 
     free(gammaProbs);
   }
+  else
+  {
+    pll_errno = PLL_ERROR_PARAM_INVALID;
+    snprintf(pll_errmsg, 200, "Invalid GAMMA disrcretization mode (%d)", rates_mode);
+    return PLL_FAILURE;
+  }
 
   return PLL_SUCCESS;
 }

src/pll.h

@@ -151,6 +151,10 @@
 #define PLL_UTREE_SHOW_SCALER_INDEX      (1 << 3)
 #define PLL_UTREE_SHOW_PMATRIX_INDEX     (1 << 4)
 
+/* GAMMA discretization modes */
+#define PLL_GAMMA_RATES_MEAN             0
+#define PLL_GAMMA_RATES_MEDIAN           1
+
 /* structures and data types */
 
 typedef struct pll_partition
@@ -603,7 +607,8 @@ PLL_EXPORT int pll_compute_likelihood_derivatives(pll_partition_t * partition,
 
 PLL_EXPORT int pll_compute_gamma_cats(double alpha,
                                       unsigned int categories,
-                                      double * output_rates);
+                                      double * output_rates,
+                                      int rates_mode);
 
 /* functions in output.c */