Add in support of multi-host

[ChrisJones687]

implement support for multi-host for up to n hosts. This should account for host susceptibility and transmissibility.

[ChrisJones687]

should be a table of hosts with hosts being rows and having the following columns:

• hosts
• susceptibility (how easily the host is infected/infested)
• competency (pathogen reproduction success if on this species or pest fecundity if feeding on this host)
• mortality (does the species experience mortality due to being infected/infested) add in mortality susceptible hosts #131

[ChrisJones687]

SOD
Host | Competency | Susceptibility | Mortality
Oak | 0 | 0.7 | 0.05
Tanoak | 0.7 | 1 | 0.1
Bay Laurel | 1 | 0.6 | 0

ECFF
Host | Competency | Susceptibility | Mortality
Cherry | 1 | 1 | 0
Honeysuckle | 0.5 | 0.8 | 0

SLF
Host | Competency | Susceptibility | Mortality
tree of heaven | 1 | 0.7 | 0.05
Black Walnut | 0.15 | 1 | 0.1
Grapes | 0 | 0.5 | 0

[wenzeslaus]

Pick-and-disperse-to-host code for MultiHostPool

Pseudocode implementing two-step process of adding a disperser to a cell. First, a target host is randomly selected based on the number of susceptibles. Then, the disperser is send to the selected host.

template<typename Container>
auto pick_item_by_probability(
    Container& container, const std::vector<double>& probabilities)
    -> decltype(container.begin())
{
    // Naive implementation which always returns the first element.
    return container.begin();
}

template<typename IntegerRaster, typename FloatRaster, typename RasterIndex>
class MultiHostPool
{
    std::vector<HostPool<IntegerRaster, FloatRaster, RasterIndex>> hosts;

    template<typename Generator>
    int disperser_to(RasterIndex row, RasterIndex col, Generator& generator)
    {
        std::vector<double> probabilities;
        for (const auto& host : hosts) {
            probabilities.push_back(host.susceptible_at(row, col));
        }
        auto item = pick_item_by_probability(hosts, probabilities);
        return item->disperser_to(i, j, generator);
    }
};

[wenzeslaus]

Host properties versus a table about host

Instead of dividing the configuration between properties of hosts, actions, and pest, use one PestHostUse table which know not only the numbers for each host, but also the interactions between hosts (and potentially between host and pest). This is the second option for implementing better handling of mortality rate and time lag in #198.

[wenzeslaus]

Spatial constant table without interaction

    class PestHostUseTable
    {
    public:
        void add_host(
            const HostPool* host,
            double competency,
            double susceptibility,
            double mortality)
        {
            records[host] = Record{competency, susceptibility, mortality};
        }
        double reproductive_rate(const HostPool* host) const
        {
            return records[host].competency;
        }

    private:
        struct Record
        {
            double competency;
            double susceptibility;
            double mortality;
        };
        std::map<const HostPool*, Record> records;
    };

Usage in int dispersers_from(RasterIndex row, RasterIndex col, Generator& generator) const:

         if (infected_at(row, col) <= 0)
             return 0;
