Merge pull request #34 from Quetzal-framework/develop

doc: add examples and documentation for spatial graphs

docs/4-tutorials.md

@@ -29,9 +29,9 @@
   - @subpage Phylogenetic Trees
   - @subpage Phylogenetic Networks
   - @subpage spatial_graphs
-    - @subpage Construction from spatial grids
-    - @subpage Setting Vertex/Edge information at construction time
+    - @subpage spatial_graph_construction
     - @subpage dispersal_kernels
+    - @subpage Setting Vertex/Edge information at construction time
     - @subpage Growth Expressions
     - @subpage Memory Management
 
@@ -891,6 +891,10 @@ There are different ways to build such graph.
 
 ## Graph from raster
 
+The provided code snippet contains the logic for constructing a spatial graph from a `quetzal::geography::raster`. 
+
+This involves identifying the cells of the raster to map `location_type` values to the vertices, and defining edges based on the desired connectivity.
+
 **Input**
 
 @include{lineno} geography_graph_2.cpp
@@ -901,6 +905,10 @@ There are different ways to build such graph.
 
 ## Graph from landscape
 
+The provided code snippet contains the logic for constructing a spatial graph from a `quetzal::geography::landscape`. 
+
+This involves identifying the cells of the landscape to map `location_type` values to the vertices, and defining edges based on the desired connectivity.
+
 **Input**
 
 @include{lineno} geography_graph_3.cpp
@@ -911,6 +919,10 @@ There are different ways to build such graph.
 
 ## Graph from abstract grid
 
+The provided code snippet contains the logic for constructing a spatial graph from a user-defined abstract grid.
+
+This involves identifying the height and width of the space to define linearly-indexed vertices, and defining edges based on the desired connectivity.
+
 **Input**
 
 @include{lineno} geography_graph_4.cpp
@@ -921,6 +933,9 @@ There are different ways to build such graph.
 
 ## From scratch
 
+The provided code snippet contains the logic for creating a sparse versus dense graph from scratch. This involves 
+dynamically adding vertices and edges based on user-defined conditions, resulting in a graph that represents the desired spatial relationships.
+
 ### Sparse Graph
 
 **Input**
