Use scaled natural kernel for overpopulation (#135)

Instead of (mis)using antropogenic kernel as a long distance kernel for pest overpopulation movement, scale (multiply) the scale (distance) parameter of natural kernel using a coefficient.

This adds one new config variable to modify the scale parameter of radial and deterministic kernels. It is meant as modification, so it is called coefficient [of scale] (`leaving_scale_coefficient`) and it defaults to 1. "Coefficient of scale" is not ideal, but "scale of scale" would not be great either.

Additionally, this moves most of the kernel creation to separate documented functions. Although the duplication between the functions would be unnecessary difficult to address, the similarities between the three kernels (natural, anthro, overflow) are now more clear than in the original intermingled code and easier to keep in sync if needed. The kernels are now created only when needed (i.e., with limited scope and potentially not for every step) since they are really not only local to the run step function, but needed only for one function call.

The new functions are protected and the existing attributes are now protected instead of private because for practical reasons, similarly to other places, we consider inheritance a case when user of a class wants to deal with internals of the object.

Finally, when config is not properly used, the new behavior is that it now consistently fails while before it was sometimes passing. In the config, we now set more defaults to ensure consistent results (that is throwing an exception when not set because the default and wrong 0 is used).

Additionally, invalid radial kernel message code is simplified and the deterministic kernel now documents the use of the dispersal reference and how invalid arguments are handled.

CHANGELOG.md

@@ -18,9 +18,16 @@ this repository.
 
 ### Added
 
-- Pest pest movement based on overpopulation (Vaclav Petras)
+- Pest movement based on overpopulation (Vaclav Petras)
   * When cell contains too many pests, pests leave and move to a different cell.
 
+### Changed
+- Model class internal attributes and functions are now protected instead of private
+  to allow derived classes to access them for greater flexibility.
+- Model class kernels are now created and returned by protected functions.
+- Config class has now more defaults and subsequent setup now consistently fails when
+  required values were not set.
+
 ## [1.0.2] - 2020-10-09
 
 - [Patch release of rpops](https://github.com/ncsu-landscape-dynamics/rpops/releases/tag/v1.0.2) (no changes in pops-core)

include/pops/config.hpp

@@ -57,15 +57,15 @@ class Config
     int latency_period_steps;
     // Kernels
     std::string natural_kernel_type;
-    double natural_scale;
+    double natural_scale{0};
     std::string natural_direction;
-    double natural_kappa;
+    double natural_kappa{0};
     bool use_anthropogenic_kernel{false};
-    double percent_natural_dispersal;
+    double percent_natural_dispersal{1};
     std::string anthro_kernel_type;
-    double anthro_scale;
+    double anthro_scale{0};
     std::string anthro_direction;
-    double anthro_kappa;
+    double anthro_kappa{0};
     double shape{1.0};
     // Treatments
     bool use_treatments{false};
@@ -89,6 +89,7 @@ class Config
     bool use_overpopulation_movements{false};
     double overpopulation_percentage{0};
     double leaving_percentage{0};
+    double leaving_scale_coefficient{1};
 
     void create_schedules()
     {

include/pops/deterministic_kernel.hpp

@@ -92,6 +92,25 @@ class DeterministicDispersalKernel
     double north_south_resolution;
 
 public:
+    /**
+     * @brief DeterministicDispersalKernel constructor
+     *
+     * When an unsupported (invalid) kernel type is passed as *dispersal_kernel*,
+     * std::invalid_argument is thrown when the function call operator is used,
+     * i.e., it is accepted by the constructor. This is to allow for constuction
+     * of invalid, but unused kernels which are created as a part of larger kernels.
+     *
+     * The reference *dispersers* is later used to obtain the current counts of
+     * dispersers.
+     *
+     * @param dispersal_kernel Type of kernel to be used
+     * @param dispersers Reference to a dispersers raster
+     * @param dispersal_percentage Percentage used to compute the kernel size
+     * @param ew_res East-west resolution
+     * @param ns_res North-south resolution
+     * @param distance_scale Scale parameter for the kernels
+     * @param shape Shape parameter for the kernels
+     */
     DeterministicDispersalKernel(
         DispersalKernelType dispersal_kernel,
         const IntegerRaster& dispersers,

include/pops/model.hpp

@@ -39,7 +39,7 @@ namespace pops {
 template<typename IntegerRaster, typename FloatRaster, typename RasterIndex>
 class Model
 {
-private:
+protected:
     Config config_;
     DispersalKernelType natural_kernel;
     DispersalKernelType anthro_kernel;
@@ -49,6 +49,120 @@ class Model
     Simulation<IntegerRaster, FloatRaster, RasterIndex> simulation_;
     unsigned last_index{0};
 
+    /**
+     * @brief Create natural kernel
+     *
+     * Kernel parameters are taken from the configuration.
+     *
+     * @param dispersers The disperser raster (reference, for deterministic kernel)
+     * @return Created kernel
+     */
+    SwitchDispersalKernel<IntegerRaster>
+    create_natural_kernel(const IntegerRaster& dispersers)
+    {
+        RadialDispersalKernel<IntegerRaster> radial_kernel(
+            config_.ew_res,
+            config_.ns_res,
+            natural_kernel,
+            config_.natural_scale,
+            direction_from_string(config_.natural_direction),
+            config_.natural_kappa,
+            config_.shape);
+        DeterministicDispersalKernel<IntegerRaster> deterministic_kernel(
+            natural_kernel,
+            dispersers,
+            config_.dispersal_percentage,
+            config_.ew_res,
+            config_.ns_res,
+            config_.natural_scale,
+            config_.shape);
+        SwitchDispersalKernel<IntegerRaster> selectable_kernel(
+            natural_kernel,
+            radial_kernel,
+            deterministic_kernel,
+            uniform_kernel,
+            natural_neighbor_kernel,
+            config_.deterministic);
+        return selectable_kernel;
+    }
+
+    /**
+     * @brief Create overpopulation movement kernel
+     *
+     * Same as the natural kernel. The natural kernel parameters are used,
+     * but the scale for radial and deterministic kernel is multiplied
+     * by the leaving scale coefficient.
+     *
+     * @param dispersers The disperser raster (reference, for deterministic kernel)
+     * @return Created kernel
+     */
+    SwitchDispersalKernel<IntegerRaster>
+    create_overpopulation_movement_kernel(const IntegerRaster& dispersers)
+    {
+        RadialDispersalKernel<IntegerRaster> radial_kernel(
+            config_.ew_res,
+            config_.ns_res,
+            natural_kernel,
+            config_.natural_scale * config_.leaving_scale_coefficient,
+            direction_from_string(config_.natural_direction),
+            config_.natural_kappa,
+            config_.shape);
+        DeterministicDispersalKernel<IntegerRaster> deterministic_kernel(
+            natural_kernel,
+            dispersers,
+            config_.dispersal_percentage,
+            config_.ew_res,
+            config_.ns_res,
+            config_.natural_scale * config_.leaving_scale_coefficient,
+            config_.shape);
+        SwitchDispersalKernel<IntegerRaster> selectable_kernel(
+            natural_kernel,
+            radial_kernel,
+            deterministic_kernel,
+            uniform_kernel,
+            natural_neighbor_kernel,
+            config_.deterministic);
+        return selectable_kernel;
+    }
+
+    /**
+     * @brief Create anthropogenic kernel
+     *
+     * Same structure as the natural kernel, but the parameters are for anthropogenic
+     * kernel when available.
+     *
+     * @param dispersers The disperser raster (reference, for deterministic kernel)
+     * @return Created kernel
+     */
+    SwitchDispersalKernel<IntegerRaster>
+    create_anthro_kernel(const IntegerRaster& dispersers)
+    {
+        RadialDispersalKernel<IntegerRaster> radial_kernel(
+            config_.ew_res,
+            config_.ns_res,
+            anthro_kernel,
+            config_.anthro_scale,
+            direction_from_string(config_.anthro_direction),
+            config_.anthro_kappa,
+            config_.shape);
+        DeterministicDispersalKernel<IntegerRaster> deterministic_kernel(
+            anthro_kernel,
+            dispersers,
+            config_.dispersal_percentage,
+            config_.ew_res,
+            config_.ns_res,
+            config_.anthro_scale,
+            config_.shape);
+        SwitchDispersalKernel<IntegerRaster> selectable_kernel(
+            anthro_kernel,
+            radial_kernel,
+            deterministic_kernel,
+            uniform_kernel,
+            anthro_neighbor_kernel,
+            config_.deterministic);
+        return selectable_kernel;
+    }
+
 public:
     Model(const Config& config)
         : config_(config),
@@ -133,57 +247,6 @@ class Model
         const std::vector<std::vector<int>> movements,
         const std::vector<std::vector<int>>& suitable_cells)
     {
-        RadialDispersalKernel<IntegerRaster> natural_radial_kernel(
-            config_.ew_res,
-            config_.ns_res,
-            natural_kernel,
-            config_.natural_scale,
-            direction_from_string(config_.natural_direction),
-            config_.natural_kappa,
-            config_.shape);
-        RadialDispersalKernel<IntegerRaster> anthro_radial_kernel(
-            config_.ew_res,
-            config_.ns_res,
-            anthro_kernel,
-            config_.anthro_scale,
-            direction_from_string(config_.anthro_direction),
-            config_.anthro_kappa,
-            config_.shape);
-        DeterministicDispersalKernel<IntegerRaster> natural_deterministic_kernel(
-            natural_kernel,
-            dispersers,
-            config_.dispersal_percentage,
-            config_.ew_res,
-            config_.ns_res,
-            config_.natural_scale,
-            config_.shape);
-        DeterministicDispersalKernel<IntegerRaster> anthro_deterministic_kernel(
-            anthro_kernel,
-            dispersers,
-            config_.dispersal_percentage,
-            config_.ew_res,
-            config_.ns_res,
-            config_.anthro_scale,
-            config_.shape);
-        SwitchDispersalKernel<IntegerRaster> natural_selectable_kernel(
-            natural_kernel,
-            natural_radial_kernel,
-            natural_deterministic_kernel,
-            uniform_kernel,
-            natural_neighbor_kernel,
-            config_.deterministic);
-        SwitchDispersalKernel<IntegerRaster> anthro_selectable_kernel(
-            anthro_kernel,
-            anthro_radial_kernel,
-            anthro_deterministic_kernel,
-            uniform_kernel,
-            anthro_neighbor_kernel,
-            config_.deterministic);
-        DispersalKernel<IntegerRaster> dispersal_kernel(
-            natural_selectable_kernel,
-            anthro_selectable_kernel,
-            config_.use_anthropogenic_kernel,
-            config_.percent_natural_dispersal);
         int mortality_simulation_year =
             simulation_step_to_action_step(config_.mortality_schedule(), step);
         // removal of dispersers due to lethal temperatures
@@ -207,6 +270,14 @@ class Model
                 config_.reproductive_rate,
                 suitable_cells);
 
+            DispersalKernel<IntegerRaster> dispersal_kernel(
+                create_natural_kernel(dispersers),
+                create_anthro_kernel(dispersers),
+                config_.use_anthropogenic_kernel,
+                config_.percent_natural_dispersal);
+            auto overpopulation_kernel =
+                create_overpopulation_movement_kernel(dispersers);
+
             simulation_.disperse_and_infect(
                 step,
                 dispersers,
@@ -227,7 +298,7 @@ class Model
                     infected,
                     total_populations,
                     outside_dispersers,
-                    anthro_selectable_kernel,
+                    overpopulation_kernel,
                     suitable_cells,
                     config_.overpopulation_percentage,
                     config_.leaving_percentage);

include/pops/radial_kernel.hpp

@@ -68,9 +68,7 @@ inline Direction direction_from_string(const std::string& text)
     }
     catch (const std::out_of_range&) {
         throw std::invalid_argument(
-            "direction_from_string: Invalid"
-            " value '"
-            + text + "' provided");
+            "direction_from_string: Invalid value '" + text + "' provided");
     }
 }
 

tests/test_overpopulation_movements.cpp

@@ -123,16 +123,13 @@ int test_model()
     config.model_type = "SI";
     config.natural_kernel_type = "cauchy";
     config.natural_scale = 0.1;
-    // We are not using anthropo kernel in standard disperser spread, but it is used by
-    // overpopulation movements.
-    config.use_anthropogenic_kernel = false;
-    config.anthro_kernel_type = "cauchy";
-    config.anthro_scale = 0.1;
+    config.anthro_scale = config.natural_scale;  // Unused, but we need to set it.
     config.ew_res = 30;
     config.ns_res = 30;
     config.use_overpopulation_movements = true;
     config.overpopulation_percentage = 0.5;
     config.leaving_percentage = 0.75;
+    config.leaving_scale_coefficient = 1;
     config.create_schedules();
     // More reference data
     auto leaving = infected(0, 0) * config.leaving_percentage;