-        double lambda =
-            environment_.influence_reproductive_rate_at(row, col, reproductive_rate_);
+        double lambda = environment_.influence_reproductive_rate_at(
+            row, col, pest_host_use_table_.reproductive_rate(this));
         int dispersers_from_cell = 0;
         if (dispersers_stochasticity_) {
             std::poisson_distribution<int> distribution(lambda);

[wenzeslaus]

Spatial variable competency based on multiple hosts

For each host:

H1	H2	H3	C
.5	.2	.1	5
.3	.2	.1	2
.5	.4	.5	8

    class PestHostUseTable
    {
    public:
        void add_host_competencies(
            const HostPool* host,
            std::vector<std::vector<double>>
                host_percentages,  // number of subsequent combinations x number of
                                   // hosts
            std::vector<double> competencies  // number of combinations
        )
        {
            // The types assumed here don't fit with the attributes.
            competency_for_host[host] = {host_percentages, competencies};
        }
        double reproductive_rate_at(
            RasterIndex row, RasterIndex col, const HostPool* host) const
        {
            auto percentages = environment_.host_distribution_at(
                row, col);  // does not consider non-hosts
            // Better would be to go over all the rows every time and keep track
            // of the highest competency which fulfilled the criteria.
            for (auto table_row : competency_for_host[host].ordered_from_worse()) {
                bool row_filled = true;
                for (auto percentage : percentages) {
                    if (percentage < table_row[i]) {
                        row_filled = false;
                        break;
                    }
                }
                if (row_filled) {
                    continue;
                }
            }
            return table_row.previous_or_this().reproductive_rate();
        }

    private:
        using TableRow =
            std::pair<std::vector<double>, double>;  // size == number of hosts
        using TableForHostAndProperty = std::vector<TableRow>;  // size variable
        std::map<const HostPool*, std::vector<TableForHostAndProperty>>
            competency_for_host;
    };

[wenzeslaus]

Depends on #200.

[wenzeslaus]

Multi-host presence-absence table

Competency table

H1	H2	H3	C
1	0	0	.5
1	0	1	.2
0	0	1	.8

Lookup of competency

    double
    competency_at(RasterIndex row, RasterIndex col, const HostPool* host) const
    {
        auto presence_absence =
            environment_.host_presence_at(row, col);
        // Better would be to go over all the rows every time and keep track
        // of the highest competency which fulfilled the criteria.
        double competency = 0;
        for (const auto& table_row : competency_for_host_table) {
            // probably faster if we just wait for the iteration over the whole thing
            if (!table_row.host_data[host_to_index(host)])
                continue;
            // Row is considered only if all required hosts are present.
            bool condition_fulfilled = true;
            for (int i = 0; i < presence_absence.size(); ++i) {
                if (presence_absence[i] && table_row.host_data[i]) {
                }
                else {
                    condition_fulfilled = false;
                    break;
                }
            }
            // Pick always the highest competency.
            if (condition_fulfilled && table_row.competency > competency)
                competency = table_row.competency;
        }
        return competency;
    }

[wenzeslaus]

Incoming disperser event in multi-host pool

template<typename IntegerRaster, typename FloatRaster, typename RasterIndex>
class MultiHostPool
{
    std::vector<HostPool<IntegerRaster, FloatRaster, RasterIndex>> hosts;

    template<typename Generator>
    int disperser_to(RasterIndex row, RasterIndex col, Generator& generator)
    {
        std::vector<HostPool*, double> probabilities;
        double total_s_score = 0;
        for (auto& host_pool : host_pools) {
            double s_for_item =
                host_pool.establishment_probability_at(row, col)  // accounts for W, N
                * host_pest_use.susceptibility(host_pool);
            probabilities.push_back({host_pool, s_for_item});
            total_s_score += s_for_item;  // we should make sure this is <=1
        }
        if (pest_or_pathogen == "pest") {
            if (do_establishment_test(total_s_score)) {  // this is now only in host
                auto& host = pick_host_by_probability(probabilities);
                host->add_disperser_at(row, col);  // simply increases the counts
            }
        }
        else if (pest_or_pathogen == "pathogen") {
            auto& host = pick_host_by_probability(probabilities);
            host->disperser_to(i, j, generator);  // with establishment test
        }
    }
};
class HostPool
{
    // ...
    // This is the current code:
    double establishment_probability_at(RasterIndex row, RasterIndex col) const
    {
        double probability_of_establishment =
            (double)(susceptible_(row, col))
            / environment_.total_population_at(row, col);
        return environment_.influence_probability_of_establishment_at(
            row, col, probability_of_establishment);
    }
    void add_disperser_at(RasterIndex row, RasterIndex col)
    {
        susceptible_(row, col) -= 1;
        if (model_type_ == ModelType::SusceptibleInfected) {
            infected_(row, col) += 1;
            mortality_tracker_vector_.back()(row, col) += 1;
        }
        else if (model_type_ == ModelType::SusceptibleExposedInfected) {
            exposed_.back()(row, col) += 1;
            total_exposed_(row, col) += 1;
        }
    }
    int disperser_to(RasterIndex row, RasterIndex col, Generator& generator)
    {
        std::uniform_real_distribution<double> distribution_uniform(0.0, 1.0);
        if (susceptible_(row, col) > 0) {
            double probability_of_establishment =
                establishment_probability_at(row, col);
            double establishment_tester = 1 - deterministic_establishment_probability_;
            if (establishment_stochasticity_)
                establishment_tester = distribution_uniform(generator);
            if (establishment_tester < probability_of_establishment) {
                add_disperser_at(row, col);
                return true;
            }
        }
        return false;
    }
};