@@ -944,7 +959,7 @@ There are different ways to build such graph.
 ---
 
 [//]: # (----------------------------------------------------------------------)
-@page dispersal_kernels Dispersal Kernels
+@page dispersal_kernels Distance-based dispersal Kernels
 @tableofcontents
 
 ## Distance-based kernel
@@ -969,5 +984,16 @@ Quetzal incorporates various kernel types available in the `quetzal::demography:
 
 ---
 
-## Resistance-based Dispersal
+[//]: # (----------------------------------------------------------------------)
+@page spatial_graph_information Embedding Vertex and/or Edge information
+@tableofcontents
+
+**Input**
+
+@include{lineno} geography_graph_7.cpp
+
+**Output**
+
+@include{lineno} geography_graph_7.txt
 
+---

example/expressive_3.cpp

@@ -22,7 +22,7 @@ int main()
     std::iota(times.begin(), times.end(), 2001);
 
     // Makes sure the rasters are aligned
-    auto landscape = landscape_type::from_files({{"suit", file1}, {"DEM", file2}}, times);
+    auto landscape = landscape_type::from_file({{"suit", file1}, {"DEM", file2}}, times);
 
     // lightweight functor returning empty optionals where NA
     auto s_view = landscape["suit"].to_view();

example/expressive_4.cpp

@@ -22,7 +22,7 @@ int main()
     std::iota(times.begin(), times.end(), 2001);
 
     // Makes sure the rasters are aligned
-    auto landscape = landscape_type::from_files({{"suit", file1}, {"DEM", file2}}, times);
+    auto landscape = landscape_type::from_file({{"suit", file1}, {"DEM", file2}}, times);
 
     // lightweight functor returning empty optionals where NA
     auto s_view = landscape["suit"].to_view();

example/geography_graph_1.cpp

@@ -17,7 +17,7 @@ int main()
 
     // 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);
+    auto land = landscape_type::from_file({{"bio1", file1}, {"bio12", file2}}, times);
 
     // Our graph will not store any useful information
     using vertex_info = geo::no_property;

example/geography_graph_3.cpp

@@ -16,7 +16,7 @@ int main()
 
     // 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);
+    auto land = landscape_type::from_file({{"bio1", file1}, {"bio12", file2}}, times);
 
     // Our graph will not store any useful information
     using vertex_info = geo::no_property;

example/geography_graph_7.cpp

@@ -0,0 +1,37 @@
+#include "quetzal/geography.hpp"
+
+namespace geo = quetzal::geography;
+
+struct MyVertexInfo {
+    std::string label = "NA";
+    std::vector<double> pop_data_chunk;
+};
+
+struct MyEdgeInfo {
+    double distance;
+    MyEdgeInfo() = default;
+    // required by from_grid
+    MyEdgeInfo(auto u, auto v, auto const& land){ 
+        distance = land.to_lonlat(u).great_circle_distance_to( land.to_lonlat(v) );
+    }
+};
+
+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);
+
+    using landscape_type = geo::raster<>;
+    auto land = geo::raster<>::from_file(file, times);
+    auto graph = geo::from_grid(land, MyVertexInfo(), MyEdgeInfo(), geo::connect_fully(), geo::isotropy(), geo::mirror());
+
+    std::cout << "Graph has " << graph.num_vertices() << " vertices, " << graph.num_edges() << " edges." << std::endl;
+
+    for( auto const& e : graph.edges() ){
+        std::cout << graph.source(e) << " <-> " << graph.target(e) << " : " << graph[e].distance << std::endl;
+    }
+
+}

example/geography_graph_complete_1.cpp

@@ -19,7 +19,7 @@ int main()
 
     // 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);
+    auto land = landscape_type::from_file({{"bio1", file1}, {"bio12", file2}}, times);
 
     // Our graph will not store any useful information
     using vertex_info = geo::no_property;
@@ -32,7 +32,7 @@ int main()
     // Define a helper function to compute distance on Earth between two location_descriptors
     auto sphere_distance = [&](auto point1, auto point2)
     {
-        return land.to_latlon(point1)->great_circle_distance_to(land.to_latlon(point2).value());
+        return land.to_latlon(point1).great_circle_distance_to(land.to_latlon(point2));
     };
 
     // Define the type of distance-based kernel to use

example/geography_landscape_1.cpp

@@ -15,7 +15,7 @@ int main()
 
     // Initialize the landscape: for each var a key and a file, for all a time series.
     using landscape_type = quetzal::geography::landscape<>;
-    auto env = quetzal::geography::landscape<>::from_files({{"bio1", file1}, {"bio12", file2}}, times);
+    auto env = quetzal::geography::landscape<>::from_file({{"bio1", file1}, {"bio12", file2}}, times);
     std::cout << env << std::endl;
 
     // We indeed recorded 2 variables: bio1 and bio12

example/geography_raster_1.cpp

@@ -70,8 +70,8 @@ int main()
     assert(point_1 == point_2);
 
     // Defines a lambda expression for checking extent
-    auto check = [&bio1](auto x) {
-        std::cout << "Point " << x << " is " << (bio1.contains(x) ? " " : "not ") << "in bio1 extent" << std::endl;
+    auto check = [&](auto x) {
+        std::cout << "Point " << x << " is" << (bio1.contains(x) ? " " : " not ") << "in bio1 extent" << std::endl;
     };
 
     check(point_1);

src/include/quetzal/geography/graph/detail/bound_policy.hpp

@@ -56,7 +56,7 @@ class sink
         int sink = num_land_vertices;
         assert( graph.num_vertices() == sink + 1 );
         if ( s < sink )
-            detail::edge_construction<edge_property>::delegate(s, sink, graph); // one-way ticket :(
+            detail::edge_construction<edge_property>::delegate(s, sink, graph, grid); // one-way ticket :(
     }
 };
 
@@ -100,7 +100,7 @@ class torus
         {
             symmetricIndex = s - width + 1;
         }
-        detail::edge_construction<edge_property>::delegate(s, symmetricIndex, graph);
+        detail::edge_construction<edge_property>::delegate(s, symmetricIndex, graph, grid);
     }
 };
 

src/include/quetzal/geography/graph/detail/concepts.hpp

@@ -29,16 +29,16 @@ concept is_any_of = (std::same_as<T, U> || ...);
 template<typename T>
 struct edge_construction
 {
-    static inline constexpr void delegate(auto s, auto t, auto& graph)
+    static inline constexpr void delegate(auto s, auto t, auto& graph, auto const& grid)
     {
-        graph.add_edge(s, t, T(s,t));
+        graph.add_edge(s, t, T(s,t, grid));
     }
 };
 
 template<>
 struct edge_construction<boost::no_property>
 {
-    static inline constexpr void delegate(auto s, auto t, auto& graph)
+    static inline constexpr void delegate(auto s, auto t, auto& graph, [[maybe_unused]] auto const& grid)
     {
         graph.add_edge(s, t);
     }

src/include/quetzal/geography/graph/detail/vicinity.hpp

@@ -70,9 +70,9 @@ struct connect_fully
         {
             for (auto t = s + 1; t < num_land_vertices; ++t)
             {
-                detail::edge_construction<edge_property>::delegate(s, t, graph); // (s -> v) == (s <- v) if undirected
+                detail::edge_construction<edge_property>::delegate(s, t, graph, grid); // (s -> v) == (s <- v) if undirected
                 if constexpr (std::same_as<directed_category, anisotropy>)
-                    detail::edge_construction<edge_property>::delegate(t, s, graph); // (s -> v) != (s <- v) because directed
+                    detail::edge_construction<edge_property>::delegate(t, s, graph, grid); // (s -> v) != (s <- v) because directed
             }
 
             if( detail::is_on_border(s, width, height ) ) 
@@ -90,9 +90,9 @@ class connect_4_neighbors
     template <typename G> void connect(auto s, auto t, G &graph, auto const &grid) const
     {
         using edge_property = typename G::edge_property;
-        detail::edge_construction<edge_property>::delegate(s, t, graph);
+        detail::edge_construction<edge_property>::delegate(s, t, graph, grid);
         if constexpr (std::same_as<typename G::directed_category, anisotropy>)
-            detail::edge_construction<edge_property>::delegate(t, s, graph);
+            detail::edge_construction<edge_property>::delegate(t, s, graph, grid);
     }
 
     void connect_top_left_corner(auto s, auto &graph, auto const &grid, auto const &bound_policy) const
@@ -240,6 +240,7 @@ class connect_8_neighbors
         using vertex_property = typename G::vertex_property;
         using edge_property = typename G::edge_property;
 
+
         int width = grid.width();
         int height = grid.height();
         int num_land_vertices = width * height;
@@ -248,9 +249,9 @@ class connect_8_neighbors
         assert( t >= 0 );
         assert( t < num_land_vertices);
 
-        detail::edge_construction<edge_property>::delegate(s, t, graph);
+        detail::edge_construction<edge_property>::delegate(s, t, graph, grid);
         if constexpr (std::same_as<directed_category, anisotropy>)
-            detail::edge_construction<edge_property>::delegate(t, s, graph);
+            detail::edge_construction<edge_property>::delegate(t, s, graph, grid);
     }
 
     void connect_top_left_corner(auto s, auto &graph, auto const &grid, auto const &bound_policy) const