Add probability-based stepping to a node in network (#168)

This adds step (i.e., jump, teleport) from one node to another based on probability assigned to an edge.

The code it uses is completely different, so a new function step is now available alongside of the existing travel function.

Node to go to is selected based on edge probability in the new step function. Ignoring of nodes, implemented for the travel function, is not implemented for the step function (or not yet) because that would change the relative probabilities.

Edges have probabilities, but they don't have costs. Probability can be any range (it is relative) but it cannot be negative. Edge probabilities are determined from the input CSV.

Network type is now driven by a single string variable with three values. Using current names network_type = snap, step, travel.

Author: wenzeslaus

include/pops/anthropogenic_kernel.hpp

@@ -59,8 +59,11 @@ std::unique_ptr<KernelInterface<Generator>> create_anthro_kernel(
     else if (anthro_kernel == DispersalKernelType::Network) {
         using Kernel =
             DynamicWrapperKernel<NetworkDispersalKernel<RasterIndex>, Generator>;
+        if (config.network_type == "step")
+            return std::unique_ptr<Kernel>(new Kernel(network));
+        bool snap = config.network_type == "snap" ? true : false;
         return std::unique_ptr<Kernel>(new Kernel(
-            network, config.network_min_distance, config.network_max_distance));
+            network, config.network_min_distance, config.network_max_distance, snap));
     }
     else if (config.deterministic) {
         using Kernel = DynamicWrapperKernel<

include/pops/config.hpp

@@ -71,9 +71,9 @@ class Config
     std::string anthro_kernel_type;
     double anthro_scale{0};
     std::string anthro_direction;
+    std::string network_type;  ///< travel, snap, step
     double network_min_distance{0};
     double network_max_distance{0};
-    bool network_snap{false};  ///< Snap resulting location to closest node in network
     double anthro_kappa{0};
     double shape{1.0};
     // Treatments

include/pops/network.hpp

@@ -64,6 +64,11 @@ class EdgeGeometry : public std::vector<RasterCell>
         return (this->size() - 1) * cost_per_cell_;
     }
 
+    void set_total_cost(double value)
+    {
+        total_cost_ = value;
+    }
+
     /** Get cost per cell for the segment (edge). */
     double cost_per_cell() const
     {
@@ -72,19 +77,27 @@ class EdgeGeometry : public std::vector<RasterCell>
         return cost_per_cell_;
     }
 
-    void set_total_cost(double value)
+    void set_cost_per_cell(double value)
     {
-        total_cost_ = value;
+        cost_per_cell_ = value;
     }
 
-    void set_cost_per_cell(double value)
+    /** Get probability of an edge */
+    double probability() const
     {
-        cost_per_cell_ = value;
+        return probability_;
+    }
+
+    /** Set probability of an edge */
+    void set_probability(double value)
+    {
+        probability_ = value;
     }
 
 private:
     double cost_per_cell_ = 0;
     double total_cost_ = 0;
+    double probability_ = 0;
 };
 
 /** Constant view of a edge geometry (to iterate a segment in either direction)
@@ -177,16 +190,14 @@ class Network
      * @param bbox Bounding box of the raster grid (in real world coordinates)
      * @param ew_res East-west resolution of the raster grid
      * @param ns_res North-south resolution of the raster grid
-     * @param snap Snap result to the closest node
      */
-    Network(BBox<double> bbox, double ew_res, double ns_res, bool snap = false)
+    Network(BBox<double> bbox, double ew_res, double ns_res)
         : bbox_(bbox),
           ew_res_(ew_res),
           ns_res_(ns_res),
           max_row_(0),
           max_col_(0),
-          distance_per_cell_((ew_res + ns_res) / 2),
-          snap_(snap)
+          distance_per_cell_((ew_res + ns_res) / 2)
     {
         std::tie(max_row_, max_col_) = xy_to_row_col(bbox_.east, bbox_.south);
     }
@@ -443,14 +454,25 @@ class Network
      * All previously visited nodes are tracked and, if possible, excluded
      * from further traveling.
      *
