Changed gamma distribution to dirichlet distribution fir continuous w…

…eights

cmd/distboot.go

@@ -13,7 +13,7 @@ var distbootSeed int64
 var distbootOutput string
 var distbootnb int
 var distbootmodel string
-var distbootgamma bool
+var distbootcontinuous bool
 var distboolRemoveGaps bool
 
 // distbootCmd represents the distboot command
@@ -35,7 +35,11 @@ Available Distances:
 
 For example:
 
-goalign build distboot -m k2p -i align.fa -o mats.txt`,
+goalign build distboot -m k2p -i align.fa -o mats.txt
+
+If -c is given, then random continuous weights are associated to all sites. 
+Weights follow a Dirichlet distribution D(n;1,...,1)
+`,
 
 	Run: func(cmd *cobra.Command, args []string) {
 		rand.Seed(distbootSeed)
@@ -45,8 +49,8 @@ goalign build distboot -m k2p -i align.fa -o mats.txt`,
 		model := distance.Model(distbootmodel, distboolRemoveGaps)
 		for i := 0; i < distbootnb; i++ {
 			var weights []float64 = nil
-			if distbootgamma {
-				weights = distance.BuildWeights(align)
+			if distbootcontinuous {
+				weights = distance.BuildWeightsDirichlet(align)
 				distMatrix := distance.DistMatrix(align, weights, model, rootcpus)
 				writeDistBootMatrix(distMatrix, f)
 			} else {
@@ -65,7 +69,7 @@ func init() {
 	distbootCmd.PersistentFlags().StringVarP(&distbootOutput, "output", "o", "stdout", "Distance matrices output file")
 	distbootCmd.PersistentFlags().StringVarP(&distbootmodel, "model", "m", "k2p", "Model for distance computation")
 	distbootCmd.PersistentFlags().IntVarP(&distbootnb, "nboot", "n", 1, "Number of bootstrap replicates to build")
-	distbootCmd.PersistentFlags().BoolVarP(&distbootgamma, "gamma", "g", false, "Bootstraps are done by weighting alignment using gamma generated weights")
+	distbootCmd.PersistentFlags().BoolVarP(&distbootcontinuous, "continuous", "c", false, "Bootstraps are done by weighting alignment with continuous weights (dirichlet)")
 	distbootCmd.PersistentFlags().BoolVarP(&distboolRemoveGaps, "rm-gaps", "r", false, "Do not take into account positions containing >=1 gaps")
 }
 

cmd/tntweightboot.go

@@ -12,8 +12,12 @@ import (
 // tntweightbootCmd represents the tntweightboot command
 var tntweightbootCmd = &cobra.Command{
 	Use:   "tntweightboot",
-	Short: "",
-	Long:  ``,
+	Short: "Generates a data file to give as TNT input, with continuous weight for sites",
+	Long: `Generates a data file to give as TNT input, with continuous weight for sites
+
+Weights follow a Dirichlet distribution D(n;1,...,1)
+
+`,
 	Run: func(cmd *cobra.Command, args []string) {
 		rand.Seed(weightbootSeed)
 
@@ -29,7 +33,7 @@ var tntweightbootCmd = &cobra.Command{
 				f = openWriteFile(weightbootOutput)
 			}
 			var weights []float64 = nil
-			weights = distance.BuildWeights(al)
+			weights = distance.BuildWeightsDirichlet(al)
 			minw := -1.0
 			maxw := -1.0
 			for _, w := range weights {

cmd/weightboot.go

@@ -15,10 +15,12 @@ var weightbootnb int
 // weightbootCmd represents the weightboot command
 var weightbootCmd = &cobra.Command{
 	Use:   "weightboot",
-	Short: "generate weights for all positions of the input alignment",
-	Long: `generate weights for all positions of the input alignment
+	Short: "generate continous weights for all positions of the input alignment",
+	Long: `generate continous weights for all positions of the input alignment
 
-If the input alignment contains several alignments, will process the first one only
+If the input alignment contains several alignments, will process the first one only.
+
+Weights follow a Dirichlet distribtion D(n;1,...,1)
 
 `,
 	Run: func(cmd *cobra.Command, args []string) {
@@ -29,7 +31,7 @@ If the input alignment contains several alignments, will process the first one o
 		f := openWriteFile(weightbootOutput)
 		for i := 0; i < weightbootnb; i++ {
 			var weights []float64 = nil
-			weights = distance.BuildWeights(al)
+			weights = distance.BuildWeightsDirichlet(al)
 			for i, w := range weights {
 				if i > 0 {
 					f.WriteString("\t")

distance/distance.go

@@ -53,8 +53,8 @@ func Model(modelType string, removegaps bool) DistModel {
 	return model
 }
 
-/* Return a normalized vector of weights */
-func BuildWeights(al align.Alignment) []float64 {
+/* Return a normalized vector of weights following a Gamma distribution*/
+func BuildWeightsGamma(al align.Alignment) []float64 {
 	outweights := make([]float64, al.Length())
 	// Parameters of the binomial
 	n := float64(al.Length())
@@ -76,6 +76,18 @@ func BuildWeights(al align.Alignment) []float64 {
 	return outweights
 }
 
+/* Returns a vector of weights following a Dirichlet distribution D(n ; 1,...,1)
+   with n alignment length
+*/
+func BuildWeightsDirichlet(al align.Alignment) []float64 {
+	alpha := make([]float64, al.Length(), al.Length())
+	for i := 0; i < al.Length(); i++ {
+		alpha[i] = 1
+	}
+	outweights := stats.Dirichlet(float64(al.Length()), alpha...)
+	return outweights
+}
+
 /* Compute a matrix distance, with weights associated to each alignment positions */
 /* If weights == nil, then all weights are considered 1 */
 func DistMatrix(al align.Alignment, weights []float64, model DistModel, cpus int) [][]float64 {

stats/dirichlet.go

@@ -9,7 +9,8 @@ import (
 // Dirichlet returns a set of random numbers from dirichlet distribution
 // See https://en.wikipedia.org/wiki/Dirichlet_distribution#Random_number_generation
 // For further information
-func Dirichlet(alpha ...float64) []float64 {
+// Factor: multiplying factor to apply to all values : Normally it is 1
+func Dirichlet(factor float64, alpha ...float64) []float64 {
 	if len(alpha) <= 2 {
 		io.ExitWithMessage(errors.New(fmt.Sprintf("Alpha parameter vector for Dirichlet sample contains less than 2 values")))
 	}
@@ -23,7 +24,7 @@ func Dirichlet(alpha ...float64) []float64 {
 		sum += sample[i]
 	}
 	for i, _ := range alpha {
-		sample[i] = sample[i] / sum
+		sample[i] = factor * sample[i] / sum
 	}
 	return sample
 }