Implemented beta diversity

Author: mtomasini

docs/source/overview.md

@@ -54,4 +54,12 @@ The Gini-Simpson index is simply the complement of the Simpson index, that is it
 
 ```{math}
 G = 1 - \sum_{i = 1}^{K} p^2_i = 1 - S.
-```
+```
+
+#### {math}`\beta`-diversity
+{math}`\beta`-diversity is a ratio between diersity over the whole metapopulation and the diversity per subpopulation. We implemented the diversity index suggested by Whittaker (1960; see also Wilson and Shmida, 1984):
+
+```{math}
+\beta_W = \frac{K_m}{\bar{K}_s} - 1
+```
+where {math}`K_m` is the number of sets of traits in the whole metapopulation, and {math}`\bar{K}_s` is the average number of sets of traits in each subpopulation.

exploration.ipynb

@@ -214,7 +214,7 @@ def metapopulation_shannon_diversity(self) -> float:
         return shannon_diversity_index
 
 
-    def metapopulation_test_sets(self) -> int:
+    def metapopulation_count_sets(self) -> int:
         """Calculates number of unique sets of traits in the whole metapopulation.
 
         Returns:
@@ -278,6 +278,23 @@ def metapopulation_gini_diversity(self) -> float:
 
         return gini_diversity_index
 
+    
+    def whittaker_beta_diversity(self) -> float:
+        """
+        Calculate beta diversity as described by Whittaker (1960).
+
+        Returns:
+            float: Whittaker beta diversity.
+        """
+        metapopulation_sets_number = self.metapopulation_count_sets()
+        subpopulation_set_counts = []
+        for subpopulation in self.subpopulations:
+            subpopulation_set_counts.append(subpopulation.count_traits_sets())
+
+        beta_diversity = metapopulation_sets_number / np.mean(subpopulation_set_counts) - 1
+
+        return beta_diversity
+
 
     def count_origin_id_spread(self) -> np.ndarray:
         """

metapypulation/metapopulation.py

@@ -145,7 +145,7 @@ def run_single_replicate(self, replicate_id: int) -> None:
                 shannon.append(np.mean(metapopulation.shannon_diversity_per_subpopulation()))
                 simpson.append(np.mean(metapopulation.simpson_diversity_per_subpopulation()))
                 gini.append(np.mean(metapopulation.gini_diversity_per_subpopulation()))
-                metapop_counts.append(metapopulation.metapopulation_test_sets())
+                metapop_counts.append(metapopulation.metapopulation_count_sets())
                 metapop_shannon.append(metapopulation.metapopulation_shannon_diversity())
                 metapop_simpson.append(metapopulation.metapopulation_simpson_diversity())
                 metapop_gini.append(metapopulation.metapopulation_gini_diversity())

metapypulation/simulation.py

@@ -55,7 +55,7 @@ def plot_comparison_gini(dataset_1, dataset_2, title, legend_1, legend_2, output
     plt.savefig(output_file)
 
 
-def metapopulation_plot_comparison(dataset_1, dataset_2, title, legend_1, legend_2, output_file, what_measure, dataset_3 = None, legend_3 = None):
+def metapopulation_plot_comparison(dataset_1, dataset_2, title, legend_1, legend_2, output_file, what_measure, number_of_pulses = 5, length_of_pulses = 1, settling_period = 99,  dataset_3 = None, legend_3 = None):
     if what_measure == 'set_counts':
         dataset_1_global = pd.read_csv(f"{dataset_1}_metapop_set_counts.csv", index_col=0)
         dataset_2_global = pd.read_csv(f"{dataset_2}_metapop_set_counts.csv", index_col=0)
@@ -68,30 +68,47 @@ def metapopulation_plot_comparison(dataset_1, dataset_2, title, legend_1, legend
     elif what_measure == 'gini':
         dataset_1_global = pd.read_csv(f"{dataset_1}_metapop_gini.csv", index_col=0)
         dataset_2_global = pd.read_csv(f"{dataset_2}_metapop_gini.csv", index_col=0)
+    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)
     else:
         print("You need to decide what measure you will plot: 'set_counts', 'shannon', 'simpson' or 'gini'?")
 
-    plt.plot(dataset_1_global.mean(axis=1), color = 'xkcd:blue')
-    plt.plot(dataset_2_global.mean(axis=1), color = 'xkcd:puce')
+    line1, = plt.plot(dataset_1_global.mean(axis=1), color = 'xkcd:blue', label=f"{legend_1}")
+    line2, = plt.plot(dataset_2_global.mean(axis=1), color = 'xkcd:puce', label=f"{legend_2}")
 
     plt.fill_between(dataset_1_global.index, dataset_1_global.mean(axis=1) - dataset_1_global.std(axis=1), dataset_1_global.mean(axis=1) + dataset_1_global.std(axis=1), color='xkcd:blue', alpha=0.3)
     plt.fill_between(dataset_2_global.index, dataset_2_global.mean(axis=1) - dataset_2_global.std(axis=1), dataset_2_global.mean(axis=1) + dataset_2_global.std(axis=1), color='xkcd:puce', alpha=0.3)
 
     plt.axvline(500, color="black", linestyle='--', ymax=1)
 
     if dataset_3 is not None:
-        # TODO FIX HERE
-        dataset_3_global = pd.read_csv(f"{dataset_3}_metapop_set_counts.csv", index_col=0)
-        plt.plot(dataset_3_global.mean(axis=1), color = 'xkcd:light purple')
+        if what_measure == 'set_counts':
+            dataset_3_global = pd.read_csv(f"{dataset_3}_metapop_set_counts.csv", index_col=0)
+        elif what_measure == 'shannon':
+            dataset_3_global == pd.read_csv(f"{dataset_3}_metapop_shannon.csv", index_col=0)
+        elif what_measure == 'simpson':
+            dataset_3_global = pd.read_csv(f"{dataset_3}_metapop_simpson.csv", index_col=0)
+        elif what_measure == 'gini':
+            dataset_3_global = pd.read_csv(f"{dataset_3}_metapop_gini.csv", index_col=0)
+        elif what_measure == 'beta':
+            dataset_3_global = pd.read_csv(f"{dataset_3}_beta_diversity.csv", index_col=0)
+        else:
+            print("You need to decide what measure you will plot: 'set_counts', 'shannon', 'simpson' or 'gini'?")
+
+        line3, = plt.plot(dataset_3_global.mean(axis=1), color = 'xkcd:light purple', label=f"{legend_3}")
         plt.fill_between(dataset_3_global.index, dataset_3_global.mean(axis=1) - dataset_3_global.std(axis=1), dataset_3_global.mean(axis=1) + dataset_3_global.std(axis=1), color='xkcd:light purple', alpha=0.3)
 
-
     if dataset_3 is not None:
-        plt.legend([f"{legend_1}", f"{legend_2}", f"{legend_3}"])
+        plt.legend(loc=3, handles=[line1, line2, line3])
     else:
         plt.legend([f"{legend_1}", f"{legend_2}"])
 
-    plt.ylabel("Number of feature sets")
+    plt.axvline(500, color="black", linestyle='--', ymax=1)
+    for i in range(1, number_of_pulses + 1):
+        plt.axvline(500 + i*(length_of_pulses + settling_period), color="dimgrey", linestyle='--', ymax=1)
+
+    plt.ylabel(r"Whittaker $\beta$ Diversity")
     plt.xlabel("Generations (x100)")
     plt.grid(linestyle=':', linewidth=0.5)
     plt.title(title)

plotting.py

@@ -10,7 +10,7 @@
 measure_timing = 100
 
 total_population = 800
-replicates = 10
+replicates = 50
 save_output = True
 interactions = ['neutral_interaction'] # 'axelrod_interaction'
 number_of_subpopulations = 8
@@ -26,7 +26,7 @@
 
 count = 0
 
-title = "singleConnection"
+title = "frontConnection"
 migration_config = np.genfromtxt(f"./configs/maritime_configs/{title}.csv", delimiter=",")
 
 for interaction in interactions:
@@ -36,6 +36,7 @@
         subpop_gini_df = pd.DataFrame()
         metapop_set_counts_df = pd.DataFrame()
         metapop_gini_df = pd.DataFrame()
+        beta_diversity_df = pd.DataFrame()
 
         carrying_capacity = int(np.ceil(total_population / number_of_subpopulations)) # [283, 39, 39, 39]# 
 
@@ -53,6 +54,7 @@
             gini = []
             metapop_counts = []
             metapop_gini = []
+            beta_diversity = []
 
             t = 0
             for gen in range(burn_in):
@@ -67,8 +69,9 @@
                 if t%measure_timing == 0:
                     set_counts.append(np.mean(metapopulation.traits_sets_per_subpopulation()))
                     gini.append(np.mean(metapopulation.gini_diversity_per_subpopulation()))
-                    metapop_counts.append(metapopulation.metapopulation_test_sets())
+                    metapop_counts.append(metapopulation.metapopulation_count_sets())
                     metapop_gini.append(metapopulation.metapopulation_gini_diversity())
+                    beta_diversity.append(metapopulation.whittaker_beta_diversity())
 
                 t += 1
 
@@ -95,15 +98,17 @@
                     if t%measure_timing == 0:
                         set_counts.append(np.mean(metapopulation.traits_sets_per_subpopulation()))
                         gini.append(np.mean(metapopulation.gini_diversity_per_subpopulation()))
-                        metapop_counts.append(metapopulation.metapopulation_test_sets())
+                        metapop_counts.append(metapopulation.metapopulation_count_sets())
                         metapop_gini.append(metapopulation.metapopulation_gini_diversity())
+                        beta_diversity.append(metapopulation.whittaker_beta_diversity())
 
                     t += 1
 
             subpop_set_counts_df = pd.concat([subpop_set_counts_df, pd.Series(set_counts, name=replicate_id)], axis=1)
             subpop_gini_df = pd.concat([subpop_gini_df, pd.Series(gini, name=replicate_id)], axis=1)
             metapop_set_counts_df = pd.concat([metapop_set_counts_df, pd.Series(metapop_counts, name=replicate_id)], axis=1)
             metapop_gini_df = pd.concat([metapop_gini_df, pd.Series(metapop_gini, name=replicate_id)], axis=1)
+            beta_diversity_df = pd.concat([beta_diversity_df, pd.Series(beta_diversity, name=replicate_id)], axis=1)
 
             if verbose:
                 print(f"Generation reached: {t}")
@@ -116,10 +121,11 @@
             # print(metapopulation.count_origin_id_spread())
 
         if save_output:
-            subpop_set_counts_df.to_csv(f"./Outputs/Kiel2025/{title}/{interaction}_{rate_of_migration}_{number_of_pulses}pulses_{pulse_length}gen_subpop_set_counts.csv", sep=",")
-            subpop_gini_df.to_csv(f"./Outputs/Kiel2025/{title}/{interaction}_{rate_of_migration}_{number_of_pulses}pulses_{pulse_length}gen_subpop_gini.csv", sep=",")
-            metapop_set_counts_df.to_csv(f"./Outputs/Kiel2025/{title}/{interaction}_{rate_of_migration}_{number_of_pulses}pulses_{pulse_length}gen_metapop_set_counts.csv", sep=",")
-            metapop_gini_df.to_csv(f"./Outputs/Kiel2025/{title}/{interaction}_{rate_of_migration}_{number_of_pulses}pulses_{pulse_length}gen_metapop_gini.csv", sep=",")
+            # subpop_set_counts_df.to_csv(f"./Outputs/Kiel2025/{title}/{interaction}_{rate_of_migration}_{number_of_pulses}pulses_{pulse_length}gen_subpop_set_counts.csv", sep=",")
+            # subpop_gini_df.to_csv(f"./Outputs/Kiel2025/{title}/{interaction}_{rate_of_migration}_{number_of_pulses}pulses_{pulse_length}gen_subpop_gini.csv", sep=",")
+            # metapop_set_counts_df.to_csv(f"./Outputs/Kiel2025/{title}/{interaction}_{rate_of_migration}_{number_of_pulses}pulses_{pulse_length}gen_metapop_set_counts.csv", sep=",")
+            # metapop_gini_df.to_csv(f"./Outputs/Kiel2025/{title}/{interaction}_{rate_of_migration}_{number_of_pulses}pulses_{pulse_length}gen_metapop_gini.csv", sep=",")
+            beta_diversity_df.to_csv(f"./Outputs/Kiel2025/{title}/{interaction}_{rate_of_migration}_{number_of_pulses}pulses_{pulse_length}gen_beta_diversity.csv", sep=",")
 
         # subpop_gini_df.to_csv("test.csv")
 

script_Kiel2025_diversity.py

@@ -65,7 +65,7 @@
             counts_pop_6.append(metapop.subpopulations[5].count_traits_sets())
             counts_pop_7.append(metapop.subpopulations[6].count_traits_sets())
             counts_pop_8.append(metapop.subpopulations[7].count_traits_sets())
-            counts_metapop.append(metapop.metapopulation_test_sets())
+            counts_metapop.append(metapop.metapopulation_count_sets())
             subpops_with_mutation.append(0)
         # if t%50000 == 0:
         #     print(f"Gen {t}!")
@@ -86,7 +86,7 @@
             counts_pop_6.append(metapop.subpopulations[5].count_traits_sets())
             counts_pop_7.append(metapop.subpopulations[6].count_traits_sets())
             counts_pop_8.append(metapop.subpopulations[7].count_traits_sets())
-            counts_metapop.append(metapop.metapopulation_test_sets())
+            counts_metapop.append(metapop.metapopulation_count_sets())
             subpops_with_mutation.append(0)
         # if t%50000 == 0:
         #     print(f"Gen {t}!")
@@ -119,7 +119,7 @@
             counts_pop_6.append(metapop.subpopulations[5].count_traits_sets())
             counts_pop_7.append(metapop.subpopulations[6].count_traits_sets())
             counts_pop_8.append(metapop.subpopulations[7].count_traits_sets())
-            counts_metapop.append(metapop.metapopulation_test_sets())
+            counts_metapop.append(metapop.metapopulation_count_sets())
 
             feature_tests = []
             for subpop in metapop.subpopulations: