diff --git a/CMakeLists.txt b/CMakeLists.txt index ed164a9..400c11d 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -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 f730a38..4005b8f 100644 --- a/data/inputs/sample.csv +++ b/data/inputs/sample.csv @@ -1,5 +1,7 @@ -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,daily_fraction_work -12,0,365,1024,8,4,4,4,4,4,4,3,2,1,442,492,208,165,51,27,0.05,3.5,1.0,7.0,1.0,28,5.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,1.0,1.0,0.02 -12,1,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,0.001 -12,2,365,16384,8,331,370,374,371,361,415,377,270,142,442,492,208,165,51,27,0.10,3.5,1.0,7.0,1.0,28,5.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.3,0.3,0.001 -12,3,1,54,8,331,370,374,371,361,415,377,270,142,0,0,0,0,0,1,0.10,3.5,1.0,7.0,1.0,28,5.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.3,0.3,0.001 +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,daily_fraction_work,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,0.001,2,2,100,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,0.001,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,0.001,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,0.001,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,0.001,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,0.001,2,2,4,4,3,3 diff --git a/readme.md b/readme.md index e1f9281..72b9802 100644 --- a/readme.md +++ b/readme.md @@ -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/) | diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index e5c9999..aa00676 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -91,6 +91,12 @@ target_compile_definitions(${LIBRARY_NAME} PUBLIC "$<$:CM 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/exatepp_abm.cu b/src/exateppabm/exatepp_abm.cu index 92dd7ad..5a439cc 100644 --- a/src/exateppabm/exatepp_abm.cu +++ b/src/exateppabm/exatepp_abm.cu @@ -100,6 +100,15 @@ int entrypoint(int argc, char* argv[]) { // Add disease progression exateppabm::disease::SEIR::appendLayers(model); + // Generate the population of agents + // prior to the simulation creation so any mutation of the model environment is applied. this is not ideal and will need adjusting for ensembles. + // @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); + if (personPopulation == nullptr) { + throw std::runtime_error("@todo - bad population generation function."); + } + // Construct the Simulation instance from the model. flamegpu::CUDASimulation simulation(model); @@ -111,6 +120,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. @@ -119,13 +133,7 @@ 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); - if (personPopulation == nullptr) { - throw std::runtime_error("@todo - bad population generation function."); - } + // add the population to this simulation instance simulation.setPopulationData(*personPopulation); perfFile.timers.preSimulate.stop(); diff --git a/src/exateppabm/input.cu b/src/exateppabm/input.cu index 954eae6..b67608e 100644 --- a/src/exateppabm/input.cu +++ b/src/exateppabm/input.cu @@ -217,6 +217,24 @@ std::shared_ptr read(const std::filesystem::path p, c 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->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"); @@ -275,6 +293,12 @@ void print(exateppabm::input::config config) { fmt::print(" relative_transmission_occupation = {}\n", config.relative_transmission_occupation); fmt::print(" relative_transmission_random = {}\n", config.relative_transmission_occupation); fmt::print(" daily_fraction_work = {}\n", config.daily_fraction_work); + 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 22a69c9..93b9f02 100644 --- a/src/exateppabm/input.h +++ b/src/exateppabm/input.h @@ -199,9 +199,40 @@ struct config { /** * Fraction of people in work network interacted with per day (by rng sampling) * - * @todo - this probably needs using differntly with more realistic networks + * @todo - this probably needs using differently with more realistic networks */ float daily_fraction_work = 0.5f; + + /** + * Mean number of random interactions for 0-19 year olds + * Arbitrary default value + */ + double mean_random_interactions_0_19 = 2u; + /** + * Standard deviation for the number of random interactions per day for 0-19 year olds + * Default value is arbitrary + */ + double sd_random_interactions_0_19 = 2u; + /** + * Mean number of random interactions for 20-69 year olds + * Arbitrary default value + */ + double mean_random_interactions_20_69 = 4u; + /** + * Standard deviation for the number of random interactions per day for 20-69 year olds + * Default value is arbitrary + */ + double sd_random_interactions_20_69 = 4u; + /** + * Mean number of random interactions for 70+ year olds + * Arbitrary default value + */ + double mean_random_interactions_70plus = 2u; + /** + * Standard deviation for the number of random interactions per day for 70+ year olds + * Default value is arbitrary + */ + double sd_random_interactions_70plus = 2u; }; /** diff --git a/src/exateppabm/person.cu b/src/exateppabm/person.cu index 90a306a..5cae40b 100644 --- a/src/exateppabm/person.cu +++ b/src/exateppabm/person.cu @@ -2,6 +2,8 @@ #include +#include + #include "exateppabm/disease/SEIR.h" #include "exateppabm/demographics.h" @@ -220,9 +222,8 @@ FLAMEGPU_AGENT_FUNCTION(interactWorkplace, flamegpu::MessageBucket, flamegpu::Me * In serial on the host, generate a vector of interactions * * @todo - more than one random interaction per day - * @todo - prevent self-interactions * @todo - prevent duplicate interactions - * @todo - more performance CPU implementation + * @todo - more performant CPU implementation * @todo - GPU implementation of this (stable matching submodel?) */ FLAMEGPU_HOST_FUNCTION(updateRandomDailyNetworkIndices) { @@ -230,32 +231,126 @@ FLAMEGPU_HOST_FUNCTION(updateRandomDailyNetworkIndices) { auto personAgent = FLAMEGPU->agent(exateppabm::person::NAME, exateppabm::person::states::DEFAULT); flamegpu::DeviceAgentVector population = personAgent.getPopulationData(); - // @temp - Interact with ID +/1 pair wise for now, as a temp fix. - if (FLAMEGPU->getStepCounter() == 0) { - for (auto person : population) { + // @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.getID(); - flamegpu::id_t other = flamegpu::ID_NOT_SET; - if (id % 2 == 0) { - other = id - 1; - } else { - other = id + 1; + 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"); } - - person.setVariable(person::v::RANDOM_INTERACTION_COUNT, 1u); - person.setVariable(person::v::RANDOM_INTERACTION_PARTNERS, other); // @todo - multiple. + 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); - // @todo implement the following - // Get the sum of the per-agent random interaction count, so we can reserve a large enough vector + // 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.getID() != aID || bAgent.getID() != bID) { + throw std::runtime_error("Agent ID does not match expected agent ID in updateRandomDailyNetworkIndices. @todo"); + } - // Build a vector of interactions, initialised to contain the id of each agent as many times as they require + // @todo - avoid self-interactions by looking ahead and swapping? + // @todo - avoid repeated interactions by looking ahead and swapping? - // Shuffle the vector of agent id's + // 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); - // Update agent data containing today's random interactions - // @todo - support more than one interaction per agent per day. + + // 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) works locally, but doesn't seem to update the full thing so do the whole array per agent to enforce the copy? + // @todo - This should not be required, possible bug in FLAME GPU / hole in the API which I'll check separately + for (auto person : population) { + person.setVariable(person::v::RANDOM_INTERACTION_PARTNERS, person.getVariable(person::v::RANDOM_INTERACTION_PARTNERS)); + } + + // @temp + /* + for (const auto &person : population) { + flamegpu::id_t id = person.getID(); + fmt::print("{}: {} [{}, {}]\n", + id, + person.getVariable(person::v::RANDOM_INTERACTION_COUNT), + person.getVariable(person::v::RANDOM_INTERACTION_PARTNERS, 0), + person.getVariable(person::v::RANDOM_INTERACTION_PARTNERS, 1) + ); + }*/ + + population.syncChanges(); +} + +// @temp +FLAMEGPU_HOST_FUNCTION(test) { + // @temp + /* + fmt::print("! test\n"); + // Get the current population of person agents + auto personAgent = FLAMEGPU->agent(exateppabm::person::NAME, exateppabm::person::states::DEFAULT); + flamegpu::DeviceAgentVector population = personAgent.getPopulationData(); + for (const auto &person : population) { + flamegpu::id_t id = person.getID(); + fmt::print("{}: {} [{}, {}]\n", + id, + person.getVariable(person::v::RANDOM_INTERACTION_COUNT), + person.getVariable(person::v::RANDOM_INTERACTION_PARTNERS, 0), + person.getVariable(person::v::RANDOM_INTERACTION_PARTNERS, 1) + ); + } + */ } /** @@ -313,7 +408,7 @@ FLAMEGPU_AGENT_FUNCTION(interactRandomDailyNetwork, flamegpu::MessageArray, flam 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); // @todo - make this an array variable and access via the index + 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 @@ -344,6 +439,11 @@ FLAMEGPU_AGENT_FUNCTION(interactRandomDailyNetwork, flamegpu::MessageArray, flam } } + // 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; } @@ -392,9 +492,16 @@ void define(flamegpu::ModelDescription& model, const exateppabm::input::config& agent.newVariable(person::v::WORKPLACE_IDX); agent.newVariable(person::v::WORKPLACE_SIZE); - // Random interaction network variables. @tood -refactor to separate location - agent.newVariable(person::v::RANDOM_INTERACTION_PARTNERS, flamegpu::ID_NOT_SET); + // 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 @@ -518,6 +625,10 @@ void appendLayers(flamegpu::ModelDescription& model) { auto layer = model.newLayer(); layer.addHostFunction(updateRandomDailyNetworkIndices); } + { + auto layer = model.newLayer(); + layer.addHostFunction(test); + } { auto layer = model.newLayer(); layer.addAgentFunction(person::NAME, "emitRandomDailyNetworkStatus"); diff --git a/src/exateppabm/person.h b/src/exateppabm/person.h index 223d91a..84ebad3 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 */ @@ -47,8 +57,10 @@ 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_SIZE[] = "workplace_size"; -DEVICE_CONSTEXPR_STRING constexpr char RANDOM_INTERACTION_COUNT[] = "random_interaction_count"; 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 @@ -91,9 +103,9 @@ void define(flamegpu::ModelDescription& model, const exateppabm::input::config& /** * 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 949ea24..c742a3e 100644 --- a/src/exateppabm/population.cu +++ b/src/exateppabm/population.cu @@ -3,6 +3,7 @@ #include #include +#include #include #include #include @@ -30,6 +31,9 @@ std::array _infectedPerDemographic = {}; std::unique_ptr generate(flamegpu::ModelDescription& model, const exateppabm::input::config params, const bool verbose) { fmt::print("@todo - validate params inputs when generated agents (pop size, initial infected count etc)\n"); + // Get a handle on the environment. + auto env = model.Environment(); + // @todo - assert that the requested initial population is non zero. auto pop = std::make_unique(model.Agent(exateppabm::person::NAME), params.n_total); @@ -98,6 +102,12 @@ std::unique_ptr generate(flamegpu::ModelDescription& mode std::array peoplePerWorkplace = {{0, 0, 0}}; // /-------- + // 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(); + // Populate agent data, by iterating households std::uint32_t personIdx = 0; for (std::uint32_t householdIdx = 0; householdIdx < households.size(); householdIdx++) { @@ -152,11 +162,59 @@ std::unique_ptr generate(flamegpu::ModelDescription& mode // Store the assigned network in the agent data structure person.setVariable(person::v::WORKPLACE_IDX, 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-binomeal 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 popualtion 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); + // Track the sum of target interaction counts + randomInteractionCountSum += randomInteractionTarget; + // Increment the person index ++personIdx; } } + if (verbose) { + fmt::print("Random Interactions: min={}, max={}, sum={}\n", randomInteractionMin, randomInteractionMax, randomInteractionCountSum); + } + + // 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); + } + // In a separate pass over the population, set the size of each workplace network per individual // @todo - refactor this once small world networks are implemented std::uint32_t idx = 0;