From 0bdf294b2ac57eaed7ba2dbebbe7b73b32a59f8f Mon Sep 17 00:00:00 2001 From: Peter Heywood Date: Wed, 18 Dec 2024 18:07:17 +0000 Subject: [PATCH] WIP: Partial transmission file implementation --- src/CMakeLists.txt | 7 +- src/exateppabm/cli.cpp | 2 + src/exateppabm/cli.h | 4 + src/exateppabm/disease/SEIR.cu | 42 +------ src/exateppabm/disease/SEIR.h | 133 ++++++++++++++++++++-- src/exateppabm/exatepp_abm.cu | 2 +- src/exateppabm/output.cu | 132 +++++++++++++++++++-- src/exateppabm/output.h | 3 +- src/exateppabm/output/TransmissionFile.cu | 57 ++++++++++ src/exateppabm/output/TransmissionFile.h | 132 +++++++++++++++++++++ src/exateppabm/person.cu | 7 ++ src/exateppabm/person.h | 14 ++- src/exateppabm/population.cu | 14 +++ tests/CMakeLists.txt | 24 +++- 14 files changed, 507 insertions(+), 66 deletions(-) create mode 100644 src/exateppabm/output/TransmissionFile.cu create mode 100644 src/exateppabm/output/TransmissionFile.h diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 999785a..b31d6ab 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -25,9 +25,10 @@ SET(LIBRARY_SRC ${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/PerformanceFile.h ${CMAKE_CURRENT_SOURCE_DIR}/exateppabm/output/TimeSeriesFile.h + ${CMAKE_CURRENT_SOURCE_DIR}/exateppabm/output/TransmissionFile.h ${CMAKE_CURRENT_SOURCE_DIR}/exateppabm/person.h ${CMAKE_CURRENT_SOURCE_DIR}/exateppabm/population.h ${CMAKE_CURRENT_SOURCE_DIR}/exateppabm/random_interactions.h @@ -43,10 +44,10 @@ SET(LIBRARY_SRC ${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/PerformanceFile.cpp ${CMAKE_CURRENT_SOURCE_DIR}/exateppabm/output/TimeSeriesFile.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/exateppabm/output/TransmissionFile.cu ${CMAKE_CURRENT_SOURCE_DIR}/exateppabm/person.cu ${CMAKE_CURRENT_SOURCE_DIR}/exateppabm/population.cu ${CMAKE_CURRENT_SOURCE_DIR}/exateppabm/random_interactions.cu diff --git a/src/exateppabm/cli.cpp b/src/exateppabm/cli.cpp index 279f285..1305aab 100644 --- a/src/exateppabm/cli.cpp +++ b/src/exateppabm/cli.cpp @@ -15,6 +15,7 @@ void setup(CLI::App& app, std::shared_ptr params) { 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"); + app.add_flag("--transmission-file", params->transmissionFile, "Enable the transmission file. This has significant performance implications"); } void print(const exateppabm::cli::params params) { @@ -25,6 +26,7 @@ void print(const exateppabm::cli::params params) { fmt::print(" inputParamLine = {}\n", params.inputParamLine); fmt::print(" outputDir = {}\n", params.outputDir); fmt::print(" individualFile = {}\n", params.individualFile); + fmt::print(" transmissionFile = {}\n", params.transmissionFile); fmt::print("}}\n"); } diff --git a/src/exateppabm/cli.h b/src/exateppabm/cli.h index c9eb3c9..8c13096 100644 --- a/src/exateppabm/cli.h +++ b/src/exateppabm/cli.h @@ -41,6 +41,10 @@ struct params { * bool indicating if the individual file should be written */ bool individualFile = false; + /** + * bool indicating if the transmission file should be written + */ + bool transmissionFile = false; }; /** diff --git a/src/exateppabm/disease/SEIR.cu b/src/exateppabm/disease/SEIR.cu index 12c2cca..83043aa 100644 --- a/src/exateppabm/disease/SEIR.cu +++ b/src/exateppabm/disease/SEIR.cu @@ -17,9 +17,6 @@ FLAMEGPU_AGENT_FUNCTION(progressDisease, flamegpu::MessageNone, flamegpu::Messag // Get the current timestep / day std::uint32_t today = FLAMEGPU->getStepCounter(); - // Get a handle to the total_infected_per_demographic macro env property - auto totalInfectedPerDemographic = FLAMEGPU->environment.getMacroProperty("total_infected_per_demographic"); - // Get the current agents infection status auto infectionState = FLAMEGPU->getVariable(person::v::INFECTION_STATE); // Get the time they last changed state @@ -28,52 +25,23 @@ FLAMEGPU_AGENT_FUNCTION(progressDisease, flamegpu::MessageNone, flamegpu::Messag float stateDuration = FLAMEGPU->getVariable(person::v::INFECTION_STATE_DURATION); // Ready to change state if today is past the next scheduled state change bool readyToChange = today >= dayOfLastStateChange + std::ceil(stateDuration); - // Get the agent's demographic - auto demographic_idx = FLAMEGPU->getVariable(person::v::AGE_DEMOGRAPHIC); - // For each different initial state, change if required and compute the next state's duration. - if (infectionState == disease::SEIR::InfectionState::Susceptible) { - // no op - } else if (infectionState == disease::SEIR::InfectionState::Exposed) { + if (infectionState == disease::SEIR::InfectionState::Exposed) { // Exposed to Infected, if enough time has passed if (readyToChange) { - // Update the state - infectionState = disease::SEIR::InfectionState::Infected; - // Update the day - dayOfLastStateChange = today; - // Compute how long the next state will last - float mean = FLAMEGPU->environment.getProperty("mean_time_to_recovered"); - float sd = FLAMEGPU->environment.getProperty("sd_time_to_recovered"); - stateDuration = (FLAMEGPU->random.normal() * sd) + mean; - // Atomically update the number of infected individuals for the current individuals demographics when they transition into the infection state - totalInfectedPerDemographic[demographic_idx]++; + exposedToInfected(FLAMEGPU, infectionState); } } else if (infectionState == disease::SEIR::InfectionState::Infected) { // Infected to Recovered if enough time has passed if (readyToChange) { - // Update the state - infectionState = disease::SEIR::InfectionState::Recovered; - // Update the day - dayOfLastStateChange = today; - // Compute how long the next state will last - float mean = FLAMEGPU->environment.getProperty("mean_time_to_susceptible"); - float sd = FLAMEGPU->environment.getProperty("sd_time_to_susceptible"); - stateDuration = (FLAMEGPU->random.normal() * sd) + mean; + infectedToRecovered(FLAMEGPU, infectionState); } } else if (infectionState == disease::SEIR::InfectionState::Recovered) { // Recovered to Susceptible, if enough time has passed. if (readyToChange) { - infectionState = disease::SEIR::Susceptible; - dayOfLastStateChange = today; - stateDuration = 0; // susceptible doesn't have a fixed duration + recoveredToSusceptible(FLAMEGPU, infectionState); } } - - // Update global agent variables from local (in register) values. - FLAMEGPU->setVariable(person::v::INFECTION_STATE, infectionState); - FLAMEGPU->setVariable(person::v::INFECTION_STATE_CHANGE_DAY, dayOfLastStateChange); - FLAMEGPU->setVariable(person::v::INFECTION_STATE_DURATION, stateDuration); - return flamegpu::ALIVE; } @@ -86,7 +54,7 @@ void define(flamegpu::ModelDescription& model, const exateppabm::input::config& env.newMacroProperty("total_infected_per_demographic"); // Add a number of model parameters to the environment, initialised with the value from the configuration file - // @todo - not all of these feel right here / add cosntexpr strings somewhere. + // @todo - not all of these feel right here / add constexpr strings somewhere. env.newProperty("mean_time_to_infected", params.mean_time_to_infected); env.newProperty("sd_time_to_infected", params.sd_time_to_infected); env.newProperty("mean_time_to_recovered", params.mean_time_to_recovered); diff --git a/src/exateppabm/disease/SEIR.h b/src/exateppabm/disease/SEIR.h index 0693475..89860a3 100644 --- a/src/exateppabm/disease/SEIR.h +++ b/src/exateppabm/disease/SEIR.h @@ -4,13 +4,15 @@ #include "flamegpu/flamegpu.h" #include "exateppabm/input.h" #include "exateppabm/person.h" +#include "exateppabm/demographics.h" + namespace exateppabm { namespace disease { /** * Methods and variables related to SEIR modelling, with no birth or death included - * + * * This is intended to demonstrate the type of model which can be implemented and how, with more complex disease modelling following a similar pattern. */ namespace SEIR { @@ -23,9 +25,9 @@ namespace SEIR { typedef std::uint32_t InfectionStateUnderlyingType; /** - * Enum representing the different states of the model. + * Enum representing the different states of the model. * Unfortunately unable to use a enum class for strong typing without excessive casting being required when passing to/from templated FLAMEGPU 2 methods. - * + * * @note - the static_casting required might be worthwhile if multiple infection models are implemented, as non-class enums are in the global scope without any prefixing. * @note - Possibly just define __host__ __device__ methods to/from_uint8() rather than having the cast? Not sure is available on device for a host device to_underlying pre c++23 * @note - not just named "state" to reduce confusion with FLAME GPU 2 agent states. @@ -40,7 +42,7 @@ enum InfectionState : InfectionStateUnderlyingType { /** * Add required parts to a FLAME GPU 2 model for this implementation. - * + * * @note - this is likely to be refactored several times * @param model the FLAME GPU Model description to mutate * @param parameters the model parameters @@ -56,9 +58,9 @@ void appendLayers(flamegpu::ModelDescription& model); /** * Device utility function to get an individuals current infection status from global agent memory - * + * * Templated for due to the templated DeviceAPI object - * * + * * * @param FLAMEGPU flame gpu device API object * @return individuals current infection state */ @@ -70,11 +72,11 @@ FLAMEGPU_DEVICE_FUNCTION disease::SEIR::InfectionStateUnderlyingType getCurrentI /** * Device utility function for when an individual is exposed, moving from susceptible to exposed - * + * * Templated for due to the templated DeviceAPI object - * + * * @param FLAMEGPU flamegpu device api object - * @param current infection status for the individual to be mutated in-place + * @param infectionStatus current infection status for the individual to be mutated in-place */ template FLAMEGPU_DEVICE_FUNCTION void susceptibleToExposed(flamegpu::DeviceAPI* FLAMEGPU, disease::SEIR::InfectionStateUnderlyingType& infectionStatus) { @@ -92,6 +94,119 @@ FLAMEGPU_DEVICE_FUNCTION void susceptibleToExposed(flamegpu::DeviceAPItemplate setVariable(person::v::TIME_EXPOSED, FLAMEGPU->getStepCounter()); +} + + +/** + * Function for updating agent data when moved from the Exposed to Infected state + * + * Templated for due to the templated DeviceAPI object + * + * @param FLAMEGPU flamegpu device api object + * @param infectionStatus current infection status for the individual to be mutated in-place + */ +template +FLAMEGPU_DEVICE_FUNCTION void exposedToInfected(flamegpu::DeviceAPI* FLAMEGPU, disease::SEIR::InfectionStateUnderlyingType& infectionStatus) { + std::uint32_t today = FLAMEGPU->getStepCounter(); + // Get a handle to the total_infected_per_demographic macro env property + auto totalInfectedPerDemographic = FLAMEGPU->environment.template getMacroProperty("total_infected_per_demographic"); + + // Get the agent's demographic + auto demographic_idx = FLAMEGPU->template getVariable(person::v::AGE_DEMOGRAPHIC); + + // Update the state + infectionStatus = disease::SEIR::InfectionState::Infected; + FLAMEGPU->template setVariable(person::v::INFECTION_STATE, infectionStatus); + + // Update the day + FLAMEGPU->template setVariable(person::v::INFECTION_STATE_CHANGE_DAY, today); + + // Compute how long the next state will last + float mean = FLAMEGPU->environment.template getProperty("mean_time_to_recovered"); + float sd = FLAMEGPU->environment.template getProperty("sd_time_to_recovered"); + float stateDuration = (FLAMEGPU->random.template normal() * sd) + mean; + FLAMEGPU->template setVariable(person::v::INFECTION_STATE_DURATION, stateDuration); + + // Atomically update the number of infected individuals for the current individual's demographics when they transition into the infection state + totalInfectedPerDemographic[demographic_idx]++; + + // Update the agent variable tracking when the "current" infection moved from exposed to infected + FLAMEGPU->template setVariable(person::v::TIME_INFECTED, today); +} + +/** + * Function for updating agent data when moved from the Infected to Recovered state + * + * Templated for due to the templated DeviceAPI object + * + * @param FLAMEGPU flamegpu device api object + * @param infectionStatus current infection status for the individual to be mutated in-place + */ +template +FLAMEGPU_DEVICE_FUNCTION void infectedToRecovered(flamegpu::DeviceAPI* FLAMEGPU, disease::SEIR::InfectionStateUnderlyingType& infectionStatus) { + std::uint32_t today = FLAMEGPU->getStepCounter(); + + // Update the state + infectionStatus = disease::SEIR::InfectionState::Recovered; + FLAMEGPU->template setVariable(person::v::INFECTION_STATE, infectionStatus); + + // Update the day + FLAMEGPU->template setVariable(person::v::INFECTION_STATE_CHANGE_DAY, today); + + // Compute how long the next state will last + float mean = FLAMEGPU->environment.template getProperty("mean_time_to_susceptible"); + float sd = FLAMEGPU->environment.template getProperty("sd_time_to_susceptible"); + float stateDuration = (FLAMEGPU->random.template normal() * sd) + mean; + FLAMEGPU->template setVariable(person::v::INFECTION_STATE_DURATION, stateDuration); + + // Update the agent variable tracking when the "current" infection moved from infected to recovered + FLAMEGPU->template setVariable(person::v::TIME_RECOVERED, today); +} + +/** + * Function for updating agent data when moved from the Recovered to Susceptible state + * + * Templated for due to the templated DeviceAPI object + * + * @param FLAMEGPU flamegpu device api object + * @param infectionStatus current infection status for the individual to be mutated in-place + + */ +template +FLAMEGPU_DEVICE_FUNCTION void recoveredToSusceptible(flamegpu::DeviceAPI* FLAMEGPU, disease::SEIR::InfectionStateUnderlyingType& infectionStatus) { + std::uint32_t today = FLAMEGPU->getStepCounter(); + + // Update the state + infectionStatus = disease::SEIR::Susceptible; + FLAMEGPU->template setVariable(person::v::INFECTION_STATE, infectionStatus); + + // Update the day + FLAMEGPU->template setVariable(person::v::INFECTION_STATE_CHANGE_DAY, today); + float stateDuration = 0; // susceptible doesn't have a fixed duration + FLAMEGPU->template setVariable(person::v::INFECTION_STATE_DURATION, stateDuration); + + // Update the agent variable tracking when the "current" infection moved from recovered to Susceptible + FLAMEGPU->template setVariable(person::v::TIME_SUSCEPTIBLE, today); +} + +/** + * Function for resetting "current" infection timing data, for use after it has been recorded + * + * Templated for due to the templated DeviceAPI object + * + * @param FLAMEGPU flamegpu device api object + */ +template +FLAMEGPU_DEVICE_FUNCTION void resetSEIRStateTimes(flamegpu::DeviceAPI* FLAMEGPU) { + // Update the agent variable tracking when the "current" infection moved from Recovered to Susceptible, and the data has been logged. + const std::uint32_t invalidTime = static_cast(-1); + // Don't reset susceptible + FLAMEGPU->template setVariable(person::v::TIME_EXPOSED, invalidTime); + FLAMEGPU->template setVariable(person::v::TIME_INFECTED, invalidTime); + FLAMEGPU->template setVariable(person::v::TIME_RECOVERED, invalidTime); } } // namespace SEIR diff --git a/src/exateppabm/exatepp_abm.cu b/src/exateppabm/exatepp_abm.cu index 8a9c9d3..380a780 100644 --- a/src/exateppabm/exatepp_abm.cu +++ b/src/exateppabm/exatepp_abm.cu @@ -95,7 +95,7 @@ int entrypoint(int argc, char* argv[]) { 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, cli_params->individualFile); + exateppabm::output::define(model, cli_params->outputDir, cli_params->individualFile, cli_params->transmissionFile); // 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. diff --git a/src/exateppabm/output.cu b/src/exateppabm/output.cu index 95382a8..81e9e07 100644 --- a/src/exateppabm/output.cu +++ b/src/exateppabm/output.cu @@ -3,15 +3,20 @@ #include #include #include +#include #include #include "exateppabm/output/OutputFile.h" #include "exateppabm/output/TimeSeriesFile.h" #include "exateppabm/output/PerIndividualFile.h" +#include "exateppabm/output/TransmissionFile.h" #include "exateppabm/person.h" #include "exateppabm/population.h" #include "exateppabm/demographics.h" #include "exateppabm/disease.h" +#include "exateppabm/household.h" +#include "exateppabm/workplace.h" + namespace exateppabm { @@ -30,13 +35,16 @@ std::unique_ptr _timeSeriesFile = nullptr; // Object representing the per individual file std::unique_ptr _perIndividualFile = nullptr; +// Object representing the transmission file +std::unique_ptr _transmissionFile = nullptr; + } // namespace /** * FLAME GPU init function to prepare for time series data capture throughout the simulation * @note - this is not sustainable for simulations with long step counts due to single output to disk, but is OK for ~365 entries. This would need refactoring to emit partial files every N iterations if finer grained data is recorded (or per agent data when that is implemented.) */ -FLAMEGPU_INIT_FUNCTION(output_init) { +FLAMEGPU_INIT_FUNCTION(output_timeseries_init) { // (re) initialise the time series file data structure with preallocated room for the number of steps. _timeSeriesFile->resetObservations(FLAMEGPU->getSimulationConfig().steps); // Set the initial number of infected individuals per age demographic. @todo. Possibly move generation into an init method instead and do it their instead? @@ -51,7 +59,7 @@ FLAMEGPU_INIT_FUNCTION(output_init) { * FLAME GPU step function, which executes at the end of each time step. * Collect relevant data from agents, and store in memory for output do disk at exit. */ -FLAMEGPU_STEP_FUNCTION(output_step) { +FLAMEGPU_STEP_FUNCTION(output_timeseries_step) { // get the current iteration number (0 indexed) auto step = FLAMEGPU->getStepCounter(); // Get an object in which to store time series data @@ -94,7 +102,7 @@ FLAMEGPU_STEP_FUNCTION(output_step) { _timeSeriesFile->appendObservations(observations); } -FLAMEGPU_EXIT_FUNCTION(output_exit) { +FLAMEGPU_EXIT_FUNCTION(output_timeseries_exit) { // Write the time series data to disk // Open the file handle _timeSeriesFile->open(); @@ -115,13 +123,12 @@ FLAMEGPU_EXIT_FUNCTION(output_exit_per_individual) { 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.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(); @@ -131,8 +138,58 @@ FLAMEGPU_EXIT_FUNCTION(output_exit_per_individual) { _perIndividualFile->close(); } +/** + * Exit function to collect transmission file data for agents who are currently in a non-susceptible state, for their current infection. Some values will be -1 in this case. + */ +FLAMEGPU_EXIT_FUNCTION(transmissionFileExit) { + fmt::print("@todo - collect transmissionData on Exit for all individuals not susceptible\n"); + // 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) { + // If the person is not susceptible, they are in an active infection which has not yet been logged to disk. + auto currentInfectionStatus = person.getVariable(person::v::INFECTION_STATE); + if (currentInfectionStatus != disease::SEIR::Susceptible) { + exateppabm::output::TransmissionFile::Event data = {}; + // Store relevant data about the current individual (recipient) + data.id_recipient = static_cast(person.getVariable(person::v::ID)); + data.age_group_recipient = person.getVariable(person::v::AGE_DEMOGRAPHIC); + data.house_no_recipient = person.getVariable(person::v::HOUSEHOLD_IDX); + data.occupation_network_recipient = person.getVariable(person::v::WORKPLACE_IDX); + + // Store which interaction network the event occured in. @todo - enum. + data.transmission_event_network = person.getVariable(person::v::TF_EVENT_NETWORK); + + // Store data about the source of the infection + data.id_source = static_cast(person.getVariable(person::v::TF_SOURCE_ID)); + // @todo - get the source workplace, household, age demo etc from disk, or do this once for all records on exit? + // data.age_group_source = something[data.id_source - 1]; + // data.house_no_source = something[data.id_source - 1]; + // data.occupation_network_source = something[data.id_source - 1]; + data.time_exposed_source = person.getVariable(person::v::TF_SOURCE_TIME_EXPOSED); + + // Log the current value for the duration in each disease state + data.time_exposed = person.getVariable(person::v::TIME_EXPOSED); + data.time_infected = person.getVariable(person::v::TIME_INFECTED); + data.time_recovered = person.getVariable(person::v::TIME_RECOVERED); + data.time_susceptible = person.getVariable(person::v::TIME_SUSCEPTIBLE); + + // Append this individuals final infection event data to disk + _transmissionFile->append(data); + } + } + // Write the transmission data to disk + // Open the file handle + _transmissionFile->open(); + // Write data to the opened file + _transmissionFile->write(); + // Close the file handle + _transmissionFile->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, const bool individualFile) { +void define(flamegpu::ModelDescription& model, const std::filesystem::path outputDirectory, const bool individualFile, const bool transmissionFile) { // Store the output directory for access in the FLAME GPU exit function (and init?) // This will want refactoring for ensembles _outputDirectory = outputDirectory; @@ -140,17 +197,74 @@ void define(flamegpu::ModelDescription& model, const std::filesystem::path outpu _timeSeriesFile = std::make_unique(_outputDirectory); // Add the init function to the model - model.addInitFunction(output_init); + model.addInitFunction(output_timeseries_init); // Add the step function to the model - model.addStepFunction(output_step); + model.addStepFunction(output_timeseries_step); // Add the exit function to the model - model.addExitFunction(output_exit); + model.addExitFunction(output_timeseries_exit); // optionally prepare for per individual file output if (individualFile) { _perIndividualFile = std::make_unique(_outputDirectory); model.addExitFunction(output_exit_per_individual); } + + // Optionally prepare for the transmission event file, which includes mutation of the model definition, or just set variables used to disable it in device code + flamegpu::EnvironmentDescription env = model.Environment(); + if (!transmissionFile) { + // Define an environment property marking this feature as disabled + env.newProperty("transmissionFileEnabled", false, true); + } else { + // Initialise the file-scoped transmissionFile object + _transmissionFile = std::make_unique(_outputDirectory); + + // Define an environment property marking this feature as enabled + env.newProperty("transmissionFileEnabled", true, true); + + // @todo - might be nicer to move some of this to person.cu, undecided. + // Get a handle to the person agent type + flamegpu::AgentDescription agent = model.Agent(person::NAME); + + // Add agent variables for the person agent to store information for the transmission File. These are not defined otherwise to reduce memory cost. + + // The network in which the transmission event occurred. @todo enum. 0 is home, 1 is work, 2 is random + agent.newVariable(person::v::TF_EVENT_NETWORK, 0); + + // The agent id for the source of infection + agent.newVariable(person::v::TF_SOURCE_ID, 0); + + // The time at which the source of infection was exposed. + agent.newVariable(person::v::TF_SOURCE_TIME_EXPOSED, std::numeric_limits::max()); + + // Add the exit function to the model, which adds data for all non-suceptible individuals for their current infection event + model.addExitFunction(transmissionFileExit); + } +} + +void appendTransmissionFileLayers(flamegpu::ModelDescription& model) { + // If the transmission file is enabled, add agent functions/conditions to control flow, which allows efficient collect + if (_transmissionFile != nullptr) { + // Move relevant agents to the state for transmission file generation + { + auto layer = model.newLayer(); + layer.addAgentFunction(person::NAME, "transmissionFileToState"); + } + // Host layer function which populates the transmission file data structure with info from re-susceptible individuals + { + auto layer = model.newLayer(); + layer.addAgentFunction(person::NAME, "transmissionFileRecordCompleted"); + } + // Move relevant agents from transmission file back to the default state + { + auto layer = model.newLayer(); + layer.addAgentFunction(person::NAME, "transmissionFileToState"); + } + // Sort agents back into their original location, for consistency with the same simulation when this file is not enabled + { + auto layer = model.newLayer(); + layer.addAgentFunction(person::NAME, "transmissionFileSortByID"); + } + } } } // namespace output diff --git a/src/exateppabm/output.h b/src/exateppabm/output.h index d75ac1f..e8f36aa 100644 --- a/src/exateppabm/output.h +++ b/src/exateppabm/output.h @@ -24,8 +24,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 + * @param transmissionFile if the per infection transmission file should be written or not */ -void define(flamegpu::ModelDescription& model, std::filesystem::path outputDirectory, bool individualFile); +void define(flamegpu::ModelDescription& model, std::filesystem::path outputDirectory, bool individualFile, bool transmissionFile); } // namespace output } // namespace exateppabm diff --git a/src/exateppabm/output/TransmissionFile.cu b/src/exateppabm/output/TransmissionFile.cu new file mode 100644 index 0000000..ec50732 --- /dev/null +++ b/src/exateppabm/output/TransmissionFile.cu @@ -0,0 +1,57 @@ +#include "TransmissionFile.h" + +#include + +#include + +namespace exateppabm { +namespace output { + +TransmissionFile::TransmissionFile(std::filesystem::path directory) : OutputFile(directory / TransmissionFile::DEFAULT_FILENAME) { } + +TransmissionFile::~TransmissionFile() { } + +void TransmissionFile::reset(const std::uint32_t n_total) { + // Reset the observations vector, with a reserved initial capacity + this->_events = std::vector(); + this->_events.reserve(n_total); +} + +void TransmissionFile::append(const Event event) { + // @todo = ensure _events has been initialised? + this->_events.push_back(event); +} + +bool TransmissionFile::write() { + if (!this->_handle) { + this->open(); + } + + // Print to the file handle + fmt::print(_handle, "id_recipient,age_group_recipient,house_no_recipient,occupation_network_recipient,transmission_event_network,time,id_source,age_group_source,house_no_source,occupation_network_source,time_infected_source,time_exposed,time_infected,time_recovered,time_susceptible\n"); + for (const auto& event : this->_events) { + fmt::print( + _handle, + "{},{},{},{},{},{},{},{},{},{},{},{},{},{}\n", + event.id_recipient, + event.age_group_recipient, + event.house_no_recipient, + event.occupation_network_recipient, + event.transmission_event_network, + event.id_source, + event.age_group_source, + event.house_no_source, + event.occupation_network_source, + event.time_exposed_source, + event.time_exposed, + event.time_infected, + event.time_recovered, + event.time_susceptible); + } + + fmt::print("Transmission File written to {}\n", std::filesystem::absolute(this->_filepath).c_str()); + return true; +} + +} // namespace output +} // namespace exateppabm diff --git a/src/exateppabm/output/TransmissionFile.h b/src/exateppabm/output/TransmissionFile.h new file mode 100644 index 0000000..ed96be0 --- /dev/null +++ b/src/exateppabm/output/TransmissionFile.h @@ -0,0 +1,132 @@ +#pragma once + +#include "OutputFile.h" + +#include +#include +#include + +#include "flamegpu/flamegpu.h" +#include "exateppabm/demographics.h" +#include "exateppabm/person.h" + +namespace exateppabm { +namespace output { + +/** + * Class for representing the transmission file, a csv containing information about each transmission event and subsequent progress of disease + */ +class TransmissionFile : public OutputFile { + public: + // Forward decl nested struct + struct Event; + + /** + * Constructor setting the path for file output. + * @param directory parent directory for the file to be stored in (using the default filename) + */ + explicit TransmissionFile(std::filesystem::path directory); + + /** + * Dtor + */ + ~TransmissionFile(); + + /** + * 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 data related to an infection event. This should be recorded when the infection is over (i.e. returns to susceptible) or on simulation exit, to avoid multiple entries for the same infection + * + * @param event data for an infection, either from time of return to susceptible, or the end of the simulation. + */ + void append(Event event); + + /** + * 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 Event { + /** + * ID for the recipient + */ + std::uint32_t id_recipient = 0; + /** + * Age group of the recipient + */ + exateppabm::demographics::AgeUnderlyingType age_group_recipient = 0; + /** + * Household index of the recipient + */ + std::uint32_t house_no_recipient = 0; + /** + * Work network index of the recipient + */ + std::uint32_t occupation_network_recipient = 0; + /** + * The network in which the recipient was infected, essentially home, work or random + */ + std::uint8_t transmission_event_network = 0; + /** + * ID for the source of infection + */ + std::uint32_t id_source = 0; + /** + * Age group of the source + */ + exateppabm::demographics::AgeUnderlyingType age_group_source = 0; + /** + * Household index of the source + */ + std::uint32_t house_no_source = 0; + /** + * Work network index of the source + */ + std::uint32_t occupation_network_source = 0; + /** + * When the source was exposed + */ + std::uint32_t time_exposed_source = std::numeric_limits::max(); + // @todo - migrate these into a structure in SEIR.h, if supporting multiple disease models + /** + * Time at which this infection event moved from susceptible to exposed, i.e. the time of transmission. + UINT_MAX if not reached. + */ + std::uint32_t time_exposed = std::numeric_limits::max(); + /** + * Time at which this infection event moved from exposed to infected. UINT_MAX if not reached + */ + std::uint32_t time_infected = std::numeric_limits::max(); + /** + * Time at which this infection event moved from infected to recovered. UINT_MAX if not reached + */ + std::uint32_t time_recovered = std::numeric_limits::max(); + /** + * Time at which this infection event moved from recovered to susceptible. UINT_MAX if not reached + */ + std::uint32_t time_susceptible = std::numeric_limits::max(); + }; + + private: + /** + * Default filename for output + * @todo - factor in the run index for ensembles? + */ + constexpr static char DEFAULT_FILENAME[] = "transmission_file.csv"; + /** + * Private member containing the observation data + */ + std::vector _events; +}; + +} // namespace output +} // namespace exateppabm diff --git a/src/exateppabm/person.cu b/src/exateppabm/person.cu index 03123ed..32448ca 100644 --- a/src/exateppabm/person.cu +++ b/src/exateppabm/person.cu @@ -2,6 +2,7 @@ #include +#include #include #include "exateppabm/disease/SEIR.h" @@ -40,6 +41,12 @@ 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); + // Store the time at which this agent most recently entered each disease state, using UINT_MAX as a not set value + agent.newVariable(person::v::TIME_SUSCEPTIBLE, 0); + agent.newVariable(person::v::TIME_EXPOSED, std::numeric_limits::max()); + agent.newVariable(person::v::TIME_INFECTED, std::numeric_limits::max()); + agent.newVariable(person::v::TIME_RECOVERED, std::numeric_limits::max()); + // Integer count for the number of times infected, defaults to 0 agent.newVariable(person::v::INFECTION_COUNT, 0u); diff --git a/src/exateppabm/person.h b/src/exateppabm/person.h index 65343ac..565c348 100644 --- a/src/exateppabm/person.h +++ b/src/exateppabm/person.h @@ -54,6 +54,10 @@ 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 TIME_SUSCEPTIBLE[] = "time_susceptible"; +DEVICE_CONSTEXPR_STRING constexpr char TIME_EXPOSED[] = "time_exposed"; +DEVICE_CONSTEXPR_STRING constexpr char TIME_INFECTED[] = "time_infected"; +DEVICE_CONSTEXPR_STRING constexpr char TIME_RECOVERED[] = "time_recovered"; 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"; @@ -62,8 +66,10 @@ DEVICE_CONSTEXPR_STRING constexpr char WORKPLACE_OUT_DEGREE[] = "workplace_out_d 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"; - - +// Optionally available agent variables related to transmission file. +DEVICE_CONSTEXPR_STRING constexpr char TF_EVENT_NETWORK[] = "tf_event_network"; +DEVICE_CONSTEXPR_STRING constexpr char TF_SOURCE_ID[] = "tf_source_id"; // id can be used to find workplace/household +DEVICE_CONSTEXPR_STRING constexpr char TF_SOURCE_TIME_EXPOSED[] = "tf_source_time_exposed"; } // namespace v /** @@ -84,9 +90,9 @@ void appendLayers(flamegpu::ModelDescription& model); /** * Device utility function for Increasing an individuals infection counter, for more legible agent code - * + * * Templated for due to the templated DeviceAPI object - * + * * @todo - consider moving to the disease namespace? */ template diff --git a/src/exateppabm/population.cu b/src/exateppabm/population.cu index 1c63c8d..fe92dd4 100644 --- a/src/exateppabm/population.cu +++ b/src/exateppabm/population.cu @@ -52,6 +52,8 @@ struct HostPerson { workplace::WorkplaceUnderlyingType workplaceIdx = 0; std::uint32_t workplaceOutDegree = 0; std::uint32_t randomInteractionTarget = 0; + std::uint32_t timeExposed = std::numeric_limits::max(); + std::uint32_t timeInfected = std::numeric_limits::max(); }; /** @@ -162,6 +164,14 @@ FLAMEGPU_INIT_FUNCTION(generatePopulation) { if (infectionStatus == disease::SEIR::Infected) { // person.setVariable(exateppabm::person::v::INFECTION_COUNT, 1u); hostPersonData[personIdx].infectionCount = 1u; + + // person.setVariable(exateppabm::person::v::TIME_EXPOSED, 0u); + hostPersonData[personIdx].timeExposed = 0u; + + // person.setVariable(exateppabm::person::v::TIME_INEFECTED, 0u); + hostPersonData[personIdx].timeInfected = 0u; + + // Increment the per-age demographic initial agent count. @todo refactor elsewhere? _infectedPerDemographic[age]++; } @@ -345,6 +355,10 @@ FLAMEGPU_INIT_FUNCTION(generatePopulation) { person.setVariable(person::v::WORKPLACE_OUT_DEGREE, hostPerson.workplaceOutDegree); person.setVariable(person::v::RANDOM_INTERACTION_COUNT_TARGET, hostPerson.randomInteractionTarget); + person.setVariable(person::v::TIME_EXPOSED, hostPerson.timeExposed); + person.setVariable(person::v::TIME_EXPOSED, hostPerson.timeInfected); + + // If this is a visualisation enabled build, set their x/y/z #if defined(FLAMEGPU_VISUALISATION) diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 511cff6..726e416 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -10,10 +10,30 @@ if(TARGET exatepp_abm) add_test( NAME integration.defaultArgs COMMAND exatepp_abm - WORKING_DIRECTORY ${CMAKE_BINARY_DIR} + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} ) set_tests_properties(integration.defaultArgs PROPERTIES LABELS "integration") set_tests_properties(integration.defaultArgs PROPERTIES WILL_FAIL FALSE) + + # Test that running the default simulation with optional individual file output enabled runs successfully. + # This does not check that the file is actually produced, just that the simulation completes successfully. + add_test( + NAME integration.individualFile + COMMAND exatepp_abm --individual-file -o + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} + ) + set_tests_properties(integration.individualFile PROPERTIES LABELS "integration") + set_tests_properties(integration.individualFile PROPERTIES WILL_FAIL FALSE) + + # Test that running the default simulation with optional transmission file output enabled runs successfully. + # This does not check that the file is actually produced, just that the simulation completes successfully. + add_test( + NAME integration.transmissionFile + COMMAND exatepp_abm --transmission-file + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} + ) + set_tests_properties(integration.transmissionFile PROPERTIES LABELS "integration") + set_tests_properties(integration.transmissionFile PROPERTIES WILL_FAIL FALSE) endif() ## ------ @@ -46,7 +66,7 @@ flamegpu_add_executable("${PROJECT_NAME}" "${TESTS_SRC}" "${FLAMEGPU_ROOT}" "${P target_include_directories("${PROJECT_NAME}" PRIVATE "${CMAKE_CURRENT_SOURCE_DIR}") # Link against GTest target_link_libraries("${PROJECT_NAME}" PRIVATE GTest::gtest) -# Link against the static library +# Link against the static library target_link_libraries("${PROJECT_NAME}" PRIVATE exatepp_abm_lib) # Workaround for incremental rebuilds on MSVC, where device link was not being performed. # https://github.com/FLAMEGPU/FLAMEGPU2/issues/483