From 7f7c034fbb2cab51cda34e594fab64aa7056acdc Mon Sep 17 00:00:00 2001 From: Peter Heywood Date: Fri, 13 Dec 2024 12:15:05 +0000 Subject: [PATCH] Refectoring: Additional person refactoring + initial population refactoring --- src/exateppabm/household.cu | 149 +++++++++++- src/exateppabm/household.h | 64 ++++- src/exateppabm/person.cu | 4 +- src/exateppabm/person.h | 28 --- src/exateppabm/population.cu | 226 +++--------------- src/exateppabm/population.h | 99 +------- src/exateppabm/random_interactions.cu | 35 +-- src/exateppabm/workplace.cu | 141 +++++++---- src/exateppabm/workplace.h | 44 ++++ tests/CMakeLists.txt | 3 +- .../{test_population.cu => test_household.cu} | 148 ++---------- tests/src/exateppabm/test_workplace.cu | 128 ++++++++++ 12 files changed, 539 insertions(+), 530 deletions(-) rename tests/src/exateppabm/{test_population.cu => test_household.cu} (50%) create mode 100644 tests/src/exateppabm/test_workplace.cu diff --git a/src/exateppabm/household.cu b/src/exateppabm/household.cu index ac0e1d6..5cf07af 100644 --- a/src/exateppabm/household.cu +++ b/src/exateppabm/household.cu @@ -1,13 +1,25 @@ #include "exateppabm/household.h" +#include + #include "flamegpu/flamegpu.h" +#include "fmt/core.h" #include "exateppabm/demographics.h" #include "exateppabm/disease.h" #include "exateppabm/person.h" +#include "exateppabm/util.h" namespace exateppabm { namespace household { +/** + * Namespace containing string constants related to the households message list + */ +namespace message_household_status { +constexpr char _NAME[] = "household_status"; +__device__ constexpr char ID[] = "id"; +} // namespace message_household_status + /** * Agent function for person agents to emit their public information, i.e. infection status to their household */ @@ -17,7 +29,7 @@ FLAMEGPU_AGENT_FUNCTION(emitHouseholdStatus, flamegpu::MessageNone, flamegpu::Me 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)); + FLAMEGPU->message_out.setVariable(message_household_status::ID, FLAMEGPU->getVariable(person::v::ID)); // Household index // @todo - typedef or using statement for the household index type? @@ -71,7 +83,7 @@ FLAMEGPU_AGENT_FUNCTION(interactHousehold, flamegpu::MessageBucket, flamegpu::Me // 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(person::message::household_status::ID) != id) { + if (message.getVariable(message_household_status::ID) != id) { // Ignore messages from other households if (message.getVariable(person::v::HOUSEHOLD_IDX) == householdIdx) { // Check if the other agent is infected @@ -114,19 +126,17 @@ void define(flamegpu::ModelDescription& model, const exateppabm::input::config& // Get a handle to the existing person agent type, which should have already been defined. flamegpu::AgentDescription agent = model.Agent(person::NAME); - // Household network variables. @todo - refactor to a separate network location? + // Household network variables agent.newVariable(person::v::HOUSEHOLD_IDX); agent.newVariable(person::v::HOUSEHOLD_SIZE); // Message list containing a persons current status for households (id, location, infection status) - flamegpu::MessageBucket::Description householdStatusMessage = model.newMessage(person::message::household_status::_NAME); - + flamegpu::MessageBucket::Description householdStatusMessage = model.newMessage(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); + householdStatusMessage.newVariable(message_household_status::ID); // Add the household index householdStatusMessage.newVariable(person::v::HOUSEHOLD_IDX); // Add a variable for the agent's infections status @@ -136,14 +146,14 @@ void define(flamegpu::ModelDescription& model, const exateppabm::input::config& // emit current status to the household flamegpu::AgentFunctionDescription emitHouseholdStatusDesc = agent.newFunction("emitHouseholdStatus", emitHouseholdStatus); - emitHouseholdStatusDesc.setMessageOutput(person::message::household_status::_NAME); + emitHouseholdStatusDesc.setMessageOutput(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.setMessageInput(message_household_status::_NAME); interactHouseholdDesc.setInitialState(person::states::DEFAULT); interactHouseholdDesc.setEndState(person::states::DEFAULT); } @@ -159,5 +169,126 @@ void appendLayers(flamegpu::ModelDescription& model) { } } + + +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; +} + } // namespace household } // namespace exateppabm diff --git a/src/exateppabm/household.h b/src/exateppabm/household.h index 54dafba..b2de951 100644 --- a/src/exateppabm/household.h +++ b/src/exateppabm/household.h @@ -1,7 +1,10 @@ #pragma once -#include "flamegpu/flamegpu.h" +#include + +#include "flamegpu/flamegpu.h" #include "exateppabm/input.h" +#include "exateppabm/demographics.h" namespace exateppabm { namespace household { @@ -23,5 +26,64 @@ void define(flamegpu::ModelDescription& model, const exateppabm::input::config& */ void appendLayers(flamegpu::ModelDescription& model); + + +/** + * 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); + + } // namespace household } // namespace exateppabm diff --git a/src/exateppabm/person.cu b/src/exateppabm/person.cu index 97c6731..03123ed 100644 --- a/src/exateppabm/person.cu +++ b/src/exateppabm/person.cu @@ -51,13 +51,13 @@ void define(flamegpu::ModelDescription& model, const exateppabm::input::config& // Define household related agent variables, functions, etc household::define(model, params); - + // Define workplace related agent variables, functions, etc workplace::define(model, params); // Define daily random interaction related agent variables, functions, etc random_interactions::define(model, params); - + // Define visualisation specific behaviours visualisation::define(model, params); } diff --git a/src/exateppabm/person.h b/src/exateppabm/person.h index 3246757..61a9e06 100644 --- a/src/exateppabm/person.h +++ b/src/exateppabm/person.h @@ -66,34 +66,6 @@ DEVICE_CONSTEXPR_STRING constexpr char RANDOM_INTERACTION_COUNT_TARGET[] = "rand } // namespace v -/** - * Namespace containing person-message related constants - */ -namespace message { -/** - * Namespace containing variable name constants for variables in household related messages - */ -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 - /** * Define the agent type representing a person in the simulation, mutating the model description object. * @param model flamegpu2 model description object to mutate diff --git a/src/exateppabm/population.cu b/src/exateppabm/population.cu index 686b95f..1c63c8d 100644 --- a/src/exateppabm/population.cu +++ b/src/exateppabm/population.cu @@ -13,10 +13,12 @@ #include "exateppabm/demographics.h" #include "exateppabm/disease.h" #include "exateppabm/input.h" +#include "exateppabm/household.h" #include "exateppabm/network.h" #include "exateppabm/person.h" #include "exateppabm/util.h" #include "exateppabm/visualisation.h" +#include "exateppabm/workplace.h" namespace exateppabm { namespace population { @@ -33,23 +35,23 @@ exateppabm::input::config _params = {}; bool _verbose = false; // Store the ID of each individual for each workplace, to form the node labels of the small world network -std::array, WORKPLACE_COUNT> workplaceMembers = {{}}; +std::array, workplace::WORKPLACE_COUNT> workplaceMembers = {{}}; } // namespace - -// Workaround struct representing a person used during host generation +// Struct used to build per-person data during multi-pass population generation. +// This is required due to bugs within FLAME GPU 2.0.0rc.2, which have since been fixed in by https://github.com/FLAMEGPU/FLAMEGPU2/pull/1270, so can be resolved in the future. struct HostPerson { - flamegpu::id_t id; + flamegpu::id_t id = flamegpu::ID_NOT_SET; disease::SEIR::InfectionStateUnderlyingType infectionStatus; - float infectionStateDuration; + float infectionStateDuration = 1.f; 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::uint32_t householdIdx = 0; + std::uint8_t householdSize = 0; + workplace::WorkplaceUnderlyingType workplaceIdx = 0; + std::uint32_t workplaceOutDegree = 0; + std::uint32_t randomInteractionTarget = 0; }; /** @@ -95,7 +97,7 @@ FLAMEGPU_INIT_FUNCTION(generatePopulation) { std::shuffle(infected_vector.begin(), infected_vector.end(), rng); // Get the number of individuals per house and their age demographics - auto households = generateHouseholdStructures(_params, rng, _verbose); + auto households = household::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 @@ -111,9 +113,9 @@ FLAMEGPU_INIT_FUNCTION(generatePopulation) { _infectedPerDemographic = {{0, 0, 0, 0, 0, 0, 0, 0, 0}}; // 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); + auto p_adultPerWorkNetwork = workplace::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 = {{}}; + // std::array, workplace::WORKPLACE_COUNT> workplaceMembers = {{}}; // Counter for random interaction count. This will be 2x the number of interactions std::uint64_t randomInteractionCountSum = 0u; @@ -126,7 +128,7 @@ FLAMEGPU_INIT_FUNCTION(generatePopulation) { 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++) { + for (household::HouseholdSizeType householdMemberIdx = 0; householdMemberIdx < household.size; householdMemberIdx++) { // assert that the household structure is complete assert(household.size == household.agePerPerson.size()); @@ -165,12 +167,12 @@ FLAMEGPU_INIT_FUNCTION(generatePopulation) { } // Assign the workplace based on age band - WorkplaceUnderlyingType workplaceIdx = generateWorkplaceForIndividual(age, p_adultPerWorkNetwork, rng); + auto workplaceIdx = workplace::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); + // 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. @@ -228,7 +230,7 @@ FLAMEGPU_INIT_FUNCTION(generatePopulation) { } // Now individuals have been assigned to each workplace network, we can generate the small world networks for workplace interactions - std::array meanInteractionsPerNetwork = {{ + 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, @@ -239,8 +241,8 @@ FLAMEGPU_INIT_FUNCTION(generatePopulation) { 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) { + std::array workplaceNetworks = {}; + for (workplace::WorkplaceUnderlyingType widx = 0u; widx < workplace::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); @@ -307,8 +309,8 @@ FLAMEGPU_INIT_FUNCTION(generatePopulation) { // 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; + // workplace::WorkplaceUnderlyingType workplaceIdx = person.getVariable(person::v::WORKPLACE_IDX); + auto 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(); @@ -339,7 +341,7 @@ FLAMEGPU_INIT_FUNCTION(generatePopulation) { 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_IDX, hostPerson.workplaceIdx); person.setVariable(person::v::WORKPLACE_OUT_DEGREE, hostPerson.workplaceOutDegree); person.setVariable(person::v::RANDOM_INTERACTION_COUNT_TARGET, hostPerson.randomInteractionTarget); @@ -370,11 +372,11 @@ FLAMEGPU_INIT_FUNCTION(generatePopulation) { 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(" 00-09: {}, 20_69 {}\n", createdPerDemographic[0], workplaceMembers[workplace::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::WORKPLACE_SCHOOL_0_9].size() - createdPerDemographic[demographics::Age::AGE_10_19]); + fmt::print(" 20_69 {}\n", workplaceMembers[workplace::Workplace::WORKPLACE_ADULT].size()); + fmt::print(" 20_69 {}, 70_79 {}\n", workplaceMembers[workplace::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::WORKPLACE_80_PLUS].size() - - createdPerDemographic[demographics::Age::AGE_80], createdPerDemographic[demographics::Age::AGE_80]); fmt::print("}}\n"); } } @@ -391,175 +393,5 @@ std::array getPerDemographicInitialInfec 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 } // namespace exateppabm diff --git a/src/exateppabm/population.h b/src/exateppabm/population.h index c1034f8..ffa7aad 100644 --- a/src/exateppabm/population.h +++ b/src/exateppabm/population.h @@ -7,6 +7,7 @@ #include "exateppabm/demographics.h" #include "exateppabm/person.h" #include "exateppabm/input.h" +#include "exateppabm/workplace.h" namespace exateppabm { @@ -35,104 +36,6 @@ void define(flamegpu::ModelDescription& model, const exateppabm::input::config p */ 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 } // namespace exateppabm diff --git a/src/exateppabm/random_interactions.cu b/src/exateppabm/random_interactions.cu index 53a6604..bb0b9a4 100644 --- a/src/exateppabm/random_interactions.cu +++ b/src/exateppabm/random_interactions.cu @@ -1,5 +1,7 @@ #include "exateppabm/random_interactions.h" +#include + #include "flamegpu/flamegpu.h" #include "exateppabm/demographics.h" #include "exateppabm/disease.h" @@ -9,6 +11,14 @@ namespace exateppabm { namespace random_interactions { +/** + * Namespace containing string constants related to messages lists within random interactions + * Use of __device__ constexpr char [] requires CUDA >= 11.4 + */ +namespace message_random_network_status { + constexpr char _NAME[] = "random_network_status"; +} // namespace message_random_network_status + /** * Update the per-day random daily interaction network * @@ -127,13 +137,12 @@ FLAMEGPU_HOST_FUNCTION(updateRandomDailyNetworkIndices) { } /** - * 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. + * Agent function for person agents to emit their public information + * i.e. infection status, for random daily network colleagues. * - * @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)); + // output public properties to array message, keyed by the agent ID so no need to include the ID. FLAMEGPU->message_out.setVariable(person::v:: INFECTION_STATE, FLAMEGPU->getVariable(person::v::INFECTION_STATE)); @@ -224,8 +233,8 @@ FLAMEGPU_AGENT_FUNCTION(interactRandomDailyNetwork, flamegpu::MessageArray, flam void define(flamegpu::ModelDescription& model, const exateppabm::input::config& params) { // Define related model environment properties flamegpu::EnvironmentDescription env = model.Environment(); - - // Environmental property containing the relative transmission scale factor for random interactions + + // Environmental property containing the relative transmission scale factor for random interactions env.newProperty("relative_transmission_random", params.relative_transmission_random); // Add an environmental variable containing the sum of each agents target number of random interactions. @@ -238,7 +247,7 @@ void define(flamegpu::ModelDescription& model, const exateppabm::input::config& flamegpu::AgentDescription agent = model.Agent(person::NAME); // Define person agent variables related to the random interactions - + // The target number of random interactions for this individual every day agent.newVariable(person::v::RANDOM_INTERACTION_COUNT_TARGET, 0u); @@ -250,18 +259,17 @@ void define(flamegpu::ModelDescription& model, const exateppabm::input::config& // 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); + flamegpu::MessageArray::Description randomNetworkStatusMessage = model.newMessage(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); + // No need to add the agent ID to the message, it's the same as the array index // 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 host and agent functions + // Define host and agent functions // Host function to generate todays random interaction list. // @todo - a GPU implementation will be needed for model scaling. @@ -270,16 +278,15 @@ void define(flamegpu::ModelDescription& model, const exateppabm::input::config& // emit current status for random interactions flamegpu::AgentFunctionDescription emitRandomDailyNetworkStatusDesc = agent.newFunction("emitRandomDailyNetworkStatus", emitRandomDailyNetworkStatus); - emitRandomDailyNetworkStatusDesc.setMessageOutput(person::message::random_network_status::_NAME); + emitRandomDailyNetworkStatusDesc.setMessageOutput(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.setMessageInput(message_random_network_status::_NAME); interactRandomDailyNetworkDesc.setInitialState(person::states::DEFAULT); interactRandomDailyNetworkDesc.setEndState(person::states::DEFAULT); - } void appendLayers(flamegpu::ModelDescription& model) { diff --git a/src/exateppabm/workplace.cu b/src/exateppabm/workplace.cu index 35ab19d..d7d8ceb 100644 --- a/src/exateppabm/workplace.cu +++ b/src/exateppabm/workplace.cu @@ -1,33 +1,33 @@ -#include "exateppabm/random_interactions.h" +#include "exateppabm/workplace.h" #include "flamegpu/flamegpu.h" #include "exateppabm/demographics.h" #include "exateppabm/disease.h" #include "exateppabm/person.h" +#include "exateppabm/util.h" namespace exateppabm { namespace workplace { +/** + * Namespace containing string constants related to the workplace message list + */ +namespace message_workplace_status { +constexpr char _NAME[] = "workplace_status"; +} // namespace message_workplace_status + /** * Agent function for person agents to emit their public information, i.e. infection status, to their workplace colleagues */ 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(person::v::WORKPLACE_IDX); - FLAMEGPU->message_out.setVariable(person::v::WORKPLACE_IDX, workplaceIdx); + // Output public properties required for interactions to an array message, indexed by the agent's ID. FLAMEGPU->message_out.setVariable(person::v:: INFECTION_STATE, FLAMEGPU->getVariable(person::v::INFECTION_STATE)); FLAMEGPU->message_out.setVariable(person::v::AGE_DEMOGRAPHIC, FLAMEGPU->getVariable(person::v::AGE_DEMOGRAPHIC)); - // Set the message key, the house hold idx for bucket messaging @Todo - FLAMEGPU->message_out.setIndex(id); + // Set the message key to be the agent ID for array access + FLAMEGPU->message_out.setIndex(FLAMEGPU->getVariable(person::v::ID)); return flamegpu::ALIVE; } @@ -43,7 +43,7 @@ FLAMEGPU_AGENT_FUNCTION(interactWorkplace, flamegpu::MessageArray, flamegpu::Mes // Get my ID const flamegpu::id_t id = FLAMEGPU->getVariable(person::v::ID); - // Get my workplace degree, if 0 return. + // Get my workplace degree, if 0 there will be no interactions so do return. const std::uint32_t degree = FLAMEGPU->getVariable(person::v::WORKPLACE_OUT_DEGREE); if (degree == 0) { return flamegpu::ALIVE; @@ -58,7 +58,7 @@ FLAMEGPU_AGENT_FUNCTION(interactWorkplace, flamegpu::MessageArray, flamegpu::Mes p_s2e *= FLAMEGPU->environment.getProperty("relative_transmission_occupation"); // Get my workplace/network index - auto workplaceIdx = FLAMEGPU->getVariable(person::v::WORKPLACE_IDX); + // auto workplaceIdx = FLAMEGPU->getVariable(person::v::WORKPLACE_IDX); // Get my age demographic auto demographic = FLAMEGPU->getVariable(person::v::AGE_DEMOGRAPHIC); @@ -85,26 +85,24 @@ FLAMEGPU_AGENT_FUNCTION(interactWorkplace, flamegpu::MessageArray, flamegpu::Mes 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(person::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(person::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(person::v::INFECTION_COUNT, FLAMEGPU->getVariable(person::v::INFECTION_COUNT) + 1); - break; - } + // No need to check workplace is correct, individuals only belong in one workplace + // 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(person::v::INFECTION_STATE) == disease::SEIR::InfectionState::Infected) { + // Roll a dice to determine if exposure occurred + 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(person::v::INFECTION_COUNT, FLAMEGPU->getVariable(person::v::INFECTION_COUNT) + 1); + break; } } } @@ -139,33 +137,29 @@ void define(flamegpu::ModelDescription& model, const exateppabm::input::config& // Get a handle to the existing person agent type, which should have already been defined. flamegpu::AgentDescription agent = model.Agent(person::NAME); - // Workplace network variables. @todo - refactor to a separate network location? - agent.newVariable(person::v::WORKPLACE_IDX); + // Workplace index [0, WORKPLACE_COUNT) + agent.newVariable(person::v::WORKPLACE_IDX); + // The number of workplace neighbours for this individual agent.newVariable(person::v::WORKPLACE_OUT_DEGREE); - // Message list containing a persons current status for workplaces (id, location, infection status) - flamegpu::MessageArray::Description workplaceStatusMessage = model.newMessage(person::message::workplace_status::_NAME); + flamegpu::MessageArray::Description workplaceStatusMessage = model.newMessage(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 workplaceStatusMessage.newVariable(person::v::INFECTION_STATE); - // Demographic? + // Demographic workplaceStatusMessage.newVariable(person::v::AGE_DEMOGRAPHIC); - // emit current status to the workplace + // emit current status to the workplace message list flamegpu::AgentFunctionDescription emitWorkplaceStatusDesc = agent.newFunction("emitWorkplaceStatus", emitWorkplaceStatus); - emitWorkplaceStatusDesc.setMessageOutput(person::message::workplace_status::_NAME); + emitWorkplaceStatusDesc.setMessageOutput(message_workplace_status::_NAME); emitWorkplaceStatusDesc.setInitialState(person::states::DEFAULT); emitWorkplaceStatusDesc.setEndState(person::states::DEFAULT); - // Interact with other agents in the workplace via their messages + // Interact with other agents in the workplace via their messages, based on the network graph flamegpu::AgentFunctionDescription interactWorkplaceDesc = agent.newFunction("interactWorkplace", interactWorkplace); - interactWorkplaceDesc.setMessageInput(person::message::workplace_status::_NAME); + interactWorkplaceDesc.setMessageInput(message_workplace_status::_NAME); interactWorkplaceDesc.setInitialState(person::states::DEFAULT); interactWorkplaceDesc.setEndState(person::states::DEFAULT); } @@ -181,5 +175,58 @@ void appendLayers(flamegpu::ModelDescription& model) { } } + + +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::WORKPLACE_SCHOOL_0_9] = (child_network_adults * n_0_9) / n_adult; + p_adultPerWorkNetwork[workplace::Workplace::WORKPLACE_SCHOOL_10_19] = (child_network_adults * n_10_19) / n_adult; + p_adultPerWorkNetwork[workplace::Workplace::WORKPLACE_70_79] = (elderly_network_adults * n_70_79) / n_adult; + p_adultPerWorkNetwork[workplace::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::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::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::WORKPLACE_SCHOOL_0_9; + } else if (age == demographics::Age::AGE_10_19) { + return workplace::Workplace::WORKPLACE_SCHOOL_10_19; + } else if (age == demographics::Age::AGE_70_79) { + return workplace::Workplace::WORKPLACE_70_79; + } else if (age == demographics::Age::AGE_80) { + return workplace::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 workplace } // namespace exateppabm diff --git a/src/exateppabm/workplace.h b/src/exateppabm/workplace.h index d54ee43..28b2447 100644 --- a/src/exateppabm/workplace.h +++ b/src/exateppabm/workplace.h @@ -2,6 +2,7 @@ #include "flamegpu/flamegpu.h" #include "exateppabm/input.h" +#include "exateppabm/demographics.h" namespace exateppabm { namespace workplace { @@ -23,5 +24,48 @@ void define(flamegpu::ModelDescription& model, const exateppabm::input::config& */ void appendLayers(flamegpu::ModelDescription& model); + +/** + * 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) + */ + +workplace::WorkplaceUnderlyingType generateWorkplaceForIndividual(const demographics::Age age, std::array p_adult_workplace, std::mt19937_64 & rng); + + } // namespace workplace } // namespace exateppabm diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index b49016e..511cff6 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -32,9 +32,10 @@ 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_household.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 + ${CMAKE_CURRENT_SOURCE_DIR}/src/exateppabm/test_workplace.cu ) # Define output location of binary files diff --git a/tests/src/exateppabm/test_population.cu b/tests/src/exateppabm/test_household.cu similarity index 50% rename from tests/src/exateppabm/test_population.cu rename to tests/src/exateppabm/test_household.cu index 046f7d0..d494574 100644 --- a/tests/src/exateppabm/test_population.cu +++ b/tests/src/exateppabm/test_household.cu @@ -8,12 +8,12 @@ #include "fmt/core.h" #include "exateppabm/input.h" -#include "exateppabm/population.h" +#include "exateppabm/household.h" /** * Test getReferenceMeanHouseholdSize using a number of manual, partially filled config objects. */ -TEST(TestPopulation, getReferenceMeanHouseholdSize) { +TEST(TestHousehold, getReferenceMeanHouseholdSize) { const double EPSILON = 1.0e-7; double value = 0.0; double expected = 0.0; @@ -27,7 +27,7 @@ TEST(TestPopulation, getReferenceMeanHouseholdSize) { params.household_size_5 = 1; params.household_size_6 = 1; expected = 21 / 6.0; - value = exateppabm::population::getReferenceMeanHouseholdSize(params); + value = exateppabm::household::getReferenceMeanHouseholdSize(params); ASSERT_NEAR(value, expected, EPSILON); // Just houses sizes of 1 @@ -38,7 +38,7 @@ TEST(TestPopulation, getReferenceMeanHouseholdSize) { params.household_size_5 = 0; params.household_size_6 = 0; expected = 1.0 / 1.0; - value = exateppabm::population::getReferenceMeanHouseholdSize(params); + value = exateppabm::household::getReferenceMeanHouseholdSize(params); ASSERT_NEAR(value, expected, EPSILON); // Just houses sizes of 6 @@ -49,7 +49,7 @@ TEST(TestPopulation, getReferenceMeanHouseholdSize) { params.household_size_5 = 0; params.household_size_6 = 1; expected = 6.0 / 1.0; - value = exateppabm::population::getReferenceMeanHouseholdSize(params); + value = exateppabm::household::getReferenceMeanHouseholdSize(params); ASSERT_NEAR(value, expected, EPSILON); // Just houses sizes of 1 or 6 @@ -60,7 +60,7 @@ TEST(TestPopulation, getReferenceMeanHouseholdSize) { params.household_size_5 = 0; params.household_size_6 = 1; expected = 7 / 2.0; - value = exateppabm::population::getReferenceMeanHouseholdSize(params); + value = exateppabm::household::getReferenceMeanHouseholdSize(params); ASSERT_NEAR(value, expected, EPSILON); // Arbitrary mixed values @@ -71,14 +71,14 @@ TEST(TestPopulation, getReferenceMeanHouseholdSize) { params.household_size_5 = 2; params.household_size_6 = 1; expected = 42 / 12.0; - value = exateppabm::population::getReferenceMeanHouseholdSize(params); + value = exateppabm::household::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) { +TEST(TestHousehold, getHouseholdSizeCumulativeProbabilityVector) { const double EPSILON = 1.0e-7; std::vector cumulativeP = {}; std::vector expected = {}; @@ -93,7 +93,7 @@ TEST(TestPopulation, getHouseholdSizeCumulativeProbabilityVector) { 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); + cumulativeP = exateppabm::household::getHouseholdSizeCumulativeProbabilityVector(params); // Check the size ASSERT_EQ(cumulativeP.size(), EXPECTED_LENGTH); // Check the final element is approx equal 1.0 @@ -111,7 +111,7 @@ TEST(TestPopulation, getHouseholdSizeCumulativeProbabilityVector) { 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); + cumulativeP = exateppabm::household::getHouseholdSizeCumulativeProbabilityVector(params); // Check the size ASSERT_EQ(cumulativeP.size(), EXPECTED_LENGTH); // Check the final element is approx equal 1.0 @@ -129,7 +129,7 @@ TEST(TestPopulation, getHouseholdSizeCumulativeProbabilityVector) { 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); + cumulativeP = exateppabm::household::getHouseholdSizeCumulativeProbabilityVector(params); // Check the size ASSERT_EQ(cumulativeP.size(), EXPECTED_LENGTH); // Check the final element is approx equal 1.0 @@ -145,9 +145,9 @@ TEST(TestPopulation, getHouseholdSizeCumulativeProbabilityVector) { * * */ -TEST(TestPopulation, generateHouseholdStructures) { +TEST(TestHousehold, generateHouseholdStructures) { const double EPSILON = 1.0e-7; - std::vector households; + std::vector households; std::mt19937_64 rng = std::mt19937_64(0); exateppabm::input::config params = {}; std::uint64_t totalPeople = 0; @@ -171,7 +171,7 @@ TEST(TestPopulation, generateHouseholdStructures) { params.population_70_79 = 0; params.population_80 = 0; rng.seed(params.rng_seed); - households = exateppabm::population::generateHouseholdStructures(params, rng, false); + households = exateppabm::household::generateHouseholdStructures(params, rng, false); // Should be 1 household ASSERT_EQ(households.size(), 1u); totalPeople = 0u; @@ -221,7 +221,7 @@ TEST(TestPopulation, generateHouseholdStructures) { params.population_70_79 = 1; params.population_80 = 1; rng.seed(params.rng_seed); - households = exateppabm::population::generateHouseholdStructures(params, rng, false); + households = exateppabm::household::generateHouseholdStructures(params, rng, false); // Should be between 1 and 32 households ASSERT_GE(households.size(), 1u); ASSERT_LE(households.size(), 32u); @@ -254,121 +254,3 @@ TEST(TestPopulation, generateHouseholdStructures) { // 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_workplace.cu b/tests/src/exateppabm/test_workplace.cu new file mode 100644 index 0000000..8da92a0 --- /dev/null +++ b/tests/src/exateppabm/test_workplace.cu @@ -0,0 +1,128 @@ +#include +#include +#include +#include +#include + +#include "gtest/gtest.h" + +#include "fmt/core.h" +#include "exateppabm/input.h" +#include "exateppabm/workplace.h" + +/** + * Test getAdultWorkplaceCumulativeProbabilityArray + */ + +TEST(TestWorkplace, 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::workplace::getAdultWorkplaceCumulativeProbabilityArray(child_network_adults, elderly_network_adults, n_per_age); +EXPECT_NEAR(p[exateppabm::workplace::WORKPLACE_SCHOOL_0_9], expected[exateppabm::workplace::WORKPLACE_SCHOOL_0_9], EPSILON); +EXPECT_NEAR(p[exateppabm::workplace::WORKPLACE_SCHOOL_10_19], expected[exateppabm::workplace::WORKPLACE_SCHOOL_10_19], EPSILON); +EXPECT_NEAR(p[exateppabm::workplace::WORKPLACE_ADULT], expected[exateppabm::workplace::WORKPLACE_ADULT], EPSILON); +EXPECT_NEAR(p[exateppabm::workplace::WORKPLACE_70_79], expected[exateppabm::workplace::WORKPLACE_70_79], EPSILON); +EXPECT_NEAR(p[exateppabm::workplace::WORKPLACE_80_PLUS], expected[exateppabm::workplace::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::workplace::getAdultWorkplaceCumulativeProbabilityArray(child_network_adults, elderly_network_adults, n_per_age); +EXPECT_NEAR(p[exateppabm::workplace::WORKPLACE_SCHOOL_0_9], expected[exateppabm::workplace::WORKPLACE_SCHOOL_0_9], EPSILON); +EXPECT_NEAR(p[exateppabm::workplace::WORKPLACE_SCHOOL_10_19], expected[exateppabm::workplace::WORKPLACE_SCHOOL_10_19], EPSILON); +EXPECT_NEAR(p[exateppabm::workplace::WORKPLACE_ADULT], expected[exateppabm::workplace::WORKPLACE_ADULT], EPSILON); +EXPECT_NEAR(p[exateppabm::workplace::WORKPLACE_70_79], expected[exateppabm::workplace::WORKPLACE_70_79], EPSILON); +EXPECT_NEAR(p[exateppabm::workplace::WORKPLACE_80_PLUS], expected[exateppabm::workplace::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::workplace::getAdultWorkplaceCumulativeProbabilityArray(child_network_adults, elderly_network_adults, n_per_age); +EXPECT_NEAR(p[exateppabm::workplace::WORKPLACE_SCHOOL_0_9], expected[exateppabm::workplace::WORKPLACE_SCHOOL_0_9], EPSILON); +EXPECT_NEAR(p[exateppabm::workplace::WORKPLACE_SCHOOL_10_19], expected[exateppabm::workplace::WORKPLACE_SCHOOL_10_19], EPSILON); +EXPECT_NEAR(p[exateppabm::workplace::WORKPLACE_ADULT], expected[exateppabm::workplace::WORKPLACE_ADULT], EPSILON); +EXPECT_NEAR(p[exateppabm::workplace::WORKPLACE_70_79], expected[exateppabm::workplace::WORKPLACE_70_79], EPSILON); +EXPECT_NEAR(p[exateppabm::workplace::WORKPLACE_80_PLUS], expected[exateppabm::workplace::WORKPLACE_80_PLUS], EPSILON); + +// @todo - expand test with more edge case coverage +} + +/** + * Test generateWorkplaceForIndividual + */ +TEST(TestWorkplace, 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::workplace::generateWorkplaceForIndividual(exateppabm::demographics::AGE_0_9, p_adult_workplace, rng); + EXPECT_EQ(w, exateppabm::workplace::Workplace::WORKPLACE_SCHOOL_0_9); +} +// 10-19 always in their workplace +for (std::uint64_t idx = 0; idx < n_samples; ++idx) { + auto w = exateppabm::workplace::generateWorkplaceForIndividual(exateppabm::demographics::AGE_10_19, p_adult_workplace, rng); + EXPECT_EQ(w, exateppabm::workplace::Workplace::WORKPLACE_SCHOOL_10_19); +} +// 70-79 always in their workplace +for (std::uint64_t idx = 0; idx < n_samples; ++idx) { + auto w = exateppabm::workplace::generateWorkplaceForIndividual(exateppabm::demographics::AGE_70_79, p_adult_workplace, rng); + EXPECT_EQ(w, exateppabm::workplace::Workplace::WORKPLACE_70_79); +} +// 80+ always in their workplace +for (std::uint64_t idx = 0; idx < n_samples; ++idx) { + auto w = exateppabm::workplace::generateWorkplaceForIndividual(exateppabm::demographics::AGE_80, p_adult_workplace, rng); + EXPECT_EQ(w, exateppabm::workplace::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::workplace::generateWorkplaceForIndividual(exateppabm::demographics::AGE_40_49, p_adult_workplace, rng); + // Ensure in range + EXPECT_GE(w, 0); + EXPECT_LT(w, exateppabm::workplace::WORKPLACE_COUNT); + // Increment the counter + ++workplaceAdults[w]; +} +// Check that each counter is roughly correct. +EXPECT_NEAR(workplaceAdults[exateppabm::workplace::Workplace::WORKPLACE_SCHOOL_0_9], target_adult_workplace[exateppabm::workplace::Workplace::WORKPLACE_SCHOOL_0_9], TARGET_ADULT_EPSILON); +EXPECT_NEAR(workplaceAdults[exateppabm::workplace::Workplace::WORKPLACE_SCHOOL_10_19], target_adult_workplace[exateppabm::workplace::Workplace::WORKPLACE_SCHOOL_10_19], TARGET_ADULT_EPSILON); +EXPECT_NEAR(workplaceAdults[exateppabm::workplace::Workplace::WORKPLACE_ADULT], target_adult_workplace[exateppabm::workplace::Workplace::WORKPLACE_ADULT], TARGET_ADULT_EPSILON); +EXPECT_NEAR(workplaceAdults[exateppabm::workplace::Workplace::WORKPLACE_70_79], target_adult_workplace[exateppabm::workplace::Workplace::WORKPLACE_70_79], TARGET_ADULT_EPSILON); +EXPECT_NEAR(workplaceAdults[exateppabm::workplace::Workplace::WORKPLACE_80_PLUS], target_adult_workplace[exateppabm::workplace::Workplace::WORKPLACE_80_PLUS], TARGET_ADULT_EPSILON); + +// @todo - expand test with more edge case coverage +}