From b04c5fcc5a393504c83bb1a4f06e2e81a2a288bd Mon Sep 17 00:00:00 2001 From: Amrita Goswami Date: Sat, 23 Mar 2024 13:48:30 +0000 Subject: [PATCH 01/17] 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 02/17] 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 03/17] 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 04/17] 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 05/17] 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 06/17] 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 07/17] 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 08/17] 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 09/17] 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 } From a0cd5ce1dcc4ba520a7e321b23e5a243ce2170ec Mon Sep 17 00:00:00 2001 From: Moritz Sallermann Date: Sat, 23 Mar 2024 15:16:20 +0000 Subject: [PATCH 10/17] Math: Implement bivariate gaussian copula (WIP) Co-authored-by: Amrita Goswami --- include/util/math.hpp | 94 ++++++++++++++++++++++++- meson.build | 1 + test/test_probability_distributions.cpp | 70 ++++++++++++++++++ 3 files changed, 164 insertions(+), 1 deletion(-) create mode 100644 test/test_probability_distributions.cpp diff --git a/include/util/math.hpp b/include/util/math.hpp index e6041fc..bbde5e6 100644 --- a/include/util/math.hpp +++ b/include/util/math.hpp @@ -1,4 +1,6 @@ #pragma once +#include "fmt/core.h" +#include "util/erfinv.hpp" #include #include #include @@ -129,9 +131,14 @@ class power_law_distribution template ScalarT operator()( Generator & gen ) + { + return inverse_cumulative_probability( dist( gen ) ); + } + + ScalarT inverse_cumulative_probability( ScalarT x ) { return std::pow( - ( 1.0 - std::pow( eps, ( 1.0 - gamma ) ) ) * dist( gen ) + std::pow( eps, ( 1.0 - gamma ) ), + ( 1.0 - std::pow( eps, ( 1.0 - gamma ) ) ) * x + std::pow( eps, ( 1.0 - gamma ) ), ( 1.0 / ( 1.0 - gamma ) ) ); } @@ -157,6 +164,16 @@ class truncated_normal_distribution std::normal_distribution normal_dist{}; size_t max_tries = 5000; + ScalarT inverse_cum_gauss( ScalarT y ) + { + return erfinv( 2.0 * y - 1 ) * std::sqrt( 2.0 ) * sigma + mean; + } + + ScalarT cum_gauss( ScalarT x ) + { + return 0.5 * ( 1 + std::erf( ( x - mean ) / ( sigma * std::sqrt( 2.0 ) ) ) ); + } + public: truncated_normal_distribution( ScalarT mean, ScalarT sigma, ScalarT eps ) : mean( mean ), sigma( sigma ), eps( eps ), normal_dist( std::normal_distribution( mean, sigma ) ) @@ -174,6 +191,81 @@ class truncated_normal_distribution } return eps; } + + ScalarT inverse_cumulative_probability( ScalarT y ) + { + return inverse_cum_gauss( + y * ( 1.0 - cum_gauss( eps, sigma, mean ) ) + cum_gauss( eps, sigma, mean ), sigma, mean ); + } +}; + +/** + * @brief Bivariate normal distribution + * with mean mu = [0,0] + * and covariance matrix Sigma = [[1, cov], [cov, 1]] + * |cov| < 1 is required + */ +template +class bivariate_normal_distribution +{ +private: + ScalarT covariance; + std::normal_distribution normal_dist{}; + +public: + bivariate_normal_distribution( ScalarT covariance ) : covariance( covariance ) {} + + template + std::array operator()( Generator & gen ) + { + ScalarT n1 = normal_dist( gen ); + ScalarT n2 = normal_dist( gen ); + + ScalarT r1 = n1; + ScalarT r2 = covariance * n1 + std::sqrt( 1 - covariance * covariance ); + + return { r1, r2 }; + } +}; + +template +class bivariate_gaussian_copula +{ +private: + ScalarT covariance; + bivariate_normal_distribution bivariate_normal_dist{}; + // std::normal_distribution normal_dist{}; + + // Cumulative probability function for gaussian with mean 0 and variance 1 + ScalarT cum_gauss( ScalarT x ) + { + return 0.5 * ( 1 + std::erf( ( x ) / std::sqrt( 2.0 ) ) ); + } + + dist1T dist1; + dist2T dist2; + +public: + bivariate_gaussian_copula( ScalarT covariance, dist1T dist1, dist2T dist2 ) + : covariance( covariance ), + dist1( dist1 ), + dist2( dist2 ), + bivariate_normal_dist( bivariate_normal_dist( covariance ) ) + { + } + + template + std::array operator()( Generator & gen ) + { + // 1. Draw from bivariate gaussian + auto z = bivariate_normal_dist( gen ); + // 2. Transform marginals to unit interval + std::array z_unit = { cum_gauss( z[0] ), cum_gauss( z[1] ) }; + // 3. Apply inverse transform sampling + std::array res + = { dist1.inverse_cumulative_probability( z_unit[0] ), dist2.inverse_cumulative_probability( z_unit[1] ) }; + return res; + } }; template diff --git a/meson.build b/meson.build index d5a4548..c9fd40d 100644 --- a/meson.build +++ b/meson.build @@ -34,6 +34,7 @@ tests = [ ['Test_Sampling', 'test/test_sampling.cpp'], ['Test_IO', 'test/test_io.cpp'], ['Test_Util', 'test/test_util.cpp'], + ['Test_Prob', 'test/test_probability_distributions.cpp'], ] Catch2 = dependency('Catch2', method : 'cmake', modules : ['Catch2::Catch2WithMain', 'Catch2::Catch2']) diff --git a/test/test_probability_distributions.cpp b/test/test_probability_distributions.cpp new file mode 100644 index 0000000..74fd6ab --- /dev/null +++ b/test/test_probability_distributions.cpp @@ -0,0 +1,70 @@ +#include +#include +#include + +#include "util/math.hpp" +#include +#include +#include +#include +#include +namespace fs = std::filesystem; + +template +std::ostream & operator<<( std::ostream & os, std::array const & v1 ) +{ + std::for_each( begin( v1 ), end( v1 ), [&os]( int val ) { os << val << " "; } ); + return os; +} + +// Samples the distribution n_samples times and writes results to file +template +void write_results_to_file( int N_Samples, distT dist, const std::string & filename ) +{ + auto proj_root_path = fs::current_path(); + auto file = proj_root_path / fs::path( "/test/output/" + filename ); + fs::create_directories( file ); + + auto gen = std::mt19937( 0 ); + + std::ofstream filestream( file ); + filestream << std::setprecision( 16 ); + + for( int i = 0; i < N_Samples; i++ ) + { + filestream << dist( gen ) << "\n"; + } + filestream.close(); +} + +TEST_CASE( "Test the probability distributions", "[prob]" ) +{ + write_results_to_file( 10000, Seldon::truncated_normal_distribution( 1.0, 0.5, 0.1 ), "truncated_normal.txt" ); + write_results_to_file( 10000, Seldon::power_law_distribution( 0.01, 2.1 ), "power_law.txt" ); + write_results_to_file( 10000, Seldon::bivariate_normal_distribution( 0.5 ), "bivariate_normal.txt" ); +} + +// TEST_CASE( "Test reading in the agents from a file", "[io_agents]" ) +// { +// using namespace Seldon; +// using namespace Catch::Matchers; + +// auto proj_root_path = fs::current_path(); +// auto network_file = proj_root_path / fs::path( "test/res/opinions.txt" ); + +// auto agents = Seldon::agents_from_file( network_file ); + +// std::vector opinions_expected = { 2.1127107987061544, 0.8088982488089491, -0.8802809369462433 }; +// std::vector activities_expected = { 0.044554683389757696, 0.015813166022685163, 0.015863953902810535 }; +// std::vector reluctances_expected = { 1.0, 1.0, 2.3 }; + +// REQUIRE( agents.size() == 3 ); + +// for( size_t i = 0; i < agents.size(); i++ ) +// { +// fmt::print( "{}", i ); +// REQUIRE_THAT( agents[i].data.opinion, Catch::Matchers::WithinAbs( opinions_expected[i], 1e-16 ) ); +// REQUIRE_THAT( agents[i].data.activity, Catch::Matchers::WithinAbs( activities_expected[i], 1e-16 ) ); +// REQUIRE_THAT( agents[i].data.reluctance, Catch::Matchers::WithinAbs( reluctances_expected[i], 1e-16 ) ); +// } +// } \ No newline at end of file From 961e4938b270f60e7294efd09cd2a700d4e9270d Mon Sep 17 00:00:00 2001 From: Moritz Sallermann Date: Sun, 24 Mar 2024 19:53:02 +0000 Subject: [PATCH 11/17] Test: Test the probability distributions Test is a generous word, as we actually just write out some results and see if it crashes. But we plotted the histograms vs the expected pdf in python and they match. It's hard to quantify this in a unit test. --- include/util/erfinv.hpp | 112 ++++++++++++++++++++++++ include/util/math.hpp | 20 ++--- test/test_probability_distributions.cpp | 14 +-- 3 files changed, 126 insertions(+), 20 deletions(-) create mode 100644 include/util/erfinv.hpp diff --git a/include/util/erfinv.hpp b/include/util/erfinv.hpp new file mode 100644 index 0000000..6096e54 --- /dev/null +++ b/include/util/erfinv.hpp @@ -0,0 +1,112 @@ +#pragma once +#include +#include + +namespace Seldon::Math +{ + +constexpr long double pi = 3.1415926535897932384626433832795028841971693993751L; +constexpr long double sqrt_pi = 1.7724538509055160272981674833411451827975494561224L; + +// Implementation adapted from https://github.com/lakshayg/erfinv same as used in golang math library +template +T erfinv( T x ) +{ + if( x < -1 || x > 1 ) + { + return std::numeric_limits::quiet_NaN(); + } + else if( x == 1.0 ) + { + return std::numeric_limits::infinity(); + } + else if( x == -1.0 ) + { + return -std::numeric_limits::infinity(); + } + + const T LN2 = 6.931471805599453094172321214581e-1L; + + const T A0 = 1.1975323115670912564578e0L; + const T A1 = 4.7072688112383978012285e1L; + const T A2 = 6.9706266534389598238465e2L; + const T A3 = 4.8548868893843886794648e3L; + const T A4 = 1.6235862515167575384252e4L; + const T A5 = 2.3782041382114385731252e4L; + const T A6 = 1.1819493347062294404278e4L; + const T A7 = 8.8709406962545514830200e2L; + + const T B0 = 1.0000000000000000000e0L; + const T B1 = 4.2313330701600911252e1L; + const T B2 = 6.8718700749205790830e2L; + const T B3 = 5.3941960214247511077e3L; + const T B4 = 2.1213794301586595867e4L; + const T B5 = 3.9307895800092710610e4L; + const T B6 = 2.8729085735721942674e4L; + const T B7 = 5.2264952788528545610e3L; + + const T C0 = 1.42343711074968357734e0L; + const T C1 = 4.63033784615654529590e0L; + const T C2 = 5.76949722146069140550e0L; + const T C3 = 3.64784832476320460504e0L; + const T C4 = 1.27045825245236838258e0L; + const T C5 = 2.41780725177450611770e-1L; + const T C6 = 2.27238449892691845833e-2L; + const T C7 = 7.74545014278341407640e-4L; + + const T D0 = 1.4142135623730950488016887e0L; + const T D1 = 2.9036514445419946173133295e0L; + const T D2 = 2.3707661626024532365971225e0L; + const T D3 = 9.7547832001787427186894837e-1L; + const T D4 = 2.0945065210512749128288442e-1L; + const T D5 = 2.1494160384252876777097297e-2L; + const T D6 = 7.7441459065157709165577218e-4L; + const T D7 = 1.4859850019840355905497876e-9L; + + const T E0 = 6.65790464350110377720e0L; + const T E1 = 5.46378491116411436990e0L; + const T E2 = 1.78482653991729133580e0L; + const T E3 = 2.96560571828504891230e-1L; + const T E4 = 2.65321895265761230930e-2L; + const T E5 = 1.24266094738807843860e-3L; + const T E6 = 2.71155556874348757815e-5L; + const T E7 = 2.01033439929228813265e-7L; + + const T F0 = 1.414213562373095048801689e0L; + const T F1 = 8.482908416595164588112026e-1L; + const T F2 = 1.936480946950659106176712e-1L; + const T F3 = 2.103693768272068968719679e-2L; + const T F4 = 1.112800997078859844711555e-3L; + const T F5 = 2.611088405080593625138020e-5L; + const T F6 = 2.010321207683943062279931e-7L; + const T F7 = 2.891024605872965461538222e-15L; + + T abs_x = std::abs( x ); + + T r, num, den; + + if( abs_x <= 0.85 ) + { + r = 0.180625 - 0.25 * x * x; + num = ( ( ( ( ( ( ( A7 * r + A6 ) * r + A5 ) * r + A4 ) * r + A3 ) * r + A2 ) * r + A1 ) * r + A0 ); + den = ( ( ( ( ( ( ( B7 * r + B6 ) * r + B5 ) * r + B4 ) * r + B3 ) * r + B2 ) * r + B1 ) * r + B0 ); + return x * num / den; + } + + r = std::sqrt( LN2 - std::log1p( -abs_x ) ); + if( r <= 5.0 ) + { + r = r - 1.6; + num = ( ( ( ( ( ( ( C7 * r + C6 ) * r + C5 ) * r + C4 ) * r + C3 ) * r + C2 ) * r + C1 ) * r + C0 ); + den = ( ( ( ( ( ( ( D7 * r + D6 ) * r + D5 ) * r + D4 ) * r + D3 ) * r + D2 ) * r + D1 ) * r + D0 ); + } + else + { + r = r - 5.0; + num = ( ( ( ( ( ( ( E7 * r + E6 ) * r + E5 ) * r + E4 ) * r + E3 ) * r + E2 ) * r + E1 ) * r + E0 ); + den = ( ( ( ( ( ( ( F7 * r + F6 ) * r + F5 ) * r + F4 ) * r + F3 ) * r + F2 ) * r + F1 ) * r + F0 ); + } + + return std::copysign( num / den, x ); +} +} // namespace Seldon::Math \ No newline at end of file diff --git a/include/util/math.hpp b/include/util/math.hpp index bbde5e6..d3a9ec8 100644 --- a/include/util/math.hpp +++ b/include/util/math.hpp @@ -161,12 +161,11 @@ class truncated_normal_distribution ScalarT mean{}; ScalarT sigma{}; ScalarT eps{}; - std::normal_distribution normal_dist{}; - size_t max_tries = 5000; + std::uniform_real_distribution uniform_dist{}; ScalarT inverse_cum_gauss( ScalarT y ) { - return erfinv( 2.0 * y - 1 ) * std::sqrt( 2.0 ) * sigma + mean; + return Math::erfinv( 2.0 * y - 1 ) * std::sqrt( 2.0 ) * sigma + mean; } ScalarT cum_gauss( ScalarT x ) @@ -176,26 +175,19 @@ class truncated_normal_distribution public: truncated_normal_distribution( ScalarT mean, ScalarT sigma, ScalarT eps ) - : mean( mean ), sigma( sigma ), eps( eps ), normal_dist( std::normal_distribution( mean, sigma ) ) + : mean( mean ), sigma( sigma ), eps( eps ), uniform_dist( 0, 1 ) { } template ScalarT operator()( Generator & gen ) { - for( size_t i = 0; i < max_tries; i++ ) - { - auto sample = normal_dist( gen ); - if( sample > eps ) - return sample; - } - return eps; + return inverse_cumulative_probability( uniform_dist( gen ) ); } ScalarT inverse_cumulative_probability( ScalarT y ) { - return inverse_cum_gauss( - y * ( 1.0 - cum_gauss( eps, sigma, mean ) ) + cum_gauss( eps, sigma, mean ), sigma, mean ); + return inverse_cum_gauss( y * ( 1.0 - cum_gauss( eps ) ) + cum_gauss( eps ) ); } }; @@ -222,7 +214,7 @@ class bivariate_normal_distribution ScalarT n2 = normal_dist( gen ); ScalarT r1 = n1; - ScalarT r2 = covariance * n1 + std::sqrt( 1 - covariance * covariance ); + ScalarT r2 = covariance * n1 + std::sqrt( 1.0 - covariance * covariance ) * n2; return { r1, r2 }; } diff --git a/test/test_probability_distributions.cpp b/test/test_probability_distributions.cpp index 74fd6ab..930e4d6 100644 --- a/test/test_probability_distributions.cpp +++ b/test/test_probability_distributions.cpp @@ -3,6 +3,7 @@ #include #include "util/math.hpp" +#include #include #include #include @@ -10,10 +11,10 @@ #include namespace fs = std::filesystem; -template -std::ostream & operator<<( std::ostream & os, std::array const & v1 ) +template +std::ostream & operator<<( std::ostream & os, std::array const & v1 ) { - std::for_each( begin( v1 ), end( v1 ), [&os]( int val ) { os << val << " "; } ); + std::for_each( begin( v1 ), end( v1 ), [&os]( T val ) { os << val << " "; } ); return os; } @@ -22,8 +23,9 @@ template void write_results_to_file( int N_Samples, distT dist, const std::string & filename ) { auto proj_root_path = fs::current_path(); - auto file = proj_root_path / fs::path( "/test/output/" + filename ); - fs::create_directories( file ); + auto file = proj_root_path / fs::path( "test/output_probability_distributions/" + filename ); + fmt::print( "file = {}\n", fmt::streamed( file ) ); + fs::create_directories( file.parent_path() ); auto gen = std::mt19937( 0 ); @@ -39,7 +41,7 @@ void write_results_to_file( int N_Samples, distT dist, const std::string & filen TEST_CASE( "Test the probability distributions", "[prob]" ) { - write_results_to_file( 10000, Seldon::truncated_normal_distribution( 1.0, 0.5, 0.1 ), "truncated_normal.txt" ); + write_results_to_file( 10000, Seldon::truncated_normal_distribution( 1.0, 0.5, 0.75 ), "truncated_normal.txt" ); write_results_to_file( 10000, Seldon::power_law_distribution( 0.01, 2.1 ), "power_law.txt" ); write_results_to_file( 10000, Seldon::bivariate_normal_distribution( 0.5 ), "bivariate_normal.txt" ); } From ea7c900fb54d8ec1e6c1e5e97a4604afb302031e Mon Sep 17 00:00:00 2001 From: Moritz Sallermann Date: Sun, 24 Mar 2024 20:04:35 +0000 Subject: [PATCH 12/17] Math: Add pdfs to distributions Added the probability density functions to the truncated gaussian and the power law distribution. Also renamed the cumulative distribution functions to cdf. Co-authored-by: Amrita Goswami --- include/util/math.hpp | 32 +++++++++++++++++++++++++------- 1 file changed, 25 insertions(+), 7 deletions(-) diff --git a/include/util/math.hpp b/include/util/math.hpp index d3a9ec8..4d49810 100644 --- a/include/util/math.hpp +++ b/include/util/math.hpp @@ -132,10 +132,15 @@ class power_law_distribution template ScalarT operator()( Generator & gen ) { - return inverse_cumulative_probability( dist( gen ) ); + return inverse_cdf( dist( gen ) ); } - ScalarT inverse_cumulative_probability( ScalarT x ) + ScalarT pdf( ScalarT x ) + { + return ( 1.0 - gamma ) / ( 1.0 - std::pow( eps, ( 1 - gamma ) ) * std::pow( x, ( -gamma ) ) ); + } + + ScalarT inverse_cdf( ScalarT x ) { return std::pow( ( 1.0 - std::pow( eps, ( 1.0 - gamma ) ) ) * x + std::pow( eps, ( 1.0 - gamma ) ), @@ -163,16 +168,21 @@ class truncated_normal_distribution ScalarT eps{}; std::uniform_real_distribution uniform_dist{}; - ScalarT inverse_cum_gauss( ScalarT y ) + ScalarT inverse_cdf_gauss( ScalarT y ) { return Math::erfinv( 2.0 * y - 1 ) * std::sqrt( 2.0 ) * sigma + mean; } - ScalarT cum_gauss( ScalarT x ) + ScalarT cdf_gauss( ScalarT x ) { return 0.5 * ( 1 + std::erf( ( x - mean ) / ( sigma * std::sqrt( 2.0 ) ) ) ); } + ScalarT pdf_gauss( ScalarT x ) + { + return 1.0 / ( sigma * std::sqrt( 2 * Math::pi ) ) * std::exp( -0.5 * std::pow( ( ( x - mean ) / sigma ), 2 ) ); + } + public: truncated_normal_distribution( ScalarT mean, ScalarT sigma, ScalarT eps ) : mean( mean ), sigma( sigma ), eps( eps ), uniform_dist( 0, 1 ) @@ -182,12 +192,20 @@ class truncated_normal_distribution template ScalarT operator()( Generator & gen ) { - return inverse_cumulative_probability( uniform_dist( gen ) ); + return inverse_cdf( uniform_dist( gen ) ); + } + + ScalarT inverse_cdf( ScalarT y ) + { + return inverse_cdf_gauss( y * ( 1.0 - cdf_gauss( eps ) ) + cdf_gauss( eps ) ); } - ScalarT inverse_cumulative_probability( ScalarT y ) + ScalarT pdf( ScalarT x ) { - return inverse_cum_gauss( y * ( 1.0 - cum_gauss( eps ) ) + cum_gauss( eps ) ); + if( x < eps ) + return 0.0; + else + return 1.0 / ( 1.0 - cdf_gauss( eps ) ) * pdf_gauss( x ); } }; From 376bb903247016ef0d4eb293c35bd38a312b4c1e Mon Sep 17 00:00:00 2001 From: Moritz Sallermann Date: Sun, 24 Mar 2024 20:20:33 +0000 Subject: [PATCH 13/17] Test: Gaussian copula Fixed issues with bivariate_gaussian_copula, due to previous renaming of cdf. Also put it into the same unit test and plotted the histograms. Co-authored-by: Amrita Goswami --- include/util/math.hpp | 15 ++++++------ test/test_probability_distributions.cpp | 31 ++++++------------------- 2 files changed, 14 insertions(+), 32 deletions(-) diff --git a/include/util/math.hpp b/include/util/math.hpp index 4d49810..493e6ab 100644 --- a/include/util/math.hpp +++ b/include/util/math.hpp @@ -243,11 +243,11 @@ class bivariate_gaussian_copula { private: ScalarT covariance; - bivariate_normal_distribution bivariate_normal_dist{}; + bivariate_normal_distribution biv_normal_dist{}; // std::normal_distribution normal_dist{}; // Cumulative probability function for gaussian with mean 0 and variance 1 - ScalarT cum_gauss( ScalarT x ) + ScalarT cdf_gauss( ScalarT x ) { return 0.5 * ( 1 + std::erf( ( x ) / std::sqrt( 2.0 ) ) ); } @@ -258,9 +258,9 @@ class bivariate_gaussian_copula public: bivariate_gaussian_copula( ScalarT covariance, dist1T dist1, dist2T dist2 ) : covariance( covariance ), + biv_normal_dist( bivariate_normal_distribution( covariance ) ), dist1( dist1 ), - dist2( dist2 ), - bivariate_normal_dist( bivariate_normal_dist( covariance ) ) + dist2( dist2 ) { } @@ -268,12 +268,11 @@ class bivariate_gaussian_copula std::array operator()( Generator & gen ) { // 1. Draw from bivariate gaussian - auto z = bivariate_normal_dist( gen ); + auto z = biv_normal_dist( gen ); // 2. Transform marginals to unit interval - std::array z_unit = { cum_gauss( z[0] ), cum_gauss( z[1] ) }; + std::array z_unit = { cdf_gauss( z[0] ), cdf_gauss( z[1] ) }; // 3. Apply inverse transform sampling - std::array res - = { dist1.inverse_cumulative_probability( z_unit[0] ), dist2.inverse_cumulative_probability( z_unit[1] ) }; + std::array res = { dist1.inverse_cdf( z_unit[0] ), dist2.inverse_cdf( z_unit[1] ) }; return res; } }; diff --git a/test/test_probability_distributions.cpp b/test/test_probability_distributions.cpp index 930e4d6..e72a3df 100644 --- a/test/test_probability_distributions.cpp +++ b/test/test_probability_distributions.cpp @@ -46,27 +46,10 @@ TEST_CASE( "Test the probability distributions", "[prob]" ) write_results_to_file( 10000, Seldon::bivariate_normal_distribution( 0.5 ), "bivariate_normal.txt" ); } -// TEST_CASE( "Test reading in the agents from a file", "[io_agents]" ) -// { -// using namespace Seldon; -// using namespace Catch::Matchers; - -// auto proj_root_path = fs::current_path(); -// auto network_file = proj_root_path / fs::path( "test/res/opinions.txt" ); - -// auto agents = Seldon::agents_from_file( network_file ); - -// std::vector opinions_expected = { 2.1127107987061544, 0.8088982488089491, -0.8802809369462433 }; -// std::vector activities_expected = { 0.044554683389757696, 0.015813166022685163, 0.015863953902810535 }; -// std::vector reluctances_expected = { 1.0, 1.0, 2.3 }; - -// REQUIRE( agents.size() == 3 ); - -// for( size_t i = 0; i < agents.size(); i++ ) -// { -// fmt::print( "{}", i ); -// REQUIRE_THAT( agents[i].data.opinion, Catch::Matchers::WithinAbs( opinions_expected[i], 1e-16 ) ); -// REQUIRE_THAT( agents[i].data.activity, Catch::Matchers::WithinAbs( activities_expected[i], 1e-16 ) ); -// REQUIRE_THAT( agents[i].data.reluctance, Catch::Matchers::WithinAbs( reluctances_expected[i], 1e-16 ) ); -// } -// } \ No newline at end of file +TEST_CASE( "Test bivariate gaussian copula", "[prob_copula]" ) +{ + auto dist1 = Seldon::power_law_distribution( 0.02, 2.5 ); + auto dist2 = Seldon::truncated_normal_distribution( 1.0, 0.75, 0.2 ); + auto copula = Seldon::bivariate_gaussian_copula( 0.5, dist1, dist2 ); + write_results_to_file( 10000, copula, "gaussian_copula.txt" ); +} From ade8ee976057d6e07a739256c0b7cbdadaecfa25 Mon Sep 17 00:00:00 2001 From: Moritz Sallermann Date: Sun, 24 Mar 2024 20:31:59 +0000 Subject: [PATCH 14/17] ActivityDrivenModel: Use the copula Use the copula to generate correlated activities and reluctances. Co-authored-by: Amrita Goswami --- include/connectivity.hpp | 2 +- include/models/ActivityDrivenModel.hpp | 20 ++++++++++---------- src/config_parser.cpp | 3 +-- test/test_network.cpp | 2 +- 4 files changed, 13 insertions(+), 14 deletions(-) diff --git a/include/connectivity.hpp b/include/connectivity.hpp index 1af9c8d..a9d0159 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 1b0ab83..7418e4c 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; @@ -113,6 +113,8 @@ class ActivityDrivenModelAbstract : public Model power_law_distribution<> dist_activity( eps, gamma ); truncated_normal_distribution<> dist_reluctance( reluctance_mean, reluctance_sigma, reluctance_eps ); + bivariate_gaussian_copula copula( covariance_factor, dist_activity, dist_reluctance ); + auto mean_activity = dist_activity.mean(); // Initial conditions for the opinions, initialize to [-1,1] @@ -121,18 +123,16 @@ class ActivityDrivenModelAbstract : public Model { 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; - } + auto res = copula( gen ); + network.agents[i].data.activity = res[0]; if( use_reluctances ) { - network.agents[i].data.reluctance = dist_reluctance( gen ); + network.agents[i].data.reluctance = res[1]; + } + if( mean_activities ) + { + network.agents[i].data.activity = mean_activity; } } diff --git a/src/config_parser.cpp b/src/config_parser.cpp index e0e054e..63c8456 100644 --- a/src/config_parser.cpp +++ b/src/config_parser.cpp @@ -253,8 +253,7 @@ 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 57239e7..f1d2bfa 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 } From 741b395d2912042ab0725f09286fb1c9e16280ee Mon Sep 17 00:00:00 2001 From: Amrita Goswami Date: Mon, 25 Mar 2024 16:51:46 +0000 Subject: [PATCH 15/17] Copula: Complete input parsing Co-authored-by: Moritz Sallermann --- examples/ActivityDrivenReluctance/conf.toml | 3 ++- include/models/ActivityDrivenModel.hpp | 1 + src/config_parser.cpp | 3 ++- 3 files changed, 5 insertions(+), 2 deletions(-) diff --git a/examples/ActivityDrivenReluctance/conf.toml b/examples/ActivityDrivenReluctance/conf.toml index 668891b..52b1005 100644 --- a/examples/ActivityDrivenReluctance/conf.toml +++ b/examples/ActivityDrivenReluctance/conf.toml @@ -24,8 +24,9 @@ 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_sigma = 0.15 # Width of normal distribution (before truncating) reluctance_eps = 0.01 # Minimum such that the normal distribution is truncated at this value +covariance_factor = 0.0 # 0.0 means that the reluctances and activities are uncorrelated. Should be in the range of [-1,1] [network] number_of_agents = 1000 diff --git a/include/models/ActivityDrivenModel.hpp b/include/models/ActivityDrivenModel.hpp index 7418e4c..64a1223 100644 --- a/include/models/ActivityDrivenModel.hpp +++ b/include/models/ActivityDrivenModel.hpp @@ -44,6 +44,7 @@ class ActivityDrivenModelAbstract : public Model reluctance_mean( settings.reluctance_mean ), reluctance_sigma( settings.reluctance_sigma ), reluctance_eps( settings.reluctance_eps ), + covariance_factor( settings.covariance_factor ), n_bots( settings.n_bots ), bot_m( settings.bot_m ), bot_activity( settings.bot_activity ), diff --git a/src/config_parser.cpp b/src/config_parser.cpp index 63c8456..ac48587 100644 --- a/src/config_parser.cpp +++ b/src/config_parser.cpp @@ -59,6 +59,7 @@ void parse_activity_settings( auto & model_settings, const auto & toml_model_opt 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.covariance_factor, toml_model_opt["covariance_factor"] ); 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"] ); @@ -213,7 +214,7 @@ void validate_settings( const SimulationOptions & options ) check( name_and_var( model_settings.reluctance_mean ), g_zero ); check( name_and_var( model_settings.reluctance_sigma ), g_zero ); check( name_and_var( model_settings.reluctance_eps ), g_zero ); - check( name_and_var( model_settings.covariance_factor ), geq_zero ); + check( name_and_var( model_settings.covariance_factor ), []( auto x ) { return x >= -1.0 && x<=1.0; } ); // Bot options size_t n_bots = model_settings.n_bots; auto check_bot_size = [&]( auto x ) { return x.size() >= n_bots; }; From 0ab1575f5ba6ecf53c3896276397a2dca18e6a82 Mon Sep 17 00:00:00 2001 From: Moritz Sallermann Date: Mon, 25 Mar 2024 16:57:24 +0000 Subject: [PATCH 16/17] config_parser: Formatting --- src/config_parser.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/config_parser.cpp b/src/config_parser.cpp index ac48587..a492e49 100644 --- a/src/config_parser.cpp +++ b/src/config_parser.cpp @@ -214,7 +214,7 @@ void validate_settings( const SimulationOptions & options ) check( name_and_var( model_settings.reluctance_mean ), g_zero ); check( name_and_var( model_settings.reluctance_sigma ), g_zero ); check( name_and_var( model_settings.reluctance_eps ), g_zero ); - check( name_and_var( model_settings.covariance_factor ), []( auto x ) { return x >= -1.0 && x<=1.0; } ); + check( name_and_var( model_settings.covariance_factor ), []( auto x ) { return x >= -1.0 && x <= 1.0; } ); // Bot options size_t n_bots = model_settings.n_bots; auto check_bot_size = [&]( auto x ) { return x.size() >= n_bots; }; From 78f1c92cbe8f69482239a61d49c6c866d65563c7 Mon Sep 17 00:00:00 2001 From: Amrita Goswami Date: Wed, 10 Apr 2024 08:11:45 +0000 Subject: [PATCH 17/17] DeffuantVector: Fixed bug in main.cpp The multi-dimensional DeffuantVector model has a different AgentT from the regular Deffuant model. Therefore, the model and simulation have to be built with the DeffuantVector agent. This requires some repititive logic in main.cpp. Also fixed a small typo in the TOML input file for the DeffuantVector model example. --- examples/DeffuantVector/conf.toml | 2 +- src/main.cpp | 13 +++++++++++-- 2 files changed, 12 insertions(+), 3 deletions(-) diff --git a/examples/DeffuantVector/conf.toml b/examples/DeffuantVector/conf.toml index 604a40e..1d5f0ca 100644 --- a/examples/DeffuantVector/conf.toml +++ b/examples/DeffuantVector/conf.toml @@ -5,7 +5,7 @@ model = "Deffuant" [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 prints +print_progress = true # Print the iteration time ; if not set, then does not print output_initial = true # Print the initial opinions and network file from step 0. If not set, this is true by default. start_output = 1 # Start writing out opinions and/or network files from this iteration. If not set, this is 1. diff --git a/src/main.cpp b/src/main.cpp index 3f429fa..ac57310 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -78,8 +78,17 @@ int main( int argc, char * argv[] ) } else if( simulation_options.model == Seldon::Config::Model::DeffuantModel ) { - simulation = std::make_unique>( - simulation_options, network_file, agent_file ); + auto model_settings = std::get( simulation_options.model_settings ); + if( model_settings.use_binary_vector ) + { + simulation = std::make_unique>( + simulation_options, network_file, agent_file ); + } + else + { + simulation = std::make_unique>( + simulation_options, network_file, agent_file ); + } } else {