diff --git a/docs/4-tutorials.md b/docs/4-tutorials.md index b0b4d0a6..887618c7 100644 --- a/docs/4-tutorials.md +++ b/docs/4-tutorials.md @@ -26,17 +26,14 @@ ## Demographic Histories - @subpage demographic_histories_in_quetzal -- @subpage Phylogenetic Trees -- @subpage Phylogenetic Networks -- @subpage spatial_graphs - - @subpage Construction from landscapes - - @subpage Construction from abstract grids - - @subpage SettingVertex/Edge information at construction - - @subpage dispersal_kernels - -- Examples - - @subpage Growth Expressions - - @subpage Memory Management + - @subpage Phylogenetic Trees + - @subpage Phylogenetic Networks + - @subpage spatial_graphs + - @subpage Construction from spatial grids + - @subpage Setting Vertex/Edge information at construction time + - @subpage dispersal_kernels + - @subpage Growth Expressions + - @subpage Memory Management ## Graphs @@ -830,34 +827,47 @@ Another category of models puts an emphasis on the geographic structure of these @tableofcontents @note -The objective of this section is to define the modality of dispersal of a species across a discrete landscape. Because dispersal is so tightly related to graph connectivity and performance, this will require to: -1. transform a `quetzal::geography::landscape` into a `quetzal::geography::graph` selecting the desired level of connectivity -2. update the edges information of the graph to reflect the desired rate of movement. +The objective of this section is to understand the concepts and assumptions required to define the modality of dispersal of a species across a landscape. Because dispersal is so tightly related to graph connectivity and performance, this will require to: +1. transform a `quetzal::geography::landscape` into a `quetzal::geography::graph`. +2. select the desired assumptions to control the level of connectivity in the graph. +2. understand how assumptions define the number of edges in the resulting spatial graph. ## Background The challenge of replicating a demographic process within a discrete landscape can be framed as a process on a spatial graph. In this graph, the vertices represent the cells of the landscape grid, while the edges link the areas where migration occurs. -The quantity of edges in this graph significantly influences both complexity and performance. For instance, in a landscape containing \f$n\f$ cells, considering migration between any two arbitrary locations would result in a complete graph with \f$\frac{n(n-1)}{2}\f$ edges, which could pose computational difficulties. +The quantity of edges in this graph significantly influences both complexity and performance. For instance, in a landscape containing \f$n\f$ cells, considering migration between any two arbitrary locations would result in a complete graph with \f$\frac{n(n-1)}{2}\f$ edges, which could pose computational difficulties. A number of alternative assumptions can be made to control the computational feasibility. -For this reason, constructing a spatial graph first requires to account for the modality of dispersal between the vertices (locations) of the graph. We distinguish at least 3 modalities: +All these choices are independently represented in the code and can be extended (although that may require some efforts). -* 4-neighbors: each cell of the landscape grid is connected only to its cardinal neighbors (N, E, S, W). -* 8-neighbors: each cell of the landscape grid is connected only to its cardinal and intercardinal (NE, SE, SW, NW) neighbors. -* fully connected graph: each cell of the landscape grid is connected to all others: the rate of migration between two cells will be a function of the distance or path of least resistance. +## Assumptions + +### Vicinity + +Constructing a spatial graph first requires to account for the modality of dispersal between the vertices (locations) of the graph. That is, from one source vertex, define which vertices can be targeted: this is the concept of vicinity, or neighborhood. + +We distinguish at least 3 modalities: + +* `connect_4_neighbors`: each cell of the landscape grid is connected only to its cardinal neighbors (N, E, S, W). +* `connect_8_neighbors`: each cell of the landscape grid is connected only to its cardinal (N, E, S, W) and intercardinal (NE, SE, SW, NW) neighbors. +* `connect_fully`: each cell of the landscape grid is connected to all others: the graph is complete. + +### Bounding policy Similarly, assumptions we make about the model's behavior at the bounds of the finite landscape influence the connectivity of the graph. We distinguish at least three bounding modalities: -* mirror: cells at the borders of the landscape are connected to their neighbors, without any further modification. In this case the individuals moving between cells will be *reflected* into the landscape. -* sink: cells at the borders of the landscape will be connected to a sink vertex: individuals that migrate past the boundaries of the landscape *escape* into the world and never come back. -* torus: the 2D landscape becomes a 3D torus world, as vertices on opposite borders are linked to their symmetric vertex: individuals that cross the Northern (resp. Western) border reappear on the Southern (resp. Eastern) border. +* `mirror`: cells at the borders of the landscape are connected to their neighbors, without any further modification. In this case the individuals moving between cells will be *reflected* into the landscape. +* `sink`: cells at the borders of the landscape will be connected to a sink vertex: individuals that migrate past the boundaries of the landscape *escape* into the world and never come back. +* `torus`: the 2D landscape becomes a 3D torus world, as vertices on opposite borders are linked to their symmetric vertex: individuals that cross the Northern (resp. Western) border reappear on the Southern (resp. Eastern) border. + +### Directionality Finally, a last assumption we make about the population process that impacts the graph representation is the directionality of the process: -* isotropic migration means that the edge \f$( u,v)\f$ and the edge \f$(v,u)\f$ are one and the same. -* anisotropic migration will distinguish between these two edges. +* `isotropy`: the migration process is independent of the direction of movement. Moving from vertex \f$u\f$ to \f$v\f$ follows the same rules as moving from \f$v\f$ to \f$u\f$: the edge \f$( u,v)\f$ is the same as \f$(v,u)\f$. +* `anisotropy` migration will distinguish between \f$( u,v)\f$ and \f$(v,u)\f$, doubling the number of edges in the graph. -All these choices are independently represented in the code and can be extended (although that may require some efforts): +## Example **Input** @@ -869,6 +879,70 @@ All these choices are independently represented in the code and can be extended --- +[//]: # (----------------------------------------------------------------------) +@page spatial_graph_construction Construction +@tableofcontents + +## Background + +Constructing a spatial graph from a spatial grid, also known as a raster, is a fundamental step in spatial analysis and modeling. Spatial graphs provide a powerful representation of relationships and connectivity between different locations within a geographic area. This process involves transforming the discrete and gridded structure of a spatial grid into a network of vertices and edges, where vertices represent spatial entities (such as cells or locations), and edges capture the spatial relationships and connections between these entities. + +There are different ways to build such graph. + +## Graph from raster + +**Input** + +@include{lineno} geography_graph_2.cpp + +**Output** + +@include{lineno} geography_graph_2.txt + +## Graph from landscape + +**Input** + +@include{lineno} geography_graph_3.cpp + +**Output** + +@include{lineno} geography_graph_3.txt + +## Graph from abstract grid + +**Input** + +@include{lineno} geography_graph_4.cpp + +**Output** + +@include{lineno} geography_graph_4.txt + +## From scratch + +### Sparse Graph + +**Input** + +@include{lineno} geography_graph_5.cpp + +**Output** + +@include{lineno} geography_graph_5.txt + +### Dense Graph + +**Input** + +@include{lineno} geography_graph_6.cpp + +**Output** + +@include{lineno} geography_graph_6.txt + +--- + [//]: # (----------------------------------------------------------------------) @page dispersal_kernels Dispersal Kernels @tableofcontents diff --git a/example/geography_graph_2.cpp b/example/geography_graph_2.cpp new file mode 100644 index 00000000..7975e9aa --- /dev/null +++ b/example/geography_graph_2.cpp @@ -0,0 +1,27 @@ +#include "quetzal/geography.hpp" + +#include +#include + +namespace geo = quetzal::geography; + +int main() +{ + auto file = std::filesystem::current_path() / "data/bio1.tif"; + + // The raster has 10 bands that we will assign to 2001 ... 2010. + std::vector times(10); + std::iota(times.begin(), times.end(), 2001); + + // Initialize the landscape + auto land = geo::raster<>::from_file(file, times); + + // Our graph will not store any useful information + using vertex_info = geo::no_property; + using edge_info = geo::no_property; + + auto graph = geo::from_grid(land, vertex_info(), edge_info(), geo::connect_fully(), geo::isotropy(), geo::mirror()); + + std::cout << "Graph has " << graph.num_vertices() << " vertices, " << graph.num_edges() << " edges." << std::endl; + +} diff --git a/example/geography_graph_3.cpp b/example/geography_graph_3.cpp new file mode 100644 index 00000000..1b7d353a --- /dev/null +++ b/example/geography_graph_3.cpp @@ -0,0 +1,29 @@ +#include "quetzal/geography.hpp" + +#include +#include + +namespace geo = quetzal::geography; + +int main() +{ + auto file1 = std::filesystem::current_path() / "data/bio1.tif"; + auto file2 = std::filesystem::current_path() / "data/bio12.tif"; + + // The raster have 10 bands that we will assign to 2001 ... 2010. + std::vector times(10); + std::iota(times.begin(), times.end(), 2001); + + // Initialize the landscape: for each var a key and a file, for all a time series. + using landscape_type = geo::landscape<>; + auto land = landscape_type::from_files({{"bio1", file1}, {"bio12", file2}}, times); + + // Our graph will not store any useful information + using vertex_info = geo::no_property; + using edge_info = geo::no_property; + + auto graph = geo::from_grid(land, vertex_info(), edge_info(), geo::connect_fully(), geo::isotropy(), geo::mirror()); + + std::cout << "Graph has " << graph.num_vertices() << " vertices, " << graph.num_edges() << " edges." << std::endl; + +} diff --git a/example/geography_graph_4.cpp b/example/geography_graph_4.cpp new file mode 100644 index 00000000..949affe3 --- /dev/null +++ b/example/geography_graph_4.cpp @@ -0,0 +1,27 @@ +#include "quetzal/geography.hpp" + +namespace geo = quetzal::geography; + +// User-defined concept of spatial grid +struct MySpatialGrid +{ + // it is just required to define these two functions + constexpr int width() const { return 300; } + constexpr int height() const { return 100; } +}; + +int main() +{ + + static_assert(geo::two_dimensional, "MySpatialGrid does not fulfill the two_dimensional requirements"); + MySpatialGrid land; + + // Our graph will not store any useful information + using vertex_info = geo::no_property; + using edge_info = geo::no_property; + + auto graph = geo::from_grid(land, vertex_info(), edge_info(), geo::connect_fully(), geo::isotropy(), geo::mirror()); + + std::cout << "Graph has " << graph.num_vertices() << " vertices, " << graph.num_edges() << " edges." << std::endl; + +} diff --git a/example/geography_graph_5.cpp b/example/geography_graph_5.cpp new file mode 100644 index 00000000..f7e554b7 --- /dev/null +++ b/example/geography_graph_5.cpp @@ -0,0 +1,20 @@ +#include "quetzal/geography.hpp" + +namespace geo = quetzal::geography; + +int main() +{ + using vertex_info = geo::no_property; + using edge_info = geo::no_property; + using graph_type = geo::graph; + + int num_vertices = 100; + graph_type graph( num_vertices ); + + graph.add_edge( 0, 1 ); + graph.add_edge( 0, 2 ); + graph.add_edge( 0, 3 ); + + std::cout << "Graph has " << graph.num_vertices() << " vertices, " << graph.num_edges() << " edges." << std::endl; + +} diff --git a/example/geography_graph_6.cpp b/example/geography_graph_6.cpp new file mode 100644 index 00000000..65d59559 --- /dev/null +++ b/example/geography_graph_6.cpp @@ -0,0 +1,20 @@ +#include "quetzal/geography.hpp" + +namespace geo = quetzal::geography; + +int main() +{ + using vertex_info = geo::no_property; + using edge_info = geo::no_property; + using graph_type = geo::graph; + + int num_vertices = 100; + graph_type graph( num_vertices ); + + graph.add_edge( 0, 1 ); + graph.add_edge( 0, 2 ); + graph.add_edge( 0, 3 ); + + std::cout << "Graph has " << graph.num_vertices() << " vertices, " << graph.num_edges() << " edges." << std::endl; + +} diff --git a/src/include/quetzal/geography/graph/detail/concepts.hpp b/src/include/quetzal/geography/graph/detail/concepts.hpp index 75fcd828..c26160da 100644 --- a/src/include/quetzal/geography/graph/detail/concepts.hpp +++ b/src/include/quetzal/geography/graph/detail/concepts.hpp @@ -12,36 +12,20 @@ #include "directionality.hpp" +#include +#include + #include #include namespace quetzal::geography { -template -concept two_dimensional = requires(T a) { - requires std::convertible_to; - requires std::convertible_to; - }; - namespace detail { template concept is_any_of = (std::same_as || ...); -} - -template -concept directional = detail::is_any_of; -template -concept bounding = two_dimensional and requires(T t, typename G::vertex_descriptor s, G &graph, S const &grid) { - { - t(s, graph, grid) - } -> std::same_as; - }; -namespace detail -{ - template struct edge_construction { @@ -61,4 +45,50 @@ struct edge_construction }; } + +/// @brief Concept to represent the directionality of edges in a graph +template +concept directional = detail::is_any_of; + +static_assert ( directional ); +static_assert ( directional ); + +/// @brief Represents no information carried by vertices or edges of a graph. +using no_property = boost::no_property; + +/// @brief Tag for sparse graph representation +struct sparse +{ + template + using rebind = boost::adjacency_list; +}; + +/// @brief Tag for sparse graph representation +struct dense +{ + template + using rebind = boost::adjacency_matrix>; +}; + +/// @brief Concept to represent the connectedness of a graph +template +concept connectedness = detail::is_any_of; + +static_assert ( connectedness ); +static_assert ( connectedness ); + +/// @brief Concept to represent a 2D spatial grid +template +concept two_dimensional = requires(T a) { + requires std::convertible_to; + requires std::convertible_to; + }; + +/// @brief Concept to represent the bounding policy of a spatial graph +template +concept bounding = two_dimensional and requires(T t, typename G::vertex_descriptor s, G &graph, S const &grid) { + { + t(s, graph, grid) + } -> std::same_as; + }; } // namespace quetzal::geography \ No newline at end of file diff --git a/src/include/quetzal/geography/graph/detail/vicinity.hpp b/src/include/quetzal/geography/graph/detail/vicinity.hpp index a0b21a84..2f1196f5 100644 --- a/src/include/quetzal/geography/graph/detail/vicinity.hpp +++ b/src/include/quetzal/geography/graph/detail/vicinity.hpp @@ -51,8 +51,7 @@ namespace detail struct connect_fully { - template - using graph_type = graph; + using connectedness = dense; template B> void connect(G &graph, S const &grid, B bound_policy) const @@ -61,8 +60,6 @@ struct connect_fully using vertex_property = typename G::vertex_property; using edge_property = typename G::edge_property; - static_assert(std::same_as>); - int width = grid.width(); int height = grid.height(); int num_land_vertices = width * height; @@ -87,8 +84,7 @@ struct connect_fully class connect_4_neighbors { public: - template - using graph_type = graph; + using connectedness = sparse; private: template void connect(auto s, auto t, G &graph, auto const &grid) const @@ -175,8 +171,6 @@ class connect_4_neighbors using vertex_property = typename G::vertex_property; using edge_property = typename G::edge_property; - static_assert(std::same_as>); - int width = grid.width(); int height = grid.height(); int num_land_vertices = width * height; @@ -237,8 +231,7 @@ class connect_4_neighbors class connect_8_neighbors { public: - template - using graph_type = graph; + using connectedness = sparse; private: template void connect(auto s, auto t, G &graph, auto const &grid) const @@ -247,8 +240,6 @@ class connect_8_neighbors using vertex_property = typename G::vertex_property; using edge_property = typename G::edge_property; - static_assert(std::same_as>); - int width = grid.width(); int height = grid.height(); int num_land_vertices = width * height; @@ -353,7 +344,6 @@ class connect_8_neighbors using directed_category = typename G::directed_category; using vertex_property = typename G::vertex_property; using edge_property = typename G::edge_property; - static_assert(std::same_as>); int width = grid.width(); int height = grid.height(); diff --git a/src/include/quetzal/geography/graph/from_grid.hpp b/src/include/quetzal/geography/graph/from_grid.hpp index ce16df59..f861afd3 100644 --- a/src/include/quetzal/geography/graph/from_grid.hpp +++ b/src/include/quetzal/geography/graph/from_grid.hpp @@ -30,7 +30,8 @@ template ; + using graph_type = quetzal::geography::graph; + // using graph_type = typename Vicinity::graph_type; // Graph size has to be static in dense graphs :/ graph_type graph( grid.width() * grid.height() + bounding_policy.num_extra_vertices() ) ; vicinity.connect(graph, grid, bounding_policy); diff --git a/src/include/quetzal/geography/graph/graph.hpp b/src/include/quetzal/geography/graph/graph.hpp index ba539432..37ed95cb 100644 --- a/src/include/quetzal/geography/graph/graph.hpp +++ b/src/include/quetzal/geography/graph/graph.hpp @@ -13,8 +13,6 @@ #include #include -#include -#include #include #include #include @@ -31,27 +29,18 @@ namespace quetzal::geography { -using no_property = boost::no_property; namespace detail { -template