network_kernel.hpp

Go to the documentation of this file.

/*
 * PoPS model - network dispersal kernel
 *
 * Copyright (C) 2020-2021 by the authors.
 *
 * Authors: Vaclav Petras (wenzeslaus gmail com)
 *
 * The code contained herein is licensed under the GNU General Public
 * License. You may obtain a copy of the GNU General Public License
 * Version 2 or later at the following locations:
 *
 * http://www.opensource.org/licenses/gpl-license.html
 * http://www.gnu.org/copyleft/gpl.html
 */
 
#ifndef POPS_NETWORK_KERNEL_HPP
#define POPS_NETWORK_KERNEL_HPP
 
#include "kernel_types.hpp"
#include "raster.hpp"
#include "utils.hpp"
 
#include <algorithm>
#include <set>
#include <random>
#include <string>
#include <sstream>
#include <map>
#include <cmath>
#include <vector>
 
namespace pops {
 
template<typename RasterIndex>
class Network
{
public:
    using NodeId = int;  
    using Statistics = std::map<std::string, int>;  
 
    Network(BBox<double> bbox, double ew_res, double ns_res, double speed)
        : bbox_(bbox),
          ew_res_(ew_res),
          ns_res_(ns_res),
          cell_travel_time_(((ew_res + ns_res) / 2) / speed)
    {}
 
    static Network null_network()
    {
        return Network(BBox<double>(), 0, 0, 0);
    }
 
    std::pair<RasterIndex, RasterIndex> xy_to_row_col(double x, double y) const
    {
        double col = (x - bbox_.west) / ew_res_;
        double row = (bbox_.north - y) / ns_res_;
        return {std::floor(row), std::floor(col)};
    }
 
    std::pair<RasterIndex, RasterIndex>
    xy_to_row_col(const std::string& x, const std::string& y) const
    {
        return xy_to_row_col(std::stod(x), std::stod(y));
    }
 
    bool out_of_bbox(double x, double y) const
    {
        return x > bbox_.east || x < bbox_.west || y > bbox_.north || y < bbox_.south;
    }
 
    template<typename InputStream>
    void load(
        InputStream& node_stream, InputStream& segment_stream, bool allow_empty = false)
    {
        std::set<NodeId> node_ids;
        load_nodes(node_stream, node_ids);
        if (node_ids.empty()) {
            if (allow_empty)
                return;
            else
                throw std::runtime_error("Network: No nodes within the extend");
        }
        load_segments(segment_stream, node_ids);
    }
 
    const std::set<NodeId>& get_nodes_at(RasterIndex row, RasterIndex col) const
    {
        auto it = nodes_by_row_col_.find(std::make_pair(row, col));
        if (it != nodes_by_row_col_.end())
            return it->second;
        // If not found, return an empty set (const reference to local static).
        static std::set<NodeId> empty;
        return empty;
    }
 
    template<typename Generator>
    NodeId
    get_random_node_at(RasterIndex row, RasterIndex col, Generator& generator) const
    {
        const auto& nodes = get_nodes_at(row, col);
        auto num_nodes = nodes.size();
        if (num_nodes == 1) {
            return *nodes.begin();
        }
        else if (num_nodes > 1) {
            return pick_random_node(nodes, generator);
        }
        else {
            throw std::invalid_argument("No nodes at a given row and column");
        }
    }
 
    bool has_node_at(RasterIndex row, RasterIndex col) const
    {
        // Replace by contains in C++20.
        // return nodes_by_row_col_.count(std::make_pair(row, col)) > 0;
        auto it = nodes_by_row_col_.find(std::make_pair(row, col));
        if (it == nodes_by_row_col_.end())
            return false;
        return true;
    }
 
    std::pair<RasterIndex, RasterIndex> get_node_row_col(NodeId node) const
    {
        for (const auto& item : nodes_by_row_col_) {
            if (container_contains(item.second, node))
                return item.first;
        }
        throw std::invalid_argument("No node with a given id");
    }
 
    template<typename Generator>
    std::tuple<int, int>
    travel(RasterIndex row, RasterIndex col, double time, Generator& generator) 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 (time >= 0) {
            auto next_node_id = next_node(node_id, visited_nodes, generator);
            visited_nodes.insert(node_id);
            // If there is no segment from the node, return the start cell.
            if (next_node_id == node_id)
                return std::make_tuple(row, col);
            auto segment = get_segment(node_id, next_node_id);
            // nodes may need special handling
            for (const auto& cell : segment) {
                time -= cell_travel_time_;
                if (time <= 0) {
                    return cell;
                    // Given the while condition, this subsequently ends the while loop
                    // as well.
                    // break;
                }
            }
            node_id = next_node_id;
        }
        throw std::invalid_argument("Time must be greater than or equal to zero");
    }
 
    std::vector<std::pair<NodeId, std::pair<RasterIndex, RasterIndex>>>
    get_all_nodes() const
    {
        std::vector<std::pair<NodeId, std::pair<RasterIndex, RasterIndex>>> nodes;
        nodes.reserve(nodes_by_row_col_.size());  // It will be at least this big.
        for (const auto& item : nodes_by_row_col_) {
            RasterIndex row = item.first.first;
            RasterIndex col = item.first.second;
            for (const auto node_id : item.second) {
                nodes.emplace_back(node_id, std::make_pair(row, col));
            }
        }
        return nodes;
    }
 
    Statistics collect_statistics() const
    {
        std::map<std::string, int> stats;
        std::set<NodeId> node_ids;
        for (const auto& item : nodes_by_row_col_) {
            for (const auto node_id : item.second) {
                node_ids.insert(node_id);
            }
        }
        stats["num_nodes"] = node_ids.size();
        // We store segements in both directions, so each segment is stored twice.
        stats["num_segments"] = node_matrix_.size() / 2;
        std::set<NodeId> nodes_with_segments;
        for (const auto& item : node_matrix_) {
            nodes_with_segments.insert(item.first);
        }
        stats["num_nodes_with_segments"] = nodes_with_segments.size();
        // TODO: Only count the nodes here. Output them in the dump network output.
        int num_standalone_nodes = 0;
        for (NodeId node_id : node_ids) {
            if (nodes_with_segments.find(node_id) == nodes_with_segments.end()) {
                stats["standalone_node_" + std::to_string(++num_standalone_nodes)] =
                    node_id;
            }
        }
        stats["num_standalone_nodes"] = num_standalone_nodes;
        // TODO: This can be faster. The list is a ordered set and we iterate over
        // elements above.
        const auto minmax_node_id = minmax_element(node_ids.begin(), node_ids.end());
        stats["min_node_id"] = *minmax_node_id.first;
        stats["max_node_id"] = *minmax_node_id.second;
        return stats;
    }
 
    template<typename OutputStream>
    void dump_yaml(OutputStream& stream) const
    {
        stream << "network:\n";
        stream << "  statistics:\n";
        auto stats = collect_statistics();
        for (const auto& item : stats) {
            stream << "    " << item.first << ": " << item.second << "\n";
        }
        stream << "  extent:\n";
        stream << "    north: " << bbox_.north << "\n";
        stream << "    south: " << bbox_.south << "\n";
        stream << "    east: " << bbox_.east << "\n";
        stream << "    west: " << bbox_.west << "\n";
        stream << "  resolution:\n";
        stream << "    ns: " << ns_res_ << "\n";
        stream << "    ew: " << ew_res_ << "\n";
        stream << "  raster:\n";
        RasterIndex min_row;
        RasterIndex min_col;
        RasterIndex max_row;
        RasterIndex max_col;
        std::tie(min_row, min_col) =
            xy_to_row_col(bbox_.west + ew_res_ / 2, bbox_.north - ns_res_ / 2);
        std::tie(max_row, max_col) =
            xy_to_row_col(bbox_.east - ew_res_ / 2, bbox_.south + ew_res_ / 2);
        stream << "    min_row: " << min_row << "\n";
        stream << "    min_col: " << min_col << "\n";
        stream << "    max_row: " << max_row << "\n";
        stream << "    max_col: " << max_col << "\n";
        stream << "    num_rows: " << max_row - min_row + 1 << "\n";
        stream << "    num_cols: " << max_col - min_col + 1 << "\n";
        stream << "  cost:\n";
        stream << "    cell_travel_time: " << cell_travel_time_ << "\n";
        std::set<std::pair<NodeId, NodeId>> edges;
        for (const auto& item : node_matrix_) {
            for (const auto& node_id : item.second) {
                edges.emplace(item.first, node_id);
            }
        }
        stream << "  edges:\n";
        for (const auto& item : edges) {
            stream << "    - [" << item.first << ", " << item.second << "]\n";
        }
        stream << "  nodes:\n";
        for (const auto& item : nodes_by_row_col_) {
            auto row = item.first.first;
            auto col = item.first.second;
            for (const auto& node_id : item.second) {
                stream << "    - id: " << node_id << "\n";
                stream << "      row: " << row << "\n";
                stream << "      col: " << col << "\n";
            }
        }
        stream << "  segments:\n";
        for (const auto& item : segments_by_nodes_) {
            stream << "    - start_node: " << item.first.first << "\n";
            stream << "      end_node: " << item.first.second << "\n";
            stream << "      cells: [";
            for (const auto& cell : item.second) {
                stream << "[" << cell.first << ", " << cell.second << "], ";
            }
            stream << "]\n";
        }
    }
 
protected:
    using NodeMatrix = std::map<NodeId, std::vector<NodeId>>;
    using Segment = std::vector<std::pair<RasterIndex, RasterIndex>>;
    using SegmentsByNodes = std::map<std::pair<NodeId, NodeId>, Segment>;
 
    static NodeId node_id_from_text(const std::string& text)
    {
        return std::stoi(text);
    }
 
    template<typename Container>
    class ConstEitherWayIterator
    {
    public:
        ConstEitherWayIterator(typename Container::const_iterator it)
            : is_forward_(true), forward_it_(it)
        {}
        ConstEitherWayIterator(typename Container::const_reverse_iterator it)
            : is_forward_(false), reverse_it_(it)
        {}
        ConstEitherWayIterator& operator++()
        {
            if (is_forward_)
                ++forward_it_;
            else
                ++reverse_it_;
            return *this;
        }
        const typename Container::value_type& operator*() const
        {
            if (is_forward_)
                return *forward_it_;
            return *reverse_it_;
        }
        bool operator!=(const ConstEitherWayIterator& other)
        {
            if (is_forward_)
                return forward_it_ != other.forward_it_;
            return reverse_it_ != other.reverse_it_;
        }
 
    private:
        // Can be improved by std::optional or any in C++17.
        bool is_forward_;
        typename Container::const_iterator forward_it_;
        typename Container::const_reverse_iterator reverse_it_;
    };
 
    class SegmentView
    {
    public:
        SegmentView(
            typename Segment::const_iterator first,
            typename Segment::const_iterator last)
            : begin_(first), end_(last)
        {}
        SegmentView(
            typename Segment::const_reverse_iterator first,
            typename Segment::const_reverse_iterator last)
            : begin_(first), end_(last)
        {}
        ConstEitherWayIterator<Segment> begin() const
        {
            return begin_;
        }
        ConstEitherWayIterator<Segment> end() const
        {
            return end_;
        }
 
    private:
        ConstEitherWayIterator<Segment> begin_;
        ConstEitherWayIterator<Segment> end_;
    };
 
    template<typename InputStream>
    void load_nodes(InputStream& stream, std::set<NodeId>& node_ids)
    {
        std::string line;
        while (std::getline(stream, line)) {
            std::istringstream line_stream{line};
            char delimeter{','};
            std::string node_text;
            std::getline(line_stream, node_text, delimeter);
            std::string x_coord_text;
            std::getline(line_stream, x_coord_text, delimeter);
            std::string y_coord_text;
            std::getline(line_stream, y_coord_text, delimeter);
            RasterIndex row;
            RasterIndex col;
            double x = std::stod(x_coord_text);
            double y = std::stod(y_coord_text);
            // Cut to extend
            if (out_of_bbox(x, y))
                continue;
            std::tie(row, col) = xy_to_row_col(x_coord_text, y_coord_text);
            NodeId node_id = node_id_from_text(node_text);
            if (node_id < 1) {
                std::runtime_error("Node ID must be greater than zero");
            }
            nodes_by_row_col_[std::make_pair(row, col)].insert(node_id);
            node_ids.insert(node_id);
        }
    }
 
    template<typename InputStream>
    void load_segments(InputStream& stream, const std::set<NodeId>& node_ids)
    {
        std::set<std::pair<NodeId, NodeId>> node_pairs;
        std::string line;
        while (std::getline(stream, line)) {
            std::istringstream line_stream{line};
            char delimeter{','};
            std::string node_1_text;
            std::getline(line_stream, node_1_text, delimeter);
            std::string node_2_text;
            std::getline(line_stream, node_2_text, delimeter);
            auto node_1_id = node_id_from_text(node_1_text);
            auto node_2_id = node_id_from_text(node_2_text);
            // If either end nodes of the segment is not in the extent, skip it.
            // This means that a segment is ignored even if one of the nodes and
            // significant portion of the segment is in the area of iterest.
            // TODO: Part of the segment may still be out, so that needs to be checked.
            // Replace find by contains for C++20.
            if (node_ids.find(node_1_id) == node_ids.end()
                || node_ids.find(node_2_id) == node_ids.end())
                continue;
            if (node_1_id == node_2_id) {
                std::runtime_error(
                    std::string("Segment cannot begin and end with the same node: ")
                    + node_1_text + " " + node_2_text);
            }
            std::string segment_text;
            std::getline(line_stream, segment_text, delimeter);
            // We don't know which way the nodes are ordered, so instead of checking
            // the order, we create a symmetric matrix since we allocated the memory
            // anyway.
            node_pairs.emplace(node_1_id, node_2_id);
            std::istringstream segment_stream{segment_text};
            char in_cell_delimeter{';'};
            std::string x_coord_text;
            std::string y_coord_text;
            Segment segment;
            while (std::getline(segment_stream, x_coord_text, in_cell_delimeter)
                   && std::getline(segment_stream, y_coord_text, in_cell_delimeter)) {
                // The same cell is possibly repeated if raster resolution is lower than
                // the detail ("resolution") of each segment, so we check if the last
                // added coordinate pair converted to same cell.
                auto new_point = xy_to_row_col(x_coord_text, y_coord_text);
                if (segment.empty() || segment.back() != new_point)
                    segment.emplace_back(new_point);
            }
            segments_by_nodes_.emplace(
                std::make_pair(node_1_id, node_2_id), std::move(segment));
        }
 
        for (auto node_id : node_ids) {
            std::vector<NodeId> nodes;
            for (const auto& item : node_pairs) {
                if (item.first == node_id)
                    nodes.push_back(item.second);
                else if (item.second == node_id)
                    nodes.push_back(item.first);
            }
            node_matrix_[node_id] = std::move(nodes);
        }
    }
 
