quarantine.hpp

Go to the documentation of this file.

/*
 * PoPS model - Quarantine escape computation
 *
 * Copyright (C) 2020 by the authors.
 *
 * Authors: Anna Petrasova <akratoc ncsu edu>
 *
 * 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_QUARANTINE_HPP
#define POPS_QUARANTINE_HPP
 
#include <tuple>
#include <map>
#include <vector>
#include <limits>
#include <type_traits>
#include <sstream>
#include <iomanip>
 
#include "utils.hpp"
 
namespace pops {
 
std::ostream& operator<<(std::ostream& os, const Direction& obj)
{
    os << static_cast<std::underlying_type<Direction>::type>(obj);
    return os;
}
 
typedef std::tuple<double, Direction> DistDir;
typedef std::tuple<bool, DistDir> EscapeDistDir;
typedef std::vector<EscapeDistDir> EscapeDistDirs;
 
template<typename IntegerRaster, typename RasterIndex = int>
class QuarantineEscape
{
private:
    RasterIndex width_;
    RasterIndex height_;
    // the west-east resolution of the pixel
    double west_east_resolution_;
    // the north-south resolution of the pixel
    double north_south_resolution_;
    unsigned num_steps;
    std::vector<BBoxInt> boundaries;
    // mapping between quarantine areas is from map and index
    std::map<int, int> boundary_id_idx_map;
    std::vector<EscapeDistDir> escape_dist_dirs;
 
    void quarantine_boundary(
        const IntegerRaster& quarantine_areas,
        const std::vector<std::vector<int>>& suitable_cells)
    {
        int n, s, e, w;
        int idx = 0;
        for (auto indices : suitable_cells) {
            int i = indices[0];
            int j = indices[1];
            auto value = quarantine_areas(i, j);
            if (value > 0) {
                auto search = boundary_id_idx_map.find(value);
                int bidx;
                if (search == boundary_id_idx_map.end()) {
                    boundary_id_idx_map.insert(std::make_pair(value, idx));
                    boundaries.push_back(
                        std::make_tuple(height_ - 1, 0, 0, width_ - 1));
                    bidx = idx;
                    ++idx;
                }
                else
                    bidx = search->second;
                std::tie(n, s, e, w) = boundaries.at(bidx);
                if (i < n)
                    n = i;
                if (i > s)
                    s = i;
                if (j > e)
                    e = j;
                if (j < w)
                    w = j;
                boundaries.at(bidx) = std::make_tuple(n, s, e, w);
            }
        }
    }
    DistDir
    closest_direction(RasterIndex i, RasterIndex j, const BBoxInt boundary) const
    {
        int n, s, e, w;
        int mindist = std::numeric_limits<int>::max();
        std::tie(n, s, e, w) = boundary;
        DistDir closest;
        if ((i - n) * north_south_resolution_ < mindist) {
            mindist = (i - n) * north_south_resolution_;
            closest = std::make_tuple(mindist, Direction::N);
        }
        if ((s - i) * north_south_resolution_ < mindist) {
            mindist = (s - i) * north_south_resolution_;
            closest = std::make_tuple(mindist, Direction::S);
        }
        if ((e - j) * west_east_resolution_ < mindist) {
            mindist = (e - j) * west_east_resolution_;
            closest = std::make_tuple(mindist, Direction::E);
        }
        if ((j - w) * west_east_resolution_ < mindist) {
            mindist = (j - w) * west_east_resolution_;
            closest = std::make_tuple(mindist, Direction::W);
        }
        return closest;
    }
 
public:
    QuarantineEscape(
        const IntegerRaster& quarantine_areas,
        double ew_res,
        double ns_res,
        unsigned num_steps,
        const std::vector<std::vector<int>>& suitable_cells)
        : width_(quarantine_areas.cols()),
          height_(quarantine_areas.rows()),
          west_east_resolution_(ew_res),
          north_south_resolution_(ns_res),
          num_steps(num_steps),
          escape_dist_dirs(
              num_steps,
              std::make_tuple(
                  false,
                  std::make_tuple(std::numeric_limits<double>::max(), Direction::None)))
    {
        quarantine_boundary(quarantine_areas, suitable_cells);
    }
 
    QuarantineEscape() = delete;
 
    void infection_escape_quarantine(
        const IntegerRaster& infected,
        const IntegerRaster& quarantine_areas,
        unsigned step,
        const std::vector<std::vector<int>>& suitable_cells)
    {
 
        DistDir min_dist_dir =
            std::make_tuple(std::numeric_limits<double>::max(), Direction::None);
        for (auto indices : suitable_cells) {
            int i = indices[0];
            int j = indices[1];
            if (!infected(i, j))
                continue;
            int area = quarantine_areas(i, j);
            if (area == 0) {
                escape_dist_dirs.at(step) = std::make_tuple(
                    true, std::make_tuple(std::nan(""), Direction::None));
                return;
            }
            double dist;
            Direction dir;
            int bindex = boundary_id_idx_map[area];
            std::tie(dist, dir) = closest_direction(i, j, boundaries.at(bindex));
            if (dist < std::get<0>(min_dist_dir)) {
                min_dist_dir = std::make_tuple(dist, dir);
            }
        }
        escape_dist_dirs.at(step) = std::make_tuple(false, min_dist_dir);
    }
    const EscapeDistDir escape_info(unsigned step) const
    {
        return escape_dist_dirs.at(step);
    }
 
    bool escaped(unsigned step) const
    {
        return std::get<0>(escape_dist_dirs.at(step));
    }
 
    double distance(unsigned step) const
    {
        auto dist_dir = std::get<1>(escape_dist_dirs.at(step));
        return std::get<0>(dist_dir);
    }
 
    Direction direction(unsigned step) const
    {
        auto dist_dir = std::get<1>(escape_dist_dirs.at(step));
        return std::get<1>(dist_dir);
    }
};
 
template<typename IntegerRaster>
double quarantine_escape_probability(
    const std::vector<QuarantineEscape<IntegerRaster>> escape_infos, unsigned step)
{
    bool escape;
    DistDir distdir;
    int escapes = 0;
    for (const auto& item : escape_infos) {
        std::tie(escape, distdir) = item.escape_info(step);
        escapes += escape;
    }
    return (double)escapes / escape_infos.size();
}
 
template<typename IntegerRaster>
std::vector<DistDir> distance_direction_to_quarantine(
    const std::vector<QuarantineEscape<IntegerRaster>> escape_infos, unsigned step)
{
    bool escape;
    DistDir distdir;
    std::vector<DistDir> distances_directions;
    for (const auto& item : escape_infos) {
        std::tie(escape, distdir) = item.escape_info(step);
        distances_directions.push_back(distdir);
    }
 
    return distances_directions;
}
 
template<typename IntegerRaster>
std::string write_quarantine_escape(
    const std::vector<QuarantineEscape<IntegerRaster>> escape_infos, unsigned num_steps)
{
    std::stringstream ss;
    ss << std::setprecision(1) << std::fixed;
    ss << "step,escape_probability";
    for (unsigned i = 0; i < escape_infos.size(); i++)
        ss << ",dist" << i << ",dir" << i;
    ss << "\n";
    for (unsigned step = 0; step < num_steps; step++) {
        ss << step;
        ss << "," << quarantine_escape_probability(escape_infos, step);
        std::vector<DistDir> dists_dirs =
            distance_direction_to_quarantine(escape_infos, step);
        double dist;
        Direction dir;
        for (unsigned i = 0; i < dists_dirs.size(); i++) {
            std::tie(dist, dir) = dists_dirs.at(i);
            if (std::isnan(dist))
                ss << ",,";
            else
                ss << "," << dist << "," << dir;
        }
        ss << "\n";
    }
    return ss.str();
}
}  // namespace pops
#endif  // POPS_QUARANTINE_HPP