From 622fe82558635437f1d90b2ccb9eeb9496ea5007 Mon Sep 17 00:00:00 2001 From: Peter Heywood Date: Wed, 11 Dec 2024 08:58:12 +0000 Subject: [PATCH] Refactoring, small world networks & workarounds - Finishes implementing small world networks and uses the resulting graph for interactions within the workplace - Includes several new parameters to influence workplace network generation - Now (and in prev commits) 5 workplace networks - Refactors code, including the graph data structure once again. More refactoring to come. - Expands test suite - Switches to agent generation in an init function, but workarounds are required due to bugs in flame gpu 2, of which a PR is already in place to address at least one of. - Vis adjustments --- src/exateppabm/exatepp_abm.cu | 15 +- src/exateppabm/network.cu | 61 +++---- src/exateppabm/network.h | 210 +++++++++++++++++++---- src/exateppabm/output.cu | 2 +- src/exateppabm/person.cu | 103 +++++------ src/exateppabm/person.h | 3 +- src/exateppabm/population.cu | 247 ++++++++++++++++++++------- src/exateppabm/population.h | 6 +- src/exateppabm/visualisation.cu | 27 ++- src/exateppabm/visualisation.h | 15 +- tests/src/exateppabm/test_network.cu | 132 +++++++------- tools/plot-per-individual.py | 8 +- 12 files changed, 551 insertions(+), 278 deletions(-) diff --git a/src/exateppabm/exatepp_abm.cu b/src/exateppabm/exatepp_abm.cu index 5a439cc..8a9c9d3 100644 --- a/src/exateppabm/exatepp_abm.cu +++ b/src/exateppabm/exatepp_abm.cu @@ -91,6 +91,9 @@ int entrypoint(int argc, char* argv[]) { // Define disease related variables and methods exateppabm::disease::SEIR::define(model, *config); + // Define init function for population generation + exateppabm::population::define(model, *config, cli_params->verbosity > 0); + // Add init, step and exit functions related to data collection and output. This may want refactoring when multiple output files are supported or collected data becomes more complex. exateppabm::output::define(model, cli_params->outputDir, cli_params->individualFile); @@ -100,15 +103,6 @@ int entrypoint(int argc, char* argv[]) { // Add disease progression exateppabm::disease::SEIR::appendLayers(model); - // Generate the population of agents - // prior to the simulation creation so any mutation of the model environment is applied. this is not ideal and will need adjusting for ensembles. - // @todo - this should probably be an in init function for ease of moving to a ensembles, but then cannot pass parameters in. - const std::uint64_t pop_seed = config->rng_seed; // @todo - split seeds - auto personPopulation = exateppabm::population::generate(model, *config, cli_params->verbosity > 0); - if (personPopulation == nullptr) { - throw std::runtime_error("@todo - bad population generation function."); - } - // Construct the Simulation instance from the model. flamegpu::CUDASimulation simulation(model); @@ -133,9 +127,6 @@ int entrypoint(int argc, char* argv[]) { // Set the GPU index simulation.CUDAConfig().device_id = cli_params->device; - // add the population to this simulation instance - simulation.setPopulationData(*personPopulation); - perfFile.timers.preSimulate.stop(); // Run the simulation diff --git a/src/exateppabm/network.cu b/src/exateppabm/network.cu index 2a14470..dbb3c4d 100644 --- a/src/exateppabm/network.cu +++ b/src/exateppabm/network.cu @@ -15,19 +15,18 @@ namespace exateppabm { namespace network { -UndirectedNetwork generateFullyConnectedUndirectedNetwork(std::vector nodes) { - UndirectedNetwork network = {}; - network.nodes = nodes; +UndirectedGraph generateFullyConnectedUndirectedGraph(std::vector nodes) { + UndirectedGraph network(nodes); if (nodes.size() <= 1) { return network; } std::uint32_t N = static_cast(nodes.size()); std::uint32_t E = (N * (N - 1)) / 2; - network.edges.reserve(E); + // network.edges.reserve(E); // Generate edges in ascending order for the undirected graph. for (std::uint32_t i = 0; i < N; ++i) { for (std::uint32_t j = i + 1; j < N; ++j) { - network.edges.push_back({i, j}); + network.addEdge(i, j); } } return network; @@ -35,13 +34,11 @@ UndirectedNetwork generateFullyConnectedUndirectedNetwork(std::vector nodes, std::uint32_t K, double p_rewire, std::mt19937_64 rng) { - // fmt::print("@todo - check if small world network should be directed or undirected\n"); - UndirectedNetwork network = {}; - network.nodes = nodes; +UndirectedGraph generateSmallWorldUndirectedGraph(std::vector nodes, std::uint32_t K, double p_rewire, std::mt19937_64 rng) { + UndirectedGraph network(nodes); // Return early with an empty network if 0 or 1 nodes, so no room for non-self edges - if (network.nodes.size() <= 1) { + if (network.getNumVertices() <= 1) { return network; } @@ -49,12 +46,12 @@ UndirectedNetwork generateSmallWorldUndirectedNetwork(std::vector 1"); return network; - } else if (K == network.nodes.size()) { + } else if (K == network.getNumVertices()) { // If K == Node count, the graph is fully connected, so return a fully connected graph - return generateFullyConnectedUndirectedNetwork(nodes); - } else if (K >= network.nodes.size()) { + return generateFullyConnectedUndirectedGraph(nodes); + } else if (K >= network.getNumVertices()) { // Raise an exception if K is too large - throw std::runtime_error(std::string("@todo values - small world network K(") + std::to_string(K) + ") must be less than |N| (" + std::to_string(network.nodes.size()) + ")"); + throw std::runtime_error(std::string("@todo values - small world network K(") + std::to_string(K) + ") must be less than |N| (" + std::to_string(network.getNumVertices()) + ")"); } // If K is odd, use K-1 if (K % 2 == 1) { @@ -70,19 +67,14 @@ UndirectedNetwork generateSmallWorldUndirectedNetwork(std::vector(nodes.size()); std::int32_t E = (N * K) / 2; - network.edges.reserve(E); + + // network.edges.reserve(E); for (std::int32_t i = 0; i < N; ++i) { for (std::int32_t j = 1; j <= static_cast(K / 2); ++j) { // only add the positive undirected edges std::uint32_t s = static_cast(i); std::uint32_t d = static_cast((i + j) % N); - // Ensure that the edge source node is lower than the destination, to simplify avoiding duplicate edges - if (s > d) { - std::swap(s, d); - } - network.edges.push_back({s, d}); - // network.edges.push_back({static_cast(i), static_cast((i + j) % N)}); - // network.edges.push_back({static_cast(i), static_cast((i - j + N) % N)}); + network.addEdge(s, d); } } @@ -99,18 +91,27 @@ UndirectedNetwork generateSmallWorldUndirectedNetwork(std::vector dest_dist(0, N-1); // Randomly rewire edges - for (auto& edge : network.edges) { + + // Take a copy of the network edges, to ensure we avoid iterator invalidation + // Only the current iterator should be removed in the undirected graph, so this is probably ok. + std::vector copyOfEdges(network.getEdges().begin(), network.getEdges().end()); + for (auto& edge : copyOfEdges) { if (p_dist(rng) < p_rewire) { + // If the source vertex is full, do not attempt to rewire it, it would loop forever so check the next original edge + if (network.degree(edge.source) >= network.getNumVertices() - 1) { + continue; + } + + // Repeatedly generate a new destination vertex until a new edge has been found. std::uint32_t newDest = dest_dist(rng); - // Repeat the generation until the edge is not a self-edge, nor duplicate - // Not using do while so the search for the edge is ok. - // @note - this will be a challenge to parallelise efficiently without atomics (i.e. stable matching) - // @note - this will be expensive for large networks, an edge matrix would be cheaper (at the cost of N*N memory) - // @todo - need to avoid directed edge duplicates. - while (newDest == edge.source || std::find(network.edges.begin(), network.edges.end(), UndirectedNetwork::Edge{edge.source, newDest}) != network.edges.end()) { + // While the new dest is the source, or already in the graph, try again. + while (newDest == edge.source || network.contains(edge.source, newDest)) { newDest = dest_dist(rng); } - edge.dest = newDest; + // Remove the old edge + network.removeEdge(edge.source, edge.dest); + // Add the new edge + network.addEdge(edge.source, newDest); } } } diff --git a/src/exateppabm/network.h b/src/exateppabm/network.h index c5f9c63..757a6b4 100644 --- a/src/exateppabm/network.h +++ b/src/exateppabm/network.h @@ -1,8 +1,10 @@ #pragma once #include +#include #include #include +#include #include #include "flamegpu/flamegpu.h" @@ -10,46 +12,196 @@ namespace exateppabm { namespace network { + /** - * Struct representing a network. - * - * @todo - re-use flamegpu's networks? Not sure this is viable yet. - * @todo - make it a class - * @todo - template the index type? Map indices to ID so a sparse matrix would make sense? + * Struct representing an Edge */ -struct UndirectedNetwork { - /** - * Struct representing an Edge - */ - struct Edge { - std::uint32_t source; - std::uint32_t dest; - /** - * const operator== comparison for std::find - */ - bool operator==(const Edge& other) const { - return source == other.source && dest == other.dest; - } - /** - * operator< for use in a set - */ - bool operator<(const Edge& other) const { +class Edge { + public: + std::uint32_t source; + std::uint32_t dest; + /** + * constrcutor for emplace_back support + */ + Edge(std::uint32_t source, std::uint32_t dest) : source(source), dest(dest) {} + /** + * const operator== comparison for std::find + */ + bool operator==(const Edge& other) const { + return source == other.source && dest == other.dest; + } + /** + * operator< for use in a set + */ + bool operator<(const Edge& other) const { if (source != other.source) { return source < other.source; } else { return dest < other.dest; } } - }; +}; + +/** + * class representing an undirected network, using edge lists as well suited to edge existan + * + * @todo - template the index type? Map indices to ID so a sparse matrix would make sense? + */ +class UndirectedGraph { + public: + /** + * Default ctor + */ + UndirectedGraph() { } + + /** + * Constructor, taking a vector of vertex labels + * + * @param vertexLabels vertex labels, setting the number of vertices for the network + */ + explicit UndirectedGraph(std::vector vertexLabels) : _vertexLabels(vertexLabels) { + this->_adjacencyList.resize(this->getNumVertices()); + } + + /** + * copy constructor + * + * @param other instance of UndirectedGraph to copy + */ + UndirectedGraph(const UndirectedGraph& other) : + _vertexLabels(other._vertexLabels), + _edges(other._edges), + _adjacencyList(other._adjacencyList) {} + + /** + * Get the number of nodes + * @return number of vertices/nodes + */ + std::uint32_t getNumVertices() const { return this->_vertexLabels.size(); } + /** + * Get the number of edges + * @return number of edges + */ + std::uint32_t getNumEdges() const { return this->_edges.size(); } + /** + * Check if an edge exists within the undirected graph + * @param u vertex index for one end of the edge + * @param v vertex index for the other end of the edge + * @return bool indicating if edge exists or not + */ + bool contains(std::uint32_t u, std::uint32_t v) const { + return this->_adjacencyList[u].count(v) > 0; + } + /** + * Check if an edge exists within the undirected graph + * @param edge to check for existance of + * @return bool indicating if edge exists or not + */ + bool contains(Edge edge) const { + return this->_adjacencyList[edge.source].count(edge.dest) > 0; + } + /** + * Add an edge to the undirected graph + * @param u vertex index for one end of the edge + * @param v vertex index for the other end of the edge + */ + void addEdge(std::uint32_t u, std::uint32_t v) { + if (u > this->getNumVertices() || v > this->getNumVertices()) { + throw std::runtime_error("Invalid u or v for UndirectedGraph::getNumVertices. @todo better error"); + } + // Undirected, so add the edge in both directions + this->_adjacencyList[u].insert(v); + this->_adjacencyList[v].insert(u); + // Only store the ascending ordered edge in the vector of edges + if (u < v) { + this->_edges.emplace_back(u, v); + } else { + this->_edges.emplace_back(v, u); + } + } + /** + * Add an edge to the undirected graph + * @param edge edge to add + */ + void addEdge(Edge edge) { + return this->addEdge(edge.source, edge.dest); + } + /** + * Remove an edge from the undirected graph + * @param edge edge to remove + */ + void removeEdge(Edge edge) { + this->_adjacencyList[edge.source].erase(edge.dest); + this->_adjacencyList[edge.dest].erase(edge.source); + // Remove the edge from the edges_ vector + auto it = std::remove_if(this->_edges.begin(), this->_edges.end(), [&edge](const Edge& e) { return e == edge; }); + this->_edges.erase(it, this->_edges.end()); + } + /** + * Remove an edge from the undirected graph + * @param u vertex index for one end of the edge + * @param v vertex index for the other end of the edge + */ + void removeEdge(std::uint32_t u, std::uint32_t v) { + return this->removeEdge({u, v}); + } + /** + * Get the degree for a vertex. this is in-out as undirected. + * + * @return vertex degree + */ + std::uint32_t degree(std::uint32_t v) { + return static_cast(this->_adjacencyList[v].size()); + } + + /** + * Get the vertex indices of neighbouring vertices + * @param v vertex index + * @return unordered set of neighbouring vertices to v + */ + const std::unordered_set& getNeighbours(std::uint32_t v) const { + return this->_adjacencyList[v]; + } + + /** + * Get the vector of edges + * + * @return vector of edges + */ + const std::vector& getEdges() const { return this->_edges; } + + /** + * Get the vertex labels, which are probably 1-indexed FLAME GPU id_t (0 is unset value) + * + * @return vector of vertex labels + */ + const std::vector& getVertexLabels() const { return this->_vertexLabels; } + + /** + * Get the a vertex label + * + * @param v vertex index + * @return label for vertex + */ + flamegpu::id_t getVertexLabel(std::uint32_t v) const { return this->_vertexLabels.at(v); } + + private: /** * Vector of nodes, which provides the mapping from the network-index based edges back to agent IDs */ - std::vector nodes; + std::vector _vertexLabels; + + /** * Vector of edge objects, each containing the directed source-dest pair. These are 0 indexed within this network. */ - std::vector edges; + std::vector _edges; + + /** + * Vector of unordered sets for fast lookup of edge existence. + */ + std::vector> _adjacencyList; }; /** @@ -57,19 +209,19 @@ struct UndirectedNetwork { * @param nodes vector of nodes (agent Id's) to include in this network. IDs may not be sequential and should be one indexed. * @return network struct for a fully connected network */ -UndirectedNetwork generateFullyConnectedUndirectedNetwork(std::vector nodes); +UndirectedGraph generateFullyConnectedUndirectedGraph(std::vector nodes); /** * Generate a watts strogatz small world (undirected) network for a given set of nodes (agent ID's, 1 indexed, not sequential), K & p. - * + * * If k is odd, then k-1 neighboours are created - * + * * @param nodes vector of nodes (i.e. agent id's) to include in this network. IDs may not be sequential and should be 1 indexed. * @param K the degree of each node * @param p rewire probability * @return a network struct for the generate small world network */ -UndirectedNetwork generateSmallWorldUndirectedNetwork(std::vector nodes, std::uint32_t K, double p, std::mt19937_64 rng); +UndirectedGraph generateSmallWorldUndirectedGraph(std::vector nodes, std::uint32_t K, double p, std::mt19937_64 rng); } // namespace network diff --git a/src/exateppabm/output.cu b/src/exateppabm/output.cu index 732138d..95382a8 100644 --- a/src/exateppabm/output.cu +++ b/src/exateppabm/output.cu @@ -113,7 +113,7 @@ FLAMEGPU_EXIT_FUNCTION(output_exit_per_individual) { flamegpu::DeviceAgentVector population = personAgent.getPopulationData(); for (const auto& person : population) { exateppabm::output::PerIndividualFile::Person personData = {}; - personData.id = static_cast(person.getID()); + personData.id = static_cast(person.getVariable(person::v::ID)); personData.age_group = person.getVariable(person::v::AGE_DEMOGRAPHIC); personData.occupation_network = person.getVariable(person::v::WORKPLACE_IDX); personData.house_no = person.getVariable(person::v::HOUSEHOLD_IDX); diff --git a/src/exateppabm/person.cu b/src/exateppabm/person.cu index d42d02c..b0b917c 100644 --- a/src/exateppabm/person.cu +++ b/src/exateppabm/person.cu @@ -21,7 +21,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->getID()); + FLAMEGPU->message_out.setVariable(person::message::household_status::ID, FLAMEGPU->getVariable(person::v::ID)); // Household index // @todo - typedef or using statement for the household index type? @@ -47,7 +47,7 @@ FLAMEGPU_AGENT_FUNCTION(emitHouseholdStatus, flamegpu::MessageNone, flamegpu::Me */ FLAMEGPU_AGENT_FUNCTION(interactHousehold, flamegpu::MessageBucket, flamegpu::MessageNone) { // Get my ID to avoid self messages - const flamegpu::id_t id = FLAMEGPU->getID(); + const flamegpu::id_t id = FLAMEGPU->getVariable(person::v::ID); // FLAMEGPU->getID(); // Get the probability of infection float p_s2e = FLAMEGPU->environment.getProperty("p_interaction_susceptible_to_exposed"); @@ -108,31 +108,26 @@ FLAMEGPU_AGENT_FUNCTION(interactHousehold, flamegpu::MessageBucket, flamegpu::Me return flamegpu::ALIVE; } - - /** * 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::MessageBucket) { - // workplaces of 1 don't need to do any messaging, there is no one to infect - std::uint32_t workplaceSize = FLAMEGPU->getVariable(v::WORKPLACE_SIZE); - if (workplaceSize > 1) { - // output public properties to bucket message, keyed by workplace - // Agent ID to avoid self messaging - FLAMEGPU->message_out.setVariable(person::message::workplace_status::ID, FLAMEGPU->getID()); +FLAMEGPU_AGENT_FUNCTION(emitWorkplaceStatus, flamegpu::MessageNone, flamegpu::MessageArray) { + // output public properties to bucket message, keyed by workplace + // Agent ID to avoid self messaging + flamegpu::id_t id = FLAMEGPU->getVariable(person::v::ID); + FLAMEGPU->message_out.setVariable(person::message::workplace_status::ID, id); - // workplace index - // @todo - typedef or using statement for the workplace index type? - std::uint32_t workplaceIdx = FLAMEGPU->getVariable(v::WORKPLACE_IDX); - FLAMEGPU->message_out.setVariable(v::WORKPLACE_IDX, workplaceIdx); + // workplace index + // @todo - typedef or using statement for the workplace index type? + std::uint32_t workplaceIdx = FLAMEGPU->getVariable(v::WORKPLACE_IDX); + FLAMEGPU->message_out.setVariable(v::WORKPLACE_IDX, workplaceIdx); - FLAMEGPU->message_out.setVariable(v:: - INFECTION_STATE, FLAMEGPU->getVariable(v::INFECTION_STATE)); - FLAMEGPU->message_out.setVariable(v::AGE_DEMOGRAPHIC, FLAMEGPU->getVariable(v::AGE_DEMOGRAPHIC)); + FLAMEGPU->message_out.setVariable(v:: + INFECTION_STATE, FLAMEGPU->getVariable(v::INFECTION_STATE)); + FLAMEGPU->message_out.setVariable(v::AGE_DEMOGRAPHIC, FLAMEGPU->getVariable(v::AGE_DEMOGRAPHIC)); - // Set the message key, the house hold idx for bucket messaging @Todo - FLAMEGPU->message_out.setKey(workplaceIdx); - } + // Set the message key, the house hold idx for bucket messaging @Todo + FLAMEGPU->message_out.setIndex(id); return flamegpu::ALIVE; } @@ -144,14 +139,20 @@ FLAMEGPU_AGENT_FUNCTION(emitWorkplaceStatus, flamegpu::MessageNone, flamegpu::Me * @todo - refactor this somewhere else? * @todo - add per network behaviours? */ -FLAMEGPU_AGENT_FUNCTION(interactWorkplace, flamegpu::MessageBucket, flamegpu::MessageNone) { - // Get my ID to avoid self messages - const flamegpu::id_t id = FLAMEGPU->getID(); +FLAMEGPU_AGENT_FUNCTION(interactWorkplace, flamegpu::MessageArray, flamegpu::MessageNone) { + // Get my ID + const flamegpu::id_t id = FLAMEGPU->getVariable(person::v::ID); + + // Get my workplace degree, if 0 return. + const std::uint32_t degree = FLAMEGPU->getVariable(person::v::WORKPLACE_OUT_DEGREE); + if (degree == 0) { + return flamegpu::ALIVE; + } // Get the probability of interaction within the workplace float p_daily_fraction_work = FLAMEGPU->environment.getProperty("daily_fraction_work"); - // Get the probability of infection + // Get the probability of an interaction leading to exposure float p_s2e = FLAMEGPU->environment.getProperty("p_interaction_susceptible_to_exposed"); // Scale it for workplace interactions p_s2e *= FLAMEGPU->environment.getProperty("relative_transmission_occupation"); @@ -163,23 +164,27 @@ FLAMEGPU_AGENT_FUNCTION(interactWorkplace, flamegpu::MessageBucket, flamegpu::Me auto demographic = FLAMEGPU->getVariable(v::AGE_DEMOGRAPHIC); // Get my susceptibility modifier and modify it. float relativeSusceptibility = FLAMEGPU->environment.getProperty("relative_susceptibility_per_demographic", demographic); - // Scale the probability of transmission + // Scale the probability of transmission for my age demographic p_s2e *= relativeSusceptibility; // Check if the current individual is susceptible to being infected auto infectionState = FLAMEGPU->getVariable(v::INFECTION_STATE); -// @todo - only interact with some others - // @todo - this will need to change for contact tracing, the message interaction will need to occur regardless. if (infectionState == disease::SEIR::Susceptible) { // Variable to store the duration of the exposed phase (if exposed) float stateDuration = 0.f; - - // Iterate messages from anyone within my workplace - for (const auto &message : FLAMEGPU->message_in(workplaceIdx)) { - // Ignore self messages (can't infect oneself) - if (message.getVariable(message::workplace_status::ID) != id) { + // Iterate my downstream neighbours (the graph is undirected, so no need to iterate in and out + auto workplaceGraph = FLAMEGPU->environment.getDirectedGraph("WORKPLACE_DIGRAPH"); + std::uint32_t myVertexIndex = workplaceGraph.getVertexIndex(id); + for (auto& edge : workplaceGraph.outEdges(myVertexIndex)) { + // Get the neighbour vertex index and then neighbours agent Id from the edge + std::uint32_t otherIndex = edge.getEdgeDestination(); + flamegpu::id_t otherID = workplaceGraph.getVertexID(otherIndex); + // if the ID is not self and not unset + if (otherID != id && otherID != flamegpu::ID_NOT_SET) { + // Get the message handle + const auto &message = FLAMEGPU->message_in.at(otherID); // Ignore messages from other workplaces if (message.getVariable(v::WORKPLACE_IDX) == workplaceIdx) { // roll a dice to determine if this interaction should occur this day @@ -253,7 +258,7 @@ FLAMEGPU_HOST_FUNCTION(updateRandomDailyNetworkIndices) { size_t first = 0; size_t last = 0; for (const auto &person : population) { - flamegpu::id_t id = person.getID(); + flamegpu::id_t id = person.getVariable(person::v::ID); std::uint32_t count = person.getVariable(person::v::RANDOM_INTERACTION_COUNT_TARGET); last = first + count; if (last > randomInteractionCountSum) { @@ -293,7 +298,7 @@ FLAMEGPU_HOST_FUNCTION(updateRandomDailyNetworkIndices) { auto aAgent = population[aID - 1]; auto bAgent = population[bID - 1]; // Raise an exception if the agents are out of order. This should not be the case as agent birth, death and sorting should not be included in this model. - if (aAgent.getID() != aID || bAgent.getID() != bID) { + if (aAgent.getVariable(person::v::ID) != aID || bAgent.getVariable(person::v::ID) != bID) { throw std::runtime_error("Agent ID does not match expected agent ID in updateRandomDailyNetworkIndices. @todo"); } @@ -340,14 +345,14 @@ FLAMEGPU_HOST_FUNCTION(updateRandomDailyNetworkIndices) { */ FLAMEGPU_AGENT_FUNCTION(emitRandomDailyNetworkStatus, flamegpu::MessageNone, flamegpu::MessageArray) { // output public properties to array message - // FLAMEGPU->message_out.setVariable(person::message::workplace_status::ID, FLAMEGPU->getID()); + // FLAMEGPU->message_out.setVariable(person::message::random_daily_status::ID, FLAMEGPU->getVariable(person::v::ID)); FLAMEGPU->message_out.setVariable(v:: INFECTION_STATE, FLAMEGPU->getVariable(v::INFECTION_STATE)); FLAMEGPU->message_out.setVariable(v::AGE_DEMOGRAPHIC, FLAMEGPU->getVariable(v::AGE_DEMOGRAPHIC)); // Set the message array message index to the agent's id. - FLAMEGPU->message_out.setIndex(FLAMEGPU->getID()); + FLAMEGPU->message_out.setIndex(FLAMEGPU->getVariable(person::v::ID)); return flamegpu::ALIVE; } @@ -362,7 +367,7 @@ FLAMEGPU_AGENT_FUNCTION(emitRandomDailyNetworkStatus, flamegpu::MessageNone, fla */ FLAMEGPU_AGENT_FUNCTION(interactRandomDailyNetwork, flamegpu::MessageArray, flamegpu::MessageNone) { // Get my ID to avoid self messages - const flamegpu::id_t id = FLAMEGPU->getID(); + const flamegpu::id_t id = FLAMEGPU->getVariable(person::v::ID); // Get the probability of infection float p_s2e = FLAMEGPU->environment.getProperty("p_interaction_susceptible_to_exposed"); @@ -449,6 +454,9 @@ void define(flamegpu::ModelDescription& model, const exateppabm::input::config& agent.newState(person::states::DEFAULT); // Define variables + // Custom ID, as getID cannot be relied upon to start at 1. + agent.newVariable(person::v::ID, flamegpu::ID_NOT_SET); + // disease related variables // @todo - define this in disease/ call a disease::SEIR::define_person() like method? agent.newVariable(person::v::INFECTION_STATE, disease::SEIR::Susceptible); @@ -470,15 +478,14 @@ void define(flamegpu::ModelDescription& model, const exateppabm::input::config& // Workplace environment directed graph // This single graph contains workplace information for all individuals, and is essentially 5 unconnected sub graphs. - flamegpu::EnvironmentDirectedGraphDescription workplaceDigraphDesc= env.newDirectedGraph("WORKPLACE_DIGRAPH"); - // Graph vertices - workplaceDigraphDesc.newVertexProperty("bar"); - workplaceDigraphDesc.newEdgeProperty("foo"); - + flamegpu::EnvironmentDirectedGraphDescription workplaceDigraphDesc = env.newDirectedGraph("WORKPLACE_DIGRAPH"); + // Graph vertices + // workplaceDigraphDesc.newVertexProperty("bar"); + // workplaceDigraphDesc.newEdgeProperty("foo"); // Workplace network variables. @todo - refactor to a separate network location? agent.newVariable(person::v::WORKPLACE_IDX); - agent.newVariable(person::v::WORKPLACE_SIZE); + agent.newVariable(person::v::WORKPLACE_OUT_DEGREE); // Random interaction network variables. @todo -refactor to separate location agent.newVariable(person::v::RANDOM_INTERACTION_PARTNERS, {flamegpu::ID_NOT_SET}); @@ -516,11 +523,10 @@ void define(flamegpu::ModelDescription& model, const exateppabm::input::config& householdStatusMessage.newVariable(person::v::AGE_DEMOGRAPHIC); // Message list containing a persons current status for workplaces (id, location, infection status) - flamegpu::MessageBucket::Description workplaceStatusMessage = model.newMessage(person::message::workplace_status::_NAME); + flamegpu::MessageArray::Description workplaceStatusMessage = model.newMessage(person::message::workplace_status::_NAME); - // Set the maximum bucket index to the population size, to the maximum number of workplace networks - // @todo - this will be replaced with a per-person message when improved network messaging is implemented (where individuals will communicate with their direct network) - workplaceStatusMessage.setUpperBound(exateppabm::population::WORKPLACE_COUNT); + // 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); @@ -551,6 +557,7 @@ 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.setMessageOutputOptional(true); emitHouseholdStatusDesc.setInitialState(person::states::DEFAULT); emitHouseholdStatusDesc.setEndState(person::states::DEFAULT); diff --git a/src/exateppabm/person.h b/src/exateppabm/person.h index 84ebad3..1422b2e 100644 --- a/src/exateppabm/person.h +++ b/src/exateppabm/person.h @@ -48,6 +48,7 @@ namespace v { DEVICE_CONSTEXPR_STRING constexpr char x[] = "x"; DEVICE_CONSTEXPR_STRING constexpr char y[] = "y"; DEVICE_CONSTEXPR_STRING constexpr char z[] = "z"; +DEVICE_CONSTEXPR_STRING constexpr char ID[] = "ID"; // Custom id, as FLAMEGPU getID isn't guaranteed to start at 1. DEVICE_CONSTEXPR_STRING constexpr char INFECTION_STATE[] = "infection_state"; DEVICE_CONSTEXPR_STRING constexpr char INFECTION_STATE_CHANGE_DAY[] = "infection_state_change_day"; DEVICE_CONSTEXPR_STRING constexpr char INFECTION_STATE_DURATION[] = "infection_state_duration"; @@ -56,7 +57,7 @@ DEVICE_CONSTEXPR_STRING constexpr char AGE_DEMOGRAPHIC[] = "age_demographic"; DEVICE_CONSTEXPR_STRING constexpr char HOUSEHOLD_IDX[] = "household_idx"; DEVICE_CONSTEXPR_STRING constexpr char HOUSEHOLD_SIZE[] = "household_size"; DEVICE_CONSTEXPR_STRING constexpr char WORKPLACE_IDX[] = "workplace_idx"; -DEVICE_CONSTEXPR_STRING constexpr char WORKPLACE_SIZE[] = "workplace_size"; +DEVICE_CONSTEXPR_STRING constexpr char WORKPLACE_OUT_DEGREE[] = "workplace_out_degree"; DEVICE_CONSTEXPR_STRING constexpr char RANDOM_INTERACTION_PARTNERS[] = "random_interaction_partners"; DEVICE_CONSTEXPR_STRING constexpr char RANDOM_INTERACTION_COUNT[] = "random_interaction_count"; DEVICE_CONSTEXPR_STRING constexpr char RANDOM_INTERACTION_COUNT_TARGET[] = "random_interaction_count_target"; diff --git a/src/exateppabm/population.cu b/src/exateppabm/population.cu index 7407b99..686b95f 100644 --- a/src/exateppabm/population.cu +++ b/src/exateppabm/population.cu @@ -24,43 +24,78 @@ namespace population { namespace { // File-scoped array containing the number of infected agents per demographic from population initialisation. This needs to be made accessible to a FLAME GPU Init func due to macro environment property limitations. - std::array _infectedPerDemographic = {}; +// File-scoped copy of the model parameters struct +exateppabm::input::config _params = {}; + +// File-scoped boolean for if generation should be verbose or not. +bool _verbose = false; + +// Store the ID of each individual for each workplace, to form the node labels of the small world network +std::array, WORKPLACE_COUNT> workplaceMembers = {{}}; + } // namespace -std::unique_ptr generate(flamegpu::ModelDescription& model, const exateppabm::input::config params, const bool verbose) { + +// Workaround struct representing a person used during host generation +struct HostPerson { + flamegpu::id_t id; + disease::SEIR::InfectionStateUnderlyingType infectionStatus; + float infectionStateDuration; + std::uint32_t infectionCount = 0; + demographics::AgeUnderlyingType ageDemographic; + std::uint32_t householdIdx; + std::uint8_t householdSize; + std::uint32_t workplaceIdx; + std::uint32_t randomInteractionTarget; + std::uint32_t workplaceOutDegree; +}; + +/** + * FLAME GPU init function which generates a population of person agents for the current simulation + * + * This function has been split into multiple init functions, due to errors encountered when attempting to do multiple passes over a newly created set of agents. There should be a working way to do this... + */ +FLAMEGPU_INIT_FUNCTION(generatePopulation) { fmt::print("@todo - validate params inputs when generated agents (pop size, initial infected count etc)\n"); // n_total must be less than uint max // n_total must be more than 0. // initial infected count must be more than 0, less than full pop. // Get a handle on the environment. - auto env = model.Environment(); + auto env = FLAMEGPU->environment; - // @todo - assert that the requested initial population is non zero. - auto pop = std::make_unique(model.Agent(exateppabm::person::NAME), params.n_total); + // Get a handle to the host agent api for Person agents + // auto personAgent = FLAMEGPU->agent(exateppabm::person::NAME, exateppabm::person::states::DEFAULT); // Do nothing if no one to create? - if (params.n_total == 0) { - return pop; + if (_params.n_total == 0) { + return; } + // Workaround: Generating agents in FLAME GPU 2 in init functions does not allow iterating the data structure multiple times, this appears to be a bug (which may have a fix in place, currently untested pr) + + // Instead, we must generate data on the host for the "first pass", storing it in non-flame gpu storage + // Then set it on the latter pass. + std::vector hostPersonData; + hostPersonData.resize(_params.n_total); + // seed host rng for population generation. // @todo - does this want to be a separate seed from the config file? - std::mt19937_64 rng(params.rng_seed); + std::mt19937_64 rng(FLAMEGPU->getSimulationConfig().random_seed); // Need to initialise a fixed number of individuals as infected. // This not very scalable way of doing it, is to create a vector with one element per individual in the simulation, initialised to false // set the first n_seed_infection elements to true/1 // Shuffle the vector, and query at agent creation time // RNG sampling in-loop would be more memory efficient, but harder to guarantee that exactly enough are created. This will likely be replaced anyway, so quick and dirty is fine. - std::vector infected_vector(params.n_total); - std::fill(infected_vector.begin(), infected_vector.begin() + std::min(params.n_total, params.n_seed_infection), true); + std::vector infected_vector(_params.n_total); + std::fill(infected_vector.begin(), infected_vector.begin() + std::min(_params.n_total, _params.n_seed_infection), true); std::shuffle(infected_vector.begin(), infected_vector.end(), rng); // Get the number of individuals per house and their age demographics - auto households = generateHouseholdStructures(params, rng, verbose); + auto households = generateHouseholdStructures(_params, rng, _verbose); // per demo total is not an output in time series. // Alternately, we need to initialise the exact number of each age band, not RNG, and just scale it down accordingly. Will look at in "realistic" population generation @@ -76,9 +111,9 @@ std::unique_ptr generate(flamegpu::ModelDescription& mode _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); - // 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 p_adultPerWorkNetwork = getAdultWorkplaceCumulativeProbabilityArray(_params.child_network_adults, _params.elderly_network_adults, createdPerDemographic); + // // Store the ID of each individual for each workplace, to form the node labels of the small world network + // std::array, WORKPLACE_COUNT> workplaceMembers = {{}}; // Counter for random interaction count. This will be 2x the number of interactions std::uint64_t randomInteractionCountSum = 0u; @@ -96,24 +131,35 @@ std::unique_ptr generate(flamegpu::ModelDescription& mode assert(household.size == household.agePerPerson.size()); // Get the flamegpu person object for the individual - auto person = pop->at(personIdx); + // auto person = personAgent.newAgent(); @temp disabled due to bug/workaround in place for multi-pass + + // Set the manual ID, as getID() isn't guaranteed to actually start at 1. + flamegpu::id_t personID = personIdx + 1; + // person.setVariable(person::v::ID, personID); + hostPersonData[personIdx].id = personID; // Set the individuals infection status. @todo - refactor into seir.cu? disease::SEIR::InfectionState infectionStatus = infected_vector.at(personIdx) ? disease::SEIR::InfectionState::Infected : disease::SEIR::InfectionState::Susceptible; - person.setVariable(exateppabm::person::v::INFECTION_STATE, infectionStatus); + // person.setVariable(exateppabm::person::v::INFECTION_STATE, infectionStatus); + hostPersonData[personIdx].infectionStatus = infectionStatus; // Also set the initial infection duration. @todo - stochastic. - float infectionStateDuration = infectionStatus == disease::SEIR::InfectionState::Infected ? params.mean_time_to_recovered: 0; - person.setVariable(exateppabm::person::v::INFECTION_STATE_DURATION, infectionStateDuration); + float infectionStateDuration = infectionStatus == disease::SEIR::InfectionState::Infected ? _params.mean_time_to_recovered: 0; + // person.setVariable(exateppabm::person::v::INFECTION_STATE_DURATION, infectionStateDuration); + hostPersonData[personIdx].infectionStateDuration = infectionStateDuration; // Set the individuals age and household properties demographics::Age age = household.agePerPerson[householdMemberIdx]; - person.setVariable(person::v::AGE_DEMOGRAPHIC, static_cast(age)); - person.setVariable(person::v::HOUSEHOLD_IDX, householdIdx); - person.setVariable(person::v::HOUSEHOLD_SIZE, household.size); + // person.setVariable(person::v::AGE_DEMOGRAPHIC, static_cast(age)); + hostPersonData[personIdx].ageDemographic = static_cast(age); + // person.setVariable(person::v::HOUSEHOLD_IDX, householdIdx); + hostPersonData[personIdx].householdIdx = householdIdx; + // person.setVariable(person::v::HOUSEHOLD_SIZE, household.size); + hostPersonData[personIdx].householdSize = household.size; // initialise the agents infection count if (infectionStatus == disease::SEIR::Infected) { - person.setVariable(exateppabm::person::v::INFECTION_COUNT, 1u); + // person.setVariable(exateppabm::person::v::INFECTION_COUNT, 1u); + hostPersonData[personIdx].infectionCount = 1u; // Increment the per-age demographic initial agent count. @todo refactor elsewhere? _infectedPerDemographic[age]++; } @@ -121,25 +167,25 @@ std::unique_ptr generate(flamegpu::ModelDescription& mode // Assign the workplace based on age band WorkplaceUnderlyingType workplaceIdx = generateWorkplaceForIndividual(age, p_adultPerWorkNetwork, rng); - // Store the agent's expected ID (getID returns 0 when not in an init func?) - // @todo - switch to init func and use getID(); - workplaceMembers[workplaceIdx].push_back(person.getIndex() + 1); + // 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. // @todo - optionally allow binomial distributions // @todo - decide if non non-binomial should be a mean or not, maybe allow non fixed normal dists? // @todo - refactor and test. - double meanRandomInteractions = params.mean_random_interactions_20_69; - double sdRandomInteractions = params.sd_random_interactions_20_69; + double meanRandomInteractions = _params.mean_random_interactions_20_69; + double sdRandomInteractions = _params.sd_random_interactions_20_69; if (age == demographics::Age::AGE_0_9 || age == demographics::Age::AGE_10_19) { - meanRandomInteractions = params.mean_random_interactions_0_19; - sdRandomInteractions = params.sd_random_interactions_0_19; + meanRandomInteractions = _params.mean_random_interactions_0_19; + sdRandomInteractions = _params.sd_random_interactions_0_19; } else if (age == demographics::Age::AGE_70_79 || age == demographics::Age::AGE_80) { - meanRandomInteractions = params.mean_random_interactions_70plus; - sdRandomInteractions = params.sd_random_interactions_70plus; + meanRandomInteractions = _params.mean_random_interactions_70plus; + sdRandomInteractions = _params.sd_random_interactions_70plus; } // Sample a normal distribution (of integers, so we can clamp to >= 0) @@ -147,7 +193,7 @@ std::unique_ptr generate(flamegpu::ModelDescription& mode // Sample from the distribution double randomInteractionsRaw = randomInteractionDist(rng); // Clamp to be between 0 and the population size, and cast to uint - std::uint32_t randomInteractionTarget = static_cast(std::clamp(randomInteractionsRaw, 0.0, static_cast(params.n_total))); + std::uint32_t randomInteractionTarget = static_cast(std::clamp(randomInteractionsRaw, 0.0, static_cast(_params.n_total))); // If the max was over the compile time upper limit due to flamegpu limitations, emit a warning and exit. if (randomInteractionTarget > person::MAX_RANDOM_DAILY_INTERACTIONS) { @@ -159,7 +205,8 @@ std::unique_ptr generate(flamegpu::ModelDescription& mode randomInteractionMin = std::min(randomInteractionMin, randomInteractionTarget); randomInteractionMax = std::max(randomInteractionMax, randomInteractionTarget); // Set for the agent - person.setVariable(person::v::RANDOM_INTERACTION_COUNT_TARGET, randomInteractionTarget); + // person.setVariable(person::v::RANDOM_INTERACTION_COUNT_TARGET, randomInteractionTarget); + hostPersonData[personIdx].randomInteractionTarget = randomInteractionTarget; // Track the sum of target interaction counts randomInteractionCountSum += randomInteractionTarget; @@ -168,7 +215,7 @@ std::unique_ptr generate(flamegpu::ModelDescription& mode } } - if (verbose) { + if (_verbose) { fmt::print("Random Interactions: min={}, max={}, sum={}\n", randomInteractionMin, randomInteractionMax, randomInteractionCountSum); } @@ -182,13 +229,17 @@ std::unique_ptr generate(flamegpu::ModelDescription& mode // Now individuals have been assigned to each workplace network, we can generate the small world networks for workplace interactions std::array meanInteractionsPerNetwork = {{ - params.mean_work_interactions_child / params.daily_fraction_work, - params.mean_work_interactions_child / params.daily_fraction_work, - params.mean_work_interactions_adult / params.daily_fraction_work, - params.mean_work_interactions_elderly / params.daily_fraction_work, - params.mean_work_interactions_elderly / params.daily_fraction_work}}; - - std::array workplaceNetworks = {}; + _params.mean_work_interactions_child / _params.daily_fraction_work, + _params.mean_work_interactions_child / _params.daily_fraction_work, + _params.mean_work_interactions_adult / _params.daily_fraction_work, + _params.mean_work_interactions_elderly / _params.daily_fraction_work, + _params.mean_work_interactions_elderly / _params.daily_fraction_work}}; + + // agent count + 1 elements, as 1 indexed + std::uint32_t totalVertices = _params.n_total; + // Count the total number of edges + std::uint32_t totalUndirectedEdges = 0u; + std::array workplaceNetworks = {}; for (population::WorkplaceUnderlyingType widx = 0; widx < population::WORKPLACE_COUNT; ++widx) { const std::uint64_t N = workplaceMembers[widx].size(); // Round to nearest even number. This might want to be always round up instead @todo @@ -204,44 +255,108 @@ std::unique_ptr generate(flamegpu::ModelDescription& mode // Shuffle the network, to avoid people working with others in their same household a disproportionate amount for low values of work_network_rewire std::shuffle(workplaceMembers[widx].begin(), workplaceMembers[widx].end(), rng); // Generate the small world network - workplaceNetworks[widx] = network::generateSmallWorldUndirectedNetwork(workplaceMembers[widx], meanInteractions, params.work_network_rewire, rng); + workplaceNetworks[widx] = network::generateSmallWorldUndirectedGraph(workplaceMembers[widx], meanInteractions, _params.work_network_rewire, rng); } else { - workplaceNetworks[widx] = network::generateFullyConnectedUndirectedNetwork(workplaceMembers[widx]); + workplaceNetworks[widx] = network::generateFullyConnectedUndirectedGraph(workplaceMembers[widx]); } - // @todo - construct a flame gpu 2 static graph. @todo - check if this can be done from an init function. + totalUndirectedEdges += workplaceNetworks[widx].getNumEdges(); } // Get a handle to the FLAME GPU workplace directed graph + flamegpu::HostEnvironmentDirectedGraph workplaceDigraph = FLAMEGPU->environment.getDirectedGraph("WORKPLACE_DIGRAPH"); // This will contain all workplace information, i.e. a single graph data structure containing 5 unconnected graphs - // The agent's ID will be used as the vertex keys, so will also be 1 indexed. + // The custom 1-indexed agent's ID will be used as the vertex keys // This will then allow agents to simply lookup their neighbours and // FLAME GPU's graphs are accessed by their string name, but agents can't have string properties, hence using a single graph. // Representing a undirected graph using a directed graph will consume twice as much device memory as required, but only digraphs are currently implemented in FLAME GPU 2. - + std::uint32_t totalDirectedEdges = 2 * totalUndirectedEdges; + workplaceDigraph.setVertexCount(totalVertices); + workplaceDigraph.setEdgeCount(totalDirectedEdges); + + // Set vertex properties + flamegpu::HostEnvironmentDirectedGraph::VertexMap vertices = workplaceDigraph.vertices(); + for (std::uint32_t v = 1; v <= totalVertices; ++v) { + flamegpu::HostEnvironmentDirectedGraph::VertexMap::Vertex vertex = vertices[v]; + // vertex.setProperty(person::v::ID, v); + } - // With the small world networks generated, associate them with the model for use in workplace interactions. - std::uint32_t idx = 0; - for (auto person : *pop) { + // Iterate each workplace network, adding it's both directions of each undirected edge to the graph + flamegpu::HostEnvironmentDirectedGraph::EdgeMap edges = workplaceDigraph.edges(); + for (auto& undirectedNetwork : workplaceNetworks) { + // For each undirected edge which contains indexes within the network, create the 2 directed edges between agent IDs + for (const auto & undirectedEdge : undirectedNetwork.getEdges()) { + flamegpu::id_t a = undirectedNetwork.getVertexLabel(undirectedEdge.source); + flamegpu::id_t b = undirectedNetwork.getVertexLabel(undirectedEdge.dest); + auto abEdge = edges[{a, b}]; + abEdge.setSourceDestinationVertexID(a, b); + auto baEdge = edges[{b, a}]; + baEdge.setSourceDestinationVertexID(b, a); + } + } + + // flamegpu::DeviceAgentVector population = personAgent.getPopulationData(); + // Update the population with additional data. + // This is (potentially) suboptimal in terms of host-device memcpy, but need to pass agents multiple times unfortunately. + + // flamegpu::DeviceAgentVector pop = FLAMEGPU->agent("person", "default").getPopulationData(); + // std::uint32_t personIdx = 0; + // for (auto person : pop) { + for (std::uint32_t personIdx = 0; personIdx < hostPersonData.size(); ++personIdx) { + // Get the agents id + // std::uint32_t agentId = person.getVariable(person::v::ID); + std::uint32_t agentId = hostPersonData[personIdx].id; // get the assigned workplace - std::uint32_t workplaceIdx = person.getVariable(person::v::WORKPLACE_IDX); - // Lookup the generated size - std::uint32_t workplaceSize = static_cast(workplaceMembers[workplaceIdx].size()); - // @todo - add the small world network access method to the model - // Store in the agent's data, for faster lookup. @todo - refactor to a env property if maximum workplace count is definitely known at model definition time? - // @todo - replace with the individuals in and/or out degree? - person.setVariable(person::v::WORKPLACE_SIZE, workplaceSize); - ++idx; + // std::uint32_t workplaceIdx = person.getVariable(person::v::WORKPLACE_IDX); + std::uint32_t workplaceIdx = hostPersonData[personIdx].workplaceIdx; + // For this individual in their small world network, get the (out)degree, i.e. the max number of interactions per day. + network::UndirectedGraph& workplaceGraph = workplaceNetworks[workplaceIdx]; + auto vertexLabels = workplaceGraph.getVertexLabels(); + auto it = std::find(vertexLabels.begin(), vertexLabels.end(), agentId); + if (it != vertexLabels.end()) { + std::uint32_t vertexIdx = std::distance(vertexLabels.begin(), it); + // Get the outdegree for this vertex index + std::uint32_t degree = workplaceGraph.degree(vertexIdx); + // person.setVariable(person::v::WORKPLACE_OUT_DEGREE, degree); + hostPersonData[personIdx].workplaceOutDegree = degree; + } else { + throw std::runtime_error("@todo - could not find agent in workplace"); + } + // ++personIdx; } - // If this is a visualisation enabled build, set their x/y/z + // Finally, in another pass create the actual agent instances, finishing the multi-pass init function (and split init function workaround) + auto personAgent = FLAMEGPU->agent(exateppabm::person::NAME, exateppabm::person::states::DEFAULT); + for (std::uint32_t personIdx = 0; personIdx < hostPersonData.size(); ++personIdx) { + const auto& hostPerson = hostPersonData[personIdx]; + // Generate the new agent + auto person = personAgent.newAgent(); + // Set each property on the agent. + person.setVariable(person::v::ID, hostPerson.id); + person.setVariable(exateppabm::person::v::INFECTION_STATE, hostPerson.infectionStatus); + person.setVariable(exateppabm::person::v::INFECTION_STATE_DURATION, hostPerson.infectionStateDuration); + person.setVariable(exateppabm::person::v::INFECTION_COUNT, hostPerson.infectionCount); + person.setVariable(person::v::AGE_DEMOGRAPHIC, hostPerson.ageDemographic); + person.setVariable(person::v::HOUSEHOLD_IDX, hostPerson.householdIdx); + person.setVariable(person::v::HOUSEHOLD_SIZE, hostPerson.householdSize); + person.setVariable(person::v::WORKPLACE_IDX, hostPerson.workplaceIdx); + person.setVariable(person::v::WORKPLACE_OUT_DEGREE, hostPerson.workplaceOutDegree); + person.setVariable(person::v::RANDOM_INTERACTION_COUNT_TARGET, hostPerson.randomInteractionTarget); + + + // If this is a visualisation enabled build, set their x/y/z #if defined(FLAMEGPU_VISUALISATION) - exateppabm::visualisation::initialiseAgentPopulation(model, params, pop, static_cast(households.size())); + auto [visX, visY, visZ] = exateppabm::visualisation::getAgentXYZ(static_cast(households.size()), hostPerson.householdIdx, 0); + person.setVariable(exateppabm::person::v::x, visX); + person.setVariable(exateppabm::person::v::y, visY); + person.setVariable(exateppabm::person::v::z, visZ); #endif // defined(FLAMEGPU_VISUALISATION) + } - if (verbose) { + // Do the verbose output + if (_verbose) { // Print a summary of population creation for now. - fmt::print("Created {} people with {} infected.\n", params.n_total, params.n_seed_infection); + fmt::print("Created {} people with {} infected.\n", _params.n_total, _params.n_seed_infection); fmt::print("Households: {}\n", households.size()); fmt::print("Demographics {{\n"); fmt::print(" 0- 9 = {}\n", createdPerDemographic[0]); @@ -262,8 +377,14 @@ std::unique_ptr generate(flamegpu::ModelDescription& mode fmt::print(" 20_69 {}, 80+ {}\n", workplaceMembers[Workplace::WORKPLACE_80_PLUS].size() - - createdPerDemographic[demographics::Age::AGE_80], createdPerDemographic[demographics::Age::AGE_80]); fmt::print("}}\n"); } +} - return pop; +void define(flamegpu::ModelDescription& model, const exateppabm::input::config params, const bool verbose) { + // Store passed in parameters in file-scoped variables + _params = params; + _verbose = verbose; + // Define the init function which will generate the population for the parameters struct stored in the anon namespace + model.addInitFunction(generatePopulation); } std::array getPerDemographicInitialInfectionCount() { diff --git a/src/exateppabm/population.h b/src/exateppabm/population.h index 5484d96..c1034f8 100644 --- a/src/exateppabm/population.h +++ b/src/exateppabm/population.h @@ -13,13 +13,15 @@ namespace exateppabm { namespace population { /** - * Generate a population of individual person agents for a given simulation configuration + * Define FLAME GPU init function to generate a population of agents for a simulation. + * + * This method additionally stores several values in an anonymous namespace for use during the init function (which cannot take additional parameters) * * @param model the model the population is to be associated with (to get information about the person agent type from) * @param params the model parameters struct for this simulation * @param verbose if verbose output is enabled */ -std::unique_ptr generate(flamegpu::ModelDescription& model, const exateppabm::input::config params, const bool verbose); +void define(flamegpu::ModelDescription& model, const exateppabm::input::config params, const bool verbose); /** * Get the number of agents per demographic which were initialised to be infected, for the most recent call to generate. diff --git a/src/exateppabm/visualisation.cu b/src/exateppabm/visualisation.cu index 643fbfd..876c5b8 100644 --- a/src/exateppabm/visualisation.cu +++ b/src/exateppabm/visualisation.cu @@ -3,6 +3,7 @@ #include #include +#include #include #include "flamegpu/flamegpu.h" #include "exateppabm/person.h" @@ -96,32 +97,28 @@ void join() { #endif // FLAMEGPU_VISUALISATION } -void initialiseAgentPopulation(const flamegpu::ModelDescription& model, const exateppabm::input::config config, std::unique_ptr & pop, const std::uint32_t householdCount) { +std::tuple getAgentXYZ(const std::uint32_t householdCount, const std::uint32_t householdIdx, const std::uint8_t idxWithinHousehold) { #if defined(FLAMEGPU_VISUALISATION) // Use the number of households to figure out the size of a 2D grid for visualisation purposes std::uint64_t visHouseholdGridwidth = static_cast(std::ceil(std::sqrt(static_cast(householdCount)))); // Prep a vector of integers to find the location within a household for each individual - auto visAssignedHouseholdCount = std::vector(householdCount, 0); + static auto visAssignedHouseholdCount = std::vector(householdCount, 0); // pre-calculate spatial offset per individual within household, for upto 6 individuals per household (current hardcoded upper limit) constexpr float OFFSET_SF = 0.7f; - std::array visHouseholdOffsetX = {{0.f + static std::array visHouseholdOffsetX = {{0.f , OFFSET_SF * static_cast(std::sin(0 * M_PI / 180.0)) , OFFSET_SF * static_cast(std::sin(72 * M_PI / 180.0)) , OFFSET_SF * static_cast(std::sin(144 * M_PI / 180.0)) , OFFSET_SF * static_cast(std::sin(216 * M_PI / 180.0)) , OFFSET_SF * static_cast(std::sin(288 * M_PI / 180.0))}}; - std::array visHouseholdOffsetY = {{0.f + static std::array visHouseholdOffsetY = {{0.f , OFFSET_SF * static_cast(std::cos(0 * M_PI / 180.0)) , OFFSET_SF * static_cast(std::cos(72 * M_PI / 180.0)) , OFFSET_SF * static_cast(std::cos(144 * M_PI / 180.0)) , OFFSET_SF * static_cast(std::cos(216 * M_PI / 180.0)) , OFFSET_SF * static_cast(std::cos(288 * M_PI / 180.0))}}; - // Iterate the population, setting the agent's x, y and z values based on their household and index within their household - unsigned idx = 0; - for (auto person : *pop) { - // Get the agent's household index - auto householdIdx = person.getVariable(person::v::HOUSEHOLD_IDX); + // Get teh x/y/z for the provided individual within their household // Get the center point of their household std::uint32_t visHouseholdRow = householdIdx / visHouseholdGridwidth; std::uint32_t visHouseholdCol = householdIdx % visHouseholdGridwidth; @@ -134,12 +131,12 @@ void initialiseAgentPopulation(const flamegpu::ModelDescription& model, const ex float visY = (visHouseholdRow * VIS_HOUSE_GRID_SPACING) + visHouseholdOffsetY[idxInHouse % visHouseholdOffsetY.size()]; float visZ = 0; // Set the x,y and z in agent data. These must be floats. - person.setVariable(exateppabm::person::v::x, visX); - person.setVariable(exateppabm::person::v::y, visY); - person.setVariable(exateppabm::person::v::z, visZ); - // Increment the agent index - ++idx; - } + // person.setVariable(exateppabm::person::v::x, visX); + // person.setVariable(exateppabm::person::v::y, visY); + // person.setVariable(exateppabm::person::v::z, visZ); + return std::tuple(visX, visY, visZ); +#else + return std::tuple(0.0, 0.0, 0.0); #endif // FLAMEGPU_VISUALISATION } diff --git a/src/exateppabm/visualisation.h b/src/exateppabm/visualisation.h index 0c13d7f..3e23a88 100644 --- a/src/exateppabm/visualisation.h +++ b/src/exateppabm/visualisation.h @@ -1,6 +1,7 @@ #pragma once #include +#include #include "flamegpu/flamegpu.h" #include "exateppabm/input.h" @@ -25,18 +26,16 @@ void setup(bool enable, flamegpu::ModelDescription& model, flamegpu::CUDASimulat void join(); /** - * Set visualisation specific agent variables for a given population + * Get visualistion properties for an agent, given model parameters and thier household index * - * I.e. set the x/y/z based on the their householdIdx + * This has had to be refactored to workaround a flame gpu limitation * - * This is more expensive as a separate loop, but only a one time cost and visualsiation runs are not as sensitive to performance - * - * @param model flame gpu model description object - * @param config simulation configuration parameters - * @param pop FLAME GPU person population object, with partially initialised population data * @param householdCount the number of generated households + * @param householdIdx, the index of the current household + * @param idxWithinHousehold the index within the household for the agent. + * @return tuple of agent properties */ -void initialiseAgentPopulation(const flamegpu::ModelDescription& model, const exateppabm::input::config config, std::unique_ptr & pop, const std::uint32_t householdCount); +std::tuple getAgentXYZ(const std::uint32_t householdCount, const std::uint32_t householdIdx, const std::uint8_t idxWithinHousehold); } // namespace visualisation } // namespace exateppabm diff --git a/tests/src/exateppabm/test_network.cu b/tests/src/exateppabm/test_network.cu index 072154a..2d4e1b1 100644 --- a/tests/src/exateppabm/test_network.cu +++ b/tests/src/exateppabm/test_network.cu @@ -15,18 +15,18 @@ * * @todo test more edge cases */ -TEST(TestNetwork, generateFullyConnectedUndirectedNetwork) { +TEST(TestNetwork, generateFullyConnectedUndirectedGraph) { // Generate nodes from [1, 4] constexpr std::uint32_t N = 4; std::vector nodes(N); std::iota(nodes.begin(), nodes.end(), 1); // Generate a fully connected network - exateppabm::network::UndirectedNetwork network = exateppabm::network::generateFullyConnectedUndirectedNetwork(nodes); + exateppabm::network::UndirectedGraph network = exateppabm::network::generateFullyConnectedUndirectedGraph(nodes); // Check the number of nodes and edges are correct - EXPECT_EQ(network.nodes.size(), N); - EXPECT_EQ(network.nodes.size(), nodes.size()); - EXPECT_EQ(network.edges.size(), (N * (N-1)) / 2); + EXPECT_EQ(network.getNumVertices(), N); + EXPECT_EQ(network.getNumVertices(), nodes.size()); + EXPECT_EQ(network.getNumEdges(), (N * (N-1)) / 2); // Check the edges are the expected values. } @@ -35,7 +35,7 @@ TEST(TestNetwork, generateFullyConnectedUndirectedNetwork) { * * @todo - edge case testing. */ -TEST(TestNetwork, generateSmallWorldUndirectedNetwork) { +TEST(TestNetwork, generateSmallWorldUndirectedGraph) { // Generate nodes from [1, 16] constexpr std::uint32_t N = 16; std::vector nodes(N); @@ -45,11 +45,11 @@ TEST(TestNetwork, generateSmallWorldUndirectedNetwork) { constexpr std::uint32_t K = 4; // degree 4 constexpr double p_rewire = 0.1; // probability of a rewire std::mt19937_64 rng(12); // seeded mersenne twister rng engine - exateppabm::network::UndirectedNetwork network = exateppabm::network::generateSmallWorldUndirectedNetwork(nodes, K, p_rewire, rng); + exateppabm::network::UndirectedGraph network = exateppabm::network::generateSmallWorldUndirectedGraph(nodes, K, p_rewire, rng); // Ensure there are the correct number of nodes and edges still - EXPECT_EQ(network.nodes.size(), nodes.size()); - EXPECT_EQ(network.edges.size(), (N * K) / 2); + EXPECT_EQ(network.getNumVertices(), nodes.size()); + EXPECT_EQ(network.getNumEdges(), (N * K) / 2); // @todo - validate the in degree and out degree of each node. // If graph is directed, out degree will still be K, but indegree will vary per node // If graph is undirected, out degree and in degree would be the same, but mean degree will be K. @@ -57,10 +57,10 @@ TEST(TestNetwork, generateSmallWorldUndirectedNetwork) { // @todo - actually test the network properties. Average degree etc. size_t edgeIdx = 0; - for (const auto &edge : network.edges) { + for (const auto &edge : network.getEdges()) { // Every source and ever dest should be a valid node index - EXPECT_LT(edge.source, network.nodes.size()); - EXPECT_LT(edge.dest, network.nodes.size()); + EXPECT_LT(edge.source, network.getNumVertices()); + EXPECT_LT(edge.dest, network.getNumVertices()); // No self-edges EXPECT_NE(edge.source, edge.dest); @@ -71,82 +71,82 @@ TEST(TestNetwork, generateSmallWorldUndirectedNetwork) { // Check invalid waltz strogatz parameters behave as intended // Odd degree - actually use K-1, - auto n_k_odd = exateppabm::network::generateSmallWorldUndirectedNetwork({{1, 2, 3, 4}}, 3, 0.1, rng); - EXPECT_EQ(n_k_odd.nodes.size(), 4); - EXPECT_EQ(n_k_odd.edges.size(), ((3-1) * n_k_odd.nodes.size())/2); + auto n_k_odd = exateppabm::network::generateSmallWorldUndirectedGraph({{1, 2, 3, 4}}, 3, 0.1, rng); + EXPECT_EQ(n_k_odd.getNumVertices(), 4); + EXPECT_EQ(n_k_odd.getNumEdges(), ((3-1) * n_k_odd.getNumVertices())/2); // too high of a degree - EXPECT_ANY_THROW(exateppabm::network::generateSmallWorldUndirectedNetwork({{1, 2, 3, 4}}, 5, 0.1, rng)); + EXPECT_ANY_THROW(exateppabm::network::generateSmallWorldUndirectedGraph({{1, 2, 3, 4}}, 5, 0.1, rng)); // too low probability - EXPECT_ANY_THROW(exateppabm::network::generateSmallWorldUndirectedNetwork({{1, 2, 3, 4}}, 2, -1.0, rng)); + EXPECT_ANY_THROW(exateppabm::network::generateSmallWorldUndirectedGraph({{1, 2, 3, 4}}, 2, -1.0, rng)); // too high probability - EXPECT_ANY_THROW(exateppabm::network::generateSmallWorldUndirectedNetwork({{1, 2, 3, 4}}, 2, 2.0, rng)); + EXPECT_ANY_THROW(exateppabm::network::generateSmallWorldUndirectedGraph({{1, 2, 3, 4}}, 2, 2.0, rng)); // Check the networks contain the expected number of edges for edge case node counts // 0 nodes, no edges - auto n0 = exateppabm::network::generateSmallWorldUndirectedNetwork({}, 2, 0.1, rng); - EXPECT_EQ(n0.nodes.size(), 0u); - EXPECT_EQ(n0.edges.size(), 0u); + auto n0 = exateppabm::network::generateSmallWorldUndirectedGraph({}, 2, 0.1, rng); + EXPECT_EQ(n0.getNumVertices(), 0u); + EXPECT_EQ(n0.getNumEdges(), 0u); // 1 node, no edges - auto n1 = exateppabm::network::generateSmallWorldUndirectedNetwork({{1}}, 2, 0.1, rng); - EXPECT_EQ(n1.nodes.size(), 1u); - EXPECT_EQ(n1.edges.size(), 0u); + auto n1 = exateppabm::network::generateSmallWorldUndirectedGraph({{1}}, 2, 0.1, rng); + EXPECT_EQ(n1.getNumVertices(), 1u); + EXPECT_EQ(n1.getNumEdges(), 0u); // 2 node, fully connected - auto n2 = exateppabm::network::generateSmallWorldUndirectedNetwork({{1, 2}}, 2, 0.1, rng); - EXPECT_EQ(n2.nodes.size(), 2u); - EXPECT_EQ(n2.edges.size(), 1u); + auto n2 = exateppabm::network::generateSmallWorldUndirectedGraph({{1, 2}}, 2, 0.1, rng); + EXPECT_EQ(n2.getNumVertices(), 2u); + EXPECT_EQ(n2.getNumEdges(), 1u); // 3 nodes, mean degree 2, 3 edges - auto n3 = exateppabm::network::generateSmallWorldUndirectedNetwork({{1, 2, 3}}, 2, 0.1, rng); - EXPECT_EQ(n3.nodes.size(), 3u); - EXPECT_EQ(n3.edges.size(), 3u); + auto n3 = exateppabm::network::generateSmallWorldUndirectedGraph({{1, 2, 3}}, 2, 0.1, rng); + EXPECT_EQ(n3.getNumVertices(), 3u); + EXPECT_EQ(n3.getNumEdges(), 3u); // 4 nodes, degree 2, 4 edges - auto n4_2 = exateppabm::network::generateSmallWorldUndirectedNetwork({{1, 2, 3, 4}}, 2, 0.1, rng); - EXPECT_EQ(n4_2.nodes.size(), 4u); - EXPECT_EQ(n4_2.edges.size(), 4u); + auto n4_2 = exateppabm::network::generateSmallWorldUndirectedGraph({{1, 2, 3, 4}}, 2, 0.1, rng); + EXPECT_EQ(n4_2.getNumVertices(), 4u); + EXPECT_EQ(n4_2.getNumEdges(), 4u); // 4 nodes, degree 4, fully connected, so only 6 edges - auto n4_4 = exateppabm::network::generateSmallWorldUndirectedNetwork({{1, 2, 3, 4}}, 4, 0.1, rng); - EXPECT_EQ(n4_4.nodes.size(), 4u); - EXPECT_EQ(n4_4.edges.size(), 6u); + auto n4_4 = exateppabm::network::generateSmallWorldUndirectedGraph({{1, 2, 3, 4}}, 4, 0.1, rng); + EXPECT_EQ(n4_4.getNumVertices(), 4u); + EXPECT_EQ(n4_4.getNumEdges(), 6u); // 12 nodes, degree 2 - auto n12_2 = exateppabm::network::generateSmallWorldUndirectedNetwork({{1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12}}, 2, 0.1, rng); - EXPECT_EQ(n12_2.nodes.size(), 12u); - EXPECT_EQ(n12_2.edges.size(), 12u); + auto n12_2 = exateppabm::network::generateSmallWorldUndirectedGraph({{1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12}}, 2, 0.1, rng); + EXPECT_EQ(n12_2.getNumVertices(), 12u); + EXPECT_EQ(n12_2.getNumEdges(), 12u); // 12 nodes, degree 2 (equiv to 2) - auto n12_3 = exateppabm::network::generateSmallWorldUndirectedNetwork({{1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12}}, 3, 0.1, rng); - EXPECT_EQ(n12_3.nodes.size(), 12u); - EXPECT_EQ(n12_3.edges.size(), 12u); + auto n12_3 = exateppabm::network::generateSmallWorldUndirectedGraph({{1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12}}, 3, 0.1, rng); + EXPECT_EQ(n12_3.getNumVertices(), 12u); + EXPECT_EQ(n12_3.getNumEdges(), 12u); // 12 nodes, degree 4 - auto n12_4 = exateppabm::network::generateSmallWorldUndirectedNetwork({{1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12}}, 4, 0.1, rng); - EXPECT_EQ(n12_4.nodes.size(), 12u); - EXPECT_EQ(n12_4.edges.size(), 24u); + auto n12_4 = exateppabm::network::generateSmallWorldUndirectedGraph({{1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12}}, 4, 0.1, rng); + EXPECT_EQ(n12_4.getNumVertices(), 12u); + EXPECT_EQ(n12_4.getNumEdges(), 24u); // 12 nodes, degree 5 (equiv to 5) - auto n12_5 = exateppabm::network::generateSmallWorldUndirectedNetwork({{1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12}}, 5, 0.1, rng); - EXPECT_EQ(n12_5.nodes.size(), 12u); - EXPECT_EQ(n12_5.edges.size(), 24u); + auto n12_5 = exateppabm::network::generateSmallWorldUndirectedGraph({{1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12}}, 5, 0.1, rng); + EXPECT_EQ(n12_5.getNumVertices(), 12u); + EXPECT_EQ(n12_5.getNumEdges(), 24u); // 12 nodes, degree 6 - auto n12_6 = exateppabm::network::generateSmallWorldUndirectedNetwork({{1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12}}, 6, 0.1, rng); - EXPECT_EQ(n12_6.nodes.size(), 12u); - EXPECT_EQ(n12_6.edges.size(), 36u); + auto n12_6 = exateppabm::network::generateSmallWorldUndirectedGraph({{1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12}}, 6, 0.1, rng); + EXPECT_EQ(n12_6.getNumVertices(), 12u); + EXPECT_EQ(n12_6.getNumEdges(), 36u); // 12 nodes, degree 8 - auto n12_8 = exateppabm::network::generateSmallWorldUndirectedNetwork({{1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12}}, 8, 0.1, rng); - EXPECT_EQ(n12_8.nodes.size(), 12u); - EXPECT_EQ(n12_8.edges.size(), 48u); + auto n12_8 = exateppabm::network::generateSmallWorldUndirectedGraph({{1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12}}, 8, 0.1, rng); + EXPECT_EQ(n12_8.getNumVertices(), 12u); + EXPECT_EQ(n12_8.getNumEdges(), 48u); // 12 nodes, degree 8 - auto n12_10 = exateppabm::network::generateSmallWorldUndirectedNetwork({{1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12}}, 10, 0.1, rng); - EXPECT_EQ(n12_10.nodes.size(), 12u); - EXPECT_EQ(n12_10.edges.size(), 60u); + auto n12_10 = exateppabm::network::generateSmallWorldUndirectedGraph({{1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12}}, 10, 0.1, rng); + EXPECT_EQ(n12_10.getNumVertices(), 12u); + EXPECT_EQ(n12_10.getNumEdges(), 60u); // 12 nodes, degree 8 - auto n12_12 = exateppabm::network::generateSmallWorldUndirectedNetwork({{1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12}}, 12, 0.1, rng); - EXPECT_EQ(n12_12.nodes.size(), 12u); - EXPECT_EQ(n12_12.edges.size(), 66u); // fully connected + auto n12_12 = exateppabm::network::generateSmallWorldUndirectedGraph({{1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12}}, 12, 0.1, rng); + EXPECT_EQ(n12_12.getNumVertices(), 12u); + EXPECT_EQ(n12_12.getNumEdges(), 66u); // fully connected // 8 nodes, degree 4, no rewiring - auto n8_4_0 = exateppabm::network::generateSmallWorldUndirectedNetwork({{1, 2, 3, 4, 5, 6, 7, 8}}, 4, 0.0, rng); - EXPECT_EQ(n8_4_0.nodes.size(), 8u); - EXPECT_EQ(n8_4_0.edges.size(), 16u); + auto n8_4_0 = exateppabm::network::generateSmallWorldUndirectedGraph({{1, 2, 3, 4, 5, 6, 7, 8}}, 4, 0.0, rng); + EXPECT_EQ(n8_4_0.getNumVertices(), 8u); + EXPECT_EQ(n8_4_0.getNumEdges(), 16u); // Check for valid edge, ascending internal order, ensure no dupes. - std::set uniqueEdges_n8_4_0 = {{}}; - for (const auto & edge : n8_4_0.edges) { + std::set uniqueEdges_n8_4_0 = {{}}; + for (const auto & edge : n8_4_0.getEdges()) { EXPECT_GE(edge.source, 0u); EXPECT_LT(edge.source, 8u); EXPECT_NE(edge.source, edge.dest); @@ -156,6 +156,6 @@ TEST(TestNetwork, generateSmallWorldUndirectedNetwork) { EXPECT_EQ(uniqueEdges_n8_4_0.count(edge), 0); uniqueEdges_n8_4_0.insert(edge); // with rewiring 0, the dest should always be source +1 or source +2 or when wrapping the boundary source will be dest -1 or dest -2 wrapped + 2) - EXPECT_TRUE((edge.dest == edge.source + 1) || (edge.dest == edge.source + 2) || (edge.source == (edge.dest + 1) % n8_4_0.nodes.size()) || (edge.source == (edge.dest + 2) % n8_4_0.nodes.size())); + EXPECT_TRUE((edge.dest == edge.source + 1) || (edge.dest == edge.source + 2) || (edge.source == (edge.dest + 1) % n8_4_0.getNumVertices()) || (edge.source == (edge.dest + 2) % n8_4_0.getNumVertices())); } } diff --git a/tools/plot-per-individual.py b/tools/plot-per-individual.py index 9aa5a54..5c9dedb 100755 --- a/tools/plot-per-individual.py +++ b/tools/plot-per-individual.py @@ -36,9 +36,11 @@ def read_inputs(args): # Add string version of occupation network group OCCUPATION_NETWORK_MAP = { - 0: "Children", - 1: "Adult", - 2: "Elderly", + 0: "0-9", + 1: "10-19", + 2: "20-69", + 3: "70-79", + 4: "80+" } df["occupation_network_str"] = pd.Categorical(df["occupation_network"].map(OCCUPATION_NETWORK_MAP), list(OCCUPATION_NETWORK_MAP.values())) return df