Skip to content

Commit

Permalink
Use multi-host version of pool and action from core (#70)
Browse files Browse the repository at this point in the history 
Use PoPS Core with host pool and action objects from the upcoming version 3 of PoPS Core.

Add test for mortality with rate 0.

---------

Co-authored-by: Anna Petrasova <[email protected]>
  • Loading branch information
@wenzeslaus @petrasovaa
wenzeslaus and petrasovaa authored Feb 20, 2024
1 parent 5b612d9 commit fbf87f0
Show file tree
Hide file tree
Showing 3 changed files with 176 additions and 47 deletions.
148 changes: 103 additions & 45 deletions main.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,7 @@
#include "graster.hpp"

#include "pops/model.hpp"
#include "pops/model_type.hpp"
#include "pops/date.hpp"
#include "pops/raster.hpp"
#include "pops/kernel.hpp"
Expand All @@ -26,6 +27,9 @@
#include "pops/statistics.hpp"
#include "pops/scheduling.hpp"
#include "pops/quarantine.hpp"
#include "pops/host_pool.hpp"
#include "pops/pest_pool.hpp"
#include "pops/multi_host_pool.hpp"

extern "C" {
#include <grass/gis.h>
Expand Down Expand Up @@ -56,9 +60,6 @@ using std::isnan;

using namespace pops;

// TODO: for backwards compatibility, update eventually
typedef Simulation<Img, DImg> Sporulation;

#define DIM 1

void fatal_option_required_for_value(struct Option* required, struct Option* given)
Expand Down Expand Up @@ -743,14 +744,17 @@ int main(int argc, char* argv[])
opt.dispersers_output = G_define_standard_option(G_OPT_R_OUTPUT);
opt.dispersers_output->key = "dispersers_output";
opt.dispersers_output->label = _("Output raster of disperses");
opt.dispersers_output->description = _("Dispersers are accumulated over all steps and stochastic runs");
opt.dispersers_output->description =
_("Dispersers are accumulated over all steps and stochastic runs");
opt.dispersers_output->required = NO;
opt.dispersers_output->guisection = _("Output");

opt.established_dispersers_output = G_define_standard_option(G_OPT_R_OUTPUT);
opt.established_dispersers_output->key = "established_dispersers_output";
opt.established_dispersers_output->label = _("Output raster of established disperses");
opt.established_dispersers_output->description = _("Dispersers are accumulated over all steps and stochastic runs");
opt.established_dispersers_output->label =
_("Output raster of established disperses");
opt.established_dispersers_output->description =
_("Dispersers are accumulated over all steps and stochastic runs");
opt.established_dispersers_output->required = NO;
opt.established_dispersers_output->guisection = _("Output");

Expand Down Expand Up @@ -1043,6 +1047,11 @@ int main(int argc, char* argv[])
? std::stoi(opt.mortality_frequency_n->answer)
: 0;
}
config.create_pest_host_table_from_parameters(1);
std::vector<std::vector<double>> competency_table_data;
competency_table_data.push_back({1, 1});
competency_table_data.push_back({0, 0});
config.read_competency_table(competency_table_data);

if (opt.survival_rate_month->answer)
config.survival_rate_month = std::stoi(opt.survival_rate_month->answer);
Expand Down Expand Up @@ -1206,6 +1215,8 @@ int main(int argc, char* argv[])
std::stod(opt.dispersers_to_soils->answer);
}

using SpreadModel = Model<Img, DImg, DImg::IndexType>;

