Added Bray-Curtis dissimilarity

Author: mtomasini

docs/source/overview.md

@@ -71,3 +71,14 @@ The fixation index is a measure of population differentiation in the presence of
 ```{math}
 \frac{\pi_{\text{Between}} - \pi_{\text{Within}}}{\pi_{\text{Between}}} .
 ```
+
+This estimate is normally too naif in "real life" as it is biased by for example very different population sizes, but it can be used in these idealized conditions, especially as the population sizes of sub-populations is quite regular. 
+
+
+#### Bray-Curtis dissimilarity
+The Bray-Curtis dissimilarity (Bray and Curtis, 1957) is a measure of how how dissimilar two subpopulations are in terms of their species richness. In the context of our model, we define a "species" as a specific set of features. What this implementation does, then, is to count how many individuals possess a certain set of features within a sub-population and compare this to a second subpopulation. In addition to the dissimilarity index between subpopulations, I implemented it between sets of population, which is useful for example to compare two groups of subpopulations as in the case of continental models. The index is implemented with the package `distancia` and is defined as:
+
+```{math}
+D_{\text{BC}} = \frac{\sum_{i=1}^{n}|x_i - y_i|}{\sum_{i=1}^{n}|x_i + y_i|}
+```
+where, for a species $i$, $x_i$ and $y_i$ are the number of individuals of that species found in subpopulations $X$ and $Y$.

exploration.ipynb

@@ -4,6 +4,7 @@
 """
 
 from collections.abc import Set, Iterator
+from distancia import BrayCurtis
 from itertools import pairwise, permutations
 import numpy as np
 import pandas as pd
@@ -355,6 +356,91 @@ def fixation_index(self, subpop_id_1: int, subpop_id_2: int) -> float:
         fixation_index = float((np.mean(pairwise_differences_between) - np.mean(pairwise_differences_within)) / np.mean(pairwise_differences_between))
 
         return fixation_index
+
+
+    def bray_curtis_by_subpopulation_pair(self, subpop_id_1: int, subpop_id_2: int) -> int:
+        """
+        Calculate the Bray-Curtis index of dissimilarity between any two subpopulations. The index is 0 for completely equal subpopulations, 
+        and 1 for completely different ones.
+
+        Args:
+            subpop_id_1 (int): id of the first subpopulation.
+            subpop_id_2 (int): id of the second subpopulation.
+
+        Returns:
+            bray_curtis (int): Bray-Curtis dissimilarity between two subpopulations
+        """
+        individuals_subpop_1 = self.subpopulations[subpop_id_1].population
+        individuals_subpop_2 = self.subpopulations[subpop_id_2].population
+
+        number_of_features = self.number_of_features
+        traits_pop_1 = np.zeros((self.subpopulations[subpop_id_1].get_population_size(), number_of_features))
+        traits_pop_2 = np.zeros((self.subpopulations[subpop_id_2].get_population_size(), number_of_features))
+
+        i = 0
+        for individual in individuals_subpop_1:
+            traits_pop_1[i] = (individual.features)
+            i += 1
+        j = 0
+        for individual in individuals_subpop_2:
+            traits_pop_2[j] = (individual.features)
+            j += 1
+
+#        print(traits_pop_1)
+        traits_pop_1_sorted = np.sort(traits_pop_1)
+        traits_pop_2_sorted = np.sort(traits_pop_2)
+        unique_traits_pop_1, counts_pop_1 = np.unique(traits_pop_1_sorted, axis = 0, return_counts = True)
+        unique_traits_pop_2, counts_pop_2 = np.unique(traits_pop_2_sorted, axis = 0, return_counts = True)
+
+        # Create an instance of the BrayCurti class
+        bray_curtis_instance = BrayCurtis()
+
+        # Calculate the Bray-Curtis distance between the two samples
+        bray_curtis_distance = bray_curtis_instance.calculate(counts_pop_1, counts_pop_2)
+
+        return bray_curtis_distance
+
+    def bray_curtis_by_sets_of_subpopulations(self, group_of_subpop_id_1: List[int], group_of_subpop_id_2: List[int]) -> int:
+        """
+        Calculate the Bray-Curtis index of dissimilarity between any two sets of subpopulations. The index is 0 for completely equal subpopulations, 
+        and 1 for completely different ones.
+
+        Args:
+            group_of_subpop_id_1 (List[int]): list of ids of the first set of subpopulations.
+            group_of_subpop_id_2 (List[int]): list of ids of the second set of subpopulations.
+
+        Returns:
+            bray_curtis (int): Bray-Curtis dissimilarity between two sets of subpopulations
+        """
+
+        traits_pop_1 = []
+        traits_pop_2 = []
+        
+        for subpop_id in group_of_subpop_id_1:
+            individuals = self.subpopulations[subpop_id].population
+            for individual in individuals:
+                traits_pop_1.append(individual.features)
+        traits_pop_1 = np.array(traits_pop_1)
+        
+        for subpop_id in group_of_subpop_id_2:
+            individuals = self.subpopulations[subpop_id].population
+            for individual in individuals:
+                traits_pop_2.append(individual.features)
+        traits_pop_2 = np.array(traits_pop_2)
+
+        traits_pop_1_sorted = np.sort(traits_pop_1)
+        traits_pop_2_sorted = np.sort(traits_pop_2)
+        unique_traits_pop_1, counts_pop_1 = np.unique(traits_pop_1_sorted, axis = 0, return_counts = True)
+        unique_traits_pop_2, counts_pop_2 = np.unique(traits_pop_2_sorted, axis = 0, return_counts = True)
+
+        # Create an instance of the BrayCurti class
+        bray_curtis_instance = BrayCurtis()
+
+        # Calculate the Bray-Curtis distance between the two samples
+        bray_curtis_distance = bray_curtis_instance.calculate(counts_pop_1, counts_pop_2)
+
+        return bray_curtis_distance
+
 
 class SubpopulationIterator(object):
     """

metapypulation/metapopulation.py

@@ -71,6 +71,9 @@ def metapopulation_plot_comparison(dataset_1, dataset_2, title, legend_1, legend
     elif what_measure == 'beta':
         dataset_1_global = pd.read_csv(f"{dataset_1}_beta_diversity.csv", index_col=0)
         dataset_2_global = pd.read_csv(f"{dataset_2}_beta_diversity.csv", index_col=0)
+    elif what_measure == 'bray-curtis-0-4':
+        dataset_1_global = pd.read_csv(f"{dataset_1}_bray-curtis_0-4.csv", index_col=0)
+        dataset_2_global = pd.read_csv(f"{dataset_2}_bray-curtis_0-4.csv", index_col=0)
     else:
         raise ValueError("You need to decide what measure you will plot: 'set_counts', 'shannon', 'simpson' or 'gini'?")