Merge pull request #42 from Quetzal-framework/develop

doc: local process with niche model and 3 species competition/predation

docs/4-tutorials.md

@@ -1029,7 +1029,7 @@ In this example, we store a time-series vector at each vertex of the graph to mo
 
 Stochastic models introduce randomness and variability into population growth processes, acknowledging the inherent uncertainty in ecological systems. These models incorporate probabilistic elements, such as random fluctuations in birth and death rates, environmental variability, and demographic stochasticity. Stochastic models are particularly useful for capturing the effects of environmental variability and demographic stochasticity on population dynamics, allowing researchers to assess the likelihood of rare events and population extinctions.
 
-We build up on the previous example by adding some stochasticity in the process. We store a time-series vector at each vertex of the graph to monitor the population dynamics. After initializing one vertex with a non-zero value, we iterate through time to update the series. It's important to note that, since we defined no migration along the graph's edges, some time-series remain at zero.
+We build up on the previous example by adding some stochasticity in the process. We store a time-series vector at each vertex of the graph to monitor the population dynamics. After initializing some vertices with a non-zero value, we iterate through time to update the series. It's important to note that, since we defined no migration along the graph's edges, some time-series remain at zero.
 
 **Input**
 
@@ -1045,7 +1045,9 @@ We build up on the previous example by adding some stochasticity in the process.
 
 Environmental niche models focus on how environmental factors influence the distribution and abundance of species within specific habitats. These models integrate ecological niche theory and environmental data to predict species' occurrence probabilities and population densities across spatial gradients. By quantifying the relationships between environmental variables (e.g., temperature, precipitation, habitat quality) and species' ecological requirements, environmental niche models help elucidate how environmental factors shape local population growth and species distributions.
 
-In this example, we store a time-series vector at each vertex of the graph to monitor the population dynamics. After initializing one vertex with a non-zero value, we iterate through time to update the series conditionally to the local contemporary environmental data. It's important to note that, since we defined no migration along the graph's edges, all time-series except one remain at zero.
+In this example, we store a time-series vector at each vertex of the graph to monitor the population dynamics. After initializing some vertices with a non-zero value, we iterate through time to update the series conditionally to the local contemporary environmental data. It's important to note that, since we defined no migration along the graph's edges, some time-series remain at zero.
+
+Since the data structures embedded in the vertices (and edges) of the graph are arbitrary, we could also store the local \f$K\f$ and \f$r\f$ time series alongside the population size. This can be achieved by defining a simple `vertex_info` class that contains three `std::vector` members. We will apply this idea in the following example for simulating a three-species competition model.
 
 **Input**
 
@@ -1061,6 +1063,31 @@ In this example, we store a time-series vector at each vertex of the graph to mo
 
 Species competition models explore how interspecific interactions, such as competition for resources or interactions with predators, influence population dynamics and species coexistence within local communities. These models typically incorporate competition coefficients, representing the strength of interspecific interactions, into population equations to simulate the effects of species interactions on population growth rates and carrying capacities. Species competition models provide insights into the mechanisms driving species coexistence, competitive exclusion, and community dynamics within ecological communities.
 
+A discrete-time model for a three-species predator-prey competition can be complex, but here's a basic version using a system of difference equations. Let's denote the populations of the three species at time \f$t\f$ as 
+\f$𝑃_1(t)\f$, \f$P_2(t)\f$ and \f$P_3(t)\f$. Assume that \f$P_1\f$ is a prey species, \f$P_2\f$ is a predator of \f$P_1\f$ and \f$P_3\f$ is a competing species that preys on \f$P_1\f$ and competes with \f$P_2\f$.
+
+\[
+P_1(t+1) = P_1(t) + r_1 P_1(t) \left(1 - \frac{P_1(t)}{K_1}\right) - a_{12} P_1(t) P_2(t) - a_{13} P_1(t) P_3(t)
+\]
+
+\[
+P_2(t+1) = P_2(t) + r_2 P_2(t) \left(1 - \frac{P_2(t)}{K_2}\right) + b_{21} P_1(t) P_2(t) - c_{23} P_2(t) P_3(t)
+\]
+
+\[
+P_3(t+1) = P_3(t) + r_3 P_3(t) \left(1 - \frac{P_3(t)}{K_3}\right) + b_{31} P_1(t) P_3(t) - c_{32} P_2(t) P_3(t)
+\]
+
+Where:
+- \( r_1, r_2, r_3 \) are the intrinsic growth rates of species 1, 2, and 3, respectively.
+- \( K_1, K_2, K_3 \) are the carrying capacities for species 1, 2, and 3, respectively.
+- \( a_{12} \) is the predation rate of species 2 on species 1.
+- \( a_{13} \) is the predation rate of species 3 on species 1.
+- \( b_{21} \) is the benefit species 2 gets from preying on species 1.
+- \( c_{23} \) is the competition rate between species 2 and 3.
+- \( b_{31} \) is the benefit species 3 gets from preying on species 1.
+- \( c_{32} \) is the competition rate between species 3 and 2.
+
 **Input**
 
 @include{lineno} geography_graph_local_process_4.cpp

example/geography_graph_local_process_3.cpp