// treatments
if (get_num_answers(opt.treatments) != get_num_answers(opt.treatment_date)
&& get_num_answers(opt.treatment_date)
Expand All @@ -1220,7 +1231,8 @@ int main(int argc, char* argv[])
TreatmentApplication treatment_app = TreatmentApplication::Ratio;
if (opt.treatment_app->answer)
treatment_app = treatment_app_enum_from_string(opt.treatment_app->answer);
Treatments<Img, DImg> treatments(config.scheduler());
Treatments<SpreadModel::StandardSingleHostPool, DImg> treatments(
config.scheduler());
config.use_treatments = false;
if (opt.treatments->answers) {
for (int i_t = 0; opt.treatment_date->answers[i_t]; i_t++) {
Expand All @@ -1234,15 +1246,17 @@ int main(int argc, char* argv[])
}
}

// build the Sporulation object
std::vector<Model<Img, DImg, int>> models;
// build the model object
std::vector<SpreadModel> models;
std::vector<Img> dispersers;
std::vector<Img> established_dispersers;
std::vector<Img> sus_species_rasts(num_runs, S_species_rast);
std::vector<Img> inf_species_rasts(num_runs, I_species_rast);
std::vector<Img> total_species_rasts(num_runs, species_rast);
std::vector<Img> resistant_rasts(num_runs, Img(S_species_rast, 0));
std::vector<Img> total_exposed_rasts(num_runs, Img(S_species_rast, 0));
std::vector<SpreadModel::StandardMultiHostPool> multi_host_pools;
std::vector<PestPool<Img, DImg, int>> pest_pools;

// We always create at least one exposed for simplicity, but we
// could also just leave it empty.
Expand All @@ -1265,6 +1279,8 @@ int main(int argc, char* argv[])
Img accumulated_dead(Img(S_species_rast, 0));

models.reserve(num_runs);
multi_host_pools.reserve(num_runs);
pest_pools.reserve(num_runs);
dispersers.reserve(num_runs);
established_dispersers.reserve(num_runs);
auto tmp_seeds = config.random_seeds;
Expand All @@ -1285,9 +1301,12 @@ int main(int argc, char* argv[])
catch (const std::invalid_argument& error) {
G_fatal_error(_("Model configuration is invalid: %s"), error.what());
}
dispersers.emplace_back(I_species_rast.rows(), I_species_rast.cols());
// For the model itself, this does not need to be initialized to 0 because the
// model only sets and then looks at suitable cells. However, the output is
// taken from all cells, so even the untouched cells need a value.
dispersers.emplace_back(I_species_rast.rows(), I_species_rast.cols(), 0);
established_dispersers.emplace_back(
I_species_rast.rows(), I_species_rast.cols());
I_species_rast.rows(), I_species_rast.cols(), 0);
}
std::vector<Img> dispersers_rasts(num_runs, Img(S_species_rast, 0));
std::vector<Img> established_dispersers_rasts(num_runs, Img(S_species_rast, 0));
Expand All @@ -1304,28 +1323,75 @@ int main(int argc, char* argv[])
}
std::vector<std::vector<std::tuple<int, int>>> outside_spores(num_runs);

// spread rate initialization
// One host pool for each run (not multi-host case).
std::vector<std::unique_ptr<SpreadModel::StandardSingleHostPool>> host_pools;
host_pools.reserve(num_runs);
std::vector<
std::unique_ptr<pops::PestHostTable<SpreadModel::StandardSingleHostPool>>>
pest_host_tables;
pest_host_tables.reserve(num_runs);
std::vector<
std::unique_ptr<pops::CompetencyTable<SpreadModel::StandardSingleHostPool>>>
competency_tables;
competency_tables.reserve(num_runs);

for (unsigned run = 0; run < num_runs; ++run) {
pest_host_tables.emplace_back(
new pops::PestHostTable<SpreadModel::StandardSingleHostPool>(
config, models[run].environment()));
competency_tables.emplace_back(
new pops::CompetencyTable<SpreadModel::StandardSingleHostPool>(
config, models[run].environment()));

host_pools.emplace_back(new SpreadModel::StandardSingleHostPool(
model_type_from_string(config.model_type),
sus_species_rasts[run],
exposed_vectors[run],
config.latency_period_steps,
inf_species_rasts[run],
total_exposed_rasts[run],
resistant_rasts[run],
mortality_tracker_vector[run],
dead_in_current_year[run],
total_species_rasts[run],
models[run].environment(),
config.generate_stochasticity,
config.reproductive_rate,
config.establishment_stochasticity,
config.establishment_probability,
config.rows,
config.cols,
suitable_cells));
std::vector<SpreadModel::StandardSingleHostPool*> tmp = {host_pools[run].get()};
multi_host_pools.emplace_back(tmp, config);
multi_host_pools[run].set_pest_host_table(*pest_host_tables[run]);
multi_host_pools[run].set_competency_table(*competency_tables[run]);

std::vector<SpreadRate<Img>> spread_rates(
pest_pools.emplace_back(
dispersers[run], established_dispersers[run], outside_spores[run]);
}
std::vector<SpreadRateAction<SpreadModel::StandardMultiHostPool, int>> spread_rates(
num_runs,
SpreadRate<Img>(
I_species_rast,
window.ew_res,
window.ns_res,
config.use_spreadrates ? config.rate_num_steps() : 0,
suitable_cells));
SpreadRateAction<SpreadModel::StandardMultiHostPool, int>(
multi_host_pools[0],
config.rows,
config.cols,
config.ew_res,
config.ns_res,
config.use_spreadrates ? config.rate_num_steps() : 0));
// Quarantine escape tracking
Img quarantine_rast(S_species_rast, 0);
if (config.use_quarantine)
quarantine_rast = raster_from_grass_integer(opt.quarantine->answer);
std::vector<QuarantineEscape<Img>> escape_infos(
std::vector<QuarantineEscapeAction<Img>> escape_infos(
num_runs,
QuarantineEscape<Img>(
QuarantineEscapeAction<Img>(
quarantine_rast,
window.ew_res,
window.ns_res,
config.use_quarantine ? config.quarantine_num_steps() : 0,
config.quarantine_directions));

// Unused movements
std::vector<std::vector<int>> movements;

Expand Down Expand Up @@ -1397,32 +1463,23 @@ int main(int argc, char* argv[])
}
models[run].run_step(
step,
inf_species_rasts[run],
sus_species_rasts[run],
multi_host_pools[run],
pest_pools[run],
lvtree_rast,
total_species_rasts[run],
dispersers[run],
established_dispersers[run],
total_exposed_rasts[run],
exposed_vectors[run],
mortality_tracker_vector[run],
dead_in_current_year[run],
treatments,
actual_temperatures,
survival_rates,
treatments,
resistant_rasts[run],
outside_spores[run],
spread_rates[run],
escape_infos[run],
quarantine_rast,
movements,
Network<Img::IndexType>::null_network(),
suitable_cells);
Network<Img::IndexType>::null_network());
++weather_step;
if (opt.dispersers_output->answer)
dispersers_rasts[run] += dispersers[run];
if (opt.established_dispersers_output->answer)
established_dispersers_rasts[run] += established_dispersers[run];
established_dispersers_rasts[run] +=
established_dispersers[run];
}
}

