diff --git a/CMakeLists.txt b/CMakeLists.txt index d21cc62..400c11d 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -2,7 +2,7 @@ cmake_minimum_required(VERSION 3.18...3.25 FATAL_ERROR) # Optionally set the version of flamegpu which should be used, ideally a tag (i.e. `v2.0.0-rc`) or branch name, or potentially a commit hash. -set(FLAMEGPU_VERSION "v2.0.0-rc.1" CACHE STRING "FLAMEGPU/FLAMEGPU2 git branch or tag to use") +set(FLAMEGPU_VERSION "v2.0.0-rc.2" CACHE STRING "FLAMEGPU/FLAMEGPU2 git branch or tag to use") # If the above version is a hash instead, also set FLAMEGPU_VERSION_ALLOW_HASH to ON # set(FLAMEGPU_VERSION_ALLOW_HASH "ON") @@ -29,6 +29,11 @@ option(BUILD_TESTING "Build the testing tree." OFF) # Option to enable google test test discovery cmake_dependent_option(ENABLE_GTEST_DISCOVER "Enable GTEST_DISCOVER for more detailed ctest output without -VV. This dramatically increases test suite runtime to CUDA context initialisation." OFF "BUILD_TESTING" OFF) +# CMAke Cache option for the maximum number of per-agent random interactions. +# If the value is changed here, CMakeCache.txt / the build dir will need deleting +# Or reconfigure with -DEXATEPP_ABM_MAX_RANDOM_DAILY_INTERACTIONS= +set(EXATEPP_ABM_MAX_RANDOM_DAILY_INTERACTIONS "20" CACHE STRING "Maximum number of random daily interactions per agent") + # Include common rules from the FLAMEGPU/FLAMEGPU2 repositories CMake include(${FLAMEGPU_ROOT}/cmake/common.cmake) diff --git a/data/inputs/sample.csv b/data/inputs/sample.csv index ee38e8b..41072ed 100644 --- a/data/inputs/sample.csv +++ b/data/inputs/sample.csv @@ -1,4 +1,7 @@ -rng_seed,param_id,duration,n_total,population_0_9,population_10_19,population_20_29,population_30_39,population_40_49,population_50_59,population_60_69,population_70_79,population_80,n_seed_infection,p_interaction_susceptible_to_exposed,mean_time_to_infected,sd_time_to_infected,mean_time_to_recovered,sd_time_to_recovered,mean_time_to_susceptible,sd_time_to_susceptible,relative_susceptibility_0_9,relative_susceptibility_10_19,relative_susceptibility_20_29,relative_susceptibility_30_39,relative_susceptibility_40_49,relative_susceptibility_50_59,relative_susceptibility_60_69,relative_susceptibility_70_79,relative_susceptibility_80 -12,0,365,16384,4,4,4,4,4,4,3,2,1,8,0.10,3.5,1.0,7.0,1.0,60,5.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0 -12,0,365,16384,4,4,4,4,4,4,4,4,4,8,0.10,3.5,1.0,7.0,1.0,180,5.0,0.0,0.0,0.0,0.01,0.05,0.01,0.05,1.0,5.0 - +rng_seed,param_id,duration,n_total,n_seed_infection,population_0_9,population_10_19,population_20_29,population_30_39,population_40_49,population_50_59,population_60_69,population_70_79,population_80,household_size_1,household_size_2,household_size_3,household_size_4,household_size_5,household_size_6,p_interaction_susceptible_to_exposed,mean_time_to_infected,sd_time_to_infected,mean_time_to_recovered,sd_time_to_recovered,mean_time_to_susceptible,sd_time_to_susceptible,relative_susceptibility_0_9,relative_susceptibility_10_19,relative_susceptibility_20_29,relative_susceptibility_30_39,relative_susceptibility_40_49,relative_susceptibility_50_59,relative_susceptibility_60_69,relative_susceptibility_70_79,relative_susceptibility_80,child_network_adults,elderly_network_adults,relative_transmission_household,relative_transmission_occupation,relative_transmission_random,mean_work_interactions_child,mean_work_interactions_adult,mean_work_interactions_elderly,daily_fraction_work,work_network_rewire,mean_random_interactions_0_19,sd_random_interactions_0_19,mean_random_interactions_20_69,sd_random_interactions_20_69,mean_random_interactions_70plus,sd_random_interactions_70plus +12,0,60,32,4,331,370,374,371,361,415,377,270,142,442,492,208,165,51,27,0.10,3.5,2.0,7.0,2.0,28,10.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,0.2,0.2,2.0,0.5,0.5,10,7,3,0.5,0.1,2,2,4,4,3,3 +12,1,365,32,4,331,370,374,371,361,415,377,270,142,442,492,208,165,51,27,0.10,3.5,2.0,7.0,2.0,28,10.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,0.2,0.2,2.0,0.5,0.5,10,7,3,0.5,0.1,2,2,4,4,3,3 +12,2,365,1024,8,331,370,374,371,361,415,377,270,142,442,492,208,165,51,27,0.10,3.5,2.0,7.0,2.0,28,10.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,0.2,0.2,2.0,0.5,0.5,10,7,3,0.5,0.1,2,2,4,4,3,3 +12,3,365,4096,8,331,370,374,371,361,415,377,270,142,442,492,208,165,51,27,0.10,3.5,2.0,7.0,2.0,28,10.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,0.2,0.2,2.0,0.5,0.5,10,7,3,0.5,0.1,2,2,4,4,3,3 +12,4,365,4096,8,331,370,374,371,361,415,377,270,142,442,492,208,165,51,27,0.10,3.5,2.0,7.0,2.0,28,10.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,0.2,0.2,2.0,0.5,0.5,10,7,3,0.5,0.1,2,2,4,4,3,3 +12,5,365,65536,8,331,370,374,371,361,415,377,270,142,442,492,208,165,51,27,0.10,3.5,2.0,7.0,2.0,28,10.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,0.2,0.2,2.0,0.5,0.5,10,7,3,0.5,0.1,2,2,4,4,3,3 diff --git a/readme.md b/readme.md index e8c5290..72b9802 100644 --- a/readme.md +++ b/readme.md @@ -2,7 +2,7 @@ A GPU accelerated Agent-Based Model (ABM) of infection disease spread within a population. -Implemented using [FLAMEGPU/FLAMEGPU2](https://github.com/FLAMEGPU/FLAMEGPU2), with model functionality inspired by [https://github.com/BDI-pathogens/OpenABM-Covid19](BDI-pathogens/OpenABM-Covid19). +Implemented using [FLAMEGPU/FLAMEGPU2](https://github.com/FLAMEGPU/FLAMEGPU2), with model functionality inspired by [BDI-pathogens/OpenABM-Covid19](https://github.com/BDI-pathogens/OpenABM-Covid19). ## Requirements / Dependencies @@ -51,6 +51,7 @@ cmake --build . --target exatepp_abm -j 8 | `CMAKE_CUDA_ARCHITECTURES` | `"50;60;70;80;90"` | ` Specify which CUDA [Compute Capability](https://developer.nvidia.com/cuda-gpus) Architectures to build for | | `CMAKE_BUILD_TYPE` | `Release` | CMake build configuration, setting optimisation levels etc. Choose from [`Release`, `RelWithDebInfo`, `MinSizeRel`, `Debug`] | | `BUILD_TESTING` | `OFF` | Enable / disable the test suite | +| `EXATEPP_ABM_MAX_RANDOM_DAILY_INTERACTIONS` | `20` | The maximum number of per agent random daily interactions for the build. If this is exceeded due to model parameters an error will be reported, and you must reconfigure and rebuild with a higher value. | | `FLAMEGPU_VISUALISATION` | `OFF` | If FLAME GPU's 3D interactive visualisation should be enabled. Requires OpenGL and local execution. | | `FLAMEGPU_SEATBELTS` | `ON` | Enable / Disable additional runtime checks which harm performance but increase usability | | `FLAMEGPU_SHARE_USAGE_STATISTICS` | `ON` | Enable / Disable FLAME GPU 2 telemetry which helps evidence use/impact of FLAME GPU 2. See the [FLAME GPU 2 user guide for more information](https://docs.flamegpu.com/guide/telemetry/) | @@ -121,7 +122,7 @@ Source code linting required `cpplint` to be installed and on your path at CMake Linting can be initiated via cmake, i.e from a build directory: ```bash -cmake --build . --target lint` +cmake --build . --target lint ``` ## Documentation diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 81f9c48..e9fe484 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -21,9 +21,11 @@ SET(LIBRARY_SRC ${CMAKE_CURRENT_SOURCE_DIR}/exateppabm/disease/SEIR.h ${CMAKE_CURRENT_SOURCE_DIR}/exateppabm/exatepp_abm.h ${CMAKE_CURRENT_SOURCE_DIR}/exateppabm/input.h + ${CMAKE_CURRENT_SOURCE_DIR}/exateppabm/network.h ${CMAKE_CURRENT_SOURCE_DIR}/exateppabm/output.h ${CMAKE_CURRENT_SOURCE_DIR}/exateppabm/output/OutputFile.h ${CMAKE_CURRENT_SOURCE_DIR}/exateppabm/output/PerformanceFile.h + ${CMAKE_CURRENT_SOURCE_DIR}/exateppabm/output/PerIndividualFile.h ${CMAKE_CURRENT_SOURCE_DIR}/exateppabm/output/TimeSeriesFile.h ${CMAKE_CURRENT_SOURCE_DIR}/exateppabm/person.h ${CMAKE_CURRENT_SOURCE_DIR}/exateppabm/population.h @@ -35,9 +37,12 @@ SET(LIBRARY_SRC ${CMAKE_CURRENT_SOURCE_DIR}/exateppabm/disease/SEIR.cu ${CMAKE_CURRENT_SOURCE_DIR}/exateppabm/exatepp_abm.cu ${CMAKE_CURRENT_SOURCE_DIR}/exateppabm/input.cu + ${CMAKE_CURRENT_SOURCE_DIR}/exateppabm/network.cu ${CMAKE_CURRENT_SOURCE_DIR}/exateppabm/output.cu ${CMAKE_CURRENT_SOURCE_DIR}/exateppabm/output/PerformanceFile.cpp ${CMAKE_CURRENT_SOURCE_DIR}/exateppabm/output/TimeSeriesFile.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/exateppabm/output/PerIndividualFile.cu + ${CMAKE_CURRENT_SOURCE_DIR}/exateppabm/output/TimeSeriesFile.cpp ${CMAKE_CURRENT_SOURCE_DIR}/exateppabm/person.cu ${CMAKE_CURRENT_SOURCE_DIR}/exateppabm/population.cu ${CMAKE_CURRENT_SOURCE_DIR}/exateppabm/util.cu @@ -83,10 +88,16 @@ endif() target_compile_options(${LIBRARY_NAME} PUBLIC "$<$:--expt-relaxed-constexpr>") # Forward on the cmake build type as a macro definition, for performance output files. There's prolly a cleaner way to do this. -target_compile_definitions(${LIBRARY_NAME} PRIVATE "$<$:CMAKE_BUILD_TYPE=\"Release\">") -target_compile_definitions(${LIBRARY_NAME} PRIVATE "$<$:CMAKE_BUILD_TYPE=\"RelWithDebInfo\">") -target_compile_definitions(${LIBRARY_NAME} PRIVATE "$<$:CMAKE_BUILD_TYPE=\"MinSizeRel\">") -target_compile_definitions(${LIBRARY_NAME} PRIVATE "$<$:CMAKE_BUILD_TYPE=\"Debug\">") +target_compile_definitions(${LIBRARY_NAME} PUBLIC "$<$:CMAKE_BUILD_TYPE=\"Release\">") +target_compile_definitions(${LIBRARY_NAME} PUBLIC "$<$:CMAKE_BUILD_TYPE=\"RelWithDebInfo\">") +target_compile_definitions(${LIBRARY_NAME} PUBLIC "$<$:CMAKE_BUILD_TYPE=\"MinSizeRel\">") +target_compile_definitions(${LIBRARY_NAME} PUBLIC "$<$:CMAKE_BUILD_TYPE=\"Debug\">") + +# Set the upper limit of the per agent max random interaction count +if(EXATEPP_ABM_MAX_RANDOM_DAILY_INTERACTIONS GREATER 0) + target_compile_definitions(${LIBRARY_NAME} PUBLIC "EXATEPP_ABM_MAX_RANDOM_DAILY_INTERACTIONS=${EXATEPP_ABM_MAX_RANDOM_DAILY_INTERACTIONS}") +endif() + # Add src directory to include path, publicly so that the target library inherits this dependency (for now). target_include_directories("${LIBRARY_NAME}" PUBLIC "${CMAKE_CURRENT_SOURCE_DIR}") diff --git a/src/exateppabm/cli.cpp b/src/exateppabm/cli.cpp index fd90e63..279f285 100644 --- a/src/exateppabm/cli.cpp +++ b/src/exateppabm/cli.cpp @@ -13,8 +13,8 @@ void setup(CLI::App& app, std::shared_ptr params) { app.add_flag("-v, --verbose", params->verbosity, "Verbosity of simulation output, forwarded to FLAME GPU 2"); app.add_option("-i, --input-file", params->inputParamFile, "Path to input parameters file"); app.add_option("-n, --param-number", params->inputParamLine, "The line from the parameters file to use. 1 indexed (assuming there is a header)"); - app.add_option("-o, --output-dir", params->outputDir, "Path to output directory"); + app.add_flag("--individual-file", params->individualFile, "Enable the creation of the per individual file"); } void print(const exateppabm::cli::params params) { @@ -22,7 +22,9 @@ void print(const exateppabm::cli::params params) { fmt::print(" device = {}\n", params.device); fmt::print(" verbosity = {}\n", params.verbosity); fmt::print(" inputParamFile = {}\n", params.inputParamFile); + fmt::print(" inputParamLine = {}\n", params.inputParamLine); fmt::print(" outputDir = {}\n", params.outputDir); + fmt::print(" individualFile = {}\n", params.individualFile); fmt::print("}}\n"); } diff --git a/src/exateppabm/cli.h b/src/exateppabm/cli.h index 6947d97..c9eb3c9 100644 --- a/src/exateppabm/cli.h +++ b/src/exateppabm/cli.h @@ -37,6 +37,10 @@ struct params { * Path to directory for file output */ std::string outputDir = std::filesystem::current_path(); + /** + * bool indicating if the individual file should be written + */ + bool individualFile = false; }; /** diff --git a/src/exateppabm/demographics.cu b/src/exateppabm/demographics.cu index 87191b4..6c7b209 100644 --- a/src/exateppabm/demographics.cu +++ b/src/exateppabm/demographics.cu @@ -1,9 +1,11 @@ #include "exateppabm/demographics.h" +#include #include #include "flamegpu/flamegpu.h" #include "exateppabm/input.h" +#include "exateppabm/util.h" namespace exateppabm { namespace demographics { @@ -27,5 +29,42 @@ void define(flamegpu::ModelDescription& model, const exateppabm::input::config& }, true); } +std::array getAllAgeDemographics() { + std::array all = {{ + demographics::Age::AGE_0_9, + demographics::Age::AGE_10_19, + demographics::Age::AGE_20_29, + demographics::Age::AGE_30_39, + demographics::Age::AGE_40_49, + demographics::Age::AGE_50_59, + demographics::Age::AGE_60_69, + demographics::Age::AGE_70_79, + demographics::Age::AGE_80 + }}; + return all; +} + +std::array getAgeDemographicCumulativeProbabilityArray(const exateppabm::input::config& params) { + // Prepare a probability matrix for selecting an age demographic for the agent based on the ratio from the configuration. + // @todo - this could probably be cleaned up + std::uint64_t configDemographicSum = params.population_0_9 + params.population_10_19 + params.population_20_29 + params.population_30_39 + params.population_40_49 + params.population_50_59 + params.population_60_69 + params.population_70_79 + params.population_80; + + std::array demographicProbabilties = {{ + params.population_0_9 / static_cast(configDemographicSum), + params.population_10_19 / static_cast(configDemographicSum), + params.population_20_29 / static_cast(configDemographicSum), + params.population_30_39 / static_cast(configDemographicSum), + params.population_40_49 / static_cast(configDemographicSum), + params.population_50_59 / static_cast(configDemographicSum), + params.population_60_69 / static_cast(configDemographicSum), + params.population_70_79 / static_cast(configDemographicSum), + params.population_80 / static_cast(configDemographicSum) + }}; + // Perform an inclusive scan to convert to cumulative probability + // Using a local method which supports inclusive scans in old libstc++ + exateppabm::util::inclusive_scan(demographicProbabilties.begin(), demographicProbabilties.end(), demographicProbabilties.begin()); + return demographicProbabilties; +} + } // namespace demographics } // namespace exateppabm diff --git a/src/exateppabm/demographics.h b/src/exateppabm/demographics.h index ff22c6c..e7b9931 100644 --- a/src/exateppabm/demographics.h +++ b/src/exateppabm/demographics.h @@ -41,5 +41,20 @@ enum Age : AgeUnderlyingType { */ void define(flamegpu::ModelDescription& model, const exateppabm::input::config& params); +/** + * Get an array containing one of each age demographic enum + * + * This is a workaround for the lack of reflection in c++17, used to simplify code elsewhere + * @return std::array which is the inverse of the Age enum. + */ +std::array getAllAgeDemographics(); + +/** + * Generate a cumulative probability distribution for age demographic sampling for a given simulation configuration + * @param params model configuration parameters + * @return per-age demographic cumulative probability, for sampling with a uniform distribution [0, 1) + */ +std::array getAgeDemographicCumulativeProbabilityArray(const exateppabm::input::config& params); + } // namespace demographics } // namespace exateppabm diff --git a/src/exateppabm/disease/SEIR.cu b/src/exateppabm/disease/SEIR.cu index 0866dc8..12c2cca 100644 --- a/src/exateppabm/disease/SEIR.cu +++ b/src/exateppabm/disease/SEIR.cu @@ -94,9 +94,6 @@ void define(flamegpu::ModelDescription& model, const exateppabm::input::config& env.newProperty("mean_time_to_susceptible", params.mean_time_to_susceptible); env.newProperty("sd_time_to_susceptible", params.sd_time_to_susceptible); - // Define the temporary, hardcoded spatial infection interaction radius? - // env.setProperty("INFECTION_INTERACTION_RADIUS", interactionRadius); - // Get a handle for the Person agent type flamegpu::AgentDescription person = model.Agent(exateppabm::person::NAME); diff --git a/src/exateppabm/exatepp_abm.cu b/src/exateppabm/exatepp_abm.cu index 9716f4f..8a9c9d3 100644 --- a/src/exateppabm/exatepp_abm.cu +++ b/src/exateppabm/exatepp_abm.cu @@ -83,9 +83,7 @@ int entrypoint(int argc, char* argv[]) { flamegpu::ModelDescription model("ExaTEPP ABM demonstrator"); // Add the person agent to the model description - const float env_width = std::ceil(std::sqrt(config->n_total)); - constexpr float interactionRadius = 1.5f; - exateppabm::person::define(model, *config, env_width, interactionRadius); + exateppabm::person::define(model, *config); // Define demographic related variables exateppabm::demographics::define(model, *config); @@ -93,8 +91,11 @@ int entrypoint(int argc, char* argv[]) { // Define disease related variables and methods exateppabm::disease::SEIR::define(model, *config); + // Define init function for population generation + exateppabm::population::define(model, *config, cli_params->verbosity > 0); + // Add init, step and exit functions related to data collection and output. This may want refactoring when multiple output files are supported or collected data becomes more complex. - exateppabm::output::define(model, cli_params->outputDir); + exateppabm::output::define(model, cli_params->outputDir, cli_params->individualFile); // Build the model control flow. This will want abstracting more in the future @todo // @note - not using the DAG control flow due to bugs encountered in another project when splitting between compilation units. @@ -113,6 +114,11 @@ int entrypoint(int argc, char* argv[]) { // Setup simulation configuration options + // If verbosity is high enough (-vvv or more) then enable flamegpu's verbose output + if (cli_params->verbosity > 2) { + simulation.SimulationConfig().verbosity = flamegpu::Verbosity::Verbose; + } + simulation.SimulationConfig().steps = config->duration; // @todo - change this to be controlled by an exit condition? // Seed the FLAME GPU 2 RNG seed. This is independent from RNG on the host, but we only have one RNG engine available in FLAME GPU 2 currently. @@ -121,15 +127,6 @@ int entrypoint(int argc, char* argv[]) { // Set the GPU index simulation.CUDAConfig().device_id = cli_params->device; - // Generate the population of agents. - // @todo - this should probably be an in init function for ease of moving to a ensembles, but then cannot pass parameters in. - const std::uint64_t pop_seed = config->rng_seed; // @todo - split seeds - auto personPopulation = exateppabm::population::generate(model, *config, cli_params->verbosity > 0, env_width, interactionRadius); - if (personPopulation == nullptr) { - throw std::runtime_error("@todo - bad population generation function."); - } - simulation.setPopulationData(*personPopulation); - perfFile.timers.preSimulate.stop(); // Run the simulation diff --git a/src/exateppabm/input.cu b/src/exateppabm/input.cu index 6986dc2..7bc89c2 100644 --- a/src/exateppabm/input.cu +++ b/src/exateppabm/input.cu @@ -15,6 +15,8 @@ namespace exateppabm { namespace input { namespace { +constexpr char header[] = "rng_seed,param_id,duration,n_total,n_seed_infection,population_0_9,population_10_19,population_20_29,population_30_39,population_40_49,population_50_59,population_60_69,population_70_79,population_80,household_size_1,household_size_2,household_size_3,household_size_4,household_size_5,household_size_6,p_interaction_susceptible_to_exposed,mean_time_to_infected,sd_time_to_infected,mean_time_to_recovered,sd_time_to_recovered,mean_time_to_susceptible,sd_time_to_susceptible,relative_susceptibility_0_9,relative_susceptibility_10_19,relative_susceptibility_20_29,relative_susceptibility_30_39,relative_susceptibility_40_49,relative_susceptibility_50_59,relative_susceptibility_60_69,relative_susceptibility_70_79,relative_susceptibility_80,child_network_adults,elderly_network_adults,relative_transmission_household,relative_transmission_occupation,relative_transmission_random,mean_work_interactions_child,mean_work_interactions_adult,mean_work_interactions_elderly,daily_fraction_work,work_network_rewire,mean_random_interactions_0_19,sd_random_interactions_0_19,mean_random_interactions_20_69,sd_random_interactions_20_69,mean_random_interactions_70plus,sd_random_interactions_70plus"; //NOLINT + /** * Get the next value from a csv row as a specified type * @@ -24,20 +26,26 @@ bool valueFromCSVLine(std::string& line, T& value) { // trim whitespace line.erase(0, line.find_first_not_of(' ')); line.erase(line.find_last_not_of(' ') + 1); - std::istringstream iss(line); + + // If the line is empty, do nothing? + if (line.length() == 0) { + return false; + } + + // Find the first comma + size_t pos = line.find(','); + std::string valueString = pos != std::string::npos ? line.substr(0, pos) : line; + // Shorten line. + line = pos != std::string::npos ? line.substr(pos + 1) : ""; + + // Create an istring stream and attempt to parse into value + std::istringstream iss(valueString); if (iss >> value) { - char comma; - iss >> comma; - if (comma != ',') { - // If the next character is not a comma, we've reached the end - line.clear(); - return true; - } - // Remove the consumed part from the string - size_t pos = line.find(','); - line = line.substr(pos + 1); + // success return true; } else { + // failure + // @todo - warn about the value not matching the type here? outside has better info though for the error message return false; } } @@ -72,6 +80,10 @@ std::shared_ptr read(const std::filesystem::path p, c if (!std::getline(fs, line)) { throw std::runtime_error("failed to read the header line @todo nicer error message"); } + // Warn if the header is not the expected value + if (strcmp(line.c_str(), header) != 0) { + fmt::print("Warning: {} header does not match the expected value. Errors may occur, or parameters may be used incorrectly.\nExpected header:\n{}\n", p.c_str(), header); + } // Discard rows until the line number is the target line number for (int currentLine = 1; currentLine < lineNumber; currentLine++) { @@ -97,7 +109,9 @@ std::shared_ptr read(const std::filesystem::path p, c if (!valueFromCSVLine(line, c->n_total)) { throw std::runtime_error("bad value for n_total during csv parsing @todo\n"); } - + if (!valueFromCSVLine(line, c->n_seed_infection)) { + throw std::runtime_error("bad value for n_seed_infection during csv parsing @todo\n"); + } if (!valueFromCSVLine(line, c->population_0_9)) { throw std::runtime_error("bad value for population_0_9 during csv parsing @todo\n"); } @@ -125,8 +139,23 @@ std::shared_ptr read(const std::filesystem::path p, c if (!valueFromCSVLine(line, c->population_80)) { throw std::runtime_error("bad value for population_80 during csv parsing @todo\n"); } - if (!valueFromCSVLine(line, c->n_seed_infection)) { - throw std::runtime_error("bad value for n_seed_infection during csv parsing @todo\n"); + if (!valueFromCSVLine(line, c->household_size_1)) { + throw std::runtime_error("bad value for household_size_1 during csv parsing @todo\n"); + } + if (!valueFromCSVLine(line, c->household_size_2)) { + throw std::runtime_error("bad value for household_size_2 during csv parsing @todo\n"); + } + if (!valueFromCSVLine(line, c->household_size_3)) { + throw std::runtime_error("bad value for household_size_3 during csv parsing @todo\n"); + } + if (!valueFromCSVLine(line, c->household_size_4)) { + throw std::runtime_error("bad value for household_size_4 during csv parsing @todo\n"); + } + if (!valueFromCSVLine(line, c->household_size_5)) { + throw std::runtime_error("bad value for household_size_5 during csv parsing @todo\n"); + } + if (!valueFromCSVLine(line, c->household_size_6)) { + throw std::runtime_error("bad value for household_size_6 during csv parsing @todo\n"); } if (!valueFromCSVLine(line, c->p_interaction_susceptible_to_exposed)) { throw std::runtime_error("bad value for p_interaction_susceptible_to_exposed during csv parsing @todo\n"); @@ -176,6 +205,54 @@ std::shared_ptr read(const std::filesystem::path p, c if (!valueFromCSVLine(line, c->relative_susceptibility_80)) { throw std::runtime_error("bad value for relative_susceptibility_80 during csv parsing @todo\n"); } + if (!valueFromCSVLine(line, c->child_network_adults)) { + throw std::runtime_error("bad value for child_network_adults during csv parsing @todo\n"); + } + if (!valueFromCSVLine(line, c->elderly_network_adults)) { + throw std::runtime_error("bad value for elderly_network_adults during csv parsing @todo\n"); + } + if (!valueFromCSVLine(line, c->relative_transmission_household)) { + throw std::runtime_error("bad value for relative_transmission_household during csv parsing @todo\n"); + } + if (!valueFromCSVLine(line, c->relative_transmission_occupation)) { + throw std::runtime_error("bad value for relative_transmission_occupation during csv parsing @todo\n"); + } + if (!valueFromCSVLine(line, c->relative_transmission_random)) { + throw std::runtime_error("bad value for relative_transmission_random during csv parsing @todo\n"); + } + if (!valueFromCSVLine(line, c->mean_work_interactions_child)) { + throw std::runtime_error("bad value for mean_work_interactions_child during csv parsing @todo\n"); + } + if (!valueFromCSVLine(line, c->mean_work_interactions_adult)) { + throw std::runtime_error("bad value for mean_work_interactions_adult during csv parsing @todo\n"); + } + if (!valueFromCSVLine(line, c->mean_work_interactions_elderly)) { + throw std::runtime_error("bad value for mean_work_interactions_elderly during csv parsing @todo\n"); + } + if (!valueFromCSVLine(line, c->daily_fraction_work)) { + throw std::runtime_error("bad value for daily_fraction_work during csv parsing @todo\n"); + } + if (!valueFromCSVLine(line, c->work_network_rewire)) { + throw std::runtime_error("bad value for work_network_rewire during csv parsing @todo\n"); + } + if (!valueFromCSVLine(line, c->mean_random_interactions_0_19)) { + throw std::runtime_error("bad value for mean_random_interactions_0_19 during csv parsing @todo\n"); + } + if (!valueFromCSVLine(line, c->sd_random_interactions_0_19)) { + throw std::runtime_error("bad value for sd_random_interactions_0_19 during csv parsing @todo\n"); + } + if (!valueFromCSVLine(line, c->mean_random_interactions_20_69)) { + throw std::runtime_error("bad value for mean_random_interactions_20_69 during csv parsing @todo\n"); + } + if (!valueFromCSVLine(line, c->sd_random_interactions_20_69)) { + throw std::runtime_error("bad value for sd_random_interactions_20_69 during csv parsing @todo\n"); + } + if (!valueFromCSVLine(line, c->mean_random_interactions_70plus)) { + throw std::runtime_error("bad value for mean_random_interactions_70plus during csv parsing @todo\n"); + } + if (!valueFromCSVLine(line, c->sd_random_interactions_70plus)) { + throw std::runtime_error("bad value for sd_random_interactions_70plus during csv parsing @todo\n"); + } } else { throw std::runtime_error("failed to read the paramameter value line @todo nicer error message"); @@ -196,6 +273,7 @@ void print(exateppabm::input::config config) { fmt::print(" param_id = {}\n", config.param_id); fmt::print(" duration = {}\n", config.duration); fmt::print(" n_total = {}\n", config.n_total); + fmt::print(" n_seed_infection = {}\n", config.n_seed_infection); fmt::print(" population_0_9 = {}\n", config.population_0_9); fmt::print(" population_10_19 = {}\n", config.population_10_19); fmt::print(" population_20_29 = {}\n", config.population_20_29); @@ -205,7 +283,12 @@ void print(exateppabm::input::config config) { fmt::print(" population_60_69 = {}\n", config.population_60_69); fmt::print(" population_70_79 = {}\n", config.population_70_79); fmt::print(" population_80 = {}\n", config.population_80); - fmt::print(" n_seed_infection = {}\n", config.n_seed_infection); + fmt::print(" household_size_1 = {}\n", config.household_size_1); + fmt::print(" household_size_2 = {}\n", config.household_size_2); + fmt::print(" household_size_3 = {}\n", config.household_size_3); + fmt::print(" household_size_4 = {}\n", config.household_size_4); + fmt::print(" household_size_5 = {}\n", config.household_size_5); + fmt::print(" household_size_6 = {}\n", config.household_size_6); fmt::print(" p_interaction_susceptible_to_exposed = {}\n", config.p_interaction_susceptible_to_exposed); fmt::print(" mean_time_to_infected = {}\n", config.mean_time_to_infected); fmt::print(" sd_time_to_infected = {}\n", config.sd_time_to_infected); @@ -222,6 +305,22 @@ void print(exateppabm::input::config config) { fmt::print(" relative_susceptibility_60_69 = {}\n", config.relative_susceptibility_60_69); fmt::print(" relative_susceptibility_70_79 = {}\n", config.relative_susceptibility_70_79); fmt::print(" relative_susceptibility_80 = {}\n", config.relative_susceptibility_80); + fmt::print(" child_network_adults = {}\n", config.child_network_adults); + fmt::print(" elderly_network_adults = {}\n", config.elderly_network_adults); + fmt::print(" relative_transmission_household = {}\n", config.relative_transmission_household); + fmt::print(" relative_transmission_occupation = {}\n", config.relative_transmission_occupation); + fmt::print(" relative_transmission_random = {}\n", config.relative_transmission_occupation); + fmt::print(" mean_work_interactions_child = {}\n", config.mean_work_interactions_child); + fmt::print(" mean_work_interactions_adult = {}\n", config.mean_work_interactions_adult); + fmt::print(" mean_work_interactions_elderly = {}\n", config.mean_work_interactions_elderly); + fmt::print(" daily_fraction_work = {}\n", config.daily_fraction_work); + fmt::print(" work_network_rewire = {}\n", config.work_network_rewire); + fmt::print(" mean_random_interactions_0_19 = {}\n", config.mean_random_interactions_0_19); + fmt::print(" sd_random_interactions_0_19 = {}\n", config.sd_random_interactions_0_19); + fmt::print(" mean_random_interactions_20_69 = {}\n", config.mean_random_interactions_20_69); + fmt::print(" sd_random_interactions_20_69 = {}\n", config.sd_random_interactions_20_69); + fmt::print(" mean_random_interactions_70plus = {}\n", config.mean_random_interactions_70plus); + fmt::print(" sd_random_interactions_70plus = {}\n", config.sd_random_interactions_70plus); fmt::print("}}\n"); } diff --git a/src/exateppabm/input.h b/src/exateppabm/input.h index 331c85b..9f9bf45 100644 --- a/src/exateppabm/input.h +++ b/src/exateppabm/input.h @@ -27,6 +27,10 @@ struct config { * The population size for the simulation */ std::uint32_t n_total = 1024; + /** + * The number of individuals who should be infected at the start of the simulation + */ + std::uint32_t n_seed_infection = 1; /** * Reference size for the number of individuals within the 0_9 age demographic, used to compute ratios for population initialisation */ @@ -64,9 +68,29 @@ struct config { */ std::uint64_t population_80 = 1; /** - * The number of individuals who should be infected at the start of the simulation + * Reference size for the number of households with 1 person, used for household generation */ - std::uint32_t n_seed_infection = 1; + std::uint64_t household_size_1 = 1; + /** + * Reference size for the number of households with 2 people, used for household generation + */ + std::uint64_t household_size_2 = 1; + /** + * Reference size for the number of households with 3 people, used for household generation + */ + std::uint64_t household_size_3 = 1; + /** + * Reference size for the number of households with 4 people, used for household generation + */ + std::uint64_t household_size_4 = 1; + /** + * Reference size for the number of households with 5 people, used for household generation + */ + std::uint64_t household_size_5 = 1; + /** + * Reference size for the number of households with 6+ people, used for household generation. 6+ chosen as largest band based on ONS estimates https://www.ons.gov.uk/peoplepopulationandcommunity/birthsdeathsandmarriages/families/datasets/householdsbyhouseholdsizeregionsofenglandandgbconstituentcountries + */ + std::uint64_t household_size_6 = 1; /** * The probability of an interaction between an infected individual and a susceptible individual resulting in an infection, prior to any susceptibility modifier. * This is not directly based on a parameter in the reference model @@ -76,7 +100,7 @@ struct config { /** * The mean time in days from exposed to infected state * Default value is arbitrary - * @todo - might be equivalent to mean_time_to_syptoms, depending on how asym works in the reference model. + * @todo - might be equivalent to mean_time_to_symptoms, depending on how the more complex disease model works in the reference model. */ float mean_time_to_infected = 4; /** @@ -149,6 +173,90 @@ struct config { * Arbitrary default value */ float relative_susceptibility_80 = 1.0; + /** + * Proportion of adults to children in school networks (0-19) + */ + double child_network_adults = 0.2; + /** + * Proportion of adults to elderly in work networks for retired/elderly individuals (70+) + */ + double elderly_network_adults = 0.2; + /** + * Relative transmission rate for interactions within the household + * Arbitrary default value + */ + float relative_transmission_household = 2.0f; + /** + * Relative transmission rate for interactions within the occupation network + * Arbitrary default value + */ + float relative_transmission_occupation = 1.0f; + /** + * Relative transmission rate for interactions within the random daily network + * Arbitrary default value + */ + float relative_transmission_random = 1.0f; + /** + * Mean daily interactions at school (0-19) + * This parameter is used during small world network generation (k) + * Arbitrary default value + */ + double mean_work_interactions_child = 10; + /** + * Mean daily interactions at adult workplaces (20-69) + * This parameter is used during small world network generation (k) + * Arbitrary default value + */ + double mean_work_interactions_adult = 7.5; + /** + * Mean daily interactions at equivalent to workplace for older people (70+) + * This parameter is used during small world network generation (k) + * Arbitrary default value + */ + double mean_work_interactions_elderly = 3; + /** + * Fraction of people in work network interacted with per day + * + * This modifies the small world generation initial mean to create larger clusters, allowing for random sampling by this factor on each day + * + * Arbitrary default value + */ + double daily_fraction_work = 0.5; + /** + * Small world network rewire parameter (probability) for workplaces, in range [0, 1] + * Arbitrary default + */ + double work_network_rewire = 0.1; + /** + * Mean number of random interactions for 0-19 age individuals + * Arbitrary default value + */ + double mean_random_interactions_0_19 = 2u; + /** + * Standard deviation for the number of random interactions per day for 0-19 age individuals + * Default value is arbitrary + */ + double sd_random_interactions_0_19 = 2u; + /** + * Mean number of random interactions for 20-69 age individuals + * Arbitrary default value + */ + double mean_random_interactions_20_69 = 4u; + /** + * Standard deviation for the number of random interactions per day for 20-69 age individuals + * Default value is arbitrary + */ + double sd_random_interactions_20_69 = 4u; + /** + * Mean number of random interactions for 70+ age individuals + * Arbitrary default value + */ + double mean_random_interactions_70plus = 2u; + /** + * Standard deviation for the number of random interactions per day for 70+ age individuals + * Default value is arbitrary + */ + double sd_random_interactions_70plus = 2u; }; /** diff --git a/src/exateppabm/network.cu b/src/exateppabm/network.cu new file mode 100644 index 0000000..dbb3c4d --- /dev/null +++ b/src/exateppabm/network.cu @@ -0,0 +1,123 @@ +#include "exateppabm/network.h" + +#include +#include +#include +#include +#include +#include +#include + +#include "flamegpu/flamegpu.h" +#include "fmt/core.h" + +namespace exateppabm { +namespace network { + + +UndirectedGraph generateFullyConnectedUndirectedGraph(std::vector nodes) { + UndirectedGraph network(nodes); + if (nodes.size() <= 1) { + return network; + } + std::uint32_t N = static_cast(nodes.size()); + std::uint32_t E = (N * (N - 1)) / 2; + // network.edges.reserve(E); + // Generate edges in ascending order for the undirected graph. + for (std::uint32_t i = 0; i < N; ++i) { + for (std::uint32_t j = i + 1; j < N; ++j) { + network.addEdge(i, j); + } + } + return network; +} + +// @todo - should the initial lattice be randomised rather than sequential? +// @todo - should the graph be directed or undirected? Data structure is directed, so do pairs need rewiring together? +UndirectedGraph generateSmallWorldUndirectedGraph(std::vector nodes, std::uint32_t K, double p_rewire, std::mt19937_64 rng) { + UndirectedGraph network(nodes); + + // Return early with an empty network if 0 or 1 nodes, so no room for non-self edges + if (network.getNumVertices() <= 1) { + return network; + } + + // If K is 0 (or 1), no edges so do nothing + if (K <= 1) { + // throw std::runtime_error("@todo values - small world network K (" + std::to_string(K) + ") must be > 1"); + return network; + } else if (K == network.getNumVertices()) { + // If K == Node count, the graph is fully connected, so return a fully connected graph + return generateFullyConnectedUndirectedGraph(nodes); + } else if (K >= network.getNumVertices()) { + // Raise an exception if K is too large + throw std::runtime_error(std::string("@todo values - small world network K(") + std::to_string(K) + ") must be less than |N| (" + std::to_string(network.getNumVertices()) + ")"); + } + // If K is odd, use K-1 + if (K % 2 == 1) { + K = K - 1; + } + + // p_rewire must be between 0 and 1. + if (p_rewire < 0 || p_rewire > 1) { + throw std::runtime_error("generateSmallWorldNetwork p must be in [0, 1]. @todo include value"); + } + + // Initialise the edges as a lattice network, going K edges in each direction. + // Use signed integers to make modulo possible + std::int32_t N = static_cast(nodes.size()); + std::int32_t E = (N * K) / 2; + + // network.edges.reserve(E); + for (std::int32_t i = 0; i < N; ++i) { + for (std::int32_t j = 1; j <= static_cast(K / 2); ++j) { + // only add the positive undirected edges + std::uint32_t s = static_cast(i); + std::uint32_t d = static_cast((i + j) % N); + network.addEdge(s, d); + } + } + + // If the network is fully connected (enough edges for each node to be connected to every other node), rewiring is not needed, so adjust the probability to prevent infinite loops. + std::int32_t MAX_E = (N * (N - 1)) / 2; + if (E >= MAX_E) { + p_rewire = 0.0; + } + + // If p_rewire is 0, no need to loop over the edges + if (p_rewire >= 0.0) { + // Get a uniform dist [0, 1) + std::uniform_real_distribution p_dist(0.0, 1.0); + // Get a uniform integer distribution from [0, N) for generating new edge indices + std::uniform_int_distribution dest_dist(0, N-1); + // Randomly rewire edges + + // Take a copy of the network edges, to ensure we avoid iterator invalidation + // Only the current iterator should be removed in the undirected graph, so this is probably ok. + std::vector copyOfEdges(network.getEdges().begin(), network.getEdges().end()); + for (auto& edge : copyOfEdges) { + if (p_dist(rng) < p_rewire) { + // If the source vertex is full, do not attempt to rewire it, it would loop forever so check the next original edge + if (network.degree(edge.source) >= network.getNumVertices() - 1) { + continue; + } + + // Repeatedly generate a new destination vertex until a new edge has been found. + std::uint32_t newDest = dest_dist(rng); + // While the new dest is the source, or already in the graph, try again. + while (newDest == edge.source || network.contains(edge.source, newDest)) { + newDest = dest_dist(rng); + } + // Remove the old edge + network.removeEdge(edge.source, edge.dest); + // Add the new edge + network.addEdge(edge.source, newDest); + } + } + } + + return network; +} + +} // namespace network +} // namespace exateppabm diff --git a/src/exateppabm/network.h b/src/exateppabm/network.h new file mode 100644 index 0000000..757a6b4 --- /dev/null +++ b/src/exateppabm/network.h @@ -0,0 +1,228 @@ +#pragma once + +#include +#include +#include +#include +#include +#include + +#include "flamegpu/flamegpu.h" + +namespace exateppabm { +namespace network { + + +/** + * Struct representing an Edge + */ +class Edge { + public: + std::uint32_t source; + std::uint32_t dest; + /** + * constrcutor for emplace_back support + */ + Edge(std::uint32_t source, std::uint32_t dest) : source(source), dest(dest) {} + /** + * const operator== comparison for std::find + */ + bool operator==(const Edge& other) const { + return source == other.source && dest == other.dest; + } + /** + * operator< for use in a set + */ + bool operator<(const Edge& other) const { + if (source != other.source) { + return source < other.source; + } else { + return dest < other.dest; + } + } +}; + +/** + * class representing an undirected network, using edge lists as well suited to edge existan + * + * @todo - template the index type? Map indices to ID so a sparse matrix would make sense? + */ +class UndirectedGraph { + public: + /** + * Default ctor + */ + UndirectedGraph() { } + + /** + * Constructor, taking a vector of vertex labels + * + * @param vertexLabels vertex labels, setting the number of vertices for the network + */ + explicit UndirectedGraph(std::vector vertexLabels) : _vertexLabels(vertexLabels) { + this->_adjacencyList.resize(this->getNumVertices()); + } + + /** + * copy constructor + * + * @param other instance of UndirectedGraph to copy + */ + UndirectedGraph(const UndirectedGraph& other) : + _vertexLabels(other._vertexLabels), + _edges(other._edges), + _adjacencyList(other._adjacencyList) {} + + /** + * Get the number of nodes + * @return number of vertices/nodes + */ + std::uint32_t getNumVertices() const { return this->_vertexLabels.size(); } + /** + * Get the number of edges + * @return number of edges + */ + std::uint32_t getNumEdges() const { return this->_edges.size(); } + /** + * Check if an edge exists within the undirected graph + * @param u vertex index for one end of the edge + * @param v vertex index for the other end of the edge + * @return bool indicating if edge exists or not + */ + bool contains(std::uint32_t u, std::uint32_t v) const { + return this->_adjacencyList[u].count(v) > 0; + } + /** + * Check if an edge exists within the undirected graph + * @param edge to check for existance of + * @return bool indicating if edge exists or not + */ + bool contains(Edge edge) const { + return this->_adjacencyList[edge.source].count(edge.dest) > 0; + } + /** + * Add an edge to the undirected graph + * @param u vertex index for one end of the edge + * @param v vertex index for the other end of the edge + */ + void addEdge(std::uint32_t u, std::uint32_t v) { + if (u > this->getNumVertices() || v > this->getNumVertices()) { + throw std::runtime_error("Invalid u or v for UndirectedGraph::getNumVertices. @todo better error"); + } + // Undirected, so add the edge in both directions + this->_adjacencyList[u].insert(v); + this->_adjacencyList[v].insert(u); + // Only store the ascending ordered edge in the vector of edges + if (u < v) { + this->_edges.emplace_back(u, v); + } else { + this->_edges.emplace_back(v, u); + } + } + /** + * Add an edge to the undirected graph + * @param edge edge to add + */ + void addEdge(Edge edge) { + return this->addEdge(edge.source, edge.dest); + } + /** + * Remove an edge from the undirected graph + * @param edge edge to remove + */ + void removeEdge(Edge edge) { + this->_adjacencyList[edge.source].erase(edge.dest); + this->_adjacencyList[edge.dest].erase(edge.source); + // Remove the edge from the edges_ vector + auto it = std::remove_if(this->_edges.begin(), this->_edges.end(), [&edge](const Edge& e) { return e == edge; }); + this->_edges.erase(it, this->_edges.end()); + } + /** + * Remove an edge from the undirected graph + * @param u vertex index for one end of the edge + * @param v vertex index for the other end of the edge + */ + void removeEdge(std::uint32_t u, std::uint32_t v) { + return this->removeEdge({u, v}); + } + + /** + * Get the degree for a vertex. this is in-out as undirected. + * + * @return vertex degree + */ + std::uint32_t degree(std::uint32_t v) { + return static_cast(this->_adjacencyList[v].size()); + } + + /** + * Get the vertex indices of neighbouring vertices + * @param v vertex index + * @return unordered set of neighbouring vertices to v + */ + const std::unordered_set& getNeighbours(std::uint32_t v) const { + return this->_adjacencyList[v]; + } + + /** + * Get the vector of edges + * + * @return vector of edges + */ + const std::vector& getEdges() const { return this->_edges; } + + /** + * Get the vertex labels, which are probably 1-indexed FLAME GPU id_t (0 is unset value) + * + * @return vector of vertex labels + */ + const std::vector& getVertexLabels() const { return this->_vertexLabels; } + + /** + * Get the a vertex label + * + * @param v vertex index + * @return label for vertex + */ + flamegpu::id_t getVertexLabel(std::uint32_t v) const { return this->_vertexLabels.at(v); } + + private: + /** + * Vector of nodes, which provides the mapping from the network-index based edges back to agent IDs + */ + std::vector _vertexLabels; + + + /** + * Vector of edge objects, each containing the directed source-dest pair. These are 0 indexed within this network. + */ + std::vector _edges; + + /** + * Vector of unordered sets for fast lookup of edge existence. + */ + std::vector> _adjacencyList; +}; + +/** + * Generate a fully connected undirected network given a vector of node labels + * @param nodes vector of nodes (agent Id's) to include in this network. IDs may not be sequential and should be one indexed. + * @return network struct for a fully connected network + */ +UndirectedGraph generateFullyConnectedUndirectedGraph(std::vector nodes); + +/** + * Generate a watts strogatz small world (undirected) network for a given set of nodes (agent ID's, 1 indexed, not sequential), K & p. + * + * If k is odd, then k-1 neighboours are created + * + * @param nodes vector of nodes (i.e. agent id's) to include in this network. IDs may not be sequential and should be 1 indexed. + * @param K the degree of each node + * @param p rewire probability + * @return a network struct for the generate small world network + */ +UndirectedGraph generateSmallWorldUndirectedGraph(std::vector nodes, std::uint32_t K, double p, std::mt19937_64 rng); + + +} // namespace network +} // namespace exateppabm diff --git a/src/exateppabm/output.cu b/src/exateppabm/output.cu index f6167f4..95382a8 100644 --- a/src/exateppabm/output.cu +++ b/src/exateppabm/output.cu @@ -7,6 +7,7 @@ #include "exateppabm/output/OutputFile.h" #include "exateppabm/output/TimeSeriesFile.h" +#include "exateppabm/output/PerIndividualFile.h" #include "exateppabm/person.h" #include "exateppabm/population.h" #include "exateppabm/demographics.h" @@ -26,6 +27,9 @@ std::filesystem::path _outputDirectory; // Object representing the time series output file std::unique_ptr _timeSeriesFile = nullptr; +// Object representing the per individual file +std::unique_ptr _perIndividualFile = nullptr; + } // namespace /** @@ -86,13 +90,6 @@ FLAMEGPU_STEP_FUNCTION(output_step) { observations.total_infected += totalInfectedPerDemographic[i]; } - /* - // @todo - not currently used direct agent data iteration - for (const auto& person : population) { - auto infected = person.getVariable(exateppabm::person::v::INFECTION_STATE); - } - */ - // Append this steps' data to the namespace-scoped data structure _timeSeriesFile->appendObservations(observations); } @@ -107,8 +104,35 @@ FLAMEGPU_EXIT_FUNCTION(output_exit) { _timeSeriesFile->close(); } + +FLAMEGPU_EXIT_FUNCTION(output_exit_per_individual) { + // Collect per agent data + // Get a handle to the person agent host api object + auto personAgent = FLAMEGPU->agent(exateppabm::person::NAME, exateppabm::person::states::DEFAULT); + // Get a handle to the person agent population on the host + flamegpu::DeviceAgentVector population = personAgent.getPopulationData(); + for (const auto& person : population) { + exateppabm::output::PerIndividualFile::Person personData = {}; + personData.id = static_cast(person.getVariable(person::v::ID)); + personData.age_group = person.getVariable(person::v::AGE_DEMOGRAPHIC); + personData.occupation_network = person.getVariable(person::v::WORKPLACE_IDX); + personData.house_no = person.getVariable(person::v::HOUSEHOLD_IDX); + personData.infection_count = person.getVariable(person::v::INFECTION_COUNT); + _perIndividualFile->appendPerson(personData); + } + + + // Write the per individual data to disk + // Open the file handle + _perIndividualFile->open(); + // Write data to the opened file + _perIndividualFile->write(); + // Close the file handle + _perIndividualFile->close(); +} + // @todo - may need to split this due to order of execution within init/step/exit funcs, if any others exist. -void define(flamegpu::ModelDescription& model, const std::filesystem::path outputDirectory) { +void define(flamegpu::ModelDescription& model, const std::filesystem::path outputDirectory, const bool individualFile) { // Store the output directory for access in the FLAME GPU exit function (and init?) // This will want refactoring for ensembles _outputDirectory = outputDirectory; @@ -121,6 +145,12 @@ void define(flamegpu::ModelDescription& model, const std::filesystem::path outpu model.addStepFunction(output_step); // Add the exit function to the model model.addExitFunction(output_exit); + + // optionally prepare for per individual file output + if (individualFile) { + _perIndividualFile = std::make_unique(_outputDirectory); + model.addExitFunction(output_exit_per_individual); + } } } // namespace output diff --git a/src/exateppabm/output.h b/src/exateppabm/output.h index eed4af8..d75ac1f 100644 --- a/src/exateppabm/output.h +++ b/src/exateppabm/output.h @@ -9,6 +9,7 @@ #include "output/OutputFile.h" #include "output/TimeSeriesFile.h" +#include "output/PerIndividualFile.h" namespace exateppabm { namespace output { @@ -22,8 +23,9 @@ namespace output { * * @param model flamegpu2 model description object to mutate * @param outputDirectory path to directory for file output + * @param individualFile if the per individual file should be written or not */ -void define(flamegpu::ModelDescription& model, std::filesystem::path outputDirectory); +void define(flamegpu::ModelDescription& model, std::filesystem::path outputDirectory, bool individualFile); } // namespace output } // namespace exateppabm diff --git a/src/exateppabm/output/PerIndividualFile.cu b/src/exateppabm/output/PerIndividualFile.cu new file mode 100644 index 0000000..39d3e35 --- /dev/null +++ b/src/exateppabm/output/PerIndividualFile.cu @@ -0,0 +1,48 @@ +#include "PerIndividualFile.h" + +#include + +#include + +namespace exateppabm { +namespace output { + +PerIndividualFile::PerIndividualFile(std::filesystem::path directory) : OutputFile(directory / PerIndividualFile::DEFAULT_FILENAME) { } + +PerIndividualFile::~PerIndividualFile() { } + +void PerIndividualFile::reset(const std::uint32_t n_total) { + // Reset the observations vector + this->_perPerson = std::vector(); + this->_perPerson.reserve(n_total); +} + +void PerIndividualFile::appendPerson(const Person person) { + // @todo = ensure _perPerson has been initialised? + this->_perPerson.push_back(person); +} + +bool PerIndividualFile::write() { + if (!this->_handle) { + this->open(); + } + + // Print to the file handle + fmt::print(_handle, "ID,age_group,occupation_network,house_no,infection_count\n"); + for (const auto& person : _perPerson) { + fmt::print( + _handle, + "{},{},{},{},{}\n", + person.id, + person.age_group, + person.occupation_network, + person.house_no, + person.infection_count); + } + + fmt::print("Per individual data written to {}\n", std::filesystem::absolute(this->_filepath).c_str()); + return true; +} + +} // namespace output +} // namespace exateppabm diff --git a/src/exateppabm/output/PerIndividualFile.h b/src/exateppabm/output/PerIndividualFile.h new file mode 100644 index 0000000..dee22be --- /dev/null +++ b/src/exateppabm/output/PerIndividualFile.h @@ -0,0 +1,92 @@ +#pragma once + +#include "OutputFile.h" + +#include +#include + +#include "exateppabm/demographics.h" +#include "exateppabm/person.h" + +namespace exateppabm { +namespace output { + +/** + * Class for representing the per individual file + */ +class PerIndividualFile : public OutputFile { + public: + // Forward decl nested struct + struct Person; + + /** + * Constructor setting the path for file output. + * @param directory parent directory for the file to be stored in (using the default filename) + */ + explicit PerIndividualFile(std::filesystem::path directory); + + /** + * Dtor + */ + ~PerIndividualFile(); + + /** + * Reset the objects internal data structures, and pre-allocate memory for storing data for the current simulation + * @param n_total number of people in total + */ + void reset(std::uint32_t n_total); + + /** + * Append a set of observations for a single time point to the internal data structure + * + * @param observations observations from a single time point + */ + void appendPerson(Person observations); + + /** + * Write the contents of the time series file to disk at the configured path + * + * @return bool indicating success + */ + bool write(); + + /** + * Structure for data observed at a single point in time within the time series + */ + struct Person { + /** + * ID for the person + */ + std::uint32_t id = 0; + /** + * Age group + */ + exateppabm::demographics::AgeUnderlyingType age_group = 0; + /** + * Work network index + */ + std::uint32_t occupation_network = 0; + /** + * Household index + */ + std::uint32_t house_no = 0; + /** + * Cumulative number of times this individual was infected + */ + std::uint32_t infection_count = 0; + }; + + private: + /** + * Default filename for output + * @todo - factor in the run index for ensembles? + */ + constexpr static char DEFAULT_FILENAME[] = "individual_file.csv"; + /** + * Private member containing the observation data + */ + std::vector _perPerson; +}; + +} // namespace output +} // namespace exateppabm diff --git a/src/exateppabm/output/TimeSeriesFile.h b/src/exateppabm/output/TimeSeriesFile.h index ac1bd88..eb0793a 100644 --- a/src/exateppabm/output/TimeSeriesFile.h +++ b/src/exateppabm/output/TimeSeriesFile.h @@ -120,6 +120,7 @@ class TimeSeriesFile : public OutputFile { private: /** * Default filename for output + * @todo - factor in the run index for ensembles? */ constexpr static char DEFAULT_FILENAME[] = "timeseries.csv"; /** diff --git a/src/exateppabm/person.cu b/src/exateppabm/person.cu index b4d3968..b0b917c 100644 --- a/src/exateppabm/person.cu +++ b/src/exateppabm/person.cu @@ -1,42 +1,378 @@ #include "exateppabm/person.h" +#include + +#include + #include "exateppabm/disease/SEIR.h" #include "exateppabm/demographics.h" +#include "exateppabm/population.h" // @todo - replace with workpalce when abstracted + namespace exateppabm { namespace person { /** - * Agent function for person agents to emit their public information, i.e. infection status + * Agent function for person agents to emit their public information, i.e. infection status to their household + */ +FLAMEGPU_AGENT_FUNCTION(emitHouseholdStatus, flamegpu::MessageNone, flamegpu::MessageBucket) { + // Households of 1 don't need to do any messaging, there is no one to infect + std::uint8_t householdSize = FLAMEGPU->getVariable(v::HOUSEHOLD_SIZE); + if (householdSize > 1) { + // output public properties to bucket message, keyed by household + // Agent ID to avoid self messaging + FLAMEGPU->message_out.setVariable(person::message::household_status::ID, FLAMEGPU->getVariable(person::v::ID)); + + // Household index + // @todo - typedef or using statement for the household index type? + std::uint32_t householdIdx = FLAMEGPU->getVariable(v::HOUSEHOLD_IDX); + FLAMEGPU->message_out.setVariable(v::HOUSEHOLD_IDX, householdIdx); + + FLAMEGPU->message_out.setVariable(v:: + INFECTION_STATE, FLAMEGPU->getVariable(v::INFECTION_STATE)); + FLAMEGPU->message_out.setVariable(v::AGE_DEMOGRAPHIC, FLAMEGPU->getVariable(v::AGE_DEMOGRAPHIC)); + // Set the message key, the house hold idx for bucket messaging @Todo + FLAMEGPU->message_out.setKey(householdIdx); + } + return flamegpu::ALIVE; +} + +/** + * Very naive agent interaction for infection spread via house hold contact + * + * Agents iterate messages from their household members, potentially becoming infected + * + * @todo - refactor this somewhere else? + * @todo - add per network behaviours? + */ +FLAMEGPU_AGENT_FUNCTION(interactHousehold, flamegpu::MessageBucket, flamegpu::MessageNone) { + // Get my ID to avoid self messages + const flamegpu::id_t id = FLAMEGPU->getVariable(person::v::ID); // FLAMEGPU->getID(); + + // Get the probability of infection + float p_s2e = FLAMEGPU->environment.getProperty("p_interaction_susceptible_to_exposed"); + // Scale it for household interactions + p_s2e *= FLAMEGPU->environment.getProperty("relative_transmission_household"); + + // Get my household index + auto householdIdx = FLAMEGPU->getVariable(v::HOUSEHOLD_IDX); + + // Get my age demographic + auto demographic = FLAMEGPU->getVariable(v::AGE_DEMOGRAPHIC); + // Get my susceptibility modifier and modify it. + float relativeSusceptibility = FLAMEGPU->environment.getProperty("relative_susceptibility_per_demographic", demographic); + // Scale the probability of transmission + p_s2e *= relativeSusceptibility; + + // Check if the current individual is susceptible to being infected + auto infectionState = FLAMEGPU->getVariable(v::INFECTION_STATE); + + // @todo - this will need to change for contact tracing, the message interaction will need to occur regardless. + if (infectionState == disease::SEIR::Susceptible) { + // Variable to store the duration of the exposed phase (if exposed) + float stateDuration = 0.f; + + // Iterate messages from anyone within the household + for (const auto &message : FLAMEGPU->message_in(householdIdx)) { + // Ignore self messages (can't infect oneself) + if (message.getVariable(message::household_status::ID) != id) { + // Ignore messages from other households + if (message.getVariable(v::HOUSEHOLD_IDX) == householdIdx) { + // Check if the other agent is infected + if (message.getVariable(v::INFECTION_STATE) == disease::SEIR::InfectionState::Infected) { + // Roll a dice + float r = FLAMEGPU->random.uniform(); + if (r < p_s2e) { + // I have been exposed + infectionState = disease::SEIR::InfectionState::Exposed; + // Generate how long until I am infected + float mean = FLAMEGPU->environment.getProperty("mean_time_to_infected"); + float sd = FLAMEGPU->environment.getProperty("sd_time_to_infected"); + stateDuration = (FLAMEGPU->random.normal() * sd) + mean; + // @todo - for now only any exposure matters. This may want to change when quantity of exposure is important? + // Increment the infection counter for this individual + FLAMEGPU->setVariable(v::INFECTION_COUNT, FLAMEGPU->getVariable(v::INFECTION_COUNT) + 1); + break; + } + } + } + } + } + // If newly exposed, store the value in global device memory. + if (infectionState == disease::SEIR::InfectionState::Exposed) { + FLAMEGPU->setVariable(v::INFECTION_STATE, infectionState); + FLAMEGPU->setVariable(person::v::INFECTION_STATE_DURATION, stateDuration); + } + } + + return flamegpu::ALIVE; +} + +/** + * Agent function for person agents to emit their public information, i.e. infection status, to their workplace colleagues */ -FLAMEGPU_AGENT_FUNCTION(emitStatus, flamegpu::MessageNone, flamegpu::MessageSpatial2D) { - // output public properties to spatial message - FLAMEGPU->message_out.setVariable(person::message::v::STATUS_ID, FLAMEGPU->getID()); - // Location which is also used for comms - FLAMEGPU->message_out.setVariable(v::x, FLAMEGPU->getVariable(v::x)); - FLAMEGPU->message_out.setVariable(v::y, FLAMEGPU->getVariable(v::y)); - // FLAMEGPU->message_out.setVariable(v::z, FLAMEGPU->getVariable(v::z)); - - // And their +FLAMEGPU_AGENT_FUNCTION(emitWorkplaceStatus, flamegpu::MessageNone, flamegpu::MessageArray) { + // output public properties to bucket message, keyed by workplace + // Agent ID to avoid self messaging + flamegpu::id_t id = FLAMEGPU->getVariable(person::v::ID); + FLAMEGPU->message_out.setVariable(person::message::workplace_status::ID, id); + + // workplace index + // @todo - typedef or using statement for the workplace index type? + std::uint32_t workplaceIdx = FLAMEGPU->getVariable(v::WORKPLACE_IDX); + FLAMEGPU->message_out.setVariable(v::WORKPLACE_IDX, workplaceIdx); + FLAMEGPU->message_out.setVariable(v:: INFECTION_STATE, FLAMEGPU->getVariable(v::INFECTION_STATE)); FLAMEGPU->message_out.setVariable(v::AGE_DEMOGRAPHIC, FLAMEGPU->getVariable(v::AGE_DEMOGRAPHIC)); + // Set the message key, the house hold idx for bucket messaging @Todo + FLAMEGPU->message_out.setIndex(id); return flamegpu::ALIVE; } /** - * Very naive agent interaction for infection spread via "contact" + * Very naive agent interaction for infection spread via workplace contact + * + * Agents iterate messages from their workplace members, potentially becoming infected * - * Agents iterate messages from their local neighbours. - * If their neighbour was infected, RNG sample to determine if "i" become infected. + * @todo - refactor this somewhere else? + * @todo - add per network behaviours? */ -FLAMEGPU_AGENT_FUNCTION(interact, flamegpu::MessageSpatial2D, flamegpu::MessageNone) { +FLAMEGPU_AGENT_FUNCTION(interactWorkplace, flamegpu::MessageArray, flamegpu::MessageNone) { + // Get my ID + const flamegpu::id_t id = FLAMEGPU->getVariable(person::v::ID); + + // Get my workplace degree, if 0 return. + const std::uint32_t degree = FLAMEGPU->getVariable(person::v::WORKPLACE_OUT_DEGREE); + if (degree == 0) { + return flamegpu::ALIVE; + } + + // Get the probability of interaction within the workplace + float p_daily_fraction_work = FLAMEGPU->environment.getProperty("daily_fraction_work"); + + // Get the probability of an interaction leading to exposure + float p_s2e = FLAMEGPU->environment.getProperty("p_interaction_susceptible_to_exposed"); + // Scale it for workplace interactions + p_s2e *= FLAMEGPU->environment.getProperty("relative_transmission_occupation"); + + // Get my workplace/network index + auto workplaceIdx = FLAMEGPU->getVariable(v::WORKPLACE_IDX); + + // Get my age demographic + auto demographic = FLAMEGPU->getVariable(v::AGE_DEMOGRAPHIC); + // Get my susceptibility modifier and modify it. + float relativeSusceptibility = FLAMEGPU->environment.getProperty("relative_susceptibility_per_demographic", demographic); + // Scale the probability of transmission for my age demographic + p_s2e *= relativeSusceptibility; + + // Check if the current individual is susceptible to being infected + auto infectionState = FLAMEGPU->getVariable(v::INFECTION_STATE); + + // @todo - this will need to change for contact tracing, the message interaction will need to occur regardless. + if (infectionState == disease::SEIR::Susceptible) { + // Variable to store the duration of the exposed phase (if exposed) + float stateDuration = 0.f; + // Iterate my downstream neighbours (the graph is undirected, so no need to iterate in and out + auto workplaceGraph = FLAMEGPU->environment.getDirectedGraph("WORKPLACE_DIGRAPH"); + std::uint32_t myVertexIndex = workplaceGraph.getVertexIndex(id); + for (auto& edge : workplaceGraph.outEdges(myVertexIndex)) { + // Get the neighbour vertex index and then neighbours agent Id from the edge + std::uint32_t otherIndex = edge.getEdgeDestination(); + flamegpu::id_t otherID = workplaceGraph.getVertexID(otherIndex); + // if the ID is not self and not unset + if (otherID != id && otherID != flamegpu::ID_NOT_SET) { + // Get the message handle + const auto &message = FLAMEGPU->message_in.at(otherID); + // Ignore messages from other workplaces + if (message.getVariable(v::WORKPLACE_IDX) == workplaceIdx) { + // roll a dice to determine if this interaction should occur this day + if (FLAMEGPU->random.uniform() < p_daily_fraction_work) { + // Check if the other agent is infected + if (message.getVariable(v::INFECTION_STATE) == disease::SEIR::InfectionState::Infected) { + // Roll a dice + float r = FLAMEGPU->random.uniform(); + if (r < p_s2e) { + // I have been exposed + infectionState = disease::SEIR::InfectionState::Exposed; + // Generate how long until I am infected + float mean = FLAMEGPU->environment.getProperty("mean_time_to_infected"); + float sd = FLAMEGPU->environment.getProperty("sd_time_to_infected"); + stateDuration = (FLAMEGPU->random.normal() * sd) + mean; + // @todo - for now only any exposure matters. This may want to change when quantity of exposure is important? + // Increment the infection counter for this individual + FLAMEGPU->setVariable(v::INFECTION_COUNT, FLAMEGPU->getVariable(v::INFECTION_COUNT) + 1); + break; + } + } + } + } + } + } + // If newly exposed, store the value in global device memory. + if (infectionState == disease::SEIR::InfectionState::Exposed) { + FLAMEGPU->setVariable(v::INFECTION_STATE, infectionState); + FLAMEGPU->setVariable(person::v::INFECTION_STATE_DURATION, stateDuration); + } + } + + return flamegpu::ALIVE; +} + + +/** + * Update the per-day random daily interaction network + * + * Individuals have a number of daily interactions, which is the same for the duration of th simulation. + * + * In serial on the host, generate a vector of interactions + * + * @todo - more than one random interaction per day + * @todo - prevent duplicate interactions + * @todo - more performant CPU implementation + * @todo - GPU implementation of this (stable matching submodel?) + */ +FLAMEGPU_HOST_FUNCTION(updateRandomDailyNetworkIndices) { + // Get the current population of person agents + auto personAgent = FLAMEGPU->agent(exateppabm::person::NAME, exateppabm::person::states::DEFAULT); + flamegpu::DeviceAgentVector population = personAgent.getPopulationData(); + + // @todo implement the following + // Get the sum of the per-agent random interaction count, so we can reserve a large enough vector + std::uint64_t randomInteractionCountSum = FLAMEGPU->environment.getProperty("RANDOM_INTERACTION_COUNT_SUM"); + + // If there are none, return. + if (randomInteractionCountSum == 0) { + return; + } + + // Build a vector of interactions, initialised to contain the id of each agent as many times as they require + // Number is fixed for the simulation, so only allocate and initialise once via a method-static variable + static std::vector randomInteractionIdVector; + // If empty, initialise to contain each agent's ID the required number of times. + // This will only trigger a device to host copy for the first pass. + // This could be optimised away, with a namespace scoped static and initialised during initial host agent population, at the cost of readability + if (randomInteractionIdVector.empty()) { + randomInteractionIdVector.resize(randomInteractionCountSum); + size_t first = 0; + size_t last = 0; + for (const auto &person : population) { + flamegpu::id_t id = person.getVariable(person::v::ID); + std::uint32_t count = person.getVariable(person::v::RANDOM_INTERACTION_COUNT_TARGET); + last = first + count; + if (last > randomInteractionCountSum) { + throw std::runtime_error("Too many random interactions being generated. @todo better error"); + } + std::fill(randomInteractionIdVector.begin() + first, randomInteractionIdVector.begin() + last, id); + first = last; + } + } + + // Shuffle the vector of agent id's. + // @todo - adjust how this RNG is seeded, this is far from ideal. + auto rng = std::mt19937_64(FLAMEGPU->random.uniform()); + std::shuffle(randomInteractionIdVector.begin(), randomInteractionIdVector.end(), rng); + + // Update agent data containing today's random interactions + // This will trigger Host to Device copies + // Iterate the pairs of ID which form the interactions, updating both agents' data for each interaction. + // @todo - avoid self-interactions and duplicate interactions. + // The number of interaction pairs is the total number of id's divided by 2, and rounded down. + std::uint64_t interactionCount = static_cast(std::floor(randomInteractionCountSum / 2.0)); + for (std::uint64_t interactionIdx = 0; interactionIdx < interactionCount; ++interactionIdx) { + // Get the indices within the vector of agent id's. An mdspan would be nice if not c++17 + size_t aIndex = (interactionIdx * 2); + size_t bIndex = aIndex + 1; + + // Get the agent ID's + flamegpu::id_t aID = randomInteractionIdVector[aIndex]; + flamegpu::id_t bID = randomInteractionIdVector[bIndex]; + + // Assuming that agents are in-oder (which is not guaranteed!) + // @todo - this validation could be done once per iteration, not once per interaction? + // @todo switch to a sparse data structure, indexed by agent ID? Not ideal for mem access, but simpler and less validation required. + + // Get a handle to each agent from the population vector. + // Agents ID's are 1 indexed + auto aAgent = population[aID - 1]; + auto bAgent = population[bID - 1]; + // Raise an exception if the agents are out of order. This should not be the case as agent birth, death and sorting should not be included in this model. + if (aAgent.getVariable(person::v::ID) != aID || bAgent.getVariable(person::v::ID) != bID) { + throw std::runtime_error("Agent ID does not match expected agent ID in updateRandomDailyNetworkIndices. @todo"); + } + + // @todo - avoid self-interactions by looking ahead and swapping? + // @todo - avoid repeated interactions by looking ahead and swapping? + + // Add the interaction to the list of interactions agent a + std::uint32_t aInteractionIdx = aAgent.getVariable(person::v::RANDOM_INTERACTION_COUNT); + // aAgent.setVariable(person::v::RANDOM_INTERACTION_PARTNERS, aInteractionIdx, bID); + aAgent.setVariable(person::v::RANDOM_INTERACTION_PARTNERS, aInteractionIdx, bID); + aAgent.setVariable(person::v::RANDOM_INTERACTION_COUNT, aInteractionIdx + 1); + + + // Add the interaction to the list of interactions agent b + std::uint32_t bInteractionIdx = bAgent.getVariable(person::v::RANDOM_INTERACTION_COUNT); + // bAgent.setVariable(person::v::RANDOM_INTERACTION_PARTNERS, bInteractionIdx, aID); + bAgent.setVariable(person::v::RANDOM_INTERACTION_PARTNERS, bInteractionIdx, aID); + bAgent.setVariable(person::v::RANDOM_INTERACTION_COUNT, bInteractionIdx + 1); + } + + // setVariable(name, idx, value) in flamegpu 2.0.0-rc.2 and earlier contains a bug where the single element setting variant does not trigger data synchronisation on the end of the host function + // This has been fixed by https://github.com/FLAMEGPU/FLAMEGPU2/pull/1266, which will be part of flamegpu 2.0.0-rc.3 + // As pre release versions are strings, this is a more complex if statement than it would ideally be + // Can macro this out if FLAME GPU is not 2.0.0 +#if defined(FLAMEGPU_VERSION) && FLAMEGPU_VERSION == 2000000 + // then at runtime we can only trigger the workaround for flamegpu2 versions rc2 and before + if (strcmp(flamegpu::VERSION_PRERELEASE, "rc.2") == 0 + || strcmp(flamegpu::VERSION_PRERELEASE, "rc.1") == 0 + || strcmp(flamegpu::VERSION_PRERELEASE, "rc") == 0 + || strcmp(flamegpu::VERSION_PRERELEASE, "alpha.2") == 0 + || strcmp(flamegpu::VERSION_PRERELEASE, "alpha.1") == 0 + || strcmp(flamegpu::VERSION_PRERELEASE, "alpha") == 0) { + for (auto person : population) { + person.setVariable(person::v::RANDOM_INTERACTION_PARTNERS, person.getVariable(person::v::RANDOM_INTERACTION_PARTNERS)); + } + } +#endif +} + +/** + * Agent function for person agents to emit their public information, i.e. infection status, for random daily network colleagues. This is put into a bucket key'd by the agent's ID. + * + * @note the bucket message is from 1...N, rather than 0 to match ID index. + */ +FLAMEGPU_AGENT_FUNCTION(emitRandomDailyNetworkStatus, flamegpu::MessageNone, flamegpu::MessageArray) { + // output public properties to array message + // FLAMEGPU->message_out.setVariable(person::message::random_daily_status::ID, FLAMEGPU->getVariable(person::v::ID)); + + FLAMEGPU->message_out.setVariable(v:: + INFECTION_STATE, FLAMEGPU->getVariable(v::INFECTION_STATE)); + FLAMEGPU->message_out.setVariable(v::AGE_DEMOGRAPHIC, FLAMEGPU->getVariable(v::AGE_DEMOGRAPHIC)); + + // Set the message array message index to the agent's id. + FLAMEGPU->message_out.setIndex(FLAMEGPU->getVariable(person::v::ID)); + return flamegpu::ALIVE; +} + +/** + * Very naive agent interaction for infection spread via random daily network contact + * + * Agents iterate messages from others within their random daily network, potentially becoming infected + * + * @todo - refactor this somewhere else? + * @todo - add per network behaviours? + * @todo - leverage a FLAME GPU 2 network structure, when they are mutable per day. + */ +FLAMEGPU_AGENT_FUNCTION(interactRandomDailyNetwork, flamegpu::MessageArray, flamegpu::MessageNone) { // Get my ID to avoid self messages - const flamegpu::id_t id = FLAMEGPU->getID(); + const flamegpu::id_t id = FLAMEGPU->getVariable(person::v::ID); // Get the probability of infection float p_s2e = FLAMEGPU->environment.getProperty("p_interaction_susceptible_to_exposed"); + // Scale it for random daily interactions + p_s2e *= FLAMEGPU->environment.getProperty("relative_transmission_random"); // Get my age demographic auto demographic = FLAMEGPU->getVariable(v::AGE_DEMOGRAPHIC); @@ -48,19 +384,20 @@ FLAMEGPU_AGENT_FUNCTION(interact, flamegpu::MessageSpatial2D, flamegpu::MessageN // Check if the current individual is susceptible to being infected auto infectionState = FLAMEGPU->getVariable(v::INFECTION_STATE); + // @todo - this will need to change for contact tracing, the message interaction will need to occur regardless. if (infectionState == disease::SEIR::Susceptible) { - // Agent position - float agent_x = FLAMEGPU->getVariable(v::x); - float agent_y = FLAMEGPU->getVariable(v::y); - // float agent_z = FLAMEGPU->getVariable(v::z); - // Variable to store the duration of the exposed phase (if exposed) float stateDuration = 0.f; - // Iterate messages from anyone within my spatial neighbourhood (i.e. cuboid not sphere) - for (const auto &message : FLAMEGPU->message_in(agent_x, agent_y)) { - // Ignore self messages (can't infect oneself) - if (message.getVariable(message::v::STATUS_ID) != id) { + // For each interaction this agent is set to perform + const std::uint32_t randomInteractionCount = FLAMEGPU->getVariable(person::v::RANDOM_INTERACTION_COUNT); + for (std::uint32_t randomInteractionIdx = 0; randomInteractionIdx < randomInteractionCount; ++randomInteractionIdx) { + // Get the (next) ID of the interaction partner. + flamegpu::id_t otherID = FLAMEGPU->getVariable(person::v::RANDOM_INTERACTION_PARTNERS, randomInteractionIdx); + // if the ID is not self and not unset + if (otherID != id && otherID != flamegpu::ID_NOT_SET) { + // Get the message handle + const auto &message = FLAMEGPU->message_in.at(otherID); // Check if the other agent is infected if (message.getVariable(v::INFECTION_STATE) == disease::SEIR::InfectionState::Infected) { // Roll a dice @@ -73,6 +410,8 @@ FLAMEGPU_AGENT_FUNCTION(interact, flamegpu::MessageSpatial2D, flamegpu::MessageN float sd = FLAMEGPU->environment.getProperty("sd_time_to_infected"); stateDuration = (FLAMEGPU->random.normal() * sd) + mean; // @todo - for now only any exposure matters. This may want to change when quantity of exposure is important? + // Increment the infection counter for this individual + FLAMEGPU->setVariable(v::INFECTION_COUNT, FLAMEGPU->getVariable(v::INFECTION_COUNT) + 1); break; } } @@ -85,23 +424,39 @@ FLAMEGPU_AGENT_FUNCTION(interact, flamegpu::MessageSpatial2D, flamegpu::MessageN } } + // reset the agent's number of interactions to 0 in advance of the next day. + // This is expensive, but the D2H copy would get triggered anyway if attempting to update on the host anyway. + // @todo - a more performant way to do the random daily interaction network would be good. + FLAMEGPU->setVariable(person::v::RANDOM_INTERACTION_COUNT, 0u); + return flamegpu::ALIVE; } -void define(flamegpu::ModelDescription& model, const exateppabm::input::config& params, const float width, const float interactionRadius) { +void define(flamegpu::ModelDescription& model, const exateppabm::input::config& params) { // Define related model environment properties (@todo - abstract these somewhere more appropriate at a later date) flamegpu::EnvironmentDescription env = model.Environment(); - env.newProperty("INFECTION_INTERACTION_RADIUS", interactionRadius); - // Define an infection probabiltiy. @todo this should be from the config file. + + // @todo this should probably be refactored elsewhere (although used in this file currently) + // Define the base probability of being exposed if interacting with an infected individual env.newProperty("p_interaction_susceptible_to_exposed", params.p_interaction_susceptible_to_exposed); + // define the per interaction type scale factors as env var + env.newProperty("relative_transmission_household", params.relative_transmission_household); + env.newProperty("relative_transmission_occupation", params.relative_transmission_occupation); + env.newProperty("relative_transmission_random", params.relative_transmission_random); + + // env var for the fraction of people in the same work network to interact with + env.newProperty("daily_fraction_work", params.daily_fraction_work); // Define the agent type - flamegpu::AgentDescription agent = model.newAgent(exateppabm::person::NAME); + flamegpu::AgentDescription agent = model.newAgent(person::NAME); // Define states - agent.newState(exateppabm::person::states::DEFAULT); + agent.newState(person::states::DEFAULT); // Define variables + // Custom ID, as getID cannot be relied upon to start at 1. + agent.newVariable(person::v::ID, flamegpu::ID_NOT_SET); + // disease related variables // @todo - define this in disease/ call a disease::SEIR::define_person() like method? agent.newVariable(person::v::INFECTION_STATE, disease::SEIR::Susceptible); @@ -110,53 +465,168 @@ void define(flamegpu::ModelDescription& model, const exateppabm::input::config& // Time until next state change? Defaults to the simulation duration + 1. agent.newVariable(person::v::INFECTION_STATE_DURATION, params.duration + 1); + // Integer count for the number of times infected, defaults to 0 + agent.newVariable(person::v::INFECTION_COUNT, 0u); + // age demographic // @todo make this an enum, and update uses of it, but flame's templating disagrees? - agent.newVariable(exateppabm::person::v::AGE_DEMOGRAPHIC, exateppabm::demographics::Age::AGE_0_9); + agent.newVariable(person::v::AGE_DEMOGRAPHIC, exateppabm::demographics::Age::AGE_0_9); + + // Household network variables. @todo - refactor to a separate network location? + agent.newVariable(person::v::HOUSEHOLD_IDX); + agent.newVariable(person::v::HOUSEHOLD_SIZE); - // @todo - temp or vis only? - agent.newVariable(exateppabm::person::v::x); - agent.newVariable(exateppabm::person::v::y); - agent.newVariable(exateppabm::person::v::z); + // Workplace environment directed graph + // This single graph contains workplace information for all individuals, and is essentially 5 unconnected sub graphs. + flamegpu::EnvironmentDirectedGraphDescription workplaceDigraphDesc = env.newDirectedGraph("WORKPLACE_DIGRAPH"); + // Graph vertices + // workplaceDigraphDesc.newVertexProperty("bar"); + // workplaceDigraphDesc.newEdgeProperty("foo"); + + // Workplace network variables. @todo - refactor to a separate network location? + agent.newVariable(person::v::WORKPLACE_IDX); + agent.newVariable(person::v::WORKPLACE_OUT_DEGREE); + + // Random interaction network variables. @todo -refactor to separate location + agent.newVariable(person::v::RANDOM_INTERACTION_PARTNERS, {flamegpu::ID_NOT_SET}); + agent.newVariable(person::v::RANDOM_INTERACTION_COUNT, 0u); + agent.newVariable(person::v::RANDOM_INTERACTION_COUNT_TARGET, 0u); + + // Add an environmental variable containing the sum of each agents target number of random interactions. + // This is not the number of pair-wise interactions, but the number each individual would like. + // there are cases where this value would be impossible to achieve for a given population. + // I.e. if there are 2 agents, with targets of [2, 1]. This would be a target of 3, but only 1 interaction is possible + env.newProperty("RANDOM_INTERACTION_COUNT_SUM", 0u); + +#if defined(FLAMEGPU_VISUALISATION) + // @vis only + agent.newVariable(person::v::x); + agent.newVariable(person::v::y); + agent.newVariable(person::v::z); +#endif // defined(FLAMEGPU_VISUALISATION) // Define relevant messages - // Message list containing a persons current status (id, location, infection status) - flamegpu::MessageSpatial2D::Description statusMessage = model.newMessage(exateppabm::person::message::STATUS); - // Set the range and bounds for the spatial component of the status message. - // This effectively limits the range of the simulation. For now we pass these values in. - statusMessage.setRadius(env.getProperty("INFECTION_INTERACTION_RADIUS")); - statusMessage.setMin(0, 0); - statusMessage.setMax(width, width); - - // Add the id to the message. x, y, z are implicit - statusMessage.newVariable(person::message::v::STATUS_ID); + // Message list containing a persons current status for households (id, location, infection status) + flamegpu::MessageBucket::Description householdStatusMessage = model.newMessage(person::message::household_status::_NAME); + + // Set the maximum bucket index to the population size. Ideally this would be household count, but that is not known at model definition time. + // In the future this would be possible once https://github.com/FLAMEGPU/FLAMEGPU2/issues/710 is implemented + householdStatusMessage.setUpperBound(params.n_total); + + // Add the agent id to the message. + householdStatusMessage.newVariable(person::message::household_status::ID); + // Add the household index + householdStatusMessage.newVariable(person::v::HOUSEHOLD_IDX); + // Add a variable for the agent's infections status + householdStatusMessage.newVariable(person::v::INFECTION_STATE); + // Demographic? + householdStatusMessage.newVariable(person::v::AGE_DEMOGRAPHIC); + + // Message list containing a persons current status for workplaces (id, location, infection status) + flamegpu::MessageArray::Description workplaceStatusMessage = model.newMessage(person::message::workplace_status::_NAME); + + // Set the maximum message array index to the maximum expected ID (n_total + 1) + workplaceStatusMessage.setLength(params.n_total + 1); + + // Add the agent id to the message. + workplaceStatusMessage.newVariable(person::message::workplace_status::ID); + // Add the household index + workplaceStatusMessage.newVariable(person::v::WORKPLACE_IDX); // Add a variable for the agent's infections status - statusMessage.newVariable(exateppabm::person::v::INFECTION_STATE); + workplaceStatusMessage.newVariable(person::v::INFECTION_STATE); // Demographic? - statusMessage.newVariable(exateppabm::person::v::AGE_DEMOGRAPHIC); + workplaceStatusMessage.newVariable(person::v::AGE_DEMOGRAPHIC); + + // Message list containing a persons current status for random,daily interaction (id, location, infection status) + // Uses an array message with 1 per agent + flamegpu::MessageArray::Description randomNetworkStatusMessage = model.newMessage(person::message::random_network_status::_NAME); + + // ID's are 1 indexed (0 is unset) so use + 1. + // For ensembles this will need to be the largest n_total until https://github.com/FLAMEGPU/FLAMEGPU2/issues/710 is implemented + randomNetworkStatusMessage.setLength(params.n_total + 1); + + // No need to add the agent ID to the message, it's implied by the bin count + // randomNetworkStatusMessage.newVariable(person::message::random_network_status::ID); + + // Add a variable for the agent's infections status + randomNetworkStatusMessage.newVariable(person::v::INFECTION_STATE); + // Agent's demographic + randomNetworkStatusMessage.newVariable(person::v::AGE_DEMOGRAPHIC); // Define agent functions - // emit current status - flamegpu::AgentFunctionDescription emitStatusDesc = agent.newFunction("emitStatus", emitStatus); - emitStatusDesc.setMessageOutput(exateppabm::person::message::STATUS); - emitStatusDesc.setInitialState(exateppabm::person::states::DEFAULT); - emitStatusDesc.setEndState(exateppabm::person::states::DEFAULT); - - // Interact with other agents via their messages - flamegpu::AgentFunctionDescription interactDesc = agent.newFunction("interact", interact); - interactDesc.setMessageInput(exateppabm::person::message::STATUS); - emitStatusDesc.setInitialState(exateppabm::person::states::DEFAULT); - emitStatusDesc.setEndState(exateppabm::person::states::DEFAULT); + // emit current status to the household + flamegpu::AgentFunctionDescription emitHouseholdStatusDesc = agent.newFunction("emitHouseholdStatus", emitHouseholdStatus); + emitHouseholdStatusDesc.setMessageOutput(person::message::household_status::_NAME); + emitHouseholdStatusDesc.setMessageOutputOptional(true); + emitHouseholdStatusDesc.setInitialState(person::states::DEFAULT); + emitHouseholdStatusDesc.setEndState(person::states::DEFAULT); + + // Interact with other agents in the household via their messages + flamegpu::AgentFunctionDescription interactHouseholdDesc = agent.newFunction("interactHousehold", interactHousehold); + interactHouseholdDesc.setMessageInput(person::message::household_status::_NAME); + interactHouseholdDesc.setInitialState(person::states::DEFAULT); + interactHouseholdDesc.setEndState(person::states::DEFAULT); + + // emit current status to the workplace + flamegpu::AgentFunctionDescription emitWorkplaceStatusDesc = agent.newFunction("emitWorkplaceStatus", emitWorkplaceStatus); + emitWorkplaceStatusDesc.setMessageOutput(person::message::workplace_status::_NAME); + emitWorkplaceStatusDesc.setInitialState(person::states::DEFAULT); + emitWorkplaceStatusDesc.setEndState(person::states::DEFAULT); + + // Interact with other agents in the workplace via their messages + flamegpu::AgentFunctionDescription interactWorkplaceDesc = agent.newFunction("interactWorkplace", interactWorkplace); + interactWorkplaceDesc.setMessageInput(person::message::workplace_status::_NAME); + interactWorkplaceDesc.setInitialState(person::states::DEFAULT); + interactWorkplaceDesc.setEndState(person::states::DEFAULT); + + // Update the daily random interactions + // @todo - a GPU implementation will be needed for model scaling. + flamegpu::HostFunctionDescription("updateRandomDailyNetworkIndices", updateRandomDailyNetworkIndices); + + // emit current status for random interactions + flamegpu::AgentFunctionDescription emitRandomDailyNetworkStatusDesc = agent.newFunction("emitRandomDailyNetworkStatus", emitRandomDailyNetworkStatus); + emitRandomDailyNetworkStatusDesc.setMessageOutput(person::message::random_network_status::_NAME); + emitRandomDailyNetworkStatusDesc.setInitialState(person::states::DEFAULT); + emitRandomDailyNetworkStatusDesc.setEndState(person::states::DEFAULT); + + // Interact with other agents in the random interactions + flamegpu::AgentFunctionDescription interactRandomDailyNetworkDesc = agent.newFunction("interactRandomDailyNetwork", interactRandomDailyNetwork); + interactRandomDailyNetworkDesc.setMessageInput(person::message::random_network_status::_NAME); + interactRandomDailyNetworkDesc.setInitialState(person::states::DEFAULT); + interactRandomDailyNetworkDesc.setEndState(person::states::DEFAULT); } void appendLayers(flamegpu::ModelDescription& model) { + // Household interactions + { + auto layer = model.newLayer(); + layer.addAgentFunction(person::NAME, "emitHouseholdStatus"); + } + { + auto layer = model.newLayer(); + layer.addAgentFunction(person::NAME, "interactHousehold"); + } + // Workplace interactions + { + auto layer = model.newLayer(); + layer.addAgentFunction(person::NAME, "emitWorkplaceStatus"); + } + { + auto layer = model.newLayer(); + layer.addAgentFunction(person::NAME, "interactWorkplace"); + } + // Random interactions + { + auto layer = model.newLayer(); + layer.addHostFunction(updateRandomDailyNetworkIndices); + } { auto layer = model.newLayer(); - layer.addAgentFunction(exateppabm::person::NAME, "emitStatus"); + layer.addAgentFunction(person::NAME, "emitRandomDailyNetworkStatus"); } { auto layer = model.newLayer(); - layer.addAgentFunction(exateppabm::person::NAME, "interact"); + layer.addAgentFunction(person::NAME, "interactRandomDailyNetwork"); } } diff --git a/src/exateppabm/person.h b/src/exateppabm/person.h index 38a0fab..1422b2e 100644 --- a/src/exateppabm/person.h +++ b/src/exateppabm/person.h @@ -23,6 +23,16 @@ namespace person { */ constexpr char NAME[] = "person"; +/** + * Maximum number of random daily interactions, which must be known at compile time due to FLAME GPU array variable limitations (they are arrays not vectors) + * @todo - make this value overrideable via CMake. + */ +#ifdef EXATEPP_ABM_MAX_RANDOM_DAILY_INTERACTIONS +constexpr std::uint32_t MAX_RANDOM_DAILY_INTERACTIONS = EXATEPP_ABM_MAX_RANDOM_DAILY_INTERACTIONS; +#else +#error "EXATEPP_ABM_MAX_RANDOM_DAILY_INTERACTIONS is not defined" +#endif // EXATEPP_ABM_MAX_RANDOM_DAILY_INTERACTIONS + /** * Namespace containing host device constants for state names related to the Person agent type */ @@ -38,23 +48,49 @@ namespace v { DEVICE_CONSTEXPR_STRING constexpr char x[] = "x"; DEVICE_CONSTEXPR_STRING constexpr char y[] = "y"; DEVICE_CONSTEXPR_STRING constexpr char z[] = "z"; +DEVICE_CONSTEXPR_STRING constexpr char ID[] = "ID"; // Custom id, as FLAMEGPU getID isn't guaranteed to start at 1. DEVICE_CONSTEXPR_STRING constexpr char INFECTION_STATE[] = "infection_state"; DEVICE_CONSTEXPR_STRING constexpr char INFECTION_STATE_CHANGE_DAY[] = "infection_state_change_day"; DEVICE_CONSTEXPR_STRING constexpr char INFECTION_STATE_DURATION[] = "infection_state_duration"; +DEVICE_CONSTEXPR_STRING constexpr char INFECTION_COUNT[] = "infection_count"; DEVICE_CONSTEXPR_STRING constexpr char AGE_DEMOGRAPHIC[] = "age_demographic"; +DEVICE_CONSTEXPR_STRING constexpr char HOUSEHOLD_IDX[] = "household_idx"; +DEVICE_CONSTEXPR_STRING constexpr char HOUSEHOLD_SIZE[] = "household_size"; +DEVICE_CONSTEXPR_STRING constexpr char WORKPLACE_IDX[] = "workplace_idx"; +DEVICE_CONSTEXPR_STRING constexpr char WORKPLACE_OUT_DEGREE[] = "workplace_out_degree"; +DEVICE_CONSTEXPR_STRING constexpr char RANDOM_INTERACTION_PARTNERS[] = "random_interaction_partners"; +DEVICE_CONSTEXPR_STRING constexpr char RANDOM_INTERACTION_COUNT[] = "random_interaction_count"; +DEVICE_CONSTEXPR_STRING constexpr char RANDOM_INTERACTION_COUNT_TARGET[] = "random_interaction_count_target"; + + } // namespace v /** * Namespace containing person-message related constants */ namespace message { -DEVICE_CONSTEXPR_STRING constexpr char STATUS[] = "status"; /** - * Namespace containing variable name constants for variables in person related messages + * Namespace containing variable name constants for variables in household related messages */ -namespace v { -DEVICE_CONSTEXPR_STRING constexpr char STATUS_ID[] = "id"; -} // namespace v +namespace household_status { +DEVICE_CONSTEXPR_STRING constexpr char _NAME[] = "household_status"; +DEVICE_CONSTEXPR_STRING constexpr char ID[] = "id"; +} // namespace household_status +/** + * Namespace containing variable name constants for variables in workplace related messages + */ +namespace workplace_status { +DEVICE_CONSTEXPR_STRING constexpr char _NAME[] = "workplace_status"; +DEVICE_CONSTEXPR_STRING constexpr char ID[] = "id"; +} // namespace workplace_status + +/** + * Namespace containing variable name constants for variables in workplace related messages + */ +namespace random_network_status { +DEVICE_CONSTEXPR_STRING constexpr char _NAME[] = "random_network_status"; +DEVICE_CONSTEXPR_STRING constexpr char ID[] = "id"; +} // namespace random_network_status } // namespace message @@ -63,16 +99,14 @@ DEVICE_CONSTEXPR_STRING constexpr char STATUS_ID[] = "id"; * Define the agent type representing a person in the simulation, mutating the model description object. * @param model flamegpu2 model description object to mutate * @param params model parameters from parameters file - * @param width the width of the 2D space currently used for spatials comms. to be removed once networks added. - * @param interactionRadius spatial interaction radius for temporary infection spread behaviour. to be removed. */ -void define(flamegpu::ModelDescription& model, const exateppabm::input::config& params, const float width, const float interactionRadius); +void define(flamegpu::ModelDescription& model, const exateppabm::input::config& params); /** * Add person related functions to the FLAMEGPU 2 layer based control flow. - * + * * Does not use the DAG abstraction due to previously encountered bugs with split compilation units which have not yet been pinned down / resolved. - * + * * @param model flamegpu2 model description object to mutate */ void appendLayers(flamegpu::ModelDescription& model); diff --git a/src/exateppabm/population.cu b/src/exateppabm/population.cu index c519d93..686b95f 100644 --- a/src/exateppabm/population.cu +++ b/src/exateppabm/population.cu @@ -3,129 +3,361 @@ #include #include +#include #include #include #include +#include #include #include "exateppabm/demographics.h" #include "exateppabm/disease.h" -#include "exateppabm/person.h" #include "exateppabm/input.h" +#include "exateppabm/network.h" +#include "exateppabm/person.h" #include "exateppabm/util.h" +#include "exateppabm/visualisation.h" namespace exateppabm { namespace population { namespace { -// File-scoped array contianing the number of infected agents per demographic from population initialisation. This needs to be made accessible to a FLAME GPU Init func due to macro environment property limitations. +// File-scoped array containing the number of infected agents per demographic from population initialisation. This needs to be made accessible to a FLAME GPU Init func due to macro environment property limitations. +std::array _infectedPerDemographic = {}; + +// File-scoped copy of the model parameters struct +exateppabm::input::config _params = {}; + +// File-scoped boolean for if generation should be verbose or not. +bool _verbose = false; -std::array infectedPerDemographic = {}; +// Store the ID of each individual for each workplace, to form the node labels of the small world network +std::array, WORKPLACE_COUNT> workplaceMembers = {{}}; } // namespace -std::unique_ptr generate(flamegpu::ModelDescription& model, const exateppabm::input::config config, const bool verbose, const float env_width, const float interactionRadius) { - fmt::print("@todo - validate config inputs when generated agents (pop size, initial infected count etc)\n"); - // @todo - assert that the requested initial population is non zero. - auto pop = std::make_unique(model.Agent(exateppabm::person::NAME), config.n_total); +// Workaround struct representing a person used during host generation +struct HostPerson { + flamegpu::id_t id; + disease::SEIR::InfectionStateUnderlyingType infectionStatus; + float infectionStateDuration; + std::uint32_t infectionCount = 0; + demographics::AgeUnderlyingType ageDemographic; + std::uint32_t householdIdx; + std::uint8_t householdSize; + std::uint32_t workplaceIdx; + std::uint32_t randomInteractionTarget; + std::uint32_t workplaceOutDegree; +}; - std::uint64_t sq_width = static_cast(env_width); - // float expectedNeighbours = interactionRadius * interactionRadius; - // fmt::print("sq_width {} interactionRadius {} expectedNeighbours {}\n", sq_width, interactionRadius, expectedNeighbours); +/** + * FLAME GPU init function which generates a population of person agents for the current simulation + * + * This function has been split into multiple init functions, due to errors encountered when attempting to do multiple passes over a newly created set of agents. There should be a working way to do this... + */ +FLAMEGPU_INIT_FUNCTION(generatePopulation) { + fmt::print("@todo - validate params inputs when generated agents (pop size, initial infected count etc)\n"); + // n_total must be less than uint max + // n_total must be more than 0. + // initial infected count must be more than 0, less than full pop. + + // Get a handle on the environment. + auto env = FLAMEGPU->environment; + + // Get a handle to the host agent api for Person agents + // auto personAgent = FLAMEGPU->agent(exateppabm::person::NAME, exateppabm::person::states::DEFAULT); + + // Do nothing if no one to create? + if (_params.n_total == 0) { + return; + } + + // Workaround: Generating agents in FLAME GPU 2 in init functions does not allow iterating the data structure multiple times, this appears to be a bug (which may have a fix in place, currently untested pr) + + // Instead, we must generate data on the host for the "first pass", storing it in non-flame gpu storage + // Then set it on the latter pass. + std::vector hostPersonData; + hostPersonData.resize(_params.n_total); // seed host rng for population generation. // @todo - does this want to be a separate seed from the config file? - std::mt19937_64 rng(config.rng_seed); - + std::mt19937_64 rng(FLAMEGPU->getSimulationConfig().random_seed); // Need to initialise a fixed number of individuals as infected. // This not very scalable way of doing it, is to create a vector with one element per individual in the simulation, initialised to false // set the first n_seed_infection elements to true/1 // Shuffle the vector, and query at agent creation time // RNG sampling in-loop would be more memory efficient, but harder to guarantee that exactly enough are created. This will likely be replaced anyway, so quick and dirty is fine. - std::vector infected_vector(config.n_total); - std::fill(infected_vector.begin(), infected_vector.begin() + std::min(config.n_total, config.n_seed_infection), true); + std::vector infected_vector(_params.n_total); + std::fill(infected_vector.begin(), infected_vector.begin() + std::min(_params.n_total, _params.n_seed_infection), true); std::shuffle(infected_vector.begin(), infected_vector.end(), rng); - // Prepare a probability matrix for selecting an age demographic for the agent based on the ratio from the configuration. - // @todo abstract this into class/methods. - // @todo - this hardcoded 9 is a bit grim. Maybe enums can help? - std::uint64_t configDemographicSum = config.population_0_9 + config.population_10_19 + config.population_20_29 + config.population_30_39 + config.population_40_49 + config.population_50_59 + config.population_60_69 + config.population_70_79 + config.population_80; - // @todo - map might be more readable than an array (incase the underlying class enum values are ever changed to be a different order?) - std::array demographicProbabilties = {{ - config.population_0_9 / static_cast(configDemographicSum), - config.population_10_19 / static_cast(configDemographicSum), - config.population_20_29 / static_cast(configDemographicSum), - config.population_30_39 / static_cast(configDemographicSum), - config.population_40_49 / static_cast(configDemographicSum), - config.population_50_59 / static_cast(configDemographicSum), - config.population_60_69 / static_cast(configDemographicSum), - config.population_70_79 / static_cast(configDemographicSum), - config.population_80 / static_cast(configDemographicSum) - }}; - // Perform an inclusive scan to convert to cumulative probability - // Using a local method which supports inclusive scans in old libstc++ - exateppabm::util::inplace_inclusive_scan(demographicProbabilties); - // std::inclusive_scan(demographicProbabilties.begin(), demographicProbabilties.end(), demographicProbabilties.begin()); - std::array allDemographics = {{ - demographics::Age::AGE_0_9, - demographics::Age::AGE_10_19, - demographics::Age::AGE_20_29, - demographics::Age::AGE_30_39, - demographics::Age::AGE_40_49, - demographics::Age::AGE_50_59, - demographics::Age::AGE_60_69, - demographics::Age::AGE_70_79, - demographics::Age::AGE_80 - }}; + // Get the number of individuals per house and their age demographics + auto households = generateHouseholdStructures(_params, rng, _verbose); // per demo total is not an output in time series. // Alternately, we need to initialise the exact number of each age band, not RNG, and just scale it down accordingly. Will look at in "realistic" population generation - std::array createdPerDemographic = {{0, 0, 0, 0, 0, 0, 0, 0, 0}}; + std::array createdPerDemographic = {}; + for (const auto & household : households) { + for (demographics::AgeUnderlyingType i = 0; i < demographics::AGE_COUNT; i++) { + createdPerDemographic[i] += household.sizePerAge[i]; + } + } + // reset per demographic count of the number initialised agents in each infection state. - infectedPerDemographic = {{0, 0, 0, 0, 0, 0, 0, 0, 0}}; + // This is used for the initial value in time series data, without having to iterate all agent data again, but may not be thread safe (i.e. probably need to change for ensembles) + _infectedPerDemographic = {{0, 0, 0, 0, 0, 0, 0, 0, 0}}; - std::uniform_real_distribution demo_dist(0.0f, 1.0f); + // Given the household and ages are known, we can compute the workplace assignment probabilities for adults + std::array p_adultPerWorkNetwork = getAdultWorkplaceCumulativeProbabilityArray(_params.child_network_adults, _params.elderly_network_adults, createdPerDemographic); + // // Store the ID of each individual for each workplace, to form the node labels of the small world network + // std::array, WORKPLACE_COUNT> workplaceMembers = {{}}; - unsigned idx = 0; - for (auto person : *pop) { - // Infections status - disease::SEIR::InfectionState infectionStatus = infected_vector.at(idx) ? disease::SEIR::InfectionState::Infected : disease::SEIR::InfectionState::Susceptible; - person.setVariable(exateppabm::person::v::INFECTION_STATE, infectionStatus); + // Counter for random interaction count. This will be 2x the number of interactions + std::uint64_t randomInteractionCountSum = 0u; + // Max/min trackers for random interaction targets, for verbose output. + std::uint32_t randomInteractionMax = 0u; + std::uint32_t randomInteractionMin = std::numeric_limits::max(); - // Demographic - // @todo - this is a bit grim, enum class aren't as nice as hoped. - float demo_random = demo_dist(rng); - // @todo - abstract this into a method. - demographics::Age demo = demographics::Age::AGE_0_9; - for (demographics::AgeUnderlyingType i = 0; i < demographics::AGE_COUNT; i++) { - if (demo_random < demographicProbabilties[i]) { - demo = allDemographics[i]; - createdPerDemographic[i]++; - if (infectionStatus == disease::SEIR::Infected) { - infectedPerDemographic[i]++; - } - break; + // Populate agent data, by iterating households + std::uint32_t personIdx = 0; + for (std::uint32_t householdIdx = 0; householdIdx < households.size(); householdIdx++) { + auto household = households.at(householdIdx); + // For each individual in the household + for (population::HouseholdSizeType householdMemberIdx = 0; householdMemberIdx < household.size; householdMemberIdx++) { + // assert that the household structure is complete + assert(household.size == household.agePerPerson.size()); + + // Get the flamegpu person object for the individual + // auto person = personAgent.newAgent(); @temp disabled due to bug/workaround in place for multi-pass + + // Set the manual ID, as getID() isn't guaranteed to actually start at 1. + flamegpu::id_t personID = personIdx + 1; + // person.setVariable(person::v::ID, personID); + hostPersonData[personIdx].id = personID; + + // Set the individuals infection status. @todo - refactor into seir.cu? + disease::SEIR::InfectionState infectionStatus = infected_vector.at(personIdx) ? disease::SEIR::InfectionState::Infected : disease::SEIR::InfectionState::Susceptible; + // person.setVariable(exateppabm::person::v::INFECTION_STATE, infectionStatus); + hostPersonData[personIdx].infectionStatus = infectionStatus; + // Also set the initial infection duration. @todo - stochastic. + float infectionStateDuration = infectionStatus == disease::SEIR::InfectionState::Infected ? _params.mean_time_to_recovered: 0; + // person.setVariable(exateppabm::person::v::INFECTION_STATE_DURATION, infectionStateDuration); + hostPersonData[personIdx].infectionStateDuration = infectionStateDuration; + + // Set the individuals age and household properties + demographics::Age age = household.agePerPerson[householdMemberIdx]; + // person.setVariable(person::v::AGE_DEMOGRAPHIC, static_cast(age)); + hostPersonData[personIdx].ageDemographic = static_cast(age); + // person.setVariable(person::v::HOUSEHOLD_IDX, householdIdx); + hostPersonData[personIdx].householdIdx = householdIdx; + // person.setVariable(person::v::HOUSEHOLD_SIZE, household.size); + hostPersonData[personIdx].householdSize = household.size; + + // initialise the agents infection count + if (infectionStatus == disease::SEIR::Infected) { + // person.setVariable(exateppabm::person::v::INFECTION_COUNT, 1u); + hostPersonData[personIdx].infectionCount = 1u; + // Increment the per-age demographic initial agent count. @todo refactor elsewhere? + _infectedPerDemographic[age]++; } + + // Assign the workplace based on age band + WorkplaceUnderlyingType workplaceIdx = generateWorkplaceForIndividual(age, p_adultPerWorkNetwork, rng); + + // Store the agent's manually set ID, cannot rely on autogenerated to be useful indices + workplaceMembers[workplaceIdx].push_back(personID); + // Store the assigned network in the agent data structure + // person.setVariable(person::v::WORKPLACE_IDX, workplaceIdx); + hostPersonData[personIdx].workplaceIdx = workplaceIdx; + + // Generate the (target) number of random interactions this individual will be involved in per day. Not all combinations will be possible on all days, hence target. + + // @todo - optionally allow binomial distributions + // @todo - decide if non non-binomial should be a mean or not, maybe allow non fixed normal dists? + // @todo - refactor and test. + double meanRandomInteractions = _params.mean_random_interactions_20_69; + double sdRandomInteractions = _params.sd_random_interactions_20_69; + if (age == demographics::Age::AGE_0_9 || age == demographics::Age::AGE_10_19) { + meanRandomInteractions = _params.mean_random_interactions_0_19; + sdRandomInteractions = _params.sd_random_interactions_0_19; + } else if (age == demographics::Age::AGE_70_79 || age == demographics::Age::AGE_80) { + meanRandomInteractions = _params.mean_random_interactions_70plus; + sdRandomInteractions = _params.sd_random_interactions_70plus; + } + + // Sample a normal distribution (of integers, so we can clamp to >= 0) + std::normal_distribution randomInteractionDist{meanRandomInteractions, sdRandomInteractions}; + // Sample from the distribution + double randomInteractionsRaw = randomInteractionDist(rng); + // Clamp to be between 0 and the population size, and cast to uint + std::uint32_t randomInteractionTarget = static_cast(std::clamp(randomInteractionsRaw, 0.0, static_cast(_params.n_total))); + + // If the max was over the compile time upper limit due to flamegpu limitations, emit a warning and exit. + if (randomInteractionTarget > person::MAX_RANDOM_DAILY_INTERACTIONS) { + fmt::print(stderr, "Fatal Error: Random Interaction Target {} exceeds fixed limit MAX_RANDOM_DAILY_INTERACTIONS {}. Please rebuild with a higher value @todo\n", randomInteractionTarget, person::MAX_RANDOM_DAILY_INTERACTIONS); + exit(EXIT_FAILURE); + } + + // Update the min and max tracking for output to stdout + randomInteractionMin = std::min(randomInteractionMin, randomInteractionTarget); + randomInteractionMax = std::max(randomInteractionMax, randomInteractionTarget); + // Set for the agent + // person.setVariable(person::v::RANDOM_INTERACTION_COUNT_TARGET, randomInteractionTarget); + hostPersonData[personIdx].randomInteractionTarget = randomInteractionTarget; + // Track the sum of target interaction counts + randomInteractionCountSum += randomInteractionTarget; + + // Increment the person index + ++personIdx; } - person.setVariable(exateppabm::person::v::AGE_DEMOGRAPHIC, demo); + } - // Location in 3D space (temp/vis) - unsigned row = idx / sq_width; - unsigned col = idx % sq_width; - person.setVariable(exateppabm::person::v::x, col); // @todo temp - person.setVariable(exateppabm::person::v::y, row); // @todo -temp - person.setVariable(exateppabm::person::v::z, 0); // @todo -temp + if (_verbose) { + fmt::print("Random Interactions: min={}, max={}, sum={}\n", randomInteractionMin, randomInteractionMax, randomInteractionCountSum); + } - // Inc counter - ++idx; + // Set the sum of per agent target interaction counts. This is not the total number of interactions that will occur, as not all configurations are possible. + // This only works if this method is called prior to simulation construction! + env.setProperty("RANDOM_INTERACTION_COUNT_SUM", randomInteractionCountSum); + // Warn if the number target number of interactions is odd. + if (randomInteractionCountSum % 2 == 1) { + fmt::print(stderr, "Warning: Total the sum of per-agent random interactions is odd ({})\n", randomInteractionCountSum); } - if (verbose) { + // Now individuals have been assigned to each workplace network, we can generate the small world networks for workplace interactions + std::array meanInteractionsPerNetwork = {{ + _params.mean_work_interactions_child / _params.daily_fraction_work, + _params.mean_work_interactions_child / _params.daily_fraction_work, + _params.mean_work_interactions_adult / _params.daily_fraction_work, + _params.mean_work_interactions_elderly / _params.daily_fraction_work, + _params.mean_work_interactions_elderly / _params.daily_fraction_work}}; + + // agent count + 1 elements, as 1 indexed + std::uint32_t totalVertices = _params.n_total; + // Count the total number of edges + std::uint32_t totalUndirectedEdges = 0u; + std::array workplaceNetworks = {}; + for (population::WorkplaceUnderlyingType widx = 0; widx < population::WORKPLACE_COUNT; ++widx) { + const std::uint64_t N = workplaceMembers[widx].size(); + // Round to nearest even number. This might want to be always round up instead @todo + std::uint32_t meanInteractions = static_cast(std::nearbyint(meanInteractionsPerNetwork[widx] / 2.0) * 2.0); + // @todo - test / cleaner handling of different N and K values. + if (meanInteractions >= N) { + // @todo - this feels brittle. + meanInteractions = N - 1; + } + + // Generate the small world network, unless N is too small, or meanInteractions is too small + if (N > 2 && meanInteractions > 1) { + // Shuffle the network, to avoid people working with others in their same household a disproportionate amount for low values of work_network_rewire + std::shuffle(workplaceMembers[widx].begin(), workplaceMembers[widx].end(), rng); + // Generate the small world network + workplaceNetworks[widx] = network::generateSmallWorldUndirectedGraph(workplaceMembers[widx], meanInteractions, _params.work_network_rewire, rng); + } else { + workplaceNetworks[widx] = network::generateFullyConnectedUndirectedGraph(workplaceMembers[widx]); + } + + totalUndirectedEdges += workplaceNetworks[widx].getNumEdges(); + } + + // Get a handle to the FLAME GPU workplace directed graph + flamegpu::HostEnvironmentDirectedGraph workplaceDigraph = FLAMEGPU->environment.getDirectedGraph("WORKPLACE_DIGRAPH"); + // This will contain all workplace information, i.e. a single graph data structure containing 5 unconnected graphs + // The custom 1-indexed agent's ID will be used as the vertex keys + // This will then allow agents to simply lookup their neighbours and + // FLAME GPU's graphs are accessed by their string name, but agents can't have string properties, hence using a single graph. + // Representing a undirected graph using a directed graph will consume twice as much device memory as required, but only digraphs are currently implemented in FLAME GPU 2. + std::uint32_t totalDirectedEdges = 2 * totalUndirectedEdges; + workplaceDigraph.setVertexCount(totalVertices); + workplaceDigraph.setEdgeCount(totalDirectedEdges); + + // Set vertex properties + flamegpu::HostEnvironmentDirectedGraph::VertexMap vertices = workplaceDigraph.vertices(); + for (std::uint32_t v = 1; v <= totalVertices; ++v) { + flamegpu::HostEnvironmentDirectedGraph::VertexMap::Vertex vertex = vertices[v]; + // vertex.setProperty(person::v::ID, v); + } + + // Iterate each workplace network, adding it's both directions of each undirected edge to the graph + flamegpu::HostEnvironmentDirectedGraph::EdgeMap edges = workplaceDigraph.edges(); + for (auto& undirectedNetwork : workplaceNetworks) { + // For each undirected edge which contains indexes within the network, create the 2 directed edges between agent IDs + for (const auto & undirectedEdge : undirectedNetwork.getEdges()) { + flamegpu::id_t a = undirectedNetwork.getVertexLabel(undirectedEdge.source); + flamegpu::id_t b = undirectedNetwork.getVertexLabel(undirectedEdge.dest); + auto abEdge = edges[{a, b}]; + abEdge.setSourceDestinationVertexID(a, b); + auto baEdge = edges[{b, a}]; + baEdge.setSourceDestinationVertexID(b, a); + } + } + + // flamegpu::DeviceAgentVector population = personAgent.getPopulationData(); + // Update the population with additional data. + // This is (potentially) suboptimal in terms of host-device memcpy, but need to pass agents multiple times unfortunately. + + // flamegpu::DeviceAgentVector pop = FLAMEGPU->agent("person", "default").getPopulationData(); + // std::uint32_t personIdx = 0; + // for (auto person : pop) { + for (std::uint32_t personIdx = 0; personIdx < hostPersonData.size(); ++personIdx) { + // Get the agents id + // std::uint32_t agentId = person.getVariable(person::v::ID); + std::uint32_t agentId = hostPersonData[personIdx].id; + // get the assigned workplace + // std::uint32_t workplaceIdx = person.getVariable(person::v::WORKPLACE_IDX); + std::uint32_t workplaceIdx = hostPersonData[personIdx].workplaceIdx; + // For this individual in their small world network, get the (out)degree, i.e. the max number of interactions per day. + network::UndirectedGraph& workplaceGraph = workplaceNetworks[workplaceIdx]; + auto vertexLabels = workplaceGraph.getVertexLabels(); + auto it = std::find(vertexLabels.begin(), vertexLabels.end(), agentId); + if (it != vertexLabels.end()) { + std::uint32_t vertexIdx = std::distance(vertexLabels.begin(), it); + // Get the outdegree for this vertex index + std::uint32_t degree = workplaceGraph.degree(vertexIdx); + // person.setVariable(person::v::WORKPLACE_OUT_DEGREE, degree); + hostPersonData[personIdx].workplaceOutDegree = degree; + } else { + throw std::runtime_error("@todo - could not find agent in workplace"); + } + // ++personIdx; + } + + // Finally, in another pass create the actual agent instances, finishing the multi-pass init function (and split init function workaround) + auto personAgent = FLAMEGPU->agent(exateppabm::person::NAME, exateppabm::person::states::DEFAULT); + for (std::uint32_t personIdx = 0; personIdx < hostPersonData.size(); ++personIdx) { + const auto& hostPerson = hostPersonData[personIdx]; + // Generate the new agent + auto person = personAgent.newAgent(); + // Set each property on the agent. + person.setVariable(person::v::ID, hostPerson.id); + person.setVariable(exateppabm::person::v::INFECTION_STATE, hostPerson.infectionStatus); + person.setVariable(exateppabm::person::v::INFECTION_STATE_DURATION, hostPerson.infectionStateDuration); + person.setVariable(exateppabm::person::v::INFECTION_COUNT, hostPerson.infectionCount); + person.setVariable(person::v::AGE_DEMOGRAPHIC, hostPerson.ageDemographic); + person.setVariable(person::v::HOUSEHOLD_IDX, hostPerson.householdIdx); + person.setVariable(person::v::HOUSEHOLD_SIZE, hostPerson.householdSize); + person.setVariable(person::v::WORKPLACE_IDX, hostPerson.workplaceIdx); + person.setVariable(person::v::WORKPLACE_OUT_DEGREE, hostPerson.workplaceOutDegree); + person.setVariable(person::v::RANDOM_INTERACTION_COUNT_TARGET, hostPerson.randomInteractionTarget); + + + // If this is a visualisation enabled build, set their x/y/z +#if defined(FLAMEGPU_VISUALISATION) + auto [visX, visY, visZ] = exateppabm::visualisation::getAgentXYZ(static_cast(households.size()), hostPerson.householdIdx, 0); + person.setVariable(exateppabm::person::v::x, visX); + person.setVariable(exateppabm::person::v::y, visY); + person.setVariable(exateppabm::person::v::z, visZ); +#endif // defined(FLAMEGPU_VISUALISATION) + } + + // Do the verbose output + if (_verbose) { // Print a summary of population creation for now. - fmt::print("Created {} people with {} infected.\n", config.n_total, config.n_seed_infection); + fmt::print("Created {} people with {} infected.\n", _params.n_total, _params.n_seed_infection); + fmt::print("Households: {}\n", households.size()); fmt::print("Demographics {{\n"); fmt::print(" 0- 9 = {}\n", createdPerDemographic[0]); fmt::print(" 10-19 = {}\n", createdPerDemographic[1]); @@ -137,13 +369,196 @@ std::unique_ptr generate(flamegpu::ModelDescription& mode fmt::print(" 70-79 = {}\n", createdPerDemographic[7]); fmt::print(" 80+ = {}\n", createdPerDemographic[8]); fmt::print("}}\n"); + fmt::print("Workplaces {{\n"); + fmt::print(" 00-09: {}, 20_69 {}\n", createdPerDemographic[0], workplaceMembers[Workplace::WORKPLACE_SCHOOL_0_9].size() - createdPerDemographic[demographics::Age::AGE_0_9]); + fmt::print(" 10-19: {}, 20_69 {}\n", createdPerDemographic[demographics::Age::AGE_10_19], workplaceMembers[Workplace::WORKPLACE_SCHOOL_0_9].size() - createdPerDemographic[demographics::Age::AGE_10_19]); + fmt::print(" 20_69 {}\n", workplaceMembers[Workplace::WORKPLACE_ADULT].size()); + fmt::print(" 20_69 {}, 70_79 {}\n", workplaceMembers[Workplace::WORKPLACE_70_79].size() - createdPerDemographic[demographics::Age::AGE_70_79], createdPerDemographic[demographics::Age::AGE_70_79]); + fmt::print(" 20_69 {}, 80+ {}\n", workplaceMembers[Workplace::WORKPLACE_80_PLUS].size() - - createdPerDemographic[demographics::Age::AGE_80], createdPerDemographic[demographics::Age::AGE_80]); + fmt::print("}}\n"); } +} - return pop; +void define(flamegpu::ModelDescription& model, const exateppabm::input::config params, const bool verbose) { + // Store passed in parameters in file-scoped variables + _params = params; + _verbose = verbose; + // Define the init function which will generate the population for the parameters struct stored in the anon namespace + model.addInitFunction(generatePopulation); } std::array getPerDemographicInitialInfectionCount() { - return infectedPerDemographic; + return _infectedPerDemographic; +} + +double getReferenceMeanHouseholdSize(const exateppabm::input::config& params) { + std::vector countPerSize = {{params.household_size_1, params.household_size_2, params.household_size_3, params.household_size_4, params.household_size_5, params.household_size_6}}; + std::uint64_t refPeople = 0u; + std::uint64_t refHouses = 0u; + for (std::size_t idx = 0; idx < countPerSize.size(); idx++) { + refPeople += (idx + 1) * countPerSize[idx]; + refHouses += countPerSize[idx]; + } + double refMeanHouseholdSize = refPeople / static_cast(refHouses); + return refMeanHouseholdSize; +} + +std::vector getHouseholdSizeCumulativeProbabilityVector(const exateppabm::input::config& params) { + // Initialise vector with each config household size + std::vector countPerSize = {{params.household_size_1, params.household_size_2, params.household_size_3, params.household_size_4, params.household_size_5, params.household_size_6}}; + // get the sum, to find relative proportions + std::uint64_t sumConfigHouseholdSizes = exateppabm::util::reduce(countPerSize.begin(), countPerSize.end(), 0ull); + // Get the number of people in each household band for the reference size + // Find the number of people that the reference household sizes can account for + std::vector peoplePerHouseSize = countPerSize; + for (std::size_t idx = 0; idx < peoplePerHouseSize.size(); idx++) { + peoplePerHouseSize[idx] = (idx + 1) * peoplePerHouseSize[idx]; + } + std::uint64_t sumConfigPeoplePerHouseSize = exateppabm::util::reduce(peoplePerHouseSize.begin(), peoplePerHouseSize.end(), 0ull); + double configMeanPeoplePerHouseSize = sumConfigPeoplePerHouseSize / static_cast(sumConfigHouseholdSizes); + + // Build a list of household sizes, by random sampling from a uniform distribution using probabilities from the reference house size counts. + std::vector householdSizeProbability(countPerSize.size()); + for (std::size_t idx = 0; idx < householdSizeProbability.size(); idx++) { + householdSizeProbability[idx] = countPerSize[idx] / static_cast(sumConfigHouseholdSizes); + } + // Perform an inclusive scan to convert to cumulative probability + exateppabm::util::inclusive_scan(householdSizeProbability.begin(), householdSizeProbability.end(), householdSizeProbability.begin()); + + return householdSizeProbability; +} + +std::vector generateHouseholdStructures(const exateppabm::input::config params, std::mt19937_64 & rng, const bool verbose) { + /* + @todo This method will want refactoring for realistic household generation. + Current structure is: + 1. Get the cumulative probability distribution of house sizes based on reference data + 2. Generate household sizes randomly, using based on reference house size data + 3. For each house, generate the age per person within the household using probabilities based on global age demographic target data + */ + + // Get the vector of household size cumulative probability + auto householdSizeProbabilityVector = getHouseholdSizeCumulativeProbabilityVector(params); + + // Get the array of age demographic cumulative probability and reverse enum map + auto ageDemographicProbabilities = demographics::getAgeDemographicCumulativeProbabilityArray(params); + auto allAgeDemographics = demographics::getAllAgeDemographics(); + + + // Specify the rng distribution to sample from, [0, 1.0) + std::uniform_real_distribution dist(0.0, 1.0); + + // initialise a value with the number of people left to generate + std::int64_t remainingPeople = static_cast(params.n_total); + // estimate the number of houses, based on the cumulative probability distribution + double refMeanHouseSize = getReferenceMeanHouseholdSize(params); + std::uint64_t householdCountEstimate = static_cast(std::ceil(remainingPeople / refMeanHouseSize)); + + // Create a vector of house structures, and reserve enough room for the estimated number of houses + std::vector households = {}; + households.reserve(householdCountEstimate); + + // create enough households for the whole population using the uniform distribution and cumulative probability vector. Ensure the last household is not too large. + while (remainingPeople > 0) { + double r_houseSize = dist(rng); + for (std::size_t idx = 0; idx < static_cast(householdSizeProbabilityVector.size()); idx++) { + if (r_houseSize < householdSizeProbabilityVector[idx]) { + HouseholdStructure household = {}; + household.size = static_cast(idx + 1) <= remainingPeople ? static_cast(idx + 1) : remainingPeople; + household.agePerPerson.reserve(household.size); + // Generate ages for members of the household + for (HouseholdSizeType pidx = 0; pidx < household.size; ++pidx) { + float r_age = dist(rng); + // @todo - abstract this into a method. + demographics::Age age = demographics::Age::AGE_0_9; + for (demographics::AgeUnderlyingType i = 0; i < demographics::AGE_COUNT; i++) { + if (r_age < ageDemographicProbabilities[i]) { + age = allAgeDemographics[i]; + break; + } + } + household.agePerPerson.push_back(age); + household.sizePerAge[static_cast(age)]++; + } + households.push_back(household); + remainingPeople -= household.size; + break; + } + } + } + // potentially shrink the vector, in case the reservation was too large + households.shrink_to_fit(); + + if (verbose) { + // Get the count of created per house size and print it. + std::vector generatedHouseSizeDistribution(householdSizeProbabilityVector.size()); + for (const auto& household : households) { + generatedHouseSizeDistribution[household.size-1]++; + } + fmt::print("generated households per household size (total {}) {{\n", households.size()); + for (const auto & v : generatedHouseSizeDistribution) { + fmt::print(" {},\n", v); + } + fmt::print("}}\n"); + // Sum the number of people per household + std::uint64_t sumPeoplePerHouse = std::accumulate(households.begin(), households.end(), 0ull, [](std::uint64_t tot, HouseholdStructure& h) {return tot + h.size;}); + // std::uint64_t sumPeoplePerHouse = exateppabm::util::reduce(households.begin(), households.end(), 0ull); + // Check the mean still agrees. + double generatedMeanPeoplePerHouseSize = sumPeoplePerHouse / static_cast(households.size()); + fmt::print("generated mean household size {}\n", generatedMeanPeoplePerHouseSize); + } + return households; +} + +std::array getAdultWorkplaceCumulativeProbabilityArray(const double child_network_adults, const double elderly_network_adults, std::array n_per_age) { + // Adults are assigned to a work place network randomly, using a probability distribution which (for a large enough sample) will match the target ratios of adults to other workplace members, based on the child_network_adults and elderly_network_adults parameters + std::uint64_t n_0_9 = n_per_age[demographics::Age::AGE_0_9]; + std::uint64_t n_10_19 = n_per_age[demographics::Age::AGE_10_19]; + std::uint64_t n_70_79 = n_per_age[demographics::Age::AGE_70_79]; + std::uint64_t n_80_plus = n_per_age[demographics::Age::AGE_80]; + std::uint64_t n_adult = n_per_age[demographics::Age::AGE_20_29] + n_per_age[demographics::Age::AGE_30_39] + n_per_age[demographics::Age::AGE_40_49] + n_per_age[demographics::Age::AGE_50_59] + n_per_age[demographics::Age::AGE_60_69]; + + // Initialise with the target number of adults in each network + std::array p_adultPerWorkNetwork = {0}; + // p of each category is target number / number adults + p_adultPerWorkNetwork[Workplace::WORKPLACE_SCHOOL_0_9] = (child_network_adults * n_0_9) / n_adult; + p_adultPerWorkNetwork[Workplace::WORKPLACE_SCHOOL_10_19] = (child_network_adults * n_10_19) / n_adult; + p_adultPerWorkNetwork[Workplace::WORKPLACE_70_79] = (elderly_network_adults * n_70_79) / n_adult; + p_adultPerWorkNetwork[Workplace::WORKPLACE_80_PLUS] = (elderly_network_adults * n_80_plus) / n_adult; + + // p of being adult is then the remaining probability + p_adultPerWorkNetwork[Workplace::WORKPLACE_ADULT] = 1.0 - std::accumulate(p_adultPerWorkNetwork.begin(), p_adultPerWorkNetwork.end(), 0.0); + + // Then convert to cumulative probability with an inclusive scan + exateppabm::util::inclusive_scan(p_adultPerWorkNetwork.begin(), p_adultPerWorkNetwork.end(), p_adultPerWorkNetwork.begin()); + + // make sure the top bracket ends in a value in case of floating point / rounding >= 1.0 + p_adultPerWorkNetwork[Workplace::WORKPLACE_80_PLUS] = 1.0; + + return p_adultPerWorkNetwork; +} + +WorkplaceUnderlyingType generateWorkplaceForIndividual(const demographics::Age age, std::array p_adult_workplace, std::mt19937_64 & rng) { + // Children, retired and elderly are assigned a network based on their age + if (age == demographics::Age::AGE_0_9) { + return Workplace::WORKPLACE_SCHOOL_0_9; + } else if (age == demographics::Age::AGE_10_19) { + return Workplace::WORKPLACE_SCHOOL_10_19; + } else if (age == demographics::Age::AGE_70_79) { + return Workplace::WORKPLACE_70_79; + } else if (age == demographics::Age::AGE_80) { + return Workplace::WORKPLACE_80_PLUS; + } else { + // Adults randomly sample using a cumulative probability distribution computed from population counts and model parameters + std::uniform_real_distribution work_network_dist(0.0, 1.0); + float r = work_network_dist(rng); + for (std::uint32_t i = 0; i < p_adult_workplace.size(); i++) { + if (r < p_adult_workplace[i]) { + return i; + } + } + throw std::runtime_error("@todo - invalid cumulative probability distribution for p_adult_workplace?"); + } } } // namespace population diff --git a/src/exateppabm/population.h b/src/exateppabm/population.h index 5e31425..c1034f8 100644 --- a/src/exateppabm/population.h +++ b/src/exateppabm/population.h @@ -1,5 +1,8 @@ #pragma once +#include #include +#include +#include #include "flamegpu/flamegpu.h" #include "exateppabm/demographics.h" #include "exateppabm/person.h" @@ -10,15 +13,15 @@ namespace exateppabm { namespace population { /** - * Generate a population of individual person agents for a given simulation configuration + * Define FLAME GPU init function to generate a population of agents for a simulation. + * + * This method additionally stores several values in an anonymous namespace for use during the init function (which cannot take additional parameters) * * @param model the model the population is to be associated with (to get information about the person agent type from) - * @param config the model parameters struct for this simulation + * @param params the model parameters struct for this simulation * @param verbose if verbose output is enabled - * @param env_width the 2D environment width used for temporary behaviour / initialisation. To be removed once networks added. - * @param interactionRadius Temporary interaction radius for this simulation. To be removed */ -std::unique_ptr generate(flamegpu::ModelDescription& model, const exateppabm::input::config config, const bool verbose, const float env_width, const float interactionRadius); +void define(flamegpu::ModelDescription& model, const exateppabm::input::config params, const bool verbose); /** * Get the number of agents per demographic which were initialised to be infected, for the most recent call to generate. @@ -32,6 +35,103 @@ std::unique_ptr generate(flamegpu::ModelDescription& mode */ std::array getPerDemographicInitialInfectionCount(); +/** + * Type definition for the type used for number of individuals within a household. + * @todo - move this to a more sensible location? + */ +typedef std::uint8_t HouseholdSizeType; + +/** + * Structure containing data per household (total size, size of each age demographic) + * + */ +struct HouseholdStructure { + /** + * total number of individuals + */ + HouseholdSizeType size = 0; + /** + * Number of individuals in each age demographic + * @todo - should this be a map indexed by the enum? Enums in C++17 aren't the best. + */ + std::array sizePerAge = {}; + /** + * Age demographic per member of the household + */ + std::vector agePerPerson = {}; +}; + +/** + * Get the reference mean household size from the simulation parameters + * @param param simulation parameters + * @return target mean household size + */ +double getReferenceMeanHouseholdSize(const exateppabm::input::config& params); + +/** + * Generate a cumulative probability distribution for household size from reference data + * @param params model configuration parameters + * @return per-household-size cumulative probability, for sampling with a uniform distribution [0, 1) + */ +std::vector getHouseholdSizeCumulativeProbabilityVector(const exateppabm::input::config& params); + +/** + * Generate household structures, including the number of people of each age within the household + * + * This is only included in the header file for testing. + * + * @todo this should try and match target household demographic reference data in some way. I.e. don't make a house of 6*0-9 year olds + * @note this supports a maximum household size of 6, rather than allowing "6+" + * + * @param config simulation parameters structure + * @param rng RNG generator (pre-seeded) + * @param verbose if verbose output should be produced + * @return vector of per-household age demographic counts. + */ + +std::vector generateHouseholdStructures(const exateppabm::input::config params, std::mt19937_64 & rng, const bool verbose); + + +/** + * Underlying type for workplace, to avoid excessive casting of an enum class due to FLAME GPU 2 templating + */ +typedef std::uint8_t WorkplaceUnderlyingType; + +/** + * The number of workplace networks, when a fixed number of workplace networks are used + */ +constexpr std::uint8_t WORKPLACE_COUNT = 5; + +/** + * Enum for workplace type for a person. + * @todo - Find a nice way to make this an enum class (i.e. scoped) but still usable in FLAME GPU templated methods. Possibly implement to_underlying(Age)/from_underlying(Age) + */ +enum Workplace: WorkplaceUnderlyingType { + WORKPLACE_SCHOOL_0_9 = 0, + WORKPLACE_SCHOOL_10_19 = 1, + WORKPLACE_ADULT = 2, + WORKPLACE_70_79 = 3, + WORKPLACE_80_PLUS = 4 +}; + +/** + * Generate a cumulative probability distribution adults within each workplace from generated populations + * @param child_network_adults ratio of children to adults in the children networks + * @param elderly_network_adults ratio of elderly to adults in the elderly networks + * @param n_per_age number of individuals per age demographic + * @return per-household-size cumulative probability, for sampling with a uniform distribution [0, 1) + */ +std::array getAdultWorkplaceCumulativeProbabilityArray(const double child_network_adults, const double elderly_network_adults, std::array n_per_age); + + +/** + * For a given age demographic, assign a the workplace network given the age, per network adult cumulative probability & rng elements/state. + * @param age age for the individual + * @param p_adult_workplace cumulative probability for adults to be assigned to each workplace + * @param rng RNG generator (pre-seeded) + */ + +WorkplaceUnderlyingType generateWorkplaceForIndividual(const demographics::Age age, std::array p_adult_workplace, std::mt19937_64 & rng); } // namespace population diff --git a/src/exateppabm/util.h b/src/exateppabm/util.h index 57e0194..4bedde0 100644 --- a/src/exateppabm/util.h +++ b/src/exateppabm/util.h @@ -3,6 +3,8 @@ #include #include #include +#include +#include namespace exateppabm { namespace util { @@ -55,40 +57,74 @@ bool getSeatbeltsEnabled(); */ std::string getCMakeBuildType(); +/** + * Inclusive scan from [first, last) writing into d_first (and onwards), for libstdc++ which does not support c++17 (i.e. GCC 8) + * + * equivalent to std::inclusive_scan(first, last, d_first); + * + * @note - this is very naive, and will happily do out of bounds reads/writes. Using GCC >= 8 is greatly preferred + * + * @param first - input iterator for the first element + * @param last - input iterator for the last element + * @param d_first destination iterator for the last element + */ +template +void naive_inclusive_scan(InputIt first, InputIt last, OutputIt d_first) { + using T = typename std::iterator_traits::value_type; + T acc = T(); // default initialise T (the inferred type) + for (; first != last; ++first, ++d_first) { + acc += *first; + *d_first = acc; + } +} /** - * Perform an inclusive scan on an interable (i.e. std::array), without using std::inclusive_scan which although part of c++17 is not available in some older stdlib implementations (i.e. gcc 8). + * Inclusive scan from [first, last) writing into d_first (and onwards), using std::inclusive_scan if possible, else a naive implementation. + * + * equivalent to std::inclusive_scan(first, last, d_first); + * + * @param first - input iterator for the first element + * @param last - input iterator for the last element + * @param d_first destination iterator for the last element */ +template +void inclusive_scan(InputIt first, InputIt last, OutputIt d_first) { +#if defined(EXATEPP_ABM_USE_STD_INCLUSIVE_SCAN) && EXATEPP_ABM_USE_STD_INCLUSIVE_SCAN + std::inclusive_scan(first, last, d_first); +#else + naive_inclusive_scan(first, last, d_first); +#endif // defined(EXATEPP_ABM_USE_STD_INCLUSIVE_SCAN) && EXATEPP_ABM_USE_STD_INCLUSIVE_SCAN +} /** - * In-place inclusive scan, for libstdc++ which does not support c++17 (i.e. GCC 8) + * reduce elements of a container, i.e. std::reduce, for libstdc++ which does not support c++17 (i.e. GCC 8) * - * equivalent to std::inclusive_scan(container.begin(), container.end(), container.begin()); + * equivalent to std::reduce(first, last, init); * - * @param container iterable container / something with size() and operator[] + * @param first - input iterator for the first element + * @param last - input iterator for the last element + * @param init - initial value to reduce into (and inferred type) + * @return reduction (sum) if values between in [first, last) */ -template -void naive_inplace_inclusive_scan(T& container) { - if (container.size() <= 1) { - return; - } - // Naive in-place inclusive scan for libstc++8 - for (size_t i = 1; i < container.size(); ++i) { - container[i] = container[i - 1] + container[i]; +template +T naive_reduce(InputIt first, InputIt last, T init) { + for (; first != last; ++first) { + init += *first; } + return init; } /** - * In-place inclusive scan, using std::inclusive_scan if possible, else a naive implementation. + * reduce elements of a container, using std::reduce if possible, else a naive implementation. * - * equivalent to std::inclusive_scan(container.begin(), container.end(), container.begin()); + * equivalent to std::reduce(first, last, init); */ -template -void inplace_inclusive_scan(T& container) { +template +T reduce(InputIt first, InputIt last, T init) { #if defined(EXATEPP_ABM_USE_STD_INCLUSIVE_SCAN) && EXATEPP_ABM_USE_STD_INCLUSIVE_SCAN - std::inclusive_scan(container.begin(), container.end(), container.begin()); + return std::reduce(first, last, init); #else - naive_inplace_inclusive_scan(container); + return naive_reduce(first, last, init); #endif // defined(EXATEPP_ABM_USE_STD_INCLUSIVE_SCAN) && EXATEPP_ABM_USE_STD_INCLUSIVE_SCAN } diff --git a/src/exateppabm/visualisation.cu b/src/exateppabm/visualisation.cu index 07708ff..876c5b8 100644 --- a/src/exateppabm/visualisation.cu +++ b/src/exateppabm/visualisation.cu @@ -1,7 +1,10 @@ #include "exateppabm/visualisation.h" #include +#include #include +#include +#include #include "flamegpu/flamegpu.h" #include "exateppabm/person.h" #include "exateppabm/disease/SEIR.h" @@ -94,5 +97,48 @@ void join() { #endif // FLAMEGPU_VISUALISATION } +std::tuple getAgentXYZ(const std::uint32_t householdCount, const std::uint32_t householdIdx, const std::uint8_t idxWithinHousehold) { +#if defined(FLAMEGPU_VISUALISATION) + // Use the number of households to figure out the size of a 2D grid for visualisation purposes + std::uint64_t visHouseholdGridwidth = static_cast(std::ceil(std::sqrt(static_cast(householdCount)))); + // Prep a vector of integers to find the location within a household for each individual + static auto visAssignedHouseholdCount = std::vector(householdCount, 0); + // pre-calculate spatial offset per individual within household, for upto 6 individuals per household (current hardcoded upper limit) + constexpr float OFFSET_SF = 0.7f; + static std::array visHouseholdOffsetX = {{0.f + , OFFSET_SF * static_cast(std::sin(0 * M_PI / 180.0)) + , OFFSET_SF * static_cast(std::sin(72 * M_PI / 180.0)) + , OFFSET_SF * static_cast(std::sin(144 * M_PI / 180.0)) + , OFFSET_SF * static_cast(std::sin(216 * M_PI / 180.0)) + , OFFSET_SF * static_cast(std::sin(288 * M_PI / 180.0))}}; + static std::array visHouseholdOffsetY = {{0.f + , OFFSET_SF * static_cast(std::cos(0 * M_PI / 180.0)) + , OFFSET_SF * static_cast(std::cos(72 * M_PI / 180.0)) + , OFFSET_SF * static_cast(std::cos(144 * M_PI / 180.0)) + , OFFSET_SF * static_cast(std::cos(216 * M_PI / 180.0)) + , OFFSET_SF * static_cast(std::cos(288 * M_PI / 180.0))}}; + + // Get teh x/y/z for the provided individual within their household + // Get the center point of their household + std::uint32_t visHouseholdRow = householdIdx / visHouseholdGridwidth; + std::uint32_t visHouseholdCol = householdIdx % visHouseholdGridwidth; + // Get their index within their household [0-5] + std::uint8_t idxInHouse = visAssignedHouseholdCount[householdIdx]; + visAssignedHouseholdCount[householdIdx]++; + // Get their arbitrary offset, given a vector of offsets (6 potential values) + constexpr float VIS_HOUSE_GRID_SPACING = 2.5f; + float visX = (visHouseholdCol * VIS_HOUSE_GRID_SPACING) + visHouseholdOffsetX[idxInHouse % visHouseholdOffsetX.size()]; + float visY = (visHouseholdRow * VIS_HOUSE_GRID_SPACING) + visHouseholdOffsetY[idxInHouse % visHouseholdOffsetY.size()]; + float visZ = 0; + // Set the x,y and z in agent data. These must be floats. + // person.setVariable(exateppabm::person::v::x, visX); + // person.setVariable(exateppabm::person::v::y, visY); + // person.setVariable(exateppabm::person::v::z, visZ); + return std::tuple(visX, visY, visZ); +#else + return std::tuple(0.0, 0.0, 0.0); +#endif // FLAMEGPU_VISUALISATION +} + } // namespace visualisation } // namespace exateppabm diff --git a/src/exateppabm/visualisation.h b/src/exateppabm/visualisation.h index 5d5f16f..3e23a88 100644 --- a/src/exateppabm/visualisation.h +++ b/src/exateppabm/visualisation.h @@ -1,6 +1,10 @@ #pragma once +#include +#include + #include "flamegpu/flamegpu.h" +#include "exateppabm/input.h" namespace exateppabm { namespace visualisation { @@ -21,5 +25,17 @@ void setup(bool enable, flamegpu::ModelDescription& model, flamegpu::CUDASimulat */ void join(); +/** + * Get visualistion properties for an agent, given model parameters and thier household index + * + * This has had to be refactored to workaround a flame gpu limitation + * + * @param householdCount the number of generated households + * @param householdIdx, the index of the current household + * @param idxWithinHousehold the index within the household for the agent. + * @return tuple of agent properties + */ +std::tuple getAgentXYZ(const std::uint32_t householdCount, const std::uint32_t householdIdx, const std::uint8_t idxWithinHousehold); + } // namespace visualisation } // namespace exateppabm diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 2950095..b49016e 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -32,6 +32,8 @@ include(${FLAMEGPU_ROOT}/cmake/common.cmake) # List test source files set(TESTS_SRC ${CMAKE_CURRENT_SOURCE_DIR}/src/main.cu + ${CMAKE_CURRENT_SOURCE_DIR}/src/exateppabm/test_network.cu + ${CMAKE_CURRENT_SOURCE_DIR}/src/exateppabm/test_population.cu ${CMAKE_CURRENT_SOURCE_DIR}/src/exateppabm/test_util.cu ) diff --git a/tests/src/exateppabm/test_network.cu b/tests/src/exateppabm/test_network.cu new file mode 100644 index 0000000..2d4e1b1 --- /dev/null +++ b/tests/src/exateppabm/test_network.cu @@ -0,0 +1,161 @@ +#include +#include +#include +#include +#include + +#include "gtest/gtest.h" + +#include "fmt/core.h" +#include "flamegpu/flamegpu.h" +#include "exateppabm/network.h" + +/** + * Test fully connected undirected network generation + * + * @todo test more edge cases + */ +TEST(TestNetwork, generateFullyConnectedUndirectedGraph) { + // Generate nodes from [1, 4] + constexpr std::uint32_t N = 4; + std::vector nodes(N); + std::iota(nodes.begin(), nodes.end(), 1); + + // Generate a fully connected network + exateppabm::network::UndirectedGraph network = exateppabm::network::generateFullyConnectedUndirectedGraph(nodes); + // Check the number of nodes and edges are correct + EXPECT_EQ(network.getNumVertices(), N); + EXPECT_EQ(network.getNumVertices(), nodes.size()); + EXPECT_EQ(network.getNumEdges(), (N * (N-1)) / 2); + // Check the edges are the expected values. +} + +/** + * Test watts strogatz undirected network generation + * + * @todo - edge case testing. + */ +TEST(TestNetwork, generateSmallWorldUndirectedGraph) { + // Generate nodes from [1, 16] + constexpr std::uint32_t N = 16; + std::vector nodes(N); + std::iota(nodes.begin(), nodes.end(), 1); + + // Generate a small world network + constexpr std::uint32_t K = 4; // degree 4 + constexpr double p_rewire = 0.1; // probability of a rewire + std::mt19937_64 rng(12); // seeded mersenne twister rng engine + exateppabm::network::UndirectedGraph network = exateppabm::network::generateSmallWorldUndirectedGraph(nodes, K, p_rewire, rng); + + // Ensure there are the correct number of nodes and edges still + EXPECT_EQ(network.getNumVertices(), nodes.size()); + EXPECT_EQ(network.getNumEdges(), (N * K) / 2); + // @todo - validate the in degree and out degree of each node. + // If graph is directed, out degree will still be K, but indegree will vary per node + // If graph is undirected, out degree and in degree would be the same, but mean degree will be K. + // @todo - decide which is intended in this case. + + // @todo - actually test the network properties. Average degree etc. + size_t edgeIdx = 0; + for (const auto &edge : network.getEdges()) { + // Every source and ever dest should be a valid node index + EXPECT_LT(edge.source, network.getNumVertices()); + EXPECT_LT(edge.dest, network.getNumVertices()); + // No self-edges + EXPECT_NE(edge.source, edge.dest); + + // fmt::print("edges[{}] = ID {} -> {} (idx {} -> {})\n", edgeIdx, network.nodes.at(edge.source), network.nodes.at(edge.dest), edge.source, edge.dest); + ++edgeIdx; + } + // @todo - check for duplicates + + // Check invalid waltz strogatz parameters behave as intended + // Odd degree - actually use K-1, + auto n_k_odd = exateppabm::network::generateSmallWorldUndirectedGraph({{1, 2, 3, 4}}, 3, 0.1, rng); + EXPECT_EQ(n_k_odd.getNumVertices(), 4); + EXPECT_EQ(n_k_odd.getNumEdges(), ((3-1) * n_k_odd.getNumVertices())/2); + // too high of a degree + EXPECT_ANY_THROW(exateppabm::network::generateSmallWorldUndirectedGraph({{1, 2, 3, 4}}, 5, 0.1, rng)); + // too low probability + EXPECT_ANY_THROW(exateppabm::network::generateSmallWorldUndirectedGraph({{1, 2, 3, 4}}, 2, -1.0, rng)); + // too high probability + EXPECT_ANY_THROW(exateppabm::network::generateSmallWorldUndirectedGraph({{1, 2, 3, 4}}, 2, 2.0, rng)); + + // Check the networks contain the expected number of edges for edge case node counts + // 0 nodes, no edges + auto n0 = exateppabm::network::generateSmallWorldUndirectedGraph({}, 2, 0.1, rng); + EXPECT_EQ(n0.getNumVertices(), 0u); + EXPECT_EQ(n0.getNumEdges(), 0u); + // 1 node, no edges + auto n1 = exateppabm::network::generateSmallWorldUndirectedGraph({{1}}, 2, 0.1, rng); + EXPECT_EQ(n1.getNumVertices(), 1u); + EXPECT_EQ(n1.getNumEdges(), 0u); + // 2 node, fully connected + auto n2 = exateppabm::network::generateSmallWorldUndirectedGraph({{1, 2}}, 2, 0.1, rng); + EXPECT_EQ(n2.getNumVertices(), 2u); + EXPECT_EQ(n2.getNumEdges(), 1u); + // 3 nodes, mean degree 2, 3 edges + auto n3 = exateppabm::network::generateSmallWorldUndirectedGraph({{1, 2, 3}}, 2, 0.1, rng); + EXPECT_EQ(n3.getNumVertices(), 3u); + EXPECT_EQ(n3.getNumEdges(), 3u); + // 4 nodes, degree 2, 4 edges + auto n4_2 = exateppabm::network::generateSmallWorldUndirectedGraph({{1, 2, 3, 4}}, 2, 0.1, rng); + EXPECT_EQ(n4_2.getNumVertices(), 4u); + EXPECT_EQ(n4_2.getNumEdges(), 4u); + // 4 nodes, degree 4, fully connected, so only 6 edges + auto n4_4 = exateppabm::network::generateSmallWorldUndirectedGraph({{1, 2, 3, 4}}, 4, 0.1, rng); + EXPECT_EQ(n4_4.getNumVertices(), 4u); + EXPECT_EQ(n4_4.getNumEdges(), 6u); + + // 12 nodes, degree 2 + auto n12_2 = exateppabm::network::generateSmallWorldUndirectedGraph({{1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12}}, 2, 0.1, rng); + EXPECT_EQ(n12_2.getNumVertices(), 12u); + EXPECT_EQ(n12_2.getNumEdges(), 12u); + // 12 nodes, degree 2 (equiv to 2) + auto n12_3 = exateppabm::network::generateSmallWorldUndirectedGraph({{1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12}}, 3, 0.1, rng); + EXPECT_EQ(n12_3.getNumVertices(), 12u); + EXPECT_EQ(n12_3.getNumEdges(), 12u); + // 12 nodes, degree 4 + auto n12_4 = exateppabm::network::generateSmallWorldUndirectedGraph({{1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12}}, 4, 0.1, rng); + EXPECT_EQ(n12_4.getNumVertices(), 12u); + EXPECT_EQ(n12_4.getNumEdges(), 24u); + // 12 nodes, degree 5 (equiv to 5) + auto n12_5 = exateppabm::network::generateSmallWorldUndirectedGraph({{1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12}}, 5, 0.1, rng); + EXPECT_EQ(n12_5.getNumVertices(), 12u); + EXPECT_EQ(n12_5.getNumEdges(), 24u); + // 12 nodes, degree 6 + auto n12_6 = exateppabm::network::generateSmallWorldUndirectedGraph({{1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12}}, 6, 0.1, rng); + EXPECT_EQ(n12_6.getNumVertices(), 12u); + EXPECT_EQ(n12_6.getNumEdges(), 36u); + // 12 nodes, degree 8 + auto n12_8 = exateppabm::network::generateSmallWorldUndirectedGraph({{1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12}}, 8, 0.1, rng); + EXPECT_EQ(n12_8.getNumVertices(), 12u); + EXPECT_EQ(n12_8.getNumEdges(), 48u); + // 12 nodes, degree 8 + auto n12_10 = exateppabm::network::generateSmallWorldUndirectedGraph({{1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12}}, 10, 0.1, rng); + EXPECT_EQ(n12_10.getNumVertices(), 12u); + EXPECT_EQ(n12_10.getNumEdges(), 60u); + // 12 nodes, degree 8 + auto n12_12 = exateppabm::network::generateSmallWorldUndirectedGraph({{1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12}}, 12, 0.1, rng); + EXPECT_EQ(n12_12.getNumVertices(), 12u); + EXPECT_EQ(n12_12.getNumEdges(), 66u); // fully connected + + // 8 nodes, degree 4, no rewiring + auto n8_4_0 = exateppabm::network::generateSmallWorldUndirectedGraph({{1, 2, 3, 4, 5, 6, 7, 8}}, 4, 0.0, rng); + EXPECT_EQ(n8_4_0.getNumVertices(), 8u); + EXPECT_EQ(n8_4_0.getNumEdges(), 16u); + // Check for valid edge, ascending internal order, ensure no dupes. + std::set uniqueEdges_n8_4_0 = {{}}; + for (const auto & edge : n8_4_0.getEdges()) { + EXPECT_GE(edge.source, 0u); + EXPECT_LT(edge.source, 8u); + EXPECT_NE(edge.source, edge.dest); + // Ensure the source node index is less than the dest index (i.e. undirected edges are always stored with the source having the lower index) + EXPECT_LT(edge.source, edge.dest); + // Edges should be unique, so this edge should not have been seen yet. + EXPECT_EQ(uniqueEdges_n8_4_0.count(edge), 0); + uniqueEdges_n8_4_0.insert(edge); + // with rewiring 0, the dest should always be source +1 or source +2 or when wrapping the boundary source will be dest -1 or dest -2 wrapped + 2) + EXPECT_TRUE((edge.dest == edge.source + 1) || (edge.dest == edge.source + 2) || (edge.source == (edge.dest + 1) % n8_4_0.getNumVertices()) || (edge.source == (edge.dest + 2) % n8_4_0.getNumVertices())); + } +} diff --git a/tests/src/exateppabm/test_population.cu b/tests/src/exateppabm/test_population.cu new file mode 100644 index 0000000..046f7d0 --- /dev/null +++ b/tests/src/exateppabm/test_population.cu @@ -0,0 +1,374 @@ +#include +#include +#include +#include +#include + +#include "gtest/gtest.h" + +#include "fmt/core.h" +#include "exateppabm/input.h" +#include "exateppabm/population.h" + +/** + * Test getReferenceMeanHouseholdSize using a number of manual, partially filled config objects. + */ +TEST(TestPopulation, getReferenceMeanHouseholdSize) { + const double EPSILON = 1.0e-7; + double value = 0.0; + double expected = 0.0; + exateppabm::input::config params = {}; + + // balanced house sizes, 21 people in 6 houses + params.household_size_1 = 1; + params.household_size_2 = 1; + params.household_size_3 = 1; + params.household_size_4 = 1; + params.household_size_5 = 1; + params.household_size_6 = 1; + expected = 21 / 6.0; + value = exateppabm::population::getReferenceMeanHouseholdSize(params); + ASSERT_NEAR(value, expected, EPSILON); + + // Just houses sizes of 1 + params.household_size_1 = 1; + params.household_size_2 = 0; + params.household_size_3 = 0; + params.household_size_4 = 0; + params.household_size_5 = 0; + params.household_size_6 = 0; + expected = 1.0 / 1.0; + value = exateppabm::population::getReferenceMeanHouseholdSize(params); + ASSERT_NEAR(value, expected, EPSILON); + + // Just houses sizes of 6 + params.household_size_1 = 0; + params.household_size_2 = 0; + params.household_size_3 = 0; + params.household_size_4 = 0; + params.household_size_5 = 0; + params.household_size_6 = 1; + expected = 6.0 / 1.0; + value = exateppabm::population::getReferenceMeanHouseholdSize(params); + ASSERT_NEAR(value, expected, EPSILON); + + // Just houses sizes of 1 or 6 + params.household_size_1 = 1; + params.household_size_2 = 0; + params.household_size_3 = 0; + params.household_size_4 = 0; + params.household_size_5 = 0; + params.household_size_6 = 1; + expected = 7 / 2.0; + value = exateppabm::population::getReferenceMeanHouseholdSize(params); + ASSERT_NEAR(value, expected, EPSILON); + + // Arbitrary mixed values + params.household_size_1 = 1; + params.household_size_2 = 2; + params.household_size_3 = 3; + params.household_size_4 = 3; + params.household_size_5 = 2; + params.household_size_6 = 1; + expected = 42 / 12.0; + value = exateppabm::population::getReferenceMeanHouseholdSize(params); + ASSERT_NEAR(value, expected, EPSILON); +} + +/** + * Get Check the cumulative household size probability vector is as expected for sample model parameters + */ +TEST(TestPopulation, getHouseholdSizeCumulativeProbabilityVector) { + const double EPSILON = 1.0e-7; + std::vector cumulativeP = {}; + std::vector expected = {}; + constexpr std::uint32_t EXPECTED_LENGTH = 6; + exateppabm::input::config params = {}; + + // balanced house sizes equal probability of each + params.household_size_1 = 1; + params.household_size_2 = 1; + params.household_size_3 = 1; + params.household_size_4 = 1; + params.household_size_5 = 1; + params.household_size_6 = 1; + expected = {{ 1/6.0, 2/6.0, 3/6.0, 4/6.0, 5/6.0, 6/6.0 }}; + cumulativeP = exateppabm::population::getHouseholdSizeCumulativeProbabilityVector(params); + // Check the size + ASSERT_EQ(cumulativeP.size(), EXPECTED_LENGTH); + // Check the final element is approx equal 1.0 + ASSERT_NEAR(cumulativeP[5], 1.0, EPSILON); + // Check each element is as expected + for (std::size_t idx = 0; idx < cumulativeP.size(); ++idx) { + ASSERT_NEAR(cumulativeP[idx], expected[idx], EPSILON); + } + + // 1 or 6 only, balanced + params.household_size_1 = 1; + params.household_size_2 = 0; + params.household_size_3 = 0; + params.household_size_4 = 0; + params.household_size_5 = 0; + params.household_size_6 = 1; + expected = {{ 1/2.0, 1/2.0, 1/2.0, 1/2.0, 1/2.0, 2/2.0 }}; + cumulativeP = exateppabm::population::getHouseholdSizeCumulativeProbabilityVector(params); + // Check the size + ASSERT_EQ(cumulativeP.size(), EXPECTED_LENGTH); + // Check the final element is approx equal 1.0 + ASSERT_NEAR(cumulativeP[5], 1.0, EPSILON); + // Check each element is as expected + for (std::size_t idx = 0; idx < cumulativeP.size(); ++idx) { + ASSERT_NEAR(cumulativeP[idx], expected[idx], EPSILON); + } + + // imbalanced values + params.household_size_1 = 0; + params.household_size_2 = 1; + params.household_size_3 = 2; + params.household_size_4 = 3; + params.household_size_5 = 2; + params.household_size_6 = 1; + expected = {{ 0/9.0, 1/9.0, 3/9.0, 6/9.0, 8/9.0, 9/9.0 }}; + cumulativeP = exateppabm::population::getHouseholdSizeCumulativeProbabilityVector(params); + // Check the size + ASSERT_EQ(cumulativeP.size(), EXPECTED_LENGTH); + // Check the final element is approx equal 1.0 + ASSERT_NEAR(cumulativeP[5], 1.0, EPSILON); + // Check each element is as expected + for (std::size_t idx = 0; idx < cumulativeP.size(); ++idx) { + ASSERT_NEAR(cumulativeP[idx], expected[idx], EPSILON); + } +} + +/** + * Test that generateHouseholdStructures produces the expected data structures + * + * + */ +TEST(TestPopulation, generateHouseholdStructures) { + const double EPSILON = 1.0e-7; + std::vector households; + std::mt19937_64 rng = std::mt19937_64(0); + exateppabm::input::config params = {}; + std::uint64_t totalPeople = 0; + + // 1 house of 1 person in the 0-9 age category + params.rng_seed = 0; + params.n_total = 1u; + params.household_size_1 = 1; + params.household_size_2 = 0; + params.household_size_3 = 0; + params.household_size_4 = 0; + params.household_size_5 = 0; + params.household_size_6 = 0; + params.population_0_9 = 1; + params.population_10_19 = 0; + params.population_20_29 = 0; + params.population_30_39 = 0; + params.population_40_49 = 0; + params.population_50_59 = 0; + params.population_60_69 = 0; + params.population_70_79 = 0; + params.population_80 = 0; + rng.seed(params.rng_seed); + households = exateppabm::population::generateHouseholdStructures(params, rng, false); + // Should be 1 household + ASSERT_EQ(households.size(), 1u); + totalPeople = 0u; + for (const auto & household : households) { + totalPeople += household.size; + // Each household should have a household size between 1 and 6 inclusive (in this case 1) + ASSERT_EQ(household.size, 1u); + // The total of the house size age distributions should be the same as the house size + std::uint64_t sizePerAgeTotal = 0u; + for (const auto & c : household.sizePerAge) { + sizePerAgeTotal += c; + } + ASSERT_EQ(sizePerAgeTotal, household.size); + // sizePerAge should be the same length as the number of age demographics + ASSERT_EQ(household.sizePerAge.size(), exateppabm::demographics::AGE_COUNT); + // There should also be a matching number of per-person ages + ASSERT_EQ(household.agePerPerson.size(), household.size); + // There should only be a person in the 0-9 age category + ASSERT_EQ(household.sizePerAge[0], 1u); + ASSERT_EQ(household.sizePerAge[1], 0u); + ASSERT_EQ(household.sizePerAge[2], 0u); + ASSERT_EQ(household.sizePerAge[3], 0u); + ASSERT_EQ(household.sizePerAge[4], 0u); + ASSERT_EQ(household.sizePerAge[5], 0u); + ASSERT_EQ(household.sizePerAge[6], 0u); + ASSERT_EQ(household.sizePerAge[7], 0u); + } + // Should be 1 person in total + ASSERT_EQ(totalPeople, 1u); + + // 32 people, in houses of 1, 2 or 3 people, aged 70-79 or 80 + params.rng_seed = 0; + params.n_total = 32u; + params.household_size_1 = 1; + params.household_size_2 = 1; + params.household_size_3 = 1; + params.household_size_4 = 0; + params.household_size_5 = 0; + params.household_size_6 = 0; + params.population_0_9 = 0; + params.population_10_19 = 0; + params.population_20_29 = 0; + params.population_30_39 = 0; + params.population_40_49 = 0; + params.population_50_59 = 0; + params.population_60_69 = 0; + params.population_70_79 = 1; + params.population_80 = 1; + rng.seed(params.rng_seed); + households = exateppabm::population::generateHouseholdStructures(params, rng, false); + // Should be between 1 and 32 households + ASSERT_GE(households.size(), 1u); + ASSERT_LE(households.size(), 32u); + totalPeople = 0u; + for (const auto & household : households) { + totalPeople += household.size; + // Each household should have a household size between 1 and 3 inclusive (in this case 1) + ASSERT_GE(household.size, 1u); + ASSERT_LE(household.size, 3u); + // The total of the house size age distributions should be the same as the house size + std::uint64_t sizePerAgeTotal = 0u; + for (const auto & c : household.sizePerAge) { + sizePerAgeTotal += c; + } + ASSERT_EQ(sizePerAgeTotal, household.size); + // sizePerAge should be the same length as the number of age demographics + ASSERT_EQ(household.sizePerAge.size(), exateppabm::demographics::AGE_COUNT); + // There should also be a matching number of per-person ages + ASSERT_EQ(household.agePerPerson.size(), household.size); + // There should only be people in the 70-79 and 80 age categories + ASSERT_EQ(household.sizePerAge[0], 0u); + ASSERT_EQ(household.sizePerAge[1], 0u); + ASSERT_EQ(household.sizePerAge[2], 0u); + ASSERT_EQ(household.sizePerAge[3], 0u); + ASSERT_EQ(household.sizePerAge[4], 0u); + ASSERT_EQ(household.sizePerAge[5], 0u); + ASSERT_LE(household.sizePerAge[6], household.size); + ASSERT_LE(household.sizePerAge[7], household.size); + } + // Should be 32 person in total + ASSERT_EQ(totalPeople, params.n_total); +} + + +/** + * Test getAdultWorkplaceCumulativeProbabilityArray + */ + +TEST(TestPopulation, getAdultWorkplaceCumulativeProbabilityArray) { +const double EPSILON = 1.0e-7; +double child_network_adults = 0; +double elderly_network_adults = 0; +std::array n_per_age = {}; +std::array p = {}; +std::array expected = {}; + + +// With 0 adults in the child or elderly networks, should all be assigned to the adult network, regardless of population counts +child_network_adults = 0; +elderly_network_adults = 0; +n_per_age = {{10, 10, 10, 10, 10, 10, 10, 10, 10}}; +expected = {{0.0, 0.0, 1.0, 1.0, 1.0}}; +p = exateppabm::population::getAdultWorkplaceCumulativeProbabilityArray(child_network_adults, elderly_network_adults, n_per_age); +EXPECT_NEAR(p[exateppabm::population::WORKPLACE_SCHOOL_0_9], expected[exateppabm::population::WORKPLACE_SCHOOL_0_9], EPSILON); +EXPECT_NEAR(p[exateppabm::population::WORKPLACE_SCHOOL_10_19], expected[exateppabm::population::WORKPLACE_SCHOOL_10_19], EPSILON); +EXPECT_NEAR(p[exateppabm::population::WORKPLACE_ADULT], expected[exateppabm::population::WORKPLACE_ADULT], EPSILON); +EXPECT_NEAR(p[exateppabm::population::WORKPLACE_70_79], expected[exateppabm::population::WORKPLACE_70_79], EPSILON); +EXPECT_NEAR(p[exateppabm::population::WORKPLACE_80_PLUS], expected[exateppabm::population::WORKPLACE_80_PLUS], EPSILON); + +// With 1 adults per member of work networks, and a population which will support this, generate probabilities that would lead to 20% of adults being assigned to each network +child_network_adults = 1.0; +elderly_network_adults = 1.0; +// I.e. 20 children/elderly per network and 100 adults total, 20 in each network in the end +n_per_age = {{20, 20, 20, 20, 20, 20, 20, 20, 20}}; +expected = {{0.2, 0.4, 0.6, 0.8, 1.0}}; +p = exateppabm::population::getAdultWorkplaceCumulativeProbabilityArray(child_network_adults, elderly_network_adults, n_per_age); +EXPECT_NEAR(p[exateppabm::population::WORKPLACE_SCHOOL_0_9], expected[exateppabm::population::WORKPLACE_SCHOOL_0_9], EPSILON); +EXPECT_NEAR(p[exateppabm::population::WORKPLACE_SCHOOL_10_19], expected[exateppabm::population::WORKPLACE_SCHOOL_10_19], EPSILON); +EXPECT_NEAR(p[exateppabm::population::WORKPLACE_ADULT], expected[exateppabm::population::WORKPLACE_ADULT], EPSILON); +EXPECT_NEAR(p[exateppabm::population::WORKPLACE_70_79], expected[exateppabm::population::WORKPLACE_70_79], EPSILON); +EXPECT_NEAR(p[exateppabm::population::WORKPLACE_80_PLUS], expected[exateppabm::population::WORKPLACE_80_PLUS], EPSILON); + +// With 10% in child networks, and 20% in elderly networks, and 100 adults total, generate the expected (unbalanced) values +child_network_adults = 0.1; +elderly_network_adults = 0.2; +n_per_age = {{20, 20, 0, 0, 20, 0, 0, 20, 10}}; +// 2/20 adult, 2/20 adults, 10/20 adults, 4/20 adults, 2/20 adults, but cumulative. +expected = {{0.1, 0.2, 0.7, 0.9, 1.0}}; +p = exateppabm::population::getAdultWorkplaceCumulativeProbabilityArray(child_network_adults, elderly_network_adults, n_per_age); +EXPECT_NEAR(p[exateppabm::population::WORKPLACE_SCHOOL_0_9], expected[exateppabm::population::WORKPLACE_SCHOOL_0_9], EPSILON); +EXPECT_NEAR(p[exateppabm::population::WORKPLACE_SCHOOL_10_19], expected[exateppabm::population::WORKPLACE_SCHOOL_10_19], EPSILON); +EXPECT_NEAR(p[exateppabm::population::WORKPLACE_ADULT], expected[exateppabm::population::WORKPLACE_ADULT], EPSILON); +EXPECT_NEAR(p[exateppabm::population::WORKPLACE_70_79], expected[exateppabm::population::WORKPLACE_70_79], EPSILON); +EXPECT_NEAR(p[exateppabm::population::WORKPLACE_80_PLUS], expected[exateppabm::population::WORKPLACE_80_PLUS], EPSILON); + +// @todo - expand test with more edge case coverage +} + +/** + * Test generateWorkplaceForIndividual + */ +TEST(TestPopulation, generateWorkplaceForIndividual) { +// number of individuals to generate, to do a large enough sample to make sure returned values are in-range. +constexpr std::uint64_t n_samples = 1000u; +constexpr std::uint64_t n_samples_adult = 10000u; +// rng state +std::mt19937_64 rng(0u); +// probability array +std::array p_adult_workplace = {{0.1, 0.2, 0.8, 0.9, 1.0}}; +// Expected/target counts for adults (non-cumulative) +std::array target_adult_workplace = {{ + 0.1 * n_samples_adult, + 0.1 * n_samples_adult, + 0.6 * n_samples_adult, + 0.1 * n_samples_adult, + 0.1 * n_samples_adult}}; +double TARGET_ADULT_EPSILON = 0.01 * n_samples_adult; + +// Repeatedly generate individuals of each age band from a given probability distribution, checking they are always the expected values + +// 0-9 always in their workplace +for (std::uint64_t idx = 0; idx < n_samples; ++idx) { + auto w = exateppabm::population::generateWorkplaceForIndividual(exateppabm::demographics::AGE_0_9, p_adult_workplace, rng); + EXPECT_EQ(w, exateppabm::population::Workplace::WORKPLACE_SCHOOL_0_9); +} +// 10-19 always in their workplace +for (std::uint64_t idx = 0; idx < n_samples; ++idx) { + auto w = exateppabm::population::generateWorkplaceForIndividual(exateppabm::demographics::AGE_10_19, p_adult_workplace, rng); + EXPECT_EQ(w, exateppabm::population::Workplace::WORKPLACE_SCHOOL_10_19); +} +// 70-79 always in their workplace +for (std::uint64_t idx = 0; idx < n_samples; ++idx) { + auto w = exateppabm::population::generateWorkplaceForIndividual(exateppabm::demographics::AGE_70_79, p_adult_workplace, rng); + EXPECT_EQ(w, exateppabm::population::Workplace::WORKPLACE_70_79); +} +// 80+ always in their workplace +for (std::uint64_t idx = 0; idx < n_samples; ++idx) { + auto w = exateppabm::population::generateWorkplaceForIndividual(exateppabm::demographics::AGE_80, p_adult_workplace, rng); + EXPECT_EQ(w, exateppabm::population::Workplace::WORKPLACE_80_PLUS); +} +// Adults are randomly assigned subject to the cumulative probability distribution. +std::array workplaceAdults = {0}; +for (std::uint64_t idx = 0; idx < n_samples_adult; ++idx) { + // Just check a single age demo for simlicity + auto w = exateppabm::population::generateWorkplaceForIndividual(exateppabm::demographics::AGE_40_49, p_adult_workplace, rng); + // Ensure in range + EXPECT_GE(w, 0); + EXPECT_LT(w, exateppabm::population::WORKPLACE_COUNT); + // Increment the counter + ++workplaceAdults[w]; +} +// Check that each counter is roughly correct. +EXPECT_NEAR(workplaceAdults[exateppabm::population::Workplace::WORKPLACE_SCHOOL_0_9], target_adult_workplace[exateppabm::population::Workplace::WORKPLACE_SCHOOL_0_9], TARGET_ADULT_EPSILON); +EXPECT_NEAR(workplaceAdults[exateppabm::population::Workplace::WORKPLACE_SCHOOL_10_19], target_adult_workplace[exateppabm::population::Workplace::WORKPLACE_SCHOOL_10_19], TARGET_ADULT_EPSILON); +EXPECT_NEAR(workplaceAdults[exateppabm::population::Workplace::WORKPLACE_ADULT], target_adult_workplace[exateppabm::population::Workplace::WORKPLACE_ADULT], TARGET_ADULT_EPSILON); +EXPECT_NEAR(workplaceAdults[exateppabm::population::Workplace::WORKPLACE_70_79], target_adult_workplace[exateppabm::population::Workplace::WORKPLACE_70_79], TARGET_ADULT_EPSILON); +EXPECT_NEAR(workplaceAdults[exateppabm::population::Workplace::WORKPLACE_80_PLUS], target_adult_workplace[exateppabm::population::Workplace::WORKPLACE_80_PLUS], TARGET_ADULT_EPSILON); + +// @todo - expand test with more edge case coverage +} diff --git a/tests/src/exateppabm/test_util.cu b/tests/src/exateppabm/test_util.cu index d0c1b9b..681d438 100644 --- a/tests/src/exateppabm/test_util.cu +++ b/tests/src/exateppabm/test_util.cu @@ -1,5 +1,6 @@ #include #include +#include #include #include "gtest/gtest.h" @@ -7,6 +8,79 @@ #include "fmt/core.h" #include "exateppabm/util.h" + +/** + * Check that getGPUName doesn't cause an error, and should not return the "unknown" string (although might, but in that case something is wrong with the system) + */ +TEST(TestUtil, getGPUName) { + // Method returns a string, shouldn't be unknown, as we assume there is atleast one gpu. + std::string gpuName = exateppabm::util::getGPUName(0); + ASSERT_NE(gpuName, ""); + ASSERT_NE(gpuName, "unknown"); + // No machines should include INT_MAX GPUs, so it should return "unknown" + std::string invalid = exateppabm::util::getGPUName(INT_MAX); + ASSERT_EQ(invalid, "unknown"); +} + +/** + * Test the getting of SM count method, as we cannot hardcode the expected number, we can just assume it should be non zero for device 0 + * This is a bit of a flakey test + */ +TEST(TestUtil, getGPUMultiProcessorCount) { + // Assuming the current machine has a gpu, should return a value >= 0 + auto result = exateppabm::util::getGPUMultiProcessorCount(0); + ASSERT_GE(result, 0); + // No machines should include INT_MAX GPUs, so it should return 0 + auto invalid = exateppabm::util::getGPUMultiProcessorCount(INT_MAX); + ASSERT_EQ(invalid, 0); +} + +/** + * Test that getting GPU memory capacity behaves roughly as expected, for sensible and not-sensible device ordinals. + */ +TEST(TestUtil, getGPUMemory) { + // Assuming the current machine has a gpu, should return a value >= 0 + auto result = exateppabm::util::getGPUMemory(0); + ASSERT_GE(result, 0u); + // No machines should include INT_MAX GPUs, so it should return "unknown" + auto invalid = exateppabm::util::getGPUMemory(INT_MAX); + ASSERT_EQ(invalid, 0u); +} + +/** + * Initialising the CUDA context has no direct effects, only side effects which are difficult to observe without digging into the CUDA Driver API. +* For now, just make sure the method doesn't trigger any execeptions + */ +TEST(TestUtil, initialiseCUDAContext) { + ASSERT_NO_THROW(exateppabm::util::initialiseCUDAContext(0)); +} + +/** + * Test that getting the seatbelts value behaves correctly + */ +TEST(TestUtil, getSeatbeltsEnabled) { + auto value = exateppabm::util::getSeatbeltsEnabled(); + #if !defined(FLAMEGPU_SEATBELTS) || FLAMEGPU_SEATBELTS + ASSERT_EQ(value, true); + #else // !defined(FLAMEGPU_SEATBELTS) || FLAMEGPU_SEATBELTS + ASSERT_EQ(value, false); + #endif // !defined(FLAMEGPU_SEATBELTS) || FLAMEGPU_SEATBELTS +} + +/** + * Test that getting the CMake build type value behaves correctly + */ +TEST(TestUtil, getCMakeBuildType) { + auto value = exateppabm::util::getCMakeBuildType(); +#if defined(CMAKE_BUILD_TYPE) + ASSERT_EQ(value, CMAKE_BUILD_TYPE); +#else // defined(CMAKE_BUILD_TYPE) + // If this occurs, the CMake is broken so we will trigger a failure. + ASSERT_FALSE(true); +#endif // defined(CMAKE_BUILD_TYPE) +} + + /** * Test the naive inplace inclusive_scan implementation for an array of integers */ @@ -15,7 +89,7 @@ TEST(TestUtil, naive_inplace_inclusive_scan_array_int) { std::array inout = {{1, 2, 3, 4}}; const std::array expected = {{1, 3, 6, 10}}; - exateppabm::util::naive_inplace_inclusive_scan(inout); + exateppabm::util::naive_inclusive_scan(inout.begin(), inout.end(), inout.begin()); for (std::uint32_t e = 0; e < inout.size(); e++) { EXPECT_EQ(inout[e], expected[e]); @@ -29,7 +103,7 @@ TEST(TestUtil, naive_inplace_inclusive_scan_vec_uint) { std::vector inout = {{1, 2, 3, 4}}; const std::vector expected = {{1, 3, 6, 10}}; - exateppabm::util::naive_inplace_inclusive_scan(inout); + exateppabm::util::naive_inclusive_scan(inout.begin(), inout.end(), inout.begin()); for (std::uint32_t e = 0; e < inout.size(); e++) { EXPECT_EQ(inout[e], expected[e]); @@ -43,7 +117,7 @@ TEST(TestUtil, naive_inplace_inclusive_scan_vec_double) { std::vector inout = {{1.1, 2.2, 3.3, 4.4}}; const std::vector expected = {{1.1, 3.3, 6.6, 11.0}}; - exateppabm::util::naive_inplace_inclusive_scan(inout); + exateppabm::util::naive_inclusive_scan(inout.begin(), inout.end(), inout.begin()); constexpr float EPSILON = 1.0e-6f; @@ -60,7 +134,7 @@ TEST(TestUtil, inplace_inclusive_scan_array_int) { std::array inout = {{1, 2, 3, 4}}; const std::array expected = {{1, 3, 6, 10}}; - exateppabm::util::inplace_inclusive_scan(inout); + exateppabm::util::inclusive_scan(inout.begin(), inout.end(), inout.begin()); for (std::uint32_t e = 0; e < inout.size(); e++) { EXPECT_EQ(inout[e], expected[e]); @@ -74,7 +148,7 @@ TEST(TestUtil, inplace_inclusive_scan_vec_uint) { std::vector inout = {{1, 2, 3, 4}}; const std::vector expected = {{1, 3, 6, 10}}; - exateppabm::util::inplace_inclusive_scan(inout); + exateppabm::util::inclusive_scan(inout.begin(), inout.end(), inout.begin()); for (std::uint32_t e = 0; e < inout.size(); e++) { EXPECT_EQ(inout[e], expected[e]); @@ -88,7 +162,7 @@ TEST(TestUtil, inplace_inclusive_scan_vec_double) { std::vector inout = {{1.1, 2.2, 3.3, 4.4}}; const std::vector expected = {{1.1, 3.3, 6.6, 11.0}}; - exateppabm::util::inplace_inclusive_scan(inout); + exateppabm::util::inclusive_scan(inout.begin(), inout.end(), inout.begin()); constexpr float EPSILON = 1.0e-6f; @@ -161,3 +235,119 @@ TEST(TestUtil, std_inclusive_scan_vec_double) { TEST(TestUtil, DISABLED_std_inclusive_scan_vec_double) { } #endif + + +/** + * Test the naive plus reduction an array of integers + */ +TEST(TestUtil, naive_reduce_array_int) { + constexpr std::uint32_t ELEMENTS = 4; + std::array in = {{1, 2, 3, 4}}; + const int expected = 10; + auto value = exateppabm::util::naive_reduce(in.begin(), in.end(), 0); + EXPECT_EQ(value, expected); +} + +/** + * Test the naive plus reduction a vector of unsigned integers + */ +TEST(TestUtil, naive_reduce_vec_uint) { + std::vector in = {{1, 2, 3, 4}}; + const std::uint32_t expected = 10u; + auto value = exateppabm::util::naive_reduce(in.begin(), in.end(), 0u); + EXPECT_EQ(value, expected); +} + +/** + * Test the naive plus reduction a vector of floats, using an arbitrary epsilon value which is good enough for the input values + */ +TEST(TestUtil, naive_reduce_vec_double) { + std::vector in = {{1.1, 2.2, 3.3, 4.4}}; + const double expected = 11.0; + auto value = exateppabm::util::naive_reduce(in.begin(), in.end(), 0.); + constexpr double EPSILON = 1.0e-6f; + EXPECT_NEAR(value, expected, EPSILON); +} + +/** + * Test the naive or std plus reduction an array of integers + */ +TEST(TestUtil, reduce_array_int) { + constexpr std::uint32_t ELEMENTS = 4; + std::array in = {{1, 2, 3, 4}}; + const int expected = 10; + auto value = exateppabm::util::reduce(in.begin(), in.end(), 0); + EXPECT_EQ(value, expected); +} + +/** + * Test the naive or std plus reduction a vector of unsigned integers + */ +TEST(TestUtil, reduce_vec_uint) { + std::vector in = {{1, 2, 3, 4}}; + const std::uint32_t expected = 10u; + auto value = exateppabm::util::reduce(in.begin(), in.end(), 0u); + EXPECT_EQ(value, expected); +} + +/** + * Test the naive or std plus reduction a vector of floats, using an arbitrary epsilon value which is good enough for the input values + */ +TEST(TestUtil, reduce_vec_double) { + std::vector in = {{1.1, 2.2, 3.3, 4.4}}; + const double expected = 11.0; + auto value = exateppabm::util::reduce(in.begin(), in.end(), 0.); + constexpr double EPSILON = 1.0e-6f; + EXPECT_NEAR(value, expected, EPSILON); +} + +/** + * Test the std plus reduction an array of integers + This is just to check the test case behaves too + */ +#if defined(EXATEPP_ABM_USE_STD_INCLUSIVE_SCAN) && EXATEPP_ABM_USE_STD_INCLUSIVE_SCAN +TEST(TestUtil, std_reduce_array_int) { + constexpr std::uint32_t ELEMENTS = 4; + std::array in = {{1, 2, 3, 4}}; + const int expected = 10; + auto value = std::reduce(in.begin(), in.end(), 0); + EXPECT_EQ(value, expected); +} +#else +TEST(TestUtil, DISABLED_std_reduce_array_int) { +} +#endif + + +/** + * Test the std plus reduction a vector of unsigned integers. + * This is just to check the test case behaves too + */ +#if defined(EXATEPP_ABM_USE_STD_INCLUSIVE_SCAN) && EXATEPP_ABM_USE_STD_INCLUSIVE_SCAN +TEST(TestUtil, std_reduce_vec_uint) { + std::vector in = {{1, 2, 3, 4}}; + const std::uint32_t expected = 10u; + auto value = std::reduce(in.begin(), in.end(), 0u); + EXPECT_EQ(value, expected); +} +#else +TEST(TestUtil, DISABLED_std_reduce_vec_uint) { +} +#endif + +/** + * Test the std plus reduction a vector of floats, using an arbitrary epsilon value which is good enough for the input values. + * This is just to check the test case behaves too + */ +#if defined(EXATEPP_ABM_USE_STD_INCLUSIVE_SCAN) && EXATEPP_ABM_USE_STD_INCLUSIVE_SCAN +TEST(TestUtil, std_reduce_vec_double) { + std::vector in = {{1.1, 2.2, 3.3, 4.4}}; + const double expected = 11.0; + auto value = std::reduce(in.begin(), in.end(), 0.); + constexpr double EPSILON = 1.0e-6f; + EXPECT_NEAR(value, expected, EPSILON); +} +#else +TEST(TestUtil, DISABLED_std_reduce_vec_double) { +} +#endif diff --git a/tools/plot-per-individual.py b/tools/plot-per-individual.py new file mode 100755 index 0000000..5c9dedb --- /dev/null +++ b/tools/plot-per-individual.py @@ -0,0 +1,131 @@ +#! /usr/bin/env python3 +import argparse +import pandas as pd +import seaborn as sns +import matplotlib.pyplot as plt +import pathlib + +def parse_cli(): + parser = argparse.ArgumentParser(description="Plotting script for per individual data") + parser.add_argument("csv", type=pathlib.Path, help="per individual CSV to plot") + parser.add_argument("-o", "--output", type=pathlib.Path, help="Path to output image location") + args = parser.parse_args() + return args + +def read_inputs(args): + df = pd.read_csv(args.csv) + + expected_columns = ["ID", "age_group", "occupation_network", "house_no", "infection_count"] + missing_columns = [c for c in expected_columns if c not in df.columns] + if len(missing_columns) > 0: + raise Exception(f"expected columns missing from csv: {missing_columns}") + + # Add string version of age group + AGE_GROUP_MAP = { + 0: "0-9", + 1: "10-19", + 2: "20-29", + 3: "30-39", + 4: "40-49", + 5: "50-59", + 6: "60-69", + 7: "70-79", + 8: "80+", + } + df["age_group_str"] = pd.Categorical(df["age_group"].map(AGE_GROUP_MAP), list(AGE_GROUP_MAP.values())) + + # Add string version of occupation network group + OCCUPATION_NETWORK_MAP = { + 0: "0-9", + 1: "10-19", + 2: "20-69", + 3: "70-79", + 4: "80+" + } + df["occupation_network_str"] = pd.Categorical(df["occupation_network"].map(OCCUPATION_NETWORK_MAP), list(OCCUPATION_NETWORK_MAP.values())) + return df + +def save_or_show(args): + # Either save the current plot to disk, or show it. + if (args.output): + args.output.parent.mkdir(exist_ok=True) + plt.savefig(args.output) + else: + plt.show() + + +def plot_age_histogram(args, df, ax): + ax.set_title("Population per age demographic") + sns.histplot(ax=ax, data=df, x="age_group_str", discrete=True) + ax.tick_params(axis='x', rotation=90) + +def plot_age_histogram_per_workplace(args, df, ax): + ax.set_title("Age Demographic per workplace") + sns.histplot(ax=ax, data=df, x="age_group_str", hue="occupation_network_str", discrete=True, element="step", palette="Dark2") + sns.move_legend(ax, "upper left", bbox_to_anchor=(1, 1)) + ax.tick_params(axis='x', rotation=90) + + + +def plot_household_size_household_age(args, df, ax): + ax.set_title("Household size against household age") + + # Group by "house_no" and calculate mean "age_group" and count + df_mean_age_per_house = df.groupby('house_no').agg({'age_group': ['mean']}).reset_index() + df_mean_age_per_house.columns = ["house_no", "mean_age_group"] + count_per_house = df.groupby('house_no').size().reset_index(name="house_size") + combined_df = df_mean_age_per_house + combined_df["house_size"] = count_per_house["house_size"] + + # sns.scatterplot(ax=ax, data=combined_df, x="house_size", y="mean_age_group") + # sns.violinplot(ax=ax, data=combined_df, x="house_size", y="mean_age_group", inner="quart") + sns.boxplot(ax=ax, data=combined_df, x="house_size", y="mean_age_group") + sns.stripplot(ax=ax, data=combined_df, x="house_size", y="mean_age_group") + # sns.swarmplot(ax=ax, data=combined_df, x="house_size", y="mean_age_group") + + + + +def plot_per_house_infection_count(args, df, ax): + ax.set_title("Total infection count per house") + + # Group by "house_no" and compute the sum of infections per hosue + df_grouped = df.groupby('house_no').agg({'infection_count': ['sum']}).reset_index() + df_grouped.columns = ["house_no", "sum_infection_count"] + count_per_house = df.groupby('house_no').size().reset_index(name="house_size") + df_grouped["style"] = 0 + + markers = ["x"] + sizes = [0.5] + sns.scatterplot(ax=ax, data=df_grouped, x="house_no", y="sum_infection_count", style="style", markers=markers, sizes=sizes, legend=None) + + +def plots(args, df): + sns.set_palette("muted") + sns.set_context("talk") + sns.set_style("darkgrid") + + # Get the current palette + palette = sns.color_palette(palette=None) + + fig, axes = plt.subplots(2, 2, figsize=(16, 9), sharex=False, sharey=False, layout="constrained") + + # plot different things on each axis + plot_age_histogram(args, df, axes[0, 0]) + plot_age_histogram_per_workplace(args, df, axes[0, 1]) + plot_household_size_household_age(args, df, axes[1, 0]) + plot_per_house_infection_count(args, df, axes[1, 1]) + + # Save to disk or show on screen + save_or_show(args) + +def main(): + # CLI parsing + args = parse_cli() + # Read the single input file + df = read_inputs(args) + # Do some plots + plots(args, df) + +if __name__ == "__main__": + main() \ No newline at end of file