    SegmentView get_segment(NodeId start, NodeId end) const
    {
        for (const auto& item : segments_by_nodes_) {
            const auto& key{item.first};
            if (key.first == start && key.second == end)
                return SegmentView(item.second.cbegin(), item.second.cend());
            if (key.second == start && key.first == end) {
                return SegmentView(item.second.crbegin(), item.second.crend());
            }
        }
        throw std::invalid_argument("No segment for given nodes");
    }
 
    const std::vector<NodeId>& nodes_connected_to(NodeId node) const
    {
        return node_matrix_.at(node);
    }
 
    template<typename Generator>
    NodeId
    next_node(NodeId node, const std::set<NodeId>& ignore, Generator& generator) const
    {
        // Get all candidate nodes.
        const auto& all_nodes = nodes_connected_to(node);
 
        // Resolve disconnected node and dead end cases.
        auto num_nodes = all_nodes.size();
        if (!num_nodes)
            return node;
        else if (num_nodes == 1)
            return all_nodes[0];
 
        // Filter out the ignored nodes.
        std::vector<int> nodes;
        std::back_insert_iterator<std::vector<int>> back_it(nodes);
        std::remove_copy_if(
            all_nodes.begin(), all_nodes.end(), back_it, [&ignore](NodeId id) {
                return container_contains(ignore, id);
            });
 
        // Pick a random node. Fallback to all candidate nodes if all are in ignore.
        num_nodes = nodes.size();
        if (!num_nodes)
            return pick_random_node(all_nodes, generator);
        else if (num_nodes == 1)
            return nodes[0];
        return pick_random_node(nodes, generator);
    }
 
    template<typename Container, typename Generator>
    static NodeId pick_random_node(const Container& nodes, Generator& generator)
    {
        auto num_nodes = nodes.size();  // Replace by std::size in C++17.
        std::uniform_int_distribution<size_t> dist(0, num_nodes - 1);
        auto index = dist(generator);
        // Our lists are expected to be short, so this is expected to be fast for
        // both sets and vectors.
        return *std::next(nodes.begin(), index);
    }
 
    BBox<double> bbox_;  
    double ew_res_;  
    double ns_res_;  
    double cell_travel_time_;  
 
    std::map<std::pair<RasterIndex, RasterIndex>, std::set<NodeId>> nodes_by_row_col_;
    NodeMatrix node_matrix_;  
    SegmentsByNodes segments_by_nodes_;  
};
 
template<typename RasterIndex>
class NetworkDispersalKernel
{
public:
    NetworkDispersalKernel(
        const Network<RasterIndex>& network, double min_time, double max_time)
        : network_(network), time_distribution_(min_time, max_time)
    {}
 
    template<typename Generator>
    std::tuple<int, int> operator()(Generator& generator, int row, int col)
    {
        double time = time_distribution_(generator);
        std::tie(row, col) = network_.travel(row, col, time, generator);
 
        return std::make_tuple(row, col);
    }
 
    bool is_cell_eligible(int row, int col)
    {
        return network_.has_node_at(row, col);
    }
 
    static bool supports_kernel(const DispersalKernelType type)
    {
        return type == DispersalKernelType::Network;
    }
 
protected:
    const Network<RasterIndex>& network_;
    std::uniform_real_distribution<double> time_distribution_;
};
 
}  // namespace pops
 
#endif  // POPS_NETWORK_KERNEL_HPP