Expand Down Expand Up @@ -1633,23 +1690,24 @@ int main(int argc, char* argv[])
if (opt.dispersers_output->answer) {
// write final result
raster_to_grass(
dispersers_rasts_sum,
opt.dispersers_output->answer,
"Sum of all dispersers",
interval.end_date());
dispersers_rasts_sum,
opt.dispersers_output->answer,
"Sum of all dispersers",
interval.end_date());
}
}
if (opt.established_dispersers_output->answer) {
Img established_dispersers_rasts_sum(I_species_rast.rows(), I_species_rast.cols(), 0);
Img established_dispersers_rasts_sum(
I_species_rast.rows(), I_species_rast.cols(), 0);
for (unsigned i = 0; i < num_runs; i++)
established_dispersers_rasts_sum += established_dispersers_rasts[i];
if (opt.established_dispersers_output->answer) {
// write final result
raster_to_grass(
established_dispersers_rasts_sum,
opt.established_dispersers_output->answer,
"Sum of all established dispersers",
interval.end_date());
established_dispersers_rasts_sum,
opt.established_dispersers_output->answer,
"Sum of all established dispersers",
interval.end_date());
}
}
return 0;
Expand Down
2 changes: 1 addition & 1 deletion pops-core
Submodule pops-core updated 47 files
+1 −1 .cppcheck_suppressions_cpp.txt
+11 −0 .cppcheck_suppressions_hpp.txt
+4 −3 .github/workflows/cppcheck.yml
+13 −0 .github/workflows/valgrind.yml
+11 −2 CMakeLists.txt
+4 −4 Doxyfile
+532 −0 include/pops/actions.hpp
+2 −0 include/pops/anthropogenic_kernel.hpp
+186 −0 include/pops/competency_table.hpp
+209 −1 include/pops/config.hpp
+4 −3 include/pops/date.hpp
+127 −6 include/pops/environment.hpp
+73 −0 include/pops/environment_interface.hpp
+0 −1 include/pops/gamma_kernel.hpp
+53 −31 include/pops/generator_provider.hpp
+1,186 −0 include/pops/host_pool.hpp
+52 −0 include/pops/host_pool_interface.hpp
+195 −125 include/pops/model.hpp
+63 −0 include/pops/model_type.hpp
+473 −0 include/pops/multi_host_pool.hpp
+2 −0 include/pops/natural_kernel.hpp
+0 −1 include/pops/network.hpp
+121 −0 include/pops/pest_host_table.hpp
+145 −0 include/pops/pest_pool.hpp
+14 −13 include/pops/quarantine.hpp
+437 −573 include/pops/simulation.hpp
+10 −4 include/pops/soils.hpp
+88 −92 include/pops/spread_rate.hpp
+84 −163 include/pops/treatments.hpp
+36 −2 include/pops/utils.hpp
+12 −6 tests/CMakeLists.txt
+151 −0 tests/test_competency_table.cpp
+28 −7 tests/test_environment.cpp
+3 −2 tests/test_generator_provider.cpp
+90 −102 tests/test_model.cpp
+76 −73 tests/test_movements.cpp
+1,571 −0 tests/test_multi_host_model.cpp
+5 −4 tests/test_network.cpp
+3 −6 tests/test_network_kernel.cpp
+5 −8 tests/test_overpopulation_movements.cpp
+44 −8 tests/test_quarantine.cpp
+71 −56 tests/test_simulation.cpp
+4 −14 tests/test_simulation_kernels.cpp
+92 −17 tests/test_soils.cpp
+40 −5 tests/test_spread_rate.cpp
+13 −12 tests/test_survival_rate.cpp
+307 −128 tests/test_treatments.cpp
73 changes: 72 additions & 1 deletion testsuite/test_r_pops_spread.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -146,11 +146,12 @@ def generate_random_raster(name, low, high):
cls.runModule("g.region", raster="lsat7_2002_30", res=85.5, flags="a")
cls.runModule(
"r.circle",
output="infected_patch",
output="raw_infected_patch",
coordinates=[639300, 220900],
max=100,
flags="b",
)
gs.mapcalc("infected_patch = min(raw_infected_patch, host)", superquiet=True)
cls.runModule(
"g.region",
n="n-800",
Expand Down Expand Up @@ -187,6 +188,7 @@ def tearDownClass(cls):
"zero",
"quarantine",
"infected_patch",
"raw_infected_patch",
*cls.weather_names,
*cls.weather_stddev_names,
],
Expand Down Expand Up @@ -831,6 +833,75 @@ def test_outputs_mortality(self):
raster=f"dead_{end}_12_31", reference=values, precision=0.001
)