@@ -0,0 +1,52 @@
+#include "quetzal/geography.hpp"
+#include <random>
+#include <cmath>
+
+namespace geo = quetzal::geography;
+
+int main()
+{
+    auto file = std::filesystem::current_path() / "data/bio1.tif";
+
+    // The raster have 10 bands that we will assign to 2001 ... 2010.
+    std::vector<int> times(10);
+    std::iota(times.begin(), times.end(), 2001);
+
+    // Landscape construction injecting times information
+    auto land = geo::raster<>::from_file(file, times);
+
+    // Graph construction with vertex and edge info
+    using vertex_info = std::vector<double>;
+    vertex_info N(times.size(), 0);
+
+    using edge_info = geo::no_property;
+    auto graph = geo::from_grid(land, N, edge_info(), geo::connect_fully(), geo::isotropy(), geo::mirror());
+
+    // Let's compare some time series
+    graph[0][0] = 1;
+    graph[1][0] = 1;
+    graph[2][0] = 1;
+
+    // Declare a random number generator
+    std::mt19937 rng{std::random_device{}()};
+
+    // Part of what could happen in a demographic simulator:
+    for( size_t t = 1; t < times.size(); ++t){
+        for( auto x : graph.vertices() ){
+            auto r = std::uniform_real_distribution<double>(1, 5)(rng);
+            auto K = 10. * land.at(x, t-1).value_or(10.);
+            auto P = graph[x][t-1];
+            graph[x][t] = std::max(P + r * P * (1 - P/K), 0.);
+        }
+    }
+
+    // Post treatment
+    std::cout <<  "Time series:\n"; 
+    for(auto x : graph.vertices()){
+         std::cout << "at vertex " << x << " : ";
+        for(auto value : graph[x] ){
+            std::cout << static_cast<int>(value) << " ";
+        }
+        std::cout << "\n";
+    }
+}

example/geography_graph_local_process_4.cpp

@@ -0,0 +1,91 @@
+#include "quetzal/geography.hpp"
+#include <random>
+#include <cmath>
+
+namespace geo = quetzal::geography;
+
+// Time series of the 3 species population sizes to embed along the graph
+struct vertex_info {
+    using time_series_type = std::vector<double>;
+    time_series_type P1, P2, P3;
+    // Constructor to initialize the vectors with a specified length and fill them with zeros
+    vertex_info(int length): P1(length, 0), P2(length, 0), P3(length, 0){}
+    // Needs to be default constructible
+    vertex_info() = default;
+};
+
+// Overload the << operator for the vertex_info struct
+std::ostream& operator<<(std::ostream& os, const vertex_info& vertex) {
+    auto format = [&](const auto & v, const auto& name){
+        os << "\t" << name << " : ";
+        for(const auto& val : v)
+            os << val << " ";
+        os << "\n";
+    };
+    format(vertex.P1, "P1");
+    format(vertex.P2, "P2");
+    format(vertex.P3, "P3");
+    return os;
+}
+
+int main()
+{
+    auto file = std::filesystem::current_path() / "data/bio1.tif";
+
+    // The raster have 10 bands that we will assign to 2001 ... 2010.
+    std::vector<int> times(10);
+    std::iota(times.begin(), times.end(), 2001);
+
+    // Landscape construction injecting times information
+    auto land = geo::raster<>::from_file(file, times);
+
+    // Graph construction with vertex and edge info
+    using edge_info = geo::no_property;
+    auto graph = geo::from_grid(land, vertex_info(times.size()), edge_info(), geo::connect_fully(), geo::isotropy(), geo::mirror());
+
+    // Let's initialize population sizes for a single vertex for demonstration
+    graph[0].P1[0] = 200;
+    graph[0].P2[0] = 1;
+    graph[0].P3[0] = 1;
+
+    // Declare a random number generator
+    std::mt19937 rng{std::random_device{}()};
+
+    auto random_double = [&](double min, double max){ return std::uniform_real_distribution<double>(min, max)(rng);};
+
+    double r1 = random_double(0.0, 1.0);  // Intrinsic growth rate for species 1
+    double r2 = random_double(0.0, 1.0);  // Intrinsic growth rate for species 2
+    double r3 = random_double(0.0, 1.0);  // Intrinsic growth rate for species 3
+
+    double K1 = random_double(50.0, 100.0);  // Carrying capacity for species 1
+    double K2 = random_double(50.0, 100.0);  // Carrying capacity for species 2
+    double K3 = random_double(50.0, 100.0);  // Carrying capacity for species 3
+
+    double a12 = random_double(0.0, 0.1);  // Predation rate of species 2 on species 1
+    double a13 = random_double(0.0, 0.1);  // Predation rate of species 3 on species 1
+
+    double b21 = random_double(0.0, 0.1);  // Benefit species 2 gets from preying on species 1
+    double c23 = random_double(0.0, 0.1);  // Competition rate between species 2 and 3
+
+    double b31 = random_double(0.0, 0.1);  // Benefit species 3 gets from preying on species 1
+    double c32 = random_double(0.0, 0.1);  // Competition rate between species 3 and 2
+
+    // Part of what could happen in a demographic simulator:
+    for( size_t t = 1; t < times.size(); ++t){
+        for( auto x : graph.vertices() ){
+            auto P1 = graph[x].P1[t-1];
+            auto P2 = graph[x].P2[t-1];
+            auto P3 = graph[x].P3[t-1];
+            graph[x].P1[t] = std::max(P1 + r1 * P1 * (1 - P1 / K1) - a12 * P1 * P2 - a13 * P1 * P3, 0.);
+            graph[x].P2[t] = std::max(P2 + r2 * P2 * (1 - P2 / K2) + b21 * P1 * P2 - c23 * P2 * P3, 0.);
+            graph[x].P3[t] = std::max(P3 + r3 * P3 * (1 - P3 / K3) + b31 * P1 * P3 - c32 * P2 * P3, 0.);
+        }
+    }
+
+    // Post treatment
+    std::cout <<  "Time series:\n"; 
+    for(auto x : graph.vertices()){
+         std::cout << "at vertex " << x << " : \n";
+         std::cout << graph[x] << std::endl;
+    }
+}