Implemented Gini index

Author: mtomasini

configs/maritime_configs/README.md

@@ -0,0 +1,13 @@
+### Migration configuration matrices for Kiel 2025
+In this folder we host the migration configuration matrices to run the experiments of maritime connectivity that we will present in Kiel. 
+
+We generally assume that a maritime domain is composed by two sides of a metapopulation that are constantly connected, with migration between the two sides changing in time and configuration. 
+
+![Basic configuration of a maritime metapopulation](../../img/metapopulation_configuration.png)
+
+From here, we consider different ways in which the two halves of the metapopulation can meet. Does each subpopulation on one side connect to a different subpopulation on the other side? Do they all connect to only one subpopulation on the other side? Or is there only a contact between two single subpopulations? And how long does this contact last? How long does it need to last for the migration event to have a lasting effect on the other subpopulation?
+
+Basic configurations that we will adopt:
+- island model
+- front connection (migration possible between all facing subpopulations)
+- single connection (migration possible only between two facing subpopulations)

docs/source/overview.md

@@ -12,7 +12,7 @@ Metapypulation has three main classes:
 
 In addition, I provide a [`Simulation`](https://mtomasini.github.io/MetapopulationsPython/metapypulation.html#module-metapypulation.simulation) class, which allows to run a simulation with several replicates, which outputs different measurements ([see below](#diversity-measures)). The class also provides some quick tools to plot the results of the simulation.
 
-## Cultural traits
+## The spread of cultural traits
 
 Currently, each individual's culture is represented by a set of {math}`N` features. Each feature in turn can assume one of {math}`\nu` traits. These features represent different (assumed) independent facets of culture: one could be language, burial tradition, boat building features, *etc*. So each individual is currently represented by a vector of integers.
 
@@ -26,23 +26,32 @@ While we plan on adding several different ways for individuals to interact and c
 
 ### Diversity measures
 
-Currently, there are two diversity measures implemented at the level of both the subpopulation and the whole metapopulation. The first is the Shannon diversity index, which for each feature is measured as
+Currently, we have implemented a few diversity measures at the level of both the subpopulation and the whole metapopulation (to avoid repeating "subpopulation / metapopulation", we refer to either as the "reference population" for the measure. 
 
-```{math}
-H^{\prime} = - \sum_{i = 1}^{\nu} p_i \ln p_i ,
-```
+#### Count of sets
+The first measure that we implemented is the number of unique sets of traits, {math}`K`. In code, this is done through the function `np.unique(..., return_counts = True)`. 
 
-where {math}`p_i` is the frequency of trait {math}`i` in the subpopulation / metapopulation. Then, the Shannon diversity index that we measure is the average of the index for each trait,
+#### Shannon diversity
+The second diversity index that we have implemented is the Shannon diversity index, which measures the entropy in the reference population. For {math}`K` sets of traits,
 
 ```{math}
-\bar{H} = H^{\prime} .
+H^{\prime} = - \sum_{i = 1}^{K} p_i \ln p_i ,
 ```
 
-The second diversity measure that we calculate is the number of unique sets of traits, {math}`K`. In code, this is done through the function `np.unique(..., return_counts = True)`. 
-We also implemented the (inversed) Simpson's diversity index, defined as 
+where {math}`p_i` is the frequency of a type of individual {math}`i` in the subpopulation / metapopulation. Shannon is a measure of how "easy" it is, given an individual of type `j`, to predict the set of traits of another individual in the same reference population.
+
+#### Simpson index
+We also implemented the Simpson's diversity index, defined as 
 
 ```{math}
-S = \frac{N(N-1)}{\sum_{k}n_k(n_k-1)} ,
+S = \sum_{i = 1}^{K} p^2_i .
 ```
 
-where {math}`n_k` is the number of individuals with the set of traits {math}`k`. This is implemented both at the level of the subpopulation and of the metapopulation.
+The Simpson index is bound between 0 and 1, and represents the probability that two individuals in the same reference population have the same set of traits.
+
+#### Gini-Simpson index
+The Gini-Simpson index is simply the complement of the Simpson index, that is it measures the probability that two random indeividuals in the same reference population do not share a set of traits. This is calculated as
+
+```{math}
+G = 1 - \sum_{i = 1}^{K} p^2_i = 1 - S
+```

exploration.ipynb

@@ -2090,53 +2090,61 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 2,
+   "execution_count": null,
    "metadata": {},
    "outputs": [
     {
      "name": "stdout",
      "output_type": "stream",
      "text": [
       "Gen 0!\n",
-      "The current number of sets in pop 0 is 100\n",
+      "The current number of sets in pop 0 is 99\n",
       "The current number of sets in pop 1 is 100\n",
       "The current number of sets in pop 2 is 100\n",
       "The current number of sets in pop 3 is 100\n",
       "Gen 50000!\n",
-      "The current number of sets in pop 0 is 4\n",
-      "The current number of sets in pop 1 is 7\n",
-      "The current number of sets in pop 2 is 5\n",
-      "The current number of sets in pop 3 is 8\n",
-      "[<bound method Subpopulation.count_deme_origin_id of <metapypulation.subpopulation.Subpopulation object at 0x7d0c293abe00>>\n",
-      " <bound method Subpopulation.count_deme_origin_id of <metapypulation.subpopulation.Subpopulation object at 0x7d0c2940f110>>\n",
-      " <bound method Subpopulation.count_deme_origin_id of <metapypulation.subpopulation.Subpopulation object at 0x7d0c2940f250>>\n",
-      " <bound method Subpopulation.count_deme_origin_id of <metapypulation.subpopulation.Subpopulation object at 0x7d0c2924a3f0>>]\n",
+      "The current number of sets in pop 0 is 10\n",
+      "The current number of sets in pop 1 is 4\n",
+      "The current number of sets in pop 2 is 12\n",
+      "The current number of sets in pop 3 is 2\n",
+      "[[100   0   0   0]\n",
+      " [  0 100   0   0]\n",
+      " [  0   0 100   0]\n",
+      " [  0   0   0 100]]\n",
       "Gen 100000!\n",
-      "The current number of sets in pop 0 is 7\n",
-      "The current number of sets in pop 1 is 5\n",
-      "The current number of sets in pop 2 is 11\n",
-      "The current number of sets in pop 3 is 7\n",
+      "The current number of sets in pop 0 is 3\n",
+      "The current number of sets in pop 1 is 1\n",
+      "The current number of sets in pop 2 is 8\n",
+      "The current number of sets in pop 3 is 4\n",
       "Gen 150000!\n",
-      "The current number of sets in pop 0 is 55\n",
-      "The current number of sets in pop 1 is 50\n",
-      "The current number of sets in pop 2 is 58\n",
-      "The current number of sets in pop 3 is 50\n"
-     ]
-    },
-    {
-     "ename": "KeyboardInterrupt",
-     "evalue": "",
-     "output_type": "error",
-     "traceback": [
-      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
-      "\u001b[0;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
-      "Cell \u001b[0;32mIn[2], line 52\u001b[0m\n\u001b[1;32m     45\u001b[0m             \u001b[38;5;28mprint\u001b[39m(\u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mThe current number of sets in pop \u001b[39m\u001b[38;5;132;01m{\u001b[39;00msubpopulation\u001b[38;5;241m.\u001b[39mid\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m is \u001b[39m\u001b[38;5;132;01m{\u001b[39;00msubpopulation\u001b[38;5;241m.\u001b[39mcount_traits_sets()\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m\"\u001b[39m)\n\u001b[1;32m     46\u001b[0m             \u001b[38;5;66;03m# sub_id = []\u001b[39;00m\n\u001b[1;32m     47\u001b[0m             \u001b[38;5;66;03m# for ind in subpopulation.population.individuals:\u001b[39;00m\n\u001b[1;32m     48\u001b[0m             \u001b[38;5;66;03m#     sub_id.append(ind.original_deme_id)\u001b[39;00m\n\u001b[1;32m     49\u001b[0m             \u001b[38;5;66;03m# indexes, counts = np.unique(sub_id, return_counts=True)\u001b[39;00m\n\u001b[1;32m     50\u001b[0m             \u001b[38;5;66;03m# print(f\"The current deme index present in population {subpopulation.id} are {indexes} with {counts} counts.\")\u001b[39;00m\n\u001b[0;32m---> 52\u001b[0m     \u001b[43mmetapopulation\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mmigrate\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m     53\u001b[0m     metapopulation\u001b[38;5;241m.\u001b[39mmake_interact()\n\u001b[1;32m     55\u001b[0m spread \u001b[38;5;241m=\u001b[39m metapopulation\u001b[38;5;241m.\u001b[39mcount_origin_id_spread()\n",
-      "File \u001b[0;32m~/Documents/06_Development/03_GRIDH/ArchaeologyMetapopulations/MetapopulationsPython/metapypulation/metapopulation.py:109\u001b[0m, in \u001b[0;36mMetapopulation.migrate\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m    106\u001b[0m         migration_rate \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mmigration_matrix[subpopulation\u001b[38;5;241m.\u001b[39mid][index]\n\u001b[1;32m    107\u001b[0m         receiving_subpopulation \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39msubpopulations[index]\n\u001b[0;32m--> 109\u001b[0m         \u001b[43mreceiving_subpopulation\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mreceive_migrants\u001b[49m\u001b[43m(\u001b[49m\u001b[43msubpopulation\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mmigration_rate\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    111\u001b[0m \u001b[38;5;28;01mfor\u001b[39;00m subpopulation \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39msubpopulations:\n\u001b[1;32m    112\u001b[0m     subpopulation\u001b[38;5;241m.\u001b[39mincorporate_migrants_in_population()\n",
-      "File \u001b[0;32m~/Documents/06_Development/03_GRIDH/ArchaeologyMetapopulations/MetapopulationsPython/metapypulation/subpopulation.py:90\u001b[0m, in \u001b[0;36mSubpopulation.receive_migrants\u001b[0;34m(self, giving_subpopulation, migration_rate)\u001b[0m\n\u001b[1;32m     88\u001b[0m number_of_migrants \u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39mrandom\u001b[38;5;241m.\u001b[39mbinomial(population_size, migration_rate)\n\u001b[1;32m     89\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m number_of_migrants \u001b[38;5;241m>\u001b[39m \u001b[38;5;241m0\u001b[39m:\n\u001b[0;32m---> 90\u001b[0m     individuals_to_remove \u001b[38;5;241m=\u001b[39m \u001b[43mgiving_subpopulation\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mpopulation\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43msample_and_remove\u001b[49m\u001b[43m(\u001b[49m\u001b[43mnumber_of_migrants\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m     91\u001b[0m     \u001b[38;5;28;01mfor\u001b[39;00m individual \u001b[38;5;129;01min\u001b[39;00m individuals_to_remove:\n\u001b[1;32m     92\u001b[0m         \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mincoming_migrants\u001b[38;5;241m.\u001b[39madd(individual)\n",
-      "File \u001b[0;32m~/Documents/06_Development/03_GRIDH/ArchaeologyMetapopulations/MetapopulationsPython/metapypulation/subpopulation.py:332\u001b[0m, in \u001b[0;36mSetOfIndividuals.sample_and_remove\u001b[0;34m(self, number_of_individuals)\u001b[0m\n\u001b[1;32m    322\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21msample_and_remove\u001b[39m(\u001b[38;5;28mself\u001b[39m, number_of_individuals: \u001b[38;5;28mint\u001b[39m) \u001b[38;5;241m-\u001b[39m\u001b[38;5;241m>\u001b[39m List[Individual]:\n\u001b[1;32m    323\u001b[0m \u001b[38;5;250m    \u001b[39m\u001b[38;5;124;03m\"\"\"\u001b[39;00m\n\u001b[1;32m    324\u001b[0m \u001b[38;5;124;03m    Sample an individual, remove it from the Set and return it in a list.\u001b[39;00m\n\u001b[1;32m    325\u001b[0m \n\u001b[0;32m   (...)\u001b[0m\n\u001b[1;32m    330\u001b[0m \u001b[38;5;124;03m        List[Individual]: List of all the individuals that have been sampled from the population.\u001b[39;00m\n\u001b[1;32m    331\u001b[0m \u001b[38;5;124;03m    \"\"\"\u001b[39;00m\n\u001b[0;32m--> 332\u001b[0m     \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mshuffle\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    334\u001b[0m     list_of_individuals \u001b[38;5;241m=\u001b[39m []\n\u001b[1;32m    335\u001b[0m     \u001b[38;5;28;01mfor\u001b[39;00m i \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28mrange\u001b[39m(number_of_individuals):\n",
-      "File \u001b[0;32m~/Documents/06_Development/03_GRIDH/ArchaeologyMetapopulations/MetapopulationsPython/metapypulation/subpopulation.py:320\u001b[0m, in \u001b[0;36mSetOfIndividuals.shuffle\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m    316\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mshuffle\u001b[39m(\u001b[38;5;28mself\u001b[39m) \u001b[38;5;241m-\u001b[39m\u001b[38;5;241m>\u001b[39m \u001b[38;5;28;01mNone\u001b[39;00m:\n\u001b[1;32m    317\u001b[0m \u001b[38;5;250m    \u001b[39m\u001b[38;5;124;03m\"\"\"\u001b[39;00m\n\u001b[1;32m    318\u001b[0m \u001b[38;5;124;03m    Shuffle the Set.\u001b[39;00m\n\u001b[1;32m    319\u001b[0m \u001b[38;5;124;03m    \"\"\"\u001b[39;00m\n\u001b[0;32m--> 320\u001b[0m     \u001b[43mrandom\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mshuffle\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mindividuals\u001b[49m\u001b[43m)\u001b[49m\n",
-      "File \u001b[0;32m~/Software/miniforge3/envs/abm/lib/python3.13/random.py:361\u001b[0m, in \u001b[0;36mRandom.shuffle\u001b[0;34m(self, x)\u001b[0m\n\u001b[1;32m    358\u001b[0m \u001b[38;5;28;01mfor\u001b[39;00m i \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28mreversed\u001b[39m(\u001b[38;5;28mrange\u001b[39m(\u001b[38;5;241m1\u001b[39m, \u001b[38;5;28mlen\u001b[39m(x))):\n\u001b[1;32m    359\u001b[0m     \u001b[38;5;66;03m# pick an element in x[:i+1] with which to exchange x[i]\u001b[39;00m\n\u001b[1;32m    360\u001b[0m     j \u001b[38;5;241m=\u001b[39m randbelow(i \u001b[38;5;241m+\u001b[39m \u001b[38;5;241m1\u001b[39m)\n\u001b[0;32m--> 361\u001b[0m     x[i], x[j] \u001b[38;5;241m=\u001b[39m x[j], x[i]\n",
-      "\u001b[0;31mKeyboardInterrupt\u001b[0m: "
+      "The current number of sets in pop 0 is 48\n",
+      "The current number of sets in pop 1 is 43\n",
+      "The current number of sets in pop 2 is 51\n",
+      "The current number of sets in pop 3 is 51\n",
+      "Gen 200000!\n",
+      "The current number of sets in pop 0 is 13\n",
+      "The current number of sets in pop 1 is 15\n",
+      "The current number of sets in pop 2 is 14\n",
+      "The current number of sets in pop 3 is 16\n",
+      "Gen 250000!\n",
+      "The current number of sets in pop 0 is 14\n",
+      "The current number of sets in pop 1 is 16\n",
+      "The current number of sets in pop 2 is 19\n",
+      "The current number of sets in pop 3 is 19\n",
+      "Gen 300000!\n",
+      "The current number of sets in pop 0 is 9\n",
+      "The current number of sets in pop 1 is 10\n",
+      "The current number of sets in pop 2 is 15\n",
+      "The current number of sets in pop 3 is 12\n",
+      "Gen 350000!\n",
+      "The current number of sets in pop 0 is 17\n",
+      "The current number of sets in pop 1 is 14\n",
+      "The current number of sets in pop 2 is 18\n",
+      "The current number of sets in pop 3 is 15\n",
+      "[[23 21 27 22]\n",
+      " [26 18 25 18]\n",
+      " [24 28 24 27]\n",
+      " [27 33 24 33]]\n"
      ]
     }
    ],
@@ -2158,11 +2166,11 @@
     "simulation_time = 300000\n",
     "for t in range(burn_in):\n",
     "    if t%1000 == 0:\n",
-    "        counts_pop_1.append(metapopulation.subpopulations[0].simpson_diversity())\n",
-    "        counts_pop_2.append(metapopulation.subpopulations[1].simpson_diversity())\n",
-    "        counts_pop_3.append(metapopulation.subpopulations[2].simpson_diversity())\n",
-    "        counts_pop_4.append(metapopulation.subpopulations[3].simpson_diversity())\n",
-    "        counts_metapop.append(metapopulation.metapopulation_simpson())\n",
+    "        counts_pop_1.append(metapopulation.subpopulations[0].gini_diversity())\n",
+    "        counts_pop_2.append(metapopulation.subpopulations[1].gini_diversity())\n",
+    "        counts_pop_3.append(metapopulation.subpopulations[2].gini_diversity())\n",
+    "        counts_pop_4.append(metapopulation.subpopulations[3].gini_diversity())\n",
+    "        counts_metapop.append(metapopulation.metapopulation_gini())\n",
     "    if t%50000 == 0:\n",
     "        print(f\"Gen {t}!\")\n",
     "        for subpopulation in metapopulation.subpopulations:\n",
@@ -2177,11 +2185,11 @@
     "metapopulation.migration_matrix = migrations\n",
     "for t in range(simulation_time):\n",
     "    if t%1000 == 0:\n",
-    "        counts_pop_1.append(metapopulation.subpopulations[0].simpson_diversity())\n",
-    "        counts_pop_2.append(metapopulation.subpopulations[1].simpson_diversity())\n",
-    "        counts_pop_3.append(metapopulation.subpopulations[2].simpson_diversity())\n",
-    "        counts_pop_4.append(metapopulation.subpopulations[3].simpson_diversity())\n",
-    "        counts_metapop.append(metapopulation.metapopulation_simpson())\n",
+    "        counts_pop_1.append(metapopulation.subpopulations[0].gini_diversity())\n",
+    "        counts_pop_2.append(metapopulation.subpopulations[1].gini_diversity())\n",
+    "        counts_pop_3.append(metapopulation.subpopulations[2].gini_diversity())\n",
+    "        counts_pop_4.append(metapopulation.subpopulations[3].gini_diversity())\n",
+    "        counts_metapop.append(metapopulation.metapopulation_gini())\n",
     "    if t%50000 == 0:\n",
     "        print(f\"Gen {t + burn_in}!\")\n",
     "        for subpopulation in metapopulation.subpopulations:\n",

img/metapopulation_configuration.odg

@@ -204,6 +204,7 @@ def metapopulation_test_sets(self) -> int:
 
         return len(uniques) 
 
+
     def metapopulation_simpson(self) -> float:
         """
         Calculates Simpson diversity index over the whole metapopulation.
@@ -227,6 +228,29 @@ def metapopulation_simpson(self) -> float:
         return simpson_diversity_index
 
 
+    def metapopulation_gini(self) -> float:
+        """
+        Calculates Gini-Simpson diversity index over the whole metapopulation.
+
+        Returns:
+            float: Gini-Simpson diversity index of the metapopulation.
+        """
+        number_of_features = self.number_of_features
+        traits = np.zeros((self.get_metapopulation_size(), number_of_features))
+        i = 0
+        for subpopulation in self.subpopulations:
+            for individual in subpopulation.population:
+                traits[i] = (individual.features)
+                i += 1
+            
+        uniques, counts = np.unique(traits, axis = 0, return_counts = True)
+
+        frequencies = counts / self.get_metapopulation_size()
+        gini_diversity_index = 1 - np.sum(frequencies*frequencies)
+        
+        return gini_diversity_index
+
+
     def count_origin_id_spread(self) -> np.ndarray:
         """
         Counts the spread of individuals from deme of origin by counting how many individuals from each deme are in each deme.