def test_outputs_mortality_enabled_rate_0(self):
"""Check with mortality rate 0 (values copied from basic outputs test)"""
start = "2019-01-01"
end = "2022-12-31"
self.assertModule(
"r.pops.spread",
host="host",
total_plants="max_host",
infected="infection",
average="average",
average_series="average",
single_series="single",
stddev="stddev",
stddev_series="stddev",
probability="probability",
probability_series="probability",
start_date=start,
end_date=end,
seasonality=[1, 12],
step_unit="week",
step_num_units=1,
reproductive_rate=1,
natural_dispersal_kernel="exponential",
natural_distance=50,
natural_direction="W",
natural_direction_strength=3,
anthropogenic_dispersal_kernel="cauchy",
anthropogenic_distance=1000,
anthropogenic_direction_strength=0,
percent_natural_dispersal=0.95,
random_seed=1,
runs=5,
nprocs=5,
flags="m",
mortality_rate=0,
mortality_time_lag=0,
mortality_series="dead",
mortality_frequency="yearly",
)
self.assertRasterExists("average")
self.assertRasterExists("stddev")
self.assertRasterExists("probability")
end = end[:4]
self.assertRasterExists(f"average_{end}_12_31")
self.assertRasterExists(f"probability_{end}_12_31")
self.assertRasterExists(f"single_{end}_12_31")
self.assertRasterExists(f"stddev_{end}_12_31")

ref_float = dict(datatype="DCELL")
ref_int = dict(datatype="CELL")
self.assertRasterFitsInfo(raster="average", reference=ref_float)
self.assertRasterFitsInfo(raster="stddev", reference=ref_float)
self.assertRasterFitsInfo(raster="probability", reference=ref_float)
self.assertRasterFitsInfo(raster=f"single_{end}_12_31", reference=ref_int)
self.assertRasterFitsInfo(raster=f"average_{end}_12_31", reference=ref_float)
self.assertRasterFitsInfo(
raster=f"probability_{end}_12_31", reference=ref_float
)
self.assertRasterFitsInfo(raster=f"stddev_{end}_12_31", reference=ref_float)

values = dict(null_cells=0, min=0, max=18, mean=1.777)
self.assertRasterFitsUnivar(raster="average", reference=values, precision=0.001)
values = dict(null_cells=0, min=0, max=100, mean=33.664)
self.assertRasterFitsUnivar(
raster="probability", reference=values, precision=0.001
)
values = dict(null_cells=0, min=0, max=7.547, mean=0.945)
self.assertRasterFitsUnivar(raster="stddev", reference=values, precision=0.001)

def test_outputs_mortality_many_runs(self):
"""Check mortality with many stochastic runs"""
start = "2019-01-01"
Expand Down

0 comments on commit fbf87f0

Please sign in to comment.