From b04c5fcc5a393504c83bb1a4f06e2e81a2a288bd Mon Sep 17 00:00:00 2001 From: Amrita Goswami Date: Sat, 23 Mar 2024 13:48:30 +0000 Subject: [PATCH 1/9] InertialModel: Created InertialAgent The InertialAgent has a velocity also. Co-authored-by: Moritz Sallermann --- include/agents/inertial_agent.hpp | 68 +++++++++++++++++++++++++++++++ 1 file changed, 68 insertions(+) create mode 100644 include/agents/inertial_agent.hpp diff --git a/include/agents/inertial_agent.hpp b/include/agents/inertial_agent.hpp new file mode 100644 index 0000000..469376a --- /dev/null +++ b/include/agents/inertial_agent.hpp @@ -0,0 +1,68 @@ +#pragma once + +#include "agent.hpp" +#include "agent_io.hpp" +#include + +namespace Seldon +{ + +struct InertialAgentData +{ + double opinion = 0; // x_i + double activity = 0; // alpha_i + double reluctance = 1.0; // m_i + double velocity = 0.0; // d(x_i)/dt +}; + +using InertialAgent = Agent; + +template<> +inline std::string agent_to_string( const InertialAgent & agent ) +{ + return fmt::format( + "{}, {}, {}, {}", agent.data.opinion, agent.data.velocity, agent.data.activity, agent.data.reluctance ); +} + +template<> +inline std::string opinion_to_string( const InertialAgent & agent ) +{ + return fmt::format( "{}", agent.data.opinion ); +} + +template<> +inline InertialAgent agent_from_string( const std::string & str ) +{ + InertialAgent res{}; + + auto callback = [&]( int idx_list, std::string & substr ) + { + if( idx_list == 0 ) + { + res.data.opinion = std::stod( substr ); + } + else if( idx_list == 1 ) + { + res.data.velocity = std::stod( substr ); + } + else if( idx_list == 2 ) + { + res.data.activity = std::stod( substr ); + } + else if( idx_list == 3 ) + { + res.data.reluctance = std::stod( substr ); + } + }; + + Seldon::parse_comma_separated_list( str, callback ); + + return res; +}; + +template<> +inline std::vector agent_to_string_column_names() +{ + return { "opinion", "velocity", "activity", "reluctance" }; +} +} // namespace Seldon \ No newline at end of file From d7bbd4003b1ddfb73bf9008eeb242e11f737f6b1 Mon Sep 17 00:00:00 2001 From: Amrita Goswami Date: Sat, 23 Mar 2024 14:29:47 +0000 Subject: [PATCH 2/9] InertialModel: WIP: using template specialization We need to implement a velocity Verlet algorithm for iteration(), in the inertial model. However, many functions are not needed in the ActivityDrivenModel. Therefore, we are considering inheritance of the abstract class instead Co-authored-by: Moritz Sallermann --- include/models/ActivityDrivenModel.hpp | 262 ++++++++++++++++++++++++- src/models/ActivityDrivenModel.cpp | 251 +---------------------- src/models/InertialModel.cpp | 50 +++++ 3 files changed, 306 insertions(+), 257 deletions(-) create mode 100644 src/models/InertialModel.cpp diff --git a/include/models/ActivityDrivenModel.hpp b/include/models/ActivityDrivenModel.hpp index 06d830f..7801b6f 100644 --- a/include/models/ActivityDrivenModel.hpp +++ b/include/models/ActivityDrivenModel.hpp @@ -1,10 +1,11 @@ #pragma once #include "agents/activity_agent.hpp" +#include "agents/inertial_agent.hpp" #include "config_parser.hpp" - #include "model.hpp" #include "network.hpp" +#include "network_generation.hpp" #include #include #include @@ -15,19 +16,55 @@ namespace Seldon { -class ActivityDrivenModel : public Model +template +class ActivityDrivenModelAbstract : public Model { public: - using AgentT = ActivityAgent; + using AgentT = AgentT_; using NetworkT = Network; + using WeightT = typename NetworkT::WeightT; + + ActivityDrivenModelAbstract( + const Config::ActivityDrivenSettings & settings, NetworkT & network, std::mt19937 & gen ) + : Model( settings.max_iterations ), + network( network ), + contact_prob_list( std::vector>( network.n_agents() ) ), + gen( gen ), + dt( settings.dt ), + m( settings.m ), + eps( settings.eps ), + gamma( settings.gamma ), + alpha( settings.alpha ), + homophily( settings.homophily ), + reciprocity( settings.reciprocity ), + K( settings.K ), + mean_activities( settings.mean_activities ), + mean_weights( settings.mean_weights ), + use_reluctances( settings.use_reluctances ), + reluctance_mean( settings.reluctance_mean ), + reluctance_sigma( settings.reluctance_sigma ), + reluctance_eps( settings.reluctance_eps ), + n_bots( settings.n_bots ), + bot_m( settings.bot_m ), + bot_activity( settings.bot_activity ), + bot_opinion( settings.bot_opinion ), + bot_homophily( settings.bot_homophily ) + { + get_agents_from_power_law(); - ActivityDrivenModel( const Config::ActivityDrivenSettings & settings, NetworkT & network, std::mt19937 & gen ); + if( mean_weights ) + { + auto agents_copy = network.agents; + network = NetworkGeneration::generate_fully_connected( network.n_agents() ); + network.agents = agents_copy; + } + } void iteration() override; private: NetworkT & network; - std::vector> contact_prob_list; // Probability of choosing i in 1 to m rounds + std::vector> contact_prob_list; // Probability of choosing i in 1 to m rounds // Random number generation std::mt19937 & gen; // reference to simulation Mersenne-Twister engine std::set> reciprocal_edge_buffer{}; @@ -66,7 +103,44 @@ class ActivityDrivenModel : public Model std::vector k3_buffer{}; std::vector k4_buffer{}; - void get_agents_from_power_law(); + void get_agents_from_power_law() + { + std::uniform_real_distribution<> dis_opinion( -1, 1 ); // Opinion initial values + power_law_distribution<> dist_activity( eps, gamma ); + truncated_normal_distribution<> dist_reluctance( reluctance_mean, reluctance_sigma, reluctance_eps ); + + auto mean_activity = dist_activity.mean(); + + // Initial conditions for the opinions, initialize to [-1,1] + // The activities should be drawn from a power law distribution + for( size_t i = 0; i < network.agents.size(); i++ ) + { + network.agents[i].data.opinion = dis_opinion( gen ); // Draw the opinion value + + if( !mean_activities ) + { // Draw from a power law distribution (1-gamma)/(1-eps^(1-gamma)) * a^(-gamma) + network.agents[i].data.activity = dist_activity( gen ); + } + else + { + network.agents[i].data.activity = mean_activity; + } + + if( use_reluctances ) + { + network.agents[i].data.reluctance = dist_reluctance( gen ); + } + } + + if( bot_present() ) + { + for( size_t bot_idx = 0; bot_idx < n_bots; bot_idx++ ) + { + network.agents[bot_idx].data.opinion = bot_opinion[bot_idx]; + network.agents[bot_idx].data.activity = bot_activity[bot_idx]; + } + } + } [[nodiscard]] bool bot_present() const { @@ -99,10 +173,178 @@ class ActivityDrivenModel : public Model } // The weight for contact between two agents - double homophily_weight( size_t idx_contacter, size_t idx_contacted ); - void update_network_probabilistic(); - void update_network_mean(); - void update_network(); + double homophily_weight( size_t idx_contacter, size_t idx_contacted ) + { + double homophily = this->homophily; + + if( idx_contacted == idx_contacter ) + return 0.0; + + if( bot_present() && idx_contacter < n_bots ) + homophily = this->bot_homophily[idx_contacter]; + + constexpr double tolerance = 1e-10; + auto opinion_diff + = std::abs( network.agents[idx_contacter].data.opinion - network.agents[idx_contacted].data.opinion ); + opinion_diff = std::max( tolerance, opinion_diff ); + + return std::pow( opinion_diff, -homophily ); + } + + void update_network_probabilistic() + { + network.switch_direction_flag(); + + std::uniform_real_distribution<> dis_activation( 0.0, 1.0 ); + std::uniform_real_distribution<> dis_reciprocation( 0.0, 1.0 ); + std::vector contacted_agents{}; + reciprocal_edge_buffer.clear(); // Clear the reciprocal edge buffer + for( size_t idx_agent = 0; idx_agent < network.n_agents(); idx_agent++ ) + { + // Test if the agent is activated + bool activated = dis_activation( gen ) < network.agents[idx_agent].data.activity; + + if( activated ) + { + // Implement the weight for the probability of agent `idx_agent` contacting agent `j` + // Not normalised since this is taken care of by the reservoir sampling + + int m_temp = this->m; + + if( bot_present() && idx_agent < n_bots ) + { + m_temp = bot_m[idx_agent]; + } + + reservoir_sampling_A_ExpJ( + m_temp, network.n_agents(), [&]( int j ) { return homophily_weight( idx_agent, j ); }, + contacted_agents, gen ); + + // Fill the outgoing edges into the reciprocal edge buffer + for( const auto & idx_outgoing : contacted_agents ) + { + reciprocal_edge_buffer.insert( + { idx_agent, idx_outgoing } ); // insert the edge idx_agent -> idx_outgoing + } + + // Set the *outgoing* edges + network.set_neighbours_and_weights( idx_agent, contacted_agents, 1.0 ); + } + else + { + network.set_neighbours_and_weights( idx_agent, {}, {} ); + } + } + + // Reciprocity check + for( size_t idx_agent = 0; idx_agent < network.n_agents(); idx_agent++ ) + { + // Get the outgoing edges + auto contacted_agents = network.get_neighbours( idx_agent ); + // For each outgoing edge we check if the reverse edge already exists + for( const auto & idx_outgoing : contacted_agents ) + { + // If the edge is not reciprocated + if( !reciprocal_edge_buffer.contains( { idx_outgoing, idx_agent } ) ) + { + if( dis_reciprocation( gen ) < reciprocity ) + { + network.push_back_neighbour_and_weight( idx_outgoing, idx_agent, 1.0 ); + } + } + } + } + + network.toggle_incoming_outgoing(); // switch direction, so that we have incoming edges + } + + void update_network_mean() + { + using WeightT = NetworkT::WeightT; + std::vector weights( network.n_agents(), 0.0 ); + + // Set all weights to zero in the beginning + for( size_t idx_agent = 0; idx_agent < network.n_agents(); idx_agent++ ) + { + network.set_weights( idx_agent, weights ); + contact_prob_list[idx_agent] = weights; // set to zero + } + + auto probability_helper = []( double omega, size_t m ) + { + double p = 0; + for( size_t i = 1; i <= m; i++ ) + p += ( std::pow( -omega, i + 1 ) + omega ) / ( omega + 1 ); + return p; + }; + + for( size_t idx_agent = 0; idx_agent < network.n_agents(); idx_agent++ ) + { + // Implement the weight for the probability of agent `idx_agent` contacting agent `j` + // Not normalised since this is taken care of by the reservoir sampling + + double normalization = 0; + for( size_t k = 0; k < network.n_agents(); k++ ) + { + normalization += homophily_weight( idx_agent, k ); + } + + // Go through all the neighbours of idx_agent + // Calculate the probability of i contacting j (in 1 to m rounds, assuming + // the agent is activated + int m_temp = m; + if( bot_present() && idx_agent < n_bots ) + { + m_temp = bot_m[idx_agent]; + } + + double activity = std::max( 1.0, network.agents[idx_agent].data.activity ); + for( size_t j = 0; j < network.n_agents(); j++ ) + { + double omega = homophily_weight( idx_agent, j ) / normalization; + // We can calculate the probability of i contacting j ( i->j ) + // Update contact prob_list (outgoing) + contact_prob_list[idx_agent][j] = activity * probability_helper( omega, m_temp ); + } + } + + for( size_t idx_agent = 0; idx_agent < network.n_agents(); idx_agent++ ) + { + // Calculate the actual weights and reciprocity + for( size_t j = 0; j < network.n_agents(); j++ ) + { + double prob_contact_ij = contact_prob_list[idx_agent][j]; // outgoing probabilites + double prob_contact_ji = contact_prob_list[j][idx_agent]; + + // Set the incoming agent weight, j-i in weight list + double & win_ji = network.get_weights( j )[idx_agent]; + win_ji += prob_contact_ij; + + // Handle the reciprocity for j->i + // Update incoming weight i-j + double & win_ij = network.get_weights( idx_agent )[j]; + + // The probability of reciprocating is + win_ij += ( 1.0 - prob_contact_ji ) * reciprocity * prob_contact_ij; + } + } + } + + void update_network() + { + + if( !mean_weights ) + { + update_network_probabilistic(); + } + else + { + update_network_mean(); + } + } }; +using ActivityDrivenModel = ActivityDrivenModelAbstract; +using InertialModel = ActivityDrivenModelAbstract; + } // namespace Seldon \ No newline at end of file diff --git a/src/models/ActivityDrivenModel.cpp b/src/models/ActivityDrivenModel.cpp index 0a0ba06..bdfc7c4 100644 --- a/src/models/ActivityDrivenModel.cpp +++ b/src/models/ActivityDrivenModel.cpp @@ -1,6 +1,5 @@ #include "models/ActivityDrivenModel.hpp" #include "network.hpp" -#include "network_generation.hpp" #include "util/math.hpp" #include #include @@ -9,252 +8,8 @@ namespace Seldon { -ActivityDrivenModel::ActivityDrivenModel( - const Config::ActivityDrivenSettings & settings, NetworkT & network, std::mt19937 & gen ) - : Model( settings.max_iterations ), - network( network ), - contact_prob_list( std::vector>( network.n_agents() ) ), - gen( gen ), - dt( settings.dt ), - m( settings.m ), - eps( settings.eps ), - gamma( settings.gamma ), - alpha( settings.alpha ), - homophily( settings.homophily ), - reciprocity( settings.reciprocity ), - K( settings.K ), - mean_activities( settings.mean_activities ), - mean_weights( settings.mean_weights ), - use_reluctances( settings.use_reluctances ), - reluctance_mean( settings.reluctance_mean ), - reluctance_sigma( settings.reluctance_sigma ), - reluctance_eps( settings.reluctance_eps ), - n_bots( settings.n_bots ), - bot_m( settings.bot_m ), - bot_activity( settings.bot_activity ), - bot_opinion( settings.bot_opinion ), - bot_homophily( settings.bot_homophily ) -{ - get_agents_from_power_law(); - - if( mean_weights ) - { - auto agents_copy = network.agents; - network = NetworkGeneration::generate_fully_connected( network.n_agents() ); - network.agents = agents_copy; - } -} - -double ActivityDrivenModel::homophily_weight( size_t idx_contacter, size_t idx_contacted ) -{ - double homophily = this->homophily; - - if( idx_contacted == idx_contacter ) - return 0.0; - - if( bot_present() && idx_contacter < n_bots ) - homophily = this->bot_homophily[idx_contacter]; - - constexpr double tolerance = 1e-10; - auto opinion_diff - = std::abs( network.agents[idx_contacter].data.opinion - network.agents[idx_contacted].data.opinion ); - opinion_diff = std::max( tolerance, opinion_diff ); - - return std::pow( opinion_diff, -homophily ); -} - -void ActivityDrivenModel::get_agents_from_power_law() -{ - std::uniform_real_distribution<> dis_opinion( -1, 1 ); // Opinion initial values - power_law_distribution<> dist_activity( eps, gamma ); - truncated_normal_distribution<> dist_reluctance( reluctance_mean, reluctance_sigma, reluctance_eps ); - - auto mean_activity = dist_activity.mean(); - - // Initial conditions for the opinions, initialize to [-1,1] - // The activities should be drawn from a power law distribution - for( size_t i = 0; i < network.agents.size(); i++ ) - { - network.agents[i].data.opinion = dis_opinion( gen ); // Draw the opinion value - - if( !mean_activities ) - { // Draw from a power law distribution (1-gamma)/(1-eps^(1-gamma)) * a^(-gamma) - network.agents[i].data.activity = dist_activity( gen ); - } - else - { - network.agents[i].data.activity = mean_activity; - } - - if( use_reluctances ) - { - network.agents[i].data.reluctance = dist_reluctance( gen ); - } - } - - if( bot_present() ) - { - for( size_t bot_idx = 0; bot_idx < n_bots; bot_idx++ ) - { - network.agents[bot_idx].data.opinion = bot_opinion[bot_idx]; - network.agents[bot_idx].data.activity = bot_activity[bot_idx]; - } - } -} - -void ActivityDrivenModel::update_network_probabilistic() -{ - network.switch_direction_flag(); - - std::uniform_real_distribution<> dis_activation( 0.0, 1.0 ); - std::uniform_real_distribution<> dis_reciprocation( 0.0, 1.0 ); - std::vector contacted_agents{}; - reciprocal_edge_buffer.clear(); // Clear the reciprocal edge buffer - for( size_t idx_agent = 0; idx_agent < network.n_agents(); idx_agent++ ) - { - // Test if the agent is activated - bool activated = dis_activation( gen ) < network.agents[idx_agent].data.activity; - - if( activated ) - { - // Implement the weight for the probability of agent `idx_agent` contacting agent `j` - // Not normalised since this is taken care of by the reservoir sampling - - int m_temp = this->m; - - if( bot_present() && idx_agent < n_bots ) - { - m_temp = bot_m[idx_agent]; - } - - reservoir_sampling_A_ExpJ( - m_temp, network.n_agents(), [&]( int j ) { return homophily_weight( idx_agent, j ); }, contacted_agents, - gen ); - - // Fill the outgoing edges into the reciprocal edge buffer - for( const auto & idx_outgoing : contacted_agents ) - { - reciprocal_edge_buffer.insert( - { idx_agent, idx_outgoing } ); // insert the edge idx_agent -> idx_outgoing - } - - // Set the *outgoing* edges - network.set_neighbours_and_weights( idx_agent, contacted_agents, 1.0 ); - } - else - { - network.set_neighbours_and_weights( idx_agent, {}, {} ); - } - } - - // Reciprocity check - for( size_t idx_agent = 0; idx_agent < network.n_agents(); idx_agent++ ) - { - // Get the outgoing edges - auto contacted_agents = network.get_neighbours( idx_agent ); - // For each outgoing edge we check if the reverse edge already exists - for( const auto & idx_outgoing : contacted_agents ) - { - // If the edge is not reciprocated - if( !reciprocal_edge_buffer.contains( { idx_outgoing, idx_agent } ) ) - { - if( dis_reciprocation( gen ) < reciprocity ) - { - network.push_back_neighbour_and_weight( idx_outgoing, idx_agent, 1.0 ); - } - } - } - } - - network.toggle_incoming_outgoing(); // switch direction, so that we have incoming edges -} - -void ActivityDrivenModel::update_network_mean() -{ - using WeightT = NetworkT::WeightT; - std::vector weights( network.n_agents(), 0.0 ); - - // Set all weights to zero in the beginning - for( size_t idx_agent = 0; idx_agent < network.n_agents(); idx_agent++ ) - { - network.set_weights( idx_agent, weights ); - contact_prob_list[idx_agent] = weights; // set to zero - } - - auto probability_helper = []( double omega, size_t m ) - { - double p = 0; - for( size_t i = 1; i <= m; i++ ) - p += ( std::pow( -omega, i + 1 ) + omega ) / ( omega + 1 ); - return p; - }; - - for( size_t idx_agent = 0; idx_agent < network.n_agents(); idx_agent++ ) - { - // Implement the weight for the probability of agent `idx_agent` contacting agent `j` - // Not normalised since this is taken care of by the reservoir sampling - - double normalization = 0; - for( size_t k = 0; k < network.n_agents(); k++ ) - { - normalization += homophily_weight( idx_agent, k ); - } - - // Go through all the neighbours of idx_agent - // Calculate the probability of i contacting j (in 1 to m rounds, assuming - // the agent is activated - int m_temp = m; - if( bot_present() && idx_agent < n_bots ) - { - m_temp = bot_m[idx_agent]; - } - - double activity = std::max( 1.0, network.agents[idx_agent].data.activity ); - for( size_t j = 0; j < network.n_agents(); j++ ) - { - double omega = homophily_weight( idx_agent, j ) / normalization; - // We can calculate the probability of i contacting j ( i->j ) - // Update contact prob_list (outgoing) - contact_prob_list[idx_agent][j] = activity * probability_helper( omega, m_temp ); - } - } - - for( size_t idx_agent = 0; idx_agent < network.n_agents(); idx_agent++ ) - { - // Calculate the actual weights and reciprocity - for( size_t j = 0; j < network.n_agents(); j++ ) - { - double prob_contact_ij = contact_prob_list[idx_agent][j]; // outgoing probabilites - double prob_contact_ji = contact_prob_list[j][idx_agent]; - - // Set the incoming agent weight, j-i in weight list - double & win_ji = network.get_weights( j )[idx_agent]; - win_ji += prob_contact_ij; - - // Handle the reciprocity for j->i - // Update incoming weight i-j - double & win_ij = network.get_weights( idx_agent )[j]; - - // The probability of reciprocating is - win_ij += ( 1.0 - prob_contact_ji ) * reciprocity * prob_contact_ij; - } - } -} - -void ActivityDrivenModel::update_network() -{ - - if( !mean_weights ) - { - update_network_probabilistic(); - } - else - { - update_network_mean(); - } -} - -void ActivityDrivenModel::iteration() +template<> +void ActivityDrivenModelAbstract::iteration() { Model::iteration(); @@ -289,4 +44,6 @@ void ActivityDrivenModel::iteration() } } } + +template class ActivityDrivenModelAbstract; } // namespace Seldon \ No newline at end of file diff --git a/src/models/InertialModel.cpp b/src/models/InertialModel.cpp new file mode 100644 index 0000000..a17d3f8 --- /dev/null +++ b/src/models/InertialModel.cpp @@ -0,0 +1,50 @@ +#include "agents/inertial_agent.hpp" +#include "models/ActivityDrivenModel.hpp" +#include "network.hpp" +#include "util/math.hpp" +#include +#include +#include + +namespace Seldon +{ + +template<> +void ActivityDrivenModelAbstract::iteration() +{ + Model::iteration(); + + update_network(); + + // Integrate the ODE using 4th order Runge-Kutta + // k_1 = hf(x_n,y_n) + get_euler_slopes( k1_buffer, [this]( size_t i ) { return network.agents[i].data.opinion; } ); + // k_2 = hf(x_n+1/2h,y_n+1/2k_1) + get_euler_slopes( + k2_buffer, [this]( size_t i ) { return network.agents[i].data.opinion + 0.5 * this->k1_buffer[i]; } ); + // k_3 = hf(x_n+1/2h,y_n+1/2k_2) + get_euler_slopes( + k3_buffer, [this]( size_t i ) { return network.agents[i].data.opinion + 0.5 * this->k2_buffer[i]; } ); + // k_4 = hf(x_n+h,y_n+k_3) + get_euler_slopes( k4_buffer, [this]( size_t i ) { return network.agents[i].data.opinion + this->k3_buffer[i]; } ); + + // Update the agent opinions + for( size_t idx_agent = 0; idx_agent < network.n_agents(); ++idx_agent ) + { + // y_(n+1) = y_n+1/6k_1+1/3k_2+1/3k_3+1/6k_4+O(h^5) + network.agents[idx_agent].data.opinion + += ( k1_buffer[idx_agent] + 2 * k2_buffer[idx_agent] + 2 * k3_buffer[idx_agent] + k4_buffer[idx_agent] ) + / 6.0; + } + + if( bot_present() ) + { + for( size_t bot_idx = 0; bot_idx < n_bots; bot_idx++ ) + { + network.agents[bot_idx].data.opinion = bot_opinion[bot_idx]; + } + } +} + +template class ActivityDrivenModelAbstract; +} // namespace Seldon \ No newline at end of file From b85a7514d758c4bb3d7fc7d909963c4e8e7eb24a Mon Sep 17 00:00:00 2001 From: Amrita Goswami Date: Sat, 23 Mar 2024 14:49:52 +0000 Subject: [PATCH 3/9] InertialModel: WIP, compiles somehow We have not yet implemented the InertialModel specific behaviour, but now InertialModel inherits from ActivityDrivenModelAbstract base class. Co-authored-by: Moritz Sallermann --- include/model_factory.hpp | 1 + include/models/ActivityDrivenModel.hpp | 70 ++++++++++++++------------ include/models/InertialModel.hpp | 33 ++++++++++++ meson.build | 1 + src/models/InertialModel.cpp | 6 +-- 5 files changed, 74 insertions(+), 37 deletions(-) create mode 100644 include/models/InertialModel.hpp diff --git a/include/model_factory.hpp b/include/model_factory.hpp index d5864ac..0aaf145 100644 --- a/include/model_factory.hpp +++ b/include/model_factory.hpp @@ -3,6 +3,7 @@ #include "models/ActivityDrivenModel.hpp" #include "models/DeGroot.hpp" #include "models/DeffuantModel.hpp" +#include "models/InertialModel.hpp" #include "network.hpp" #include #include diff --git a/include/models/ActivityDrivenModel.hpp b/include/models/ActivityDrivenModel.hpp index 7801b6f..5cd4401 100644 --- a/include/models/ActivityDrivenModel.hpp +++ b/include/models/ActivityDrivenModel.hpp @@ -60,10 +60,12 @@ class ActivityDrivenModelAbstract : public Model } } - void iteration() override; + void iteration() override {}; -private: +protected: NetworkT & network; + +private: std::vector> contact_prob_list; // Probability of choosing i in 1 to m rounds // Random number generation std::mt19937 & gen; // reference to simulation Mersenne-Twister engine @@ -91,6 +93,7 @@ class ActivityDrivenModelAbstract : public Model double reluctance_eps{}; double covariance_factor{}; +protected: size_t n_bots = 0; // The first n_bots agents are bots std::vector bot_m = std::vector( 0 ); std::vector bot_activity = std::vector( 0 ); @@ -103,6 +106,7 @@ class ActivityDrivenModelAbstract : public Model std::vector k3_buffer{}; std::vector k4_buffer{}; +private: void get_agents_from_power_law() { std::uniform_real_distribution<> dis_opinion( -1, 1 ); // Opinion initial values @@ -142,36 +146,6 @@ class ActivityDrivenModelAbstract : public Model } } - [[nodiscard]] bool bot_present() const - { - return n_bots > 0; - } - - template - void get_euler_slopes( std::vector & k_buffer, Opinion_Callback opinion ) - { - // h is the timestep - size_t j_index = 0; - - k_buffer.resize( network.n_agents() ); - - for( size_t idx_agent = 0; idx_agent < network.n_agents(); ++idx_agent ) - { - auto neighbour_buffer = network.get_neighbours( idx_agent ); // Get the incoming neighbours - auto weight_buffer = network.get_weights( idx_agent ); // Get incoming weights - k_buffer[idx_agent] = -opinion( idx_agent ); - // Loop through neighbouring agents - for( size_t j = 0; j < neighbour_buffer.size(); j++ ) - { - j_index = neighbour_buffer[j]; - k_buffer[idx_agent] += 1.0 / network.agents[idx_agent].data.reluctance * K * weight_buffer[j] - * std::tanh( alpha * opinion( j_index ) ); - } - // Multiply by the timestep - k_buffer[idx_agent] *= dt; - } - } - // The weight for contact between two agents double homophily_weight( size_t idx_contacter, size_t idx_contacted ) { @@ -330,6 +304,12 @@ class ActivityDrivenModelAbstract : public Model } } +protected: + [[nodiscard]] bool bot_present() const + { + return n_bots > 0; + } + void update_network() { @@ -342,9 +322,33 @@ class ActivityDrivenModelAbstract : public Model update_network_mean(); } } + + template + void get_euler_slopes( std::vector & k_buffer, Opinion_Callback opinion ) + { + // h is the timestep + size_t j_index = 0; + + k_buffer.resize( network.n_agents() ); + + for( size_t idx_agent = 0; idx_agent < network.n_agents(); ++idx_agent ) + { + auto neighbour_buffer = network.get_neighbours( idx_agent ); // Get the incoming neighbours + auto weight_buffer = network.get_weights( idx_agent ); // Get incoming weights + k_buffer[idx_agent] = -opinion( idx_agent ); + // Loop through neighbouring agents + for( size_t j = 0; j < neighbour_buffer.size(); j++ ) + { + j_index = neighbour_buffer[j]; + k_buffer[idx_agent] += 1.0 / network.agents[idx_agent].data.reluctance * K * weight_buffer[j] + * std::tanh( alpha * opinion( j_index ) ); + } + // Multiply by the timestep + k_buffer[idx_agent] *= dt; + } + } }; using ActivityDrivenModel = ActivityDrivenModelAbstract; -using InertialModel = ActivityDrivenModelAbstract; } // namespace Seldon \ No newline at end of file diff --git a/include/models/InertialModel.hpp b/include/models/InertialModel.hpp new file mode 100644 index 0000000..89adb93 --- /dev/null +++ b/include/models/InertialModel.hpp @@ -0,0 +1,33 @@ +#pragma once + +#include "agents/activity_agent.hpp" +#include "agents/inertial_agent.hpp" +#include "config_parser.hpp" +#include "model.hpp" +#include "models/ActivityDrivenModel.hpp" +#include "network.hpp" +#include "network_generation.hpp" +#include +#include +#include +#include +#include +#include + +namespace Seldon +{ + +class InertialModel : public ActivityDrivenModelAbstract +{ +public: + using AgentT = InertialAgent; + using NetworkT = Network; + using WeightT = typename NetworkT::WeightT; + + void iteration() override; + +private: + double friction_coefficient = 1.0; +}; + +} // namespace Seldon \ No newline at end of file diff --git a/meson.build b/meson.build index e7017c3..564c1ee 100644 --- a/meson.build +++ b/meson.build @@ -12,6 +12,7 @@ sources_seldon = [ 'src/models/ActivityDrivenModel.cpp', 'src/models/DeffuantModel.cpp', 'src/models/DeffuantModelVector.cpp', + 'src/models/InertialModel.cpp', 'src/util/tomlplusplus.cpp', ] diff --git a/src/models/InertialModel.cpp b/src/models/InertialModel.cpp index a17d3f8..aaf825b 100644 --- a/src/models/InertialModel.cpp +++ b/src/models/InertialModel.cpp @@ -1,5 +1,5 @@ +#include "models/InertialModel.hpp" #include "agents/inertial_agent.hpp" -#include "models/ActivityDrivenModel.hpp" #include "network.hpp" #include "util/math.hpp" #include @@ -9,8 +9,7 @@ namespace Seldon { -template<> -void ActivityDrivenModelAbstract::iteration() +void InertialModel::iteration() { Model::iteration(); @@ -46,5 +45,4 @@ void ActivityDrivenModelAbstract::iteration() } } -template class ActivityDrivenModelAbstract; } // namespace Seldon \ No newline at end of file From ce74473eb06f23785117c91b1ef6e77090fc374c Mon Sep 17 00:00:00 2001 From: Amrita Goswami Date: Sat, 23 Mar 2024 15:54:11 +0000 Subject: [PATCH 4/9] ActivityDrivenModel: Changed get_euler_slopes We changed the definition of get_euler_slopes. Now the timestep is multiplied in the update rule. This enables us to re-use get_euler_slopes in the InertialModel. Co-authored-by: Moritz Sallermann --- include/models/ActivityDrivenModel.hpp | 6 +++--- src/models/ActivityDrivenModel.cpp | 10 ++++++---- 2 files changed, 9 insertions(+), 7 deletions(-) diff --git a/include/models/ActivityDrivenModel.hpp b/include/models/ActivityDrivenModel.hpp index 5cd4401..237bbe9 100644 --- a/include/models/ActivityDrivenModel.hpp +++ b/include/models/ActivityDrivenModel.hpp @@ -71,6 +71,7 @@ class ActivityDrivenModelAbstract : public Model std::mt19937 & gen; // reference to simulation Mersenne-Twister engine std::set> reciprocal_edge_buffer{}; +protected: // Model-specific parameters double dt{}; // Timestep for the integration of the coupled ODEs // Various free parameters @@ -93,7 +94,6 @@ class ActivityDrivenModelAbstract : public Model double reluctance_eps{}; double covariance_factor{}; -protected: size_t n_bots = 0; // The first n_bots agents are bots std::vector bot_m = std::vector( 0 ); std::vector bot_activity = std::vector( 0 ); @@ -343,8 +343,8 @@ class ActivityDrivenModelAbstract : public Model k_buffer[idx_agent] += 1.0 / network.agents[idx_agent].data.reluctance * K * weight_buffer[j] * std::tanh( alpha * opinion( j_index ) ); } - // Multiply by the timestep - k_buffer[idx_agent] *= dt; + // Here, we won't multiply by the timestep. + // Instead multiply in the update rule } } }; diff --git a/src/models/ActivityDrivenModel.cpp b/src/models/ActivityDrivenModel.cpp index bdfc7c4..0b0dd9c 100644 --- a/src/models/ActivityDrivenModel.cpp +++ b/src/models/ActivityDrivenModel.cpp @@ -20,19 +20,21 @@ void ActivityDrivenModelAbstract::iteration() get_euler_slopes( k1_buffer, [this]( size_t i ) { return network.agents[i].data.opinion; } ); // k_2 = hf(x_n+1/2h,y_n+1/2k_1) get_euler_slopes( - k2_buffer, [this]( size_t i ) { return network.agents[i].data.opinion + 0.5 * this->k1_buffer[i]; } ); + k2_buffer, [this]( size_t i ) { return network.agents[i].data.opinion + 0.5 * dt * this->k1_buffer[i]; } ); // k_3 = hf(x_n+1/2h,y_n+1/2k_2) get_euler_slopes( - k3_buffer, [this]( size_t i ) { return network.agents[i].data.opinion + 0.5 * this->k2_buffer[i]; } ); + k3_buffer, [this]( size_t i ) { return network.agents[i].data.opinion + 0.5 * dt * this->k2_buffer[i]; } ); // k_4 = hf(x_n+h,y_n+k_3) - get_euler_slopes( k4_buffer, [this]( size_t i ) { return network.agents[i].data.opinion + this->k3_buffer[i]; } ); + get_euler_slopes( + k4_buffer, [this]( size_t i ) { return network.agents[i].data.opinion + dt * this->k3_buffer[i]; } ); // Update the agent opinions for( size_t idx_agent = 0; idx_agent < network.n_agents(); ++idx_agent ) { // y_(n+1) = y_n+1/6k_1+1/3k_2+1/3k_3+1/6k_4+O(h^5) network.agents[idx_agent].data.opinion - += ( k1_buffer[idx_agent] + 2 * k2_buffer[idx_agent] + 2 * k3_buffer[idx_agent] + k4_buffer[idx_agent] ) + += dt + * ( k1_buffer[idx_agent] + 2 * k2_buffer[idx_agent] + 2 * k3_buffer[idx_agent] + k4_buffer[idx_agent] ) / 6.0; } From 1d7e8389eff142b836fb816d3d79adf63e4c4a46 Mon Sep 17 00:00:00 2001 From: Amrita Goswami Date: Sat, 23 Mar 2024 16:19:32 +0000 Subject: [PATCH 5/9] InertialModel: Implemented update rule Implemented a Velocity Verlet algorithm for the Inertial model. Co-authored-by: Moritz Sallermann --- include/models/InertialModel.hpp | 5 +++ src/models/InertialModel.cpp | 62 +++++++++++++++++++++----------- 2 files changed, 47 insertions(+), 20 deletions(-) diff --git a/include/models/InertialModel.hpp b/include/models/InertialModel.hpp index 89adb93..8f5577d 100644 --- a/include/models/InertialModel.hpp +++ b/include/models/InertialModel.hpp @@ -28,6 +28,11 @@ class InertialModel : public ActivityDrivenModelAbstract private: double friction_coefficient = 1.0; + std::vector drift_t_buffer{}; + std::vector drift_next_t_buffer{}; + + void calc_velocity(); + void calc_position(); }; } // namespace Seldon \ No newline at end of file diff --git a/src/models/InertialModel.cpp b/src/models/InertialModel.cpp index aaf825b..62a5795 100644 --- a/src/models/InertialModel.cpp +++ b/src/models/InertialModel.cpp @@ -9,32 +9,54 @@ namespace Seldon { +// X(t+dt) +// Also updates the position +void InertialModel::calc_position() +{ + // Calculating 'drift' = a(t)-friction + get_euler_slopes( drift_t_buffer, [this]( size_t i ) { return network.agents[i].data.opinion; } ); + + for( size_t idx_agent = 0; idx_agent < network.n_agents(); idx_agent++ ) + { + auto & agent_data = network.agents[idx_agent].data; + auto accleration = drift_t_buffer[idx_agent] - friction_coefficient * agent_data.velocity; + double next_position = agent_data.opinion + agent_data.velocity * dt + 0.5 * (accleration)*dt * dt; + + // Update the position to the new position + agent_data.opinion = next_position; + } +} + +// V(t+dt) +// Also updates the velocity +void InertialModel::calc_velocity() +{ + // Calculating new 'drift' + get_euler_slopes( drift_next_t_buffer, [this]( size_t i ) { return network.agents[i].data.opinion; } ); + + for( size_t idx_agent = 0; idx_agent < network.n_agents(); idx_agent++ ) + { + auto & agent_data = network.agents[idx_agent].data; + double next_velocity = agent_data.velocity + + 0.5 * dt + * ( drift_t_buffer[idx_agent] - friction_coefficient * agent_data.velocity + + drift_next_t_buffer[idx_agent] ); + next_velocity /= 1.0 + 0.5 * friction_coefficient * dt; + + // Update velocity + agent_data.velocity = next_velocity; + } +} + void InertialModel::iteration() { Model::iteration(); update_network(); - // Integrate the ODE using 4th order Runge-Kutta - // k_1 = hf(x_n,y_n) - get_euler_slopes( k1_buffer, [this]( size_t i ) { return network.agents[i].data.opinion; } ); - // k_2 = hf(x_n+1/2h,y_n+1/2k_1) - get_euler_slopes( - k2_buffer, [this]( size_t i ) { return network.agents[i].data.opinion + 0.5 * this->k1_buffer[i]; } ); - // k_3 = hf(x_n+1/2h,y_n+1/2k_2) - get_euler_slopes( - k3_buffer, [this]( size_t i ) { return network.agents[i].data.opinion + 0.5 * this->k2_buffer[i]; } ); - // k_4 = hf(x_n+h,y_n+k_3) - get_euler_slopes( k4_buffer, [this]( size_t i ) { return network.agents[i].data.opinion + this->k3_buffer[i]; } ); - - // Update the agent opinions - for( size_t idx_agent = 0; idx_agent < network.n_agents(); ++idx_agent ) - { - // y_(n+1) = y_n+1/6k_1+1/3k_2+1/3k_3+1/6k_4+O(h^5) - network.agents[idx_agent].data.opinion - += ( k1_buffer[idx_agent] + 2 * k2_buffer[idx_agent] + 2 * k3_buffer[idx_agent] + k4_buffer[idx_agent] ) - / 6.0; - } + // Use Velocity Verlet algorithm + calc_position(); + calc_velocity(); if( bot_present() ) { From dee1ae6dd0855b5cfbcabc27c0193dc69765897e Mon Sep 17 00:00:00 2001 From: Amrita Goswami Date: Sat, 23 Mar 2024 18:01:17 +0000 Subject: [PATCH 6/9] InertialModel: Init of actual instantiation Co-authored-by: Moritz Sallermann --- include/config_parser.hpp | 6 ++++++ include/models/InertialModel.hpp | 2 ++ src/models/InertialModel.cpp | 6 ++++++ 3 files changed, 14 insertions(+) diff --git a/include/config_parser.hpp b/include/config_parser.hpp index 208e27d..6099a3f 100644 --- a/include/config_parser.hpp +++ b/include/config_parser.hpp @@ -25,6 +25,7 @@ enum class Model { DeGroot, ActivityDrivenModel, + InertialModel, DeffuantModel }; @@ -81,6 +82,11 @@ struct ActivityDrivenSettings double covariance_factor = 0.0; }; +struct InertialSettings : public ActivityDrivenSettings +{ + double friction_coefficient = 1.0; +}; + struct InitialNetworkSettings { std::optional file; diff --git a/include/models/InertialModel.hpp b/include/models/InertialModel.hpp index 8f5577d..70b83ef 100644 --- a/include/models/InertialModel.hpp +++ b/include/models/InertialModel.hpp @@ -24,6 +24,8 @@ class InertialModel : public ActivityDrivenModelAbstract using NetworkT = Network; using WeightT = typename NetworkT::WeightT; + InertialModel( const Config::InertialSettings & settings, NetworkT & network, std::mt19937 & gen ); + void iteration() override; private: diff --git a/src/models/InertialModel.cpp b/src/models/InertialModel.cpp index 62a5795..0d982c2 100644 --- a/src/models/InertialModel.cpp +++ b/src/models/InertialModel.cpp @@ -9,6 +9,12 @@ namespace Seldon { +InertialModel::InertialModel( const Config::InertialSettings & settings, NetworkT & network, std::mt19937 & gen ) + : ActivityDrivenModelAbstract( settings, network, gen ), + friction_coefficient( settings.friction_coefficient ) +{ +} + // X(t+dt) // Also updates the position void InertialModel::calc_position() From ece80c64d8bb3bd2de6b35423f568c31cbc7758d Mon Sep 17 00:00:00 2001 From: Amrita Goswami Date: Sat, 23 Mar 2024 18:36:41 +0000 Subject: [PATCH 7/9] InertialModel: Completed implementation We have implemented the Inertial model (more specifically the activity driven inertial model). The example added runs and produces output. TODO: unit test the shit of this! Co-authored-by: Moritz Sallermann --- examples/ActivityDrivenInertial/conf.toml | 34 ++++ include/config_parser.hpp | 9 +- include/model_factory.hpp | 17 ++ include/models/InertialModel.hpp | 2 +- include/simulation.hpp | 4 + src/config_parser.cpp | 181 +++++++++++++--------- src/main.cpp | 6 + src/models/InertialModel.cpp | 3 +- 8 files changed, 179 insertions(+), 77 deletions(-) create mode 100644 examples/ActivityDrivenInertial/conf.toml diff --git a/examples/ActivityDrivenInertial/conf.toml b/examples/ActivityDrivenInertial/conf.toml new file mode 100644 index 0000000..383ba76 --- /dev/null +++ b/examples/ActivityDrivenInertial/conf.toml @@ -0,0 +1,34 @@ +[simulation] +model = "ActivityDrivenInertial" +# rng_seed = 120 # Leaving this empty will pick a random seed + +[io] +n_output_network = 20 # Write the network every 20 iterations +n_output_agents = 1 # Write the opinions of agents after every iteration +print_progress = true # Print the iteration time ; if not set, then does not print + +[model] +max_iterations = 500 # If not set, max iterations is infinite + +[ActivityDrivenInertial] +dt = 0.01 # Timestep for the integration of the coupled ODEs +m = 10 # Number of agents contacted, when the agent is active +eps = 0.01 # Minimum activity epsilon; a_i belongs to [epsilon,1] +gamma = 2.1 # Exponent of activity power law distribution of activities +reciprocity = 0.65 # probability that when agent i contacts j via weighted reservoir sampling, j also sends feedback to i. So every agent can have more than m incoming connections +homophily = 1.0 # aka beta. if zero, agents pick their interaction partners at random +alpha = 3.0 # Controversialness of the issue, must be greater than 0. +K = 2.0 +mean_activities = false # Use the mean value of the powerlaw distribution for the activities of all agents +mean_weights = false # Use the meanfield approximation of the network edges + +reluctances = true # Assigns a "reluctance" (m_i) to each agent. By default; false and every agent has a reluctance of 1 +reluctance_mean = 1.0 # Mean of distribution before drawing from a truncated normal distribution (default set to 1.0) +reluctance_sigma = 0.25 # Width of normal distribution (before truncating) +reluctance_eps = 0.01 # Minimum such that the normal distribution is truncated at this value + +friction_coefficient = 1.0 # Friction coefficient, making agents tend to go to rest without acceleration + +[network] +number_of_agents = 1000 +connections_per_agent = 10 diff --git a/include/config_parser.hpp b/include/config_parser.hpp index 6099a3f..326dc59 100644 --- a/include/config_parser.hpp +++ b/include/config_parser.hpp @@ -24,8 +24,8 @@ namespace Seldon::Config enum class Model { DeGroot, - ActivityDrivenModel, - InertialModel, + ActivityDrivenModel, // @TODO : no need for model here + ActivityDrivenInertial, DeffuantModel }; @@ -82,7 +82,7 @@ struct ActivityDrivenSettings double covariance_factor = 0.0; }; -struct InertialSettings : public ActivityDrivenSettings +struct ActivityDrivenInertialSettings : public ActivityDrivenSettings { double friction_coefficient = 1.0; }; @@ -96,7 +96,8 @@ struct InitialNetworkSettings struct SimulationOptions { - using ModelVariantT = std::variant; + using ModelVariantT + = std::variant; Model model; std::string model_string; int rng_seed = std::random_device()(); diff --git a/include/model_factory.hpp b/include/model_factory.hpp index 0aaf145..a87f9dd 100644 --- a/include/model_factory.hpp +++ b/include/model_factory.hpp @@ -62,6 +62,23 @@ create_model_activity_driven( Network & network, const ModelVariantT & m } } +template +inline auto create_model_activity_driven_inertial( + Network & network, const ModelVariantT & model_settings, std::mt19937 & gen ) +{ + if constexpr( std::is_same_v ) + { + auto settings = std::get( model_settings ); + auto model = std::make_unique( settings, network, gen ); + return model; + } + else + { + throw std::runtime_error( "Incompatible agent and model type!" ); + return std::unique_ptr>{}; + } +} + template inline auto create_model_deffuant( Network & network, const ModelVariantT & model_settings, std::mt19937 & gen ) { diff --git a/include/models/InertialModel.hpp b/include/models/InertialModel.hpp index 70b83ef..ba95a0b 100644 --- a/include/models/InertialModel.hpp +++ b/include/models/InertialModel.hpp @@ -24,7 +24,7 @@ class InertialModel : public ActivityDrivenModelAbstract using NetworkT = Network; using WeightT = typename NetworkT::WeightT; - InertialModel( const Config::InertialSettings & settings, NetworkT & network, std::mt19937 & gen ); + InertialModel( const Config::ActivityDrivenInertialSettings & settings, NetworkT & network, std::mt19937 & gen ); void iteration() override; diff --git a/include/simulation.hpp b/include/simulation.hpp index 4eca2fc..118b824 100644 --- a/include/simulation.hpp +++ b/include/simulation.hpp @@ -70,6 +70,10 @@ class Simulation : public SimulationInterface { model = ModelFactory::create_model_activity_driven( network, options.model_settings, gen ); } + else if( options.model == Config::Model::ActivityDrivenInertial ) + { + model = ModelFactory::create_model_activity_driven_inertial( network, options.model_settings, gen ); + } else if( options.model == Config::Model::DeffuantModel ) { auto deffuant_settings = std::get( options.model_settings ); diff --git a/src/config_parser.cpp b/src/config_parser.cpp index b65d83a..63c8456 100644 --- a/src/config_parser.cpp +++ b/src/config_parser.cpp @@ -29,9 +29,82 @@ Model model_string_to_enum( std::string_view model_string ) { return Model::DeffuantModel; } + else if( model_string == "ActivityDrivenInertial" ) + { + return Model::ActivityDrivenInertial; + } throw std::runtime_error( fmt::format( "Invalid model string {}", model_string ) ); } +void set_if_specified( auto & opt, const auto & toml_opt ) +{ + using T = typename std::remove_reference::type; + auto t_opt = toml_opt.template value(); + if( t_opt.has_value() ) + opt = t_opt.value(); +} + +// Convenience function to parse activity settings +void parse_activity_settings( auto & model_settings, const auto & toml_model_opt, const auto & tbl ) +{ + set_if_specified( model_settings.dt, toml_model_opt["dt"] ); + set_if_specified( model_settings.m, toml_model_opt["m"] ); + set_if_specified( model_settings.eps, toml_model_opt["eps"] ); + set_if_specified( model_settings.gamma, toml_model_opt["gamma"] ); + set_if_specified( model_settings.homophily, toml_model_opt["homophily"] ); + set_if_specified( model_settings.reciprocity, toml_model_opt["reciprocity"] ); + set_if_specified( model_settings.alpha, toml_model_opt["alpha"] ); + set_if_specified( model_settings.K, toml_model_opt["K"] ); + // Mean activity model options + set_if_specified( model_settings.mean_activities, toml_model_opt["mean_activities"] ); + set_if_specified( model_settings.mean_weights, toml_model_opt["mean_weights"] ); + // Reluctances + set_if_specified( model_settings.use_reluctances, toml_model_opt["reluctances"] ); + set_if_specified( model_settings.reluctance_mean, toml_model_opt["reluctance_mean"] ); + set_if_specified( model_settings.reluctance_sigma, toml_model_opt["reluctance_sigma"] ); + set_if_specified( model_settings.reluctance_eps, toml_model_opt["reluctance_eps"] ); + + model_settings.max_iterations = tbl["model"]["max_iterations"].template value(); + + // bot + set_if_specified( model_settings.n_bots, toml_model_opt["n_bots"] ); + + auto push_back_bot_array = [&]( auto toml_node, auto & options_array, auto default_value ) + { + if( toml_node.is_array() ) + { + toml::array * toml_arr = toml_node.as_array(); + + toml_arr->for_each( + [&]( auto && elem ) + { + if( elem.is_integer() ) + { + options_array.push_back( elem.as_integer()->get() ); + } + else if( elem.is_floating_point() ) + { + options_array.push_back( elem.as_floating_point()->get() ); + } + } ); + } + else + { + options_array = std::vector( model_settings.n_bots, default_value ); + } + }; + + auto bot_opinion = toml_model_opt["bot_opinion"]; + auto bot_m = toml_model_opt["bot_m"]; + auto bot_activity = toml_model_opt["bot_activity"]; + auto bot_homophily = toml_model_opt["bot_homophily"]; + + push_back_bot_array( bot_m, model_settings.bot_m, model_settings.m ); + push_back_bot_array( bot_opinion, model_settings.bot_opinion, 0.0 ); + push_back_bot_array( bot_activity, model_settings.bot_activity, 1.0 ); + push_back_bot_array( bot_homophily, model_settings.bot_homophily, model_settings.homophily ); +} + SimulationOptions parse_config_file( std::string_view config_file_path ) { SimulationOptions options; @@ -41,14 +114,6 @@ SimulationOptions parse_config_file( std::string_view config_file_path ) options.rng_seed = tbl["simulation"]["rng_seed"].value_or( int( options.rng_seed ) ); - auto set_if_specified = [&]( auto & opt, const auto & toml_opt ) - { - using T = typename std::remove_reference::type; - auto t_opt = toml_opt.template value(); - if( t_opt.has_value() ) - opt = t_opt.value(); - }; - // Parse output settings options.output_settings.n_output_network = tbl["io"]["n_output_network"].value(); options.output_settings.n_output_agents = tbl["io"]["n_output_agents"].value(); @@ -87,63 +152,15 @@ SimulationOptions parse_config_file( std::string_view config_file_path ) { auto model_settings = ActivityDrivenSettings(); - set_if_specified( model_settings.dt, tbl[options.model_string]["dt"] ); - set_if_specified( model_settings.m, tbl[options.model_string]["m"] ); - set_if_specified( model_settings.eps, tbl[options.model_string]["eps"] ); - set_if_specified( model_settings.gamma, tbl[options.model_string]["gamma"] ); - set_if_specified( model_settings.homophily, tbl[options.model_string]["homophily"] ); - set_if_specified( model_settings.reciprocity, tbl[options.model_string]["reciprocity"] ); - set_if_specified( model_settings.alpha, tbl[options.model_string]["alpha"] ); - set_if_specified( model_settings.K, tbl[options.model_string]["K"] ); - // Mean activity model options - set_if_specified( model_settings.mean_activities, tbl[options.model_string]["mean_activities"] ); - set_if_specified( model_settings.mean_weights, tbl[options.model_string]["mean_weights"] ); - // Reluctances - set_if_specified( model_settings.use_reluctances, tbl[options.model_string]["reluctances"] ); - set_if_specified( model_settings.reluctance_mean, tbl[options.model_string]["reluctance_mean"] ); - set_if_specified( model_settings.reluctance_sigma, tbl[options.model_string]["reluctance_sigma"] ); - set_if_specified( model_settings.reluctance_eps, tbl[options.model_string]["reluctance_eps"] ); - - model_settings.max_iterations = tbl["model"]["max_iterations"].value(); - - // bot - set_if_specified( model_settings.n_bots, tbl[options.model_string]["n_bots"] ); - - auto push_back_bot_array = [&]( auto toml_node, auto & options_array, auto default_value ) - { - if( toml_node.is_array() ) - { - toml::array * toml_arr = toml_node.as_array(); - - toml_arr->for_each( - [&]( auto && elem ) - { - if( elem.is_integer() ) - { - options_array.push_back( elem.as_integer()->get() ); - } - else if( elem.is_floating_point() ) - { - options_array.push_back( elem.as_floating_point()->get() ); - } - } ); - } - else - { - options_array = std::vector( model_settings.n_bots, default_value ); - } - }; - - auto bot_opinion = tbl[options.model_string]["bot_opinion"]; - auto bot_m = tbl[options.model_string]["bot_m"]; - auto bot_activity = tbl[options.model_string]["bot_activity"]; - auto bot_homophily = tbl[options.model_string]["bot_homophily"]; - - push_back_bot_array( bot_m, model_settings.bot_m, model_settings.m ); - push_back_bot_array( bot_opinion, model_settings.bot_opinion, 0.0 ); - push_back_bot_array( bot_activity, model_settings.bot_activity, 1.0 ); - push_back_bot_array( bot_homophily, model_settings.bot_homophily, model_settings.homophily ); + parse_activity_settings( model_settings, tbl[options.model_string], tbl ); + options.model_settings = model_settings; + } + else if( options.model == Model::ActivityDrivenInertial ) + { + auto model_settings = ActivityDrivenInertialSettings(); + parse_activity_settings( model_settings, tbl[options.model_string], tbl ); + set_if_specified( model_settings.friction_coefficient, tbl[options.model_string]["friction_coefficient"] ); options.model_settings = model_settings; } @@ -183,10 +200,8 @@ void validate_settings( const SimulationOptions & options ) // @TODO: Check that start_output is less than the max_iterations? check( name_and_var( options.output_settings.start_output ), g_zero ); - if( options.model == Model::ActivityDrivenModel ) + auto validate_activity = [&]( const auto & model_settings ) { - auto model_settings = std::get( options.model_settings ); - check( name_and_var( model_settings.dt ), g_zero ); check( name_and_var( model_settings.m ), geq_zero ); check( name_and_var( model_settings.eps ), g_zero ); @@ -207,6 +222,19 @@ void validate_settings( const SimulationOptions & options ) check( name_and_var( model_settings.bot_activity ), check_bot_size, bot_msg ); check( name_and_var( model_settings.bot_opinion ), check_bot_size, bot_msg ); check( name_and_var( model_settings.bot_homophily ), check_bot_size, bot_msg ); + }; + + if( options.model == Model::ActivityDrivenModel ) + { + auto model_settings = std::get( options.model_settings ); + + validate_activity( model_settings ); + } + else if( options.model == Model::ActivityDrivenInertial ) + { + auto model_settings = std::get( options.model_settings ); + check( name_and_var( model_settings.friction_coefficient ), geq_zero ); + validate_activity( model_settings ); } else if( options.model == Model::DeGroot ) { @@ -237,11 +265,8 @@ void print_settings( const SimulationOptions & options ) fmt::print( "[Model]\n" ); fmt::print( " type {}\n", options.model_string ); - // @TODO: Optionally print *all* settings to the console, including defaults that were set - if( options.model == Model::ActivityDrivenModel ) + auto print_activity_settings = [&]( auto model_settings ) { - auto model_settings = std::get( options.model_settings ); - fmt::print( " max_iterations {}\n", model_settings.max_iterations ); fmt::print( " dt {} \n", model_settings.dt ); fmt::print( " m {} \n", model_settings.m ); @@ -270,6 +295,20 @@ void print_settings( const SimulationOptions & options ) fmt::print( " reluctance_eps {}\n", model_settings.reluctance_eps ); fmt::print( " covariance_factor {}\n", model_settings.covariance_factor ); } + }; + + // @TODO: Optionally print *all* settings to the console, including defaults that were set + if( options.model == Model::ActivityDrivenModel ) + { + auto model_settings = std::get( options.model_settings ); + + print_activity_settings( model_settings ); + } + else if( options.model == Model::ActivityDrivenInertial ) + { + auto model_settings = std::get( options.model_settings ); + print_activity_settings( model_settings ); + fmt::print( " friction_coefficient {}\n", model_settings.friction_coefficient ); } else if( options.model == Model::DeGroot ) { diff --git a/src/main.cpp b/src/main.cpp index 9343730..3f429fa 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -1,6 +1,7 @@ #include "config_parser.hpp" #include "models/DeGroot.hpp" #include "models/DeffuantModel.hpp" +#include "models/InertialModel.hpp" #include "simulation.hpp" #include #include @@ -70,6 +71,11 @@ int main( int argc, char * argv[] ) simulation = std::make_unique>( simulation_options, network_file, agent_file ); } + else if( simulation_options.model == Seldon::Config::Model::ActivityDrivenInertial ) + { + simulation = std::make_unique>( + simulation_options, network_file, agent_file ); + } else if( simulation_options.model == Seldon::Config::Model::DeffuantModel ) { simulation = std::make_unique>( diff --git a/src/models/InertialModel.cpp b/src/models/InertialModel.cpp index 0d982c2..df46cc1 100644 --- a/src/models/InertialModel.cpp +++ b/src/models/InertialModel.cpp @@ -9,7 +9,8 @@ namespace Seldon { -InertialModel::InertialModel( const Config::InertialSettings & settings, NetworkT & network, std::mt19937 & gen ) +InertialModel::InertialModel( + const Config::ActivityDrivenInertialSettings & settings, NetworkT & network, std::mt19937 & gen ) : ActivityDrivenModelAbstract( settings, network, gen ), friction_coefficient( settings.friction_coefficient ) { From 6c0ff3ef5458b78f6d8185c280244dea4a6def79 Mon Sep 17 00:00:00 2001 From: Amrita Goswami Date: Sat, 23 Mar 2024 19:27:15 +0000 Subject: [PATCH 8/9] InertialModel: Implemented unit test We have implemented a unit test with 1 agent and 1 bot for the activity driven inertial model Co-authored-by: Moritz Sallermann --- meson.build | 1 + test/res/1bot_1agent_inertial.toml | 42 +++++++++++++++++ test/test_activity_inertial.cpp | 76 ++++++++++++++++++++++++++++++ 3 files changed, 119 insertions(+) create mode 100644 test/res/1bot_1agent_inertial.toml create mode 100644 test/test_activity_inertial.cpp diff --git a/meson.build b/meson.build index 564c1ee..d5a4548 100644 --- a/meson.build +++ b/meson.build @@ -27,6 +27,7 @@ tests = [ ['Test_Tarjan', 'test/test_tarjan.cpp'], ['Test_DeGroot', 'test/test_deGroot.cpp'], ['Test_Activity_Driven', 'test/test_activity.cpp'], + ['Test_Activity_Driven_Inertial', 'test/test_activity_inertial.cpp'], ['Test_Deffuant', 'test/test_deffuant.cpp'], ['Test_Network', 'test/test_network.cpp'], ['Test_Network_Generation', 'test/test_network_generation.cpp'], diff --git a/test/res/1bot_1agent_inertial.toml b/test/res/1bot_1agent_inertial.toml new file mode 100644 index 0000000..6520309 --- /dev/null +++ b/test/res/1bot_1agent_inertial.toml @@ -0,0 +1,42 @@ +[simulation] +model = "ActivityDrivenInertial" +rng_seed = 120 # Leaving this empty will pick a random seed + +[io] +# n_output_network = 1 # Write the network every 20 iterations +# n_output_agents = 1 +print_progress = false # Print the iteration time ; if not set, then does not print + +[model] +max_iterations = 1000 # If not set, max iterations is infinite + +[ActivityDrivenInertial] +dt = 0.001 # Timestep for the integration of the coupled ODEs +m = 1 # Number of agents contacted, when the agent is active +eps = 1 # Minimum activity epsilon; a_i belongs to [epsilon,1] +gamma = 2.1 # Exponent of activity power law distribution of activities +reciprocity = 1 # probability that when agent i contacts j via weighted reservoir sampling, j also sends feedback to i. So every agent can have more than m incoming connections +homophily = 0.5 # aka beta. if zero, agents pick their interaction partners at random +alpha = 1.5 # Controversialness of the issue, must be greater than 0. +K = 2.0 # Social interaction strength +mean_activities = false # Use the mean value of the powerlaw distribution for the activities of all agents +mean_weights = false # Use the meanfield approximation of the network edges + +reluctances = true # Assigns a "reluctance" (m_i) to each agent. By default; false and every agent has a reluctance of 1 +reluctance_mean = 1.5 # Mean of distribution before drawing from a truncated normal distribution (default set to 1.0) +reluctance_sigma = 0.1 # Width of normal distribution (before truncating) +reluctance_eps = 0.01 # Minimum such that the normal distribution is truncated at this value + +n_bots = 1 # The number of bots to be used; if not specified defaults to 0 (which means bots are deactivated) +# Bots are agents with fixed opinions and different parameters, the parameters are specified in the following lists +# If n_bots is smaller than the length of any of the lists, the first n_bots entries are used. If n_bots is greater the code will throw an exception. +bot_m = [1] # If not specified, defaults to `m` +bot_homophily = [0.7] # If not specified, defaults to `homophily` +bot_activity = [1.0] # If not specified, defaults to 0 +bot_opinion = [2] # The fixed opinions of the bots + +friction_coefficient = 0.5 # Friction coefficient + +[network] +number_of_agents = 2 +connections_per_agent = 1 diff --git a/test/test_activity_inertial.cpp b/test/test_activity_inertial.cpp new file mode 100644 index 0000000..60f3eef --- /dev/null +++ b/test/test_activity_inertial.cpp @@ -0,0 +1,76 @@ +#include "catch2/matchers/catch_matchers.hpp" +#include "models/InertialModel.hpp" +#include "util/math.hpp" +#include +#include +#include +#include +#include +#include +#include +#include + +namespace fs = std::filesystem; + +TEST_CASE( "Test the probabilistic inertial activity driven model with one bot and one agent", "[inertial1Bot1Agent]" ) +{ + using namespace Seldon; + using namespace Catch::Matchers; + using AgentT = InertialModel::AgentT; + + auto proj_root_path = fs::current_path(); + auto input_file = proj_root_path / fs::path( "test/res/1bot_1agent_inertial.toml" ); + + auto options = Config::parse_config_file( input_file.string() ); + Config::print_settings( options ); + + auto simulation = Simulation( options, std::nullopt, std::nullopt ); + + // We need an output path for Simulation, but we won't write anything out there + fs::path output_dir_path = proj_root_path / fs::path( "test/output_inertial" ); + + // Get the bot opinion (which won't change) + auto bot = simulation.network.agents[0]; + auto x_bot = bot.data.opinion; // Bot opinion + fmt::print( "We have set the bot opinion = {}\n", x_bot ); + + // Get the initial agent opinion + auto & agent = simulation.network.agents[1]; + auto x_0 = agent.data.opinion; // Agent opinion + fmt::print( "We have set agent x_0 = {}\n", x_0 ); + + simulation.run( output_dir_path ); + + auto model_settings = std::get( options.model_settings ); + auto K = model_settings.K; + auto alpha = model_settings.alpha; + auto iterations = model_settings.max_iterations.value(); + auto dt = model_settings.dt; + std::complex mu = model_settings.friction_coefficient; + auto time_elapsed = iterations * dt; + + // Final agent and bot opinions after the simulation run + auto x_t = agent.data.opinion; + auto x_t_bot = bot.data.opinion; + auto reluctance = agent.data.reluctance; + + fmt::print( "Agent reluctance is = {}\n", reluctance ); + + // The bot opinion should not change during the simulation + REQUIRE_THAT( x_t_bot, WithinAbs( x_bot, 1e-16 ) ); + + // C = K/m tanh (alpha*x_bot) + auto C = K / reluctance * std::tanh( alpha * x_bot ); + + fmt::print( "C = {}\n", C ); + auto a1 = 0.5 * ( -std::sqrt( mu * mu - 4.0 ) - mu ); + auto a2 = 0.5 * ( std::sqrt( mu * mu - 4.0 ) - mu ); + auto c1 = ( x_0 - C ) / ( 1.0 - a1 / a2 ); + auto c2 = -c1 * a1 / a2; + // Test that the agent opinion matches the analytical solution for an agent with a bot + // Analytical solution is + // x_t = c1 * exp(a1*t) + c2 *exp(a2*t) + C + auto x_t_analytical = c1 * std::exp( a1 * time_elapsed ) + c2 * std::exp( a2 * time_elapsed ) + C; + + REQUIRE_THAT( x_t, WithinAbs( x_t_analytical.real(), 1e-5 ) ); +} \ No newline at end of file From c28c2f32a81fd9fbd476477798c01db8e0b7cc57 Mon Sep 17 00:00:00 2001 From: Amrita Goswami Date: Sun, 24 Mar 2024 15:09:34 +0000 Subject: [PATCH 9/9] CI: Fixed small bug that caused crash on the CI We removed a superfluous `using WeightT = typename Network::WeightT;`. These are present at the top of the class Co-authored-by: Moritz Sallermann --- include/connectivity.hpp | 2 +- include/models/ActivityDrivenModel.hpp | 3 +-- src/config_parser.cpp | 3 ++- test/test_network.cpp | 2 +- 4 files changed, 5 insertions(+), 5 deletions(-) diff --git a/include/connectivity.hpp b/include/connectivity.hpp index a9d0159..1af9c8d 100644 --- a/include/connectivity.hpp +++ b/include/connectivity.hpp @@ -75,7 +75,7 @@ class TarjanConnectivityAlgo { lowest[v] = std::min( lowest[v], num[u] ); } // u not processed - } // u has been visited + } // u has been visited } // Now v has been processed diff --git a/include/models/ActivityDrivenModel.hpp b/include/models/ActivityDrivenModel.hpp index 237bbe9..1b0ab83 100644 --- a/include/models/ActivityDrivenModel.hpp +++ b/include/models/ActivityDrivenModel.hpp @@ -60,7 +60,7 @@ class ActivityDrivenModelAbstract : public Model } } - void iteration() override {}; + void iteration() override{}; protected: NetworkT & network; @@ -234,7 +234,6 @@ class ActivityDrivenModelAbstract : public Model void update_network_mean() { - using WeightT = NetworkT::WeightT; std::vector weights( network.n_agents(), 0.0 ); // Set all weights to zero in the beginning diff --git a/src/config_parser.cpp b/src/config_parser.cpp index 63c8456..e0e054e 100644 --- a/src/config_parser.cpp +++ b/src/config_parser.cpp @@ -253,7 +253,8 @@ void validate_settings( const SimulationOptions & options ) { const std::string basic_deff_msg = "The basic Deffuant model has not been implemented with multiple dimensions"; - check( name_and_var( model_settings.dim ), []( auto x ) { return x == 1; }, basic_deff_msg ); + check( + name_and_var( model_settings.dim ), []( auto x ) { return x == 1; }, basic_deff_msg ); } } } diff --git a/test/test_network.cpp b/test/test_network.cpp index f1d2bfa..57239e7 100644 --- a/test/test_network.cpp +++ b/test/test_network.cpp @@ -112,7 +112,7 @@ TEST_CASE( "Testing the network class" ) auto weight = buffer_w[i_neighbour]; std::tuple edge{ neighbour, i_agent, weight - }; // Note that i_agent and neighbour are flipped compared to before + }; // Note that i_agent and neighbour are flipped compared to before REQUIRE( old_edges.contains( edge ) ); // can we find the transposed edge? old_edges.extract( edge ); // extract the edge afterwards }