+     * The function assumes there is a node at the given *row* and *column*, i.e., that
+     * the decision to call this function was based on the caller knowing there is a
+     * node. If there is no node, an std::invalid_argument exception is thrown.
+     * If there is more than one node at the given *row* and *column*, a random node is
+     * picked and used for traveling.
+     *
+     * If *snap* is true, then results are snapped to the closest node, otherwise
+     * result can be anywhere in between the nodes based on the edge geomerty (segment).
+     *
      * @returns Final row and column pair
      */
     template<typename Generator>
     std::tuple<int, int> travel(
-        RasterIndex row, RasterIndex col, double distance, Generator& generator) const
+        RasterIndex row,
+        RasterIndex col,
+        double distance,
+        Generator& generator,
+        bool snap = false) const
     {
-        // We assume there is a node here, i.e., that we are made decision
-        // to use this kernel knowing there is a node.
         auto node_id = get_random_node_at(row, col, generator);
         std::set<NodeId> visited_nodes;
         while (distance >= 0) {
@@ -472,7 +494,7 @@ class Network
                 distance -= segment.cost();
                 continue;
             }
-            if (snap_) {
+            if (snap) {
                 if (distance < segment.cost() / 2) {
                     // Less than half snaps to the start node.
                     return segment.front();
@@ -487,6 +509,38 @@ class Network
         throw std::invalid_argument("Distance must be greater than or equal to zero");
     }
 
+    /**
+     * Step to a different node in the network from given row and column.
+     *
+     * Returns any node of the nodes connected to the start node possibly based on the
+     * edge probability if probability was assigned to the edges.
+     *
+     * If *num_steps* is greater than 1, multiple steps are perfomed and the last node
+     * is returned. In each node, the probability of picking a specific connection is
+     * either determined by the provided edge probabilities or is equal among the
+     * connections. Consequently, previously visited nodes can be visited again. In
+     * other words, for highly probable connections, most likely next step is a step
+     * back to the starting node (or generally previous node for `num_steps >= 3`).
+     *
+     * The function assumes a node is at the *row*, *col* coordinates, i.e., that this
+     * was either checked beforehand or otherwise ensured. If there is no node, an
+     * std::invalid_argument exception is thrown.
+     * If there is more than one node at the given *row* and *column*, a random node is
+     * picked and used for traveling.
+     *
+     * @returns Destination row and column pair
+     */
+    template<typename Generator>
+    std::tuple<int, int> step(
+        RasterIndex row, RasterIndex col, Generator& generator, int num_steps = 1) const
+    {
+        auto node_id = get_random_node_at(row, col, generator);
+        for (int i = 0; i < num_steps; ++i) {
+            node_id = next_probable_node(node_id, generator);
+        }
+        return get_node_row_col(node_id);
+    }
+
     /**
      * @brief Get all nodes as vector of all ids with their row and column.
      *
@@ -594,7 +648,7 @@ class Network
         stream << "    distance_per_cell: " << distance_per_cell_ << "\n";
         std::set<std::pair<NodeId, NodeId>> edges;
         for (const auto& item : node_matrix_) {
-            for (const auto& node_id : item.second) {
+            for (const auto& node_id : item.second.second) {
                 edges.emplace(item.first, node_id);
             }
         }
@@ -631,7 +685,8 @@ class Network
      * Connection in both directions is stored explicitly (see load_segments() source
      * code).
      */
-    using NodeMatrix = std::map<NodeId, std::vector<NodeId>>;
+    using NodeMatrix =
+        std::map<NodeId, std::pair<std::vector<double>, std::vector<NodeId>>>;
 
     /** Smallest component of edge geometry */
     using Cell = std::pair<RasterIndex, RasterIndex>;
@@ -680,6 +735,48 @@ class Network
         }
     }
 
+    /**
+     * @brief Convert string to probability
+     *
+     * @param text String with probability
+     * @return Probability as number
+     */
+    static double probability_from_text(const std::string& text)
+    {
+        double value{0};
+        try {
+            value = std::stod(text);
+        }
+        catch (const std::invalid_argument& err) {
+            if (string_contains(text, '"') || string_contains(text, '\'')) {
+                throw std::invalid_argument(
+                    std::string("Text for connection probabilty cannot contain quotes "
+                                "(only digits are allowed): ")
+                    + text);
+            }
+            if (text.empty()) {
+                throw std::invalid_argument(
+                    "Text for connection probabilty cannot be an empty string");
+            }
+            else {
+                throw std::invalid_argument(
+                    std::string("Text cannot be converted to connection probabilty "
+                                "(only digits are allowed): ")
+                    + text);
+            }
+        }
+        catch (const std::out_of_range& err) {
+            throw std::out_of_range(
+                std::string("Numerical value too large for connection probabilty: ")
+                + text);
+        }
+        if (value < 0) {
+            throw std::invalid_argument(
+                std::string("Probability needs to be >=0, not ") + text);
+        }
+        return value;
+    }
+
     /**
      * @brief Convert string to cost
      *
@@ -715,8 +812,10 @@ class Network
     }
 
     template<typename InputStream>
-    static bool stream_has_cost_column(InputStream& stream, char delimeter)
+    static std::pair<bool, bool> stream_has_columns(InputStream& stream, char delimeter)
     {
+        bool has_cost{false};
+        bool has_probability{false};
         // Get header to determine what is included.
         auto starting_position = stream.tellg();
         std::string line;
@@ -731,15 +830,29 @@ class Network
                 stream.seekg(starting_position);
                 break;
             }
-            if (label == "cost") {
+            if (label == "probability") {
+                if (has_cost) {
+                    // Detailed check to give a more relevant error message.
+                    throw std::runtime_error(
+                        "The cost column must be after the probability column");
+                }
                 if (column_number != 3) {
                     throw std::runtime_error(
-                        "The cost column must be the third column");
+                        "The probability column must be the third column");
                 }
-                return true;
+                has_probability = true;
+                continue;
+            }
+            if (label == "cost") {
+                if (!(column_number == 3 || column_number == 4)) {
+                    throw std::runtime_error(
+                        "The cost column must be the third or fourth column");
+                }
+                has_cost = true;
+                continue;
             }
         }
-        return false;
+        return {has_cost, has_probability};
     }
 
     /**
@@ -752,7 +865,9 @@ class Network
     void load_segments(InputStream& stream)
     {
         char delimeter{','};
-        bool has_cost = stream_has_cost_column(stream, delimeter);
+        bool has_cost{false};
+        bool has_probability{false};
+        std::tie(has_cost, has_probability) = stream_has_columns(stream, delimeter);
 
         std::string line;
         while (std::getline(stream, line)) {
@@ -775,6 +890,12 @@ class Network
             }
             Segment segment;
 
+            if (has_probability) {
+                std::string probability_text;
+                std::getline(line_stream, probability_text, delimeter);
+                double connection_probability = probability_from_text(probability_text);
+                segment.set_probability(connection_probability);
+            }
             if (has_cost) {
                 std::string cost_text;
                 std::getline(line_stream, cost_text, delimeter);
@@ -842,8 +963,12 @@ class Network
             const auto& segment{node_segment.second};
             nodes_by_row_col_[segment.front()].insert(start_node_id);
             nodes_by_row_col_[segment.back()].insert(end_node_id);
-            node_matrix_[start_node_id].push_back(end_node_id);
-            node_matrix_[end_node_id].push_back(start_node_id);
+            if (has_probability) {
+                node_matrix_[start_node_id].first.push_back(segment.probability());
+                node_matrix_[end_node_id].first.push_back(segment.probability());
+            }
+            node_matrix_[start_node_id].second.push_back(end_node_id);
+            node_matrix_[end_node_id].second.push_back(start_node_id);
         }
     }
 
@@ -881,7 +1006,7 @@ class Network
      */
     const std::vector<NodeId>& nodes_connected_to(NodeId node) const
     {
-        return node_matrix_.at(node);
+        return node_matrix_.at(node).second;
     }
 
     /**
@@ -929,13 +1054,45 @@ class Network
         return pick_random_item(nodes, generator);
     }
 
+    /**
+     * @brief Pick a probable node from the given node.
+     *
+     * If there is more than one edge leading from the given node, a random node is
+     * picked based on edge probabilities. If there are no probabilities assigned, a
+     * random node is picked. If there are no edges leading from the given node, the
+     * node itself is returned.
+     */
+    template<typename Generator>
+    NodeId next_probable_node(NodeId node, Generator& generator) const
+    {
+        // Get all candidate nodes.
+        const auto& record{node_matrix_.at(node)};
+        const auto& probabilities{record.first};
+        const auto& nodes{record.second};
+
+        // Resolve disconnected node and dead end cases.
+        auto num_nodes = nodes.size();
+        if (!num_nodes)
+            return node;
+        else if (num_nodes == 1)
+            return nodes[0];
+
+        // Pick nodes based on edge probabilities if they are available.
+        if (!probabilities.empty()) {
+            std::discrete_distribution<int> dd{
+                probabilities.begin(), probabilities.end()};
+            return nodes.at(dd(generator));
+        }
+        // Pick a connected node with equal edge probabilities.
+        return pick_random_item(nodes, generator);
+    }
+
     BBox<double> bbox_;  ///< Bounding box of the network grid in real world coordinates
     double ew_res_;  ///< East-west resolution of the grid
     double ns_res_;  ///< North-south resolution of the grid
     RasterIndex max_row_;  ///< Maximum row index in the grid
     RasterIndex max_col_;  ///< Maximum column index in the grid
     double distance_per_cell_;  ///< Distance to travel through one cell (cost)
-    bool snap_ = false;
     /** Node IDs stored by row and column (multiple nodes per cell) */
     std::map<std::pair<RasterIndex, RasterIndex>, std::set<NodeId>> nodes_by_row_col_;
     NodeMatrix node_matrix_;  ///< List of node neighbors by node ID (edges)