diff --git a/examples/ActivityDriven/conf.toml b/examples/ActivityDriven/conf.toml new file mode 100644 index 0000000..699eedf --- /dev/null +++ b/examples/ActivityDriven/conf.toml @@ -0,0 +1,20 @@ +[simulation] +model = "ActivityDriven" +# rng_seed = 120 # Leaving this empty will pick a random seed + +[model] +max_iterations = 20 # If not set, max iterations is infinite + +[ActivityDriven] +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.5 # 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 = 3.0 # Controversialness of the issue, must be greater than 0. +K = 3.0 # Social interaction strength + +[network] +number_of_agents = 1000 +connections_per_agent = 10 diff --git a/include/agent.hpp b/include/agent.hpp index c2a5200..7c41ed8 100644 --- a/include/agent.hpp +++ b/include/agent.hpp @@ -4,15 +4,17 @@ namespace Seldon { -/* Agent is a base class for actual agents with opinion type T, it needs to implement to_string and from_string*/ +/* Agent is a base class for actual agents with data type T +(which contains an opinion and perhaps some other things), +it needs to implement to_string and from_string*/ template class Agent : public AgentBase { public: - using opinion_t = T; - opinion_t opinion; + using data_t = T; + data_t data; Agent() = default; - Agent( opinion_t opinion ) : opinion( opinion ) {} + Agent( data_t data ) : data( data ) {} void from_string( const std::string & str ); @@ -25,7 +27,7 @@ class Agent : public AgentBase template<> inline void Agent::from_string( const std::string & str ) { - opinion = std::stod( str ); + data = std::stod( str ); } } // namespace Seldon \ No newline at end of file diff --git a/include/model.hpp b/include/model.hpp index 5b44e5e..4a8967f 100644 --- a/include/model.hpp +++ b/include/model.hpp @@ -19,7 +19,7 @@ class Model : public ModelBase Model( size_t n_agents ) : agents( std::vector( int( n_agents ), AgentT() ) ) {} Model( std::vector && agents ) : agents( agents ) {} - void Agents_from_File( const std::string & file ) + void agents_from_file( const std::string & file ) override { agents.clear(); @@ -39,7 +39,7 @@ class Model : public ModelBase auto line = file_contents.substr( start_of_line, end_of_line - start_of_line ); start_of_line = end_of_line + 1; // TODO: check if empty or comment - if( line.size() == 0 ) + if( line.empty() ) { break; } @@ -57,7 +57,10 @@ class Model : public ModelBase } } - void iteration() override = 0; + void iteration() override + { + n_iterations++; + }; bool finished() override { @@ -77,10 +80,4 @@ class Model : public ModelBase } }; -template -void Seldon::Model::iteration() -{ - n_iterations++; -}; - } // namespace Seldon \ No newline at end of file diff --git a/include/model_base.hpp b/include/model_base.hpp index 86159ae..7612031 100644 --- a/include/model_base.hpp +++ b/include/model_base.hpp @@ -16,7 +16,8 @@ class ModelBase virtual AgentBase * get_agent( int idx ) = 0; // Use this to get an abstract representation of the agent at idx virtual void iteration() = 0; virtual bool finished() = 0; - virtual ~ModelBase() = default; + virtual void agents_from_file( const std::string & file ) = 0; + virtual ~ModelBase() = default; }; } // namespace Seldon \ No newline at end of file diff --git a/include/models/ActivityDrivenModel.hpp b/include/models/ActivityDrivenModel.hpp new file mode 100644 index 0000000..6518fa3 --- /dev/null +++ b/include/models/ActivityDrivenModel.hpp @@ -0,0 +1,104 @@ +#pragma once + +#include "agent.hpp" +#include "model.hpp" +#include "network.hpp" +#include +#include +#include +#include +#include +#include + +namespace Seldon +{ + +struct ActivityAgentData +{ + double opinion = 0; // x_i + double activity = 0; // a_i +}; + +template<> +inline std::string Agent::to_string() const +{ + return fmt::format( "{}, {}", data.opinion, data.activity ); +}; + +template<> +inline void Agent::from_string( const std::string & str ) +{ + auto pos_comma = str.find_first_of( ',' ); + data.opinion = std::stod( str.substr( 0, pos_comma ) ); + data.activity = std::stod( str.substr( pos_comma + 1, str.size() ) ); +}; + +class ActivityAgentModel : public Model> +{ + using AgentT = Agent; + +private: + double max_opinion_diff = 0; + Network & network; + std::vector agents_current_copy; + // Random number generation + std::mt19937 & gen; // reference to simulation Mersenne-Twister engine + std::set> reciprocal_edge_buffer{}; + + // Buffers for RK4 integration + std::vector k1_buffer{}; + std::vector k2_buffer{}; + std::vector k3_buffer{}; + std::vector k4_buffer{}; + + template + void get_euler_slopes( std::vector & k_buffer, Opinion_Callback opinion ) + { + // h is the timestep + auto neighbour_buffer = std::vector(); + size_t j_index = 0; + + k_buffer.resize( network.n_agents() ); + + for( size_t idx_agent = 0; idx_agent < network.n_agents(); ++idx_agent ) + { + network.get_neighbours( idx_agent, neighbour_buffer ); // Get the incoming neighbours + 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] += K * std::tanh( alpha * opinion( j_index ) ); + } + // Multiply by the timestep + k_buffer[idx_agent] *= dt; + } + } + +public: + // Model-specific parameters + double dt = 0.01; // Timestep for the integration of the coupled ODEs + // Various free parameters + int m = 10; // Number of agents contacted, when the agent is active + double eps = 0.01; // Minimum activity epsilon; a_i belongs to [epsilon,1] + double gamma = 2.1; // Exponent of activity power law distribution of activities + double alpha = 3.0; // Controversialness of the issue, must be greater than 0. + double homophily = 0.5; // aka beta. if zero, agents pick their interaction partners at random + // Reciprocity aka r. 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 + double reciprocity = 0.5; + double K = 3.0; // Social interaction strength; K>0 + + double convergence_tol = 1e-12; // TODO: ?? + + ActivityAgentModel( int n_agents, Network & network, std::mt19937 & gen ); + + void get_agents_from_power_law(); // This needs to be called after eps and gamma have been set + + void iteration() override; + + // bool finished() overteration() override; + // bool finished() override; +}; + +} // namespace Seldon \ No newline at end of file diff --git a/include/network.hpp b/include/network.hpp index f433b63..fb602be 100644 --- a/include/network.hpp +++ b/include/network.hpp @@ -14,9 +14,21 @@ class Network Network( std::vector> && neighbour_list, std::vector> && weight_list ); std::size_t n_agents() const; + void get_neighbours( std::size_t agent_idx, std::vector & buffer ) const; void get_weights( std::size_t agent_idx, std::vector & buffer ) const; + void set_neighbours_and_weights( + std::size_t agent_idx, const std::vector & buffer_neighbours, const WeightT & weight ); + + void set_neighbours_and_weights( + std::size_t agent_idx, const std::vector & buffer_neighbours, + const std::vector & buffer_weights ); + + void push_back_neighbour_and_weight( size_t i, size_t j, WeightT w ); + + void transpose(); + private: std::vector> neighbour_list; // Neighbour list for the connections std::vector> weight_list; // List for the interaction weights of each connection diff --git a/include/util/math.hpp b/include/util/math.hpp new file mode 100644 index 0000000..a69b575 --- /dev/null +++ b/include/util/math.hpp @@ -0,0 +1,132 @@ +#pragma once +#include +#include +#include +#include +#include +#include + +namespace Seldon +{ + +// Function for getting a vector of k agents (corresponding to connections) +// drawing from n agents (without duplication) +// ignore_idx ignores the index of the agent itself, since we will later add the agent itself ourselves to prevent duplication +inline void draw_unique_k_from_n( + std::size_t ignore_idx, std::size_t k, std::size_t n, std::vector & buffer, std::mt19937 & gen ) +{ + struct SequenceGenerator + { + /* An iterator that generates a sequence of integers 2, 3, 4 ...*/ + using iterator_category = std::forward_iterator_tag; + using difference_type = std::ptrdiff_t; + using value_type = size_t; + using pointer = size_t *; // or also value_type* + using reference = size_t &; + + SequenceGenerator( const size_t i_, const size_t ignore_idx ) : i( i_ ), ignore_idx( ignore_idx ) + { + if( i == ignore_idx ) + { + i++; + } + } + size_t i; + size_t ignore_idx; + + size_t & operator*() + { + return i; + }; + bool operator==( const SequenceGenerator & it1 ) + { + return i == it1.i; + }; + SequenceGenerator & operator++() + { + i++; + if( i == ignore_idx ) + i++; + return *this; + } + }; + + buffer.resize( k ); + std::sample( SequenceGenerator( 0, ignore_idx ), SequenceGenerator( n, ignore_idx ), buffer.begin(), k, gen ); +} + +template +void reservoir_sampling_A_ExpJ( + size_t k, size_t n, WeightCallbackT weight, std::vector & buffer, std::mt19937 & mt ) +{ + std::uniform_real_distribution distribution( 0.0, 1.0 ); + + std::vector reservoir( k ); + using QueueItemT = std::pair; + + auto compare = []( const QueueItemT & item1, const QueueItemT & item2 ) { return item1.second > item2.second; }; + std::priority_queue, decltype( compare )> H; + + size_t idx = 0; + while( idx < n & H.size() < k ) + { + double r = std::pow( distribution( mt ), 1.0 / weight( idx ) ); + H.push( { idx, r } ); + idx++; + } + + auto X = std::log( distribution( mt ) ) / std::log( H.top().second ); + while( idx < n ) + { + auto w = weight( idx ); + X -= w; + if( X <= 0 ) + { + auto t = std::pow( H.top().second, w ); + auto uniform_from_t_to_one = distribution( mt ) * ( 1.0 - t ) + t; // Random number in interval [t, 1.0] + auto r = std::pow( uniform_from_t_to_one, 1.0 / w ); + H.pop(); + H.push( { idx, r } ); + X = std::log( distribution( mt ) ) / std::log( H.top().second ); + } + idx++; + } + + buffer.resize( H.size() ); + + for( size_t i = 0; i < k; i++ ) + { + buffer[i] = H.top().first; + H.pop(); + } +} + +/** + * @brief Power law distribution for random numbers. + * A continuous random distribution on the range [eps, infty) + * with p(x) ~ x^(-gamma) + * Including normalization the PDF is + * p(x) = (1-gamma)/(1-eps^(1-gamma)) * x^(-gamma) + */ +template +class power_law_distribution +{ +private: + ScalarT eps; + ScalarT gamma; + std::uniform_real_distribution dist + = std::uniform_real_distribution( 0.0, 1.0 ); // Uniform random variable for activities + +public: + power_law_distribution( ScalarT eps, ScalarT gamma ) : eps( eps ), gamma( gamma ) {} + + template + ScalarT operator()( Generator & gen ) + { + return std::pow( + ( 1.0 - std::pow( eps, ( 1.0 - gamma ) ) ) * dist( gen ) + std::pow( eps, ( 1.0 - gamma ) ), + ( 1.0 / ( 1.0 - gamma ) ) ); + } +}; + +} // namespace Seldon \ No newline at end of file diff --git a/meson.build b/meson.build index ff6d5c2..cd58a5b 100644 --- a/meson.build +++ b/meson.build @@ -10,6 +10,7 @@ sources_seldon = [ 'src/network_generation.cpp', 'src/model.cpp', 'src/models/DeGroot.cpp', + 'src/models/ActivityDrivenModel.cpp', 'src/util/tomlplusplus.cpp' ] @@ -22,7 +23,9 @@ exe = executable('seldon', sources_seldon + 'src/main.cpp', tests = [ ['Test Tarjan', 'test/test_tarjan.cpp'], - ['Test DeGroot', 'test/test_deGroot.cpp'] + ['Test DeGroot', 'test/test_deGroot.cpp'], + ['Test Network', 'test/test_network.cpp'], + ['Test Sampling', 'test/test_sampling.cpp'], ] Catch2 = dependency('Catch2', method : 'cmake', modules : ['Catch2::Catch2WithMain', 'Catch2::Catch2']) diff --git a/src/main.cpp b/src/main.cpp index 8006cac..71483c1 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -41,6 +41,9 @@ int main( int argc, char * argv[] ) fs::create_directories( output_dir_path ); // Create the output directory auto simulation = Seldon::Simulation( config_file_path.string(), network_file, agent_file ); + + fmt::print( "{}", simulation.model->get_agent( 0 )->to_string() ); + Seldon::IO::network_to_dot_file( *simulation.network, ( output_dir_path / fs::path( "network.dot" ) ).string() ); Seldon::IO::network_to_file( simulation, ( output_dir_path / fs::path( "network.txt" ) ).string() ); auto filename = fmt::format( "opinions_{}.txt", 0 ); diff --git a/src/models/ActivityDrivenModel.cpp b/src/models/ActivityDrivenModel.cpp new file mode 100644 index 0000000..eb37181 --- /dev/null +++ b/src/models/ActivityDrivenModel.cpp @@ -0,0 +1,118 @@ +#include "models/ActivityDrivenModel.hpp" +#include "util/math.hpp" +#include +#include +#include + +Seldon::ActivityAgentModel::ActivityAgentModel( int n_agents, Network & network, std::mt19937 & gen ) + : Model( n_agents ), + network( network ), + agents_current_copy( std::vector( n_agents ) ), + gen( gen ) +{ +} + +void Seldon::ActivityAgentModel::get_agents_from_power_law() +{ + std::uniform_real_distribution<> dis_opinion( -1, 1 ); // Opinion initial values + power_law_distribution<> dist_activity( eps, gamma ); + + // 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 < agents.size(); i++ ) + { + agents[i].data.opinion = dis_opinion( gen ); // Draw the opinion value + // Draw from a power law distribution (1-gamma)/(1-eps^(1-gamma)) * a^(-gamma) + agents[i].data.activity = dist_activity( gen ); + } +} + +void Seldon::ActivityAgentModel::iteration() +{ + Model::iteration(); + + std::uniform_real_distribution<> dis_activation( 0.0, 1.0 ); // Opinion initial values + std::uniform_real_distribution<> dis_reciprocation( 0.0, 1.0 ); // Opinion initial values + 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 ) < 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 + auto weight_callback = [idx_agent, this]( size_t j ) + { + if( idx_agent == j ) // The agent does not contact itself + return 0.0; + return std::pow( + std::abs( this->agents[idx_agent].data.opinion - this->agents[j].data.opinion ), -this->homophily ); + }; + + Seldon::reservoir_sampling_A_ExpJ( m, network.n_agents(), weight_callback, 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 + network.get_neighbours( idx_agent, contacted_agents ); + // 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.transpose(); // transpose the network, so that we have incoming edges + + // 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 this->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 this->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 this->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 this->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) + 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; + } +} \ No newline at end of file diff --git a/src/models/DeGroot.cpp b/src/models/DeGroot.cpp index 153d6e9..d8954b2 100644 --- a/src/models/DeGroot.cpp +++ b/src/models/DeGroot.cpp @@ -7,7 +7,7 @@ Seldon::DeGrootModel::DeGrootModel( int n_agents, Network & network ) { for( size_t i = 0; i < agents.size(); i++ ) { - agents[i].opinion = double( i ) / double( agents.size() ); + agents[i].data = double( i ) / double( agents.size() ); } } @@ -24,12 +24,12 @@ void Seldon::DeGrootModel::iteration() { network.get_neighbours( i, neighbour_buffer ); network.get_weights( i, weight_buffer ); - agents_current_copy[i].opinion = 0.0; + agents_current_copy[i].data = 0.0; for( size_t j = 0; j < neighbour_buffer.size(); j++ ) { j_index = neighbour_buffer[j]; weight = weight_buffer[j]; - agents_current_copy[i].opinion += weight * agents[j_index].opinion; + agents_current_copy[i].data += weight * agents[j_index].data; } } @@ -37,7 +37,7 @@ void Seldon::DeGrootModel::iteration() // Update the original agent opinions for( std::size_t i = 0; i < agents.size(); i++ ) { - max_opinion_diff = std::max( max_opinion_diff, std::abs( agents[i].opinion - agents_current_copy[i].opinion ) ); + max_opinion_diff = std::max( max_opinion_diff, std::abs( agents[i].data - agents_current_copy[i].data ) ); agents[i] = agents_current_copy[i]; } } diff --git a/src/network.cpp b/src/network.cpp index 721809a..d693540 100644 --- a/src/network.cpp +++ b/src/network.cpp @@ -3,6 +3,7 @@ #include #include #include +#include Seldon::Network::Network( std::vector> && neighbour_list, std::vector> && weight_list ) @@ -28,22 +29,61 @@ size_t Seldon::Network::n_agents() const void Seldon::Network::get_neighbours( std::size_t agent_idx, std::vector & buffer ) const { - // TODO: rewrite this using std::span - const size_t n_edges = neighbour_list[agent_idx].size(); - buffer.resize( n_edges ); - for( size_t i_edge = 0; i_edge < n_edges; i_edge++ ) + buffer = neighbour_list[agent_idx]; +} + +void Seldon::Network::set_neighbours_and_weights( + std::size_t agent_idx, const std::vector & buffer_neighbours, const WeightT & weight ) +{ + neighbour_list[agent_idx] = buffer_neighbours; + + weight_list[agent_idx].resize( buffer_neighbours.size() ); + for( auto w : weight_list[agent_idx] ) + { + w = weight; + } +} + +void Seldon::Network::set_neighbours_and_weights( + std::size_t agent_idx, const std::vector & buffer_neighbours, + const std::vector & buffer_weights ) +{ + if( buffer_neighbours.size() != buffer_weights.size() ) [[unlikely]] { - buffer[i_edge] = neighbour_list[agent_idx][i_edge]; + throw std::runtime_error( "Network::set_neighbours_and_weights: both buffers need to have the same length!" ); } + + neighbour_list[agent_idx] = buffer_neighbours; + weight_list[agent_idx] = buffer_weights; } -void Seldon::Network::get_weights( std::size_t agent_idx, std::vector & buffer ) const +void Seldon::Network::push_back_neighbour_and_weight( size_t i, size_t j, WeightT w ) { - // TODO: rewrite this using std::span - const size_t n_edges = weight_list[agent_idx].size(); - buffer.resize( n_edges ); - for( size_t i_edge = 0; i_edge < n_edges; i_edge++ ) + neighbour_list[i].push_back( j ); + weight_list[i].push_back( w ); +} + +void Seldon::Network::get_weights( std::size_t agent_idx, std::vector & buffer ) const +{ + buffer = weight_list[agent_idx]; +} + +void Seldon::Network::transpose() +{ + std::vector> neighbour_list_transpose( n_agents(), std::vector( 0 ) ); + std::vector> weight_list_transpose( n_agents(), std::vector( 0 ) ); + + for( size_t i_agent = 0; i_agent < n_agents(); i_agent++ ) { - buffer[i_edge] = weight_list[agent_idx][i_edge]; + for( size_t i_neighbour = 0; i_neighbour < neighbour_list[i_agent].size(); i_neighbour++ ) + { + const auto neighbour = neighbour_list[i_agent][i_neighbour]; + const auto weight = weight_list[i_agent][i_neighbour]; + neighbour_list_transpose[neighbour].push_back( i_agent ); + weight_list_transpose[neighbour].push_back( weight ); + } } + + neighbour_list = std::move( neighbour_list_transpose ); + weight_list = std::move( weight_list_transpose ); } \ No newline at end of file diff --git a/src/network_generation.cpp b/src/network_generation.cpp index a8c9a88..0580068 100644 --- a/src/network_generation.cpp +++ b/src/network_generation.cpp @@ -8,54 +8,9 @@ #include #include #include +#include #include -// Function for getting a vector of k agents (corresponding to connections) -// drawing from n agents (without duplication) -// ignore_idx ignores the index of the agent itself, since we will later add the agent itself ourselves to prevent duplication -void draw_unique_k_from_n( - std::size_t ignore_idx, std::size_t k, std::size_t n, std::vector & buffer, std::mt19937 & gen ) -{ - struct SequenceGenerator - { - /* An iterator that generates a sequence of integers 2, 3, 4 ...*/ - using iterator_category = std::forward_iterator_tag; - using difference_type = std::ptrdiff_t; - using value_type = size_t; - using pointer = size_t *; // or also value_type* - using reference = size_t &; - - SequenceGenerator( const size_t i_, const size_t ignore_idx ) : i( i_ ), ignore_idx( ignore_idx ) - { - if( i == ignore_idx ) - { - i++; - } - } - size_t i; - size_t ignore_idx; - - size_t & operator*() - { - return i; - }; - bool operator==( const SequenceGenerator & it1 ) - { - return i == it1.i; - }; - SequenceGenerator & operator++() - { - i++; - if( i == ignore_idx ) - i++; - return *this; - } - }; - - buffer.resize( k ); - std::sample( SequenceGenerator( 0, ignore_idx ), SequenceGenerator( n, ignore_idx ), buffer.begin(), k, gen ); -} - std::unique_ptr Seldon::generate_n_connections( int n_agents, int n_connections, std::mt19937 & gen ) { using WeightT = Network::WeightT; @@ -78,7 +33,7 @@ std::unique_ptr Seldon::generate_n_connections( int n_agents, i // Get the vector of sorted adjacencies, excluding i (include i later) // TODO: option for making the n_conections variable - draw_unique_k_from_n( i_agent, n_connections, n_agents, incoming_neighbour_buffer, gen ); + Seldon::draw_unique_k_from_n( i_agent, n_connections, n_agents, incoming_neighbour_buffer, gen ); incoming_neighbour_weights.resize( incoming_neighbour_buffer.size() ); for( size_t j = 0; j < incoming_neighbour_buffer.size(); ++j ) diff --git a/src/simulation.cpp b/src/simulation.cpp index f0238f1..0878403 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -1,4 +1,5 @@ #include "simulation.hpp" +#include "models/ActivityDrivenModel.hpp" #include "models/DeGroot.hpp" #include "network_generation.hpp" #include "util/tomlplusplus.hpp" @@ -9,15 +10,10 @@ #include #include -enum class ModelType : unsigned int -{ - DeGroot -}; - Seldon::Simulation::Simulation( std::string config_file, std::optional cli_network_file, std::optional cli_agent_file ) { - std::set allowed_models = { "DeGroot" }; + std::set allowed_models = { "DeGroot", "ActivityDriven" }; toml::table tbl; tbl = toml::parse_file( config_file ); @@ -48,6 +44,8 @@ Seldon::Simulation::Simulation( throw std::runtime_error( fmt::format( "Unknown model type: '{}'!", model_string ) ); } + fmt::print( "Model type: {}\n", model_string ); + // Construct the network std::optional file = cli_network_file; if( !file.has_value() ) // Check if toml file should be superceded by cli_network_file @@ -73,7 +71,6 @@ Seldon::Simulation::Simulation( // Generic model parameters std::optional max_iterations = tbl["model"]["max_iterations"].value(); - ModelType model_type; if( model_string == "DeGroot" ) { // DeGroot specific parameters @@ -82,15 +79,29 @@ Seldon::Simulation::Simulation( auto model_DeGroot = std::make_unique( n_agents, *network ); model_DeGroot->max_iterations = max_iterations; model_DeGroot->convergence_tol = convergence; + model = std::move( model_DeGroot ); + } + else if( model_string == "ActivityDriven" ) + { + auto model_activityDriven = std::make_unique( n_agents, *network, gen ); + model_activityDriven->dt = tbl["ActivityDriven"]["dt"].value_or( 0.01 ); + model_activityDriven->m = tbl["ActivityDriven"]["m"].value_or( 10 ); + model_activityDriven->eps = tbl["ActivityDriven"]["eps"].value_or( 0.01 ); + model_activityDriven->gamma = tbl["ActivityDriven"]["gamma"].value_or( 2.1 ); + model_activityDriven->homophily = tbl["ActivityDriven"]["homophily"].value_or( 0.5 ); + model_activityDriven->reciprocity = tbl["ActivityDriven"]["reciprocity"].value_or( 0.5 ); + model_activityDriven->alpha = tbl["ActivityDriven"]["alpha"].value_or( 3.0 ); + model_activityDriven->K = tbl["ActivityDriven"]["K"].value_or( 3.0 ); + + model_activityDriven->max_iterations = max_iterations; - // TODO: make this more general, ivoke on ModelBase or Model? - if( cli_agent_file.has_value() ) - { - fmt::print( "Reading agents from file {}\n", cli_agent_file.value() ); - model_DeGroot->Agents_from_File( cli_agent_file.value() ); - } + model_activityDriven->get_agents_from_power_law(); + model = std::move( model_activityDriven ); + } - model = std::move( model_DeGroot ); - model_type = ModelType::DeGroot; + if( cli_agent_file.has_value() ) + { + fmt::print( "Reading agents from file {}\n", cli_agent_file.value() ); + model->agents_from_file( cli_agent_file.value() ); } } \ No newline at end of file diff --git a/test/test_deGroot.cpp b/test/test_deGroot.cpp index 65128c0..781d70c 100644 --- a/test/test_deGroot.cpp +++ b/test/test_deGroot.cpp @@ -21,24 +21,23 @@ TEST_CASE( "Test the DeGroot Model Symmetric", "[DeGroot]" ) { 0.2, 0.8 }, }; - auto gen = std::mt19937(); auto network = Network( std::move( neighbour_list ), std::move( weight_list ) ); auto model = DeGrootModel( n_agents, network ); - model.convergence_tol = 1e-6; - model.max_iterations = 100; - model.agents[0].opinion = 0.0; - model.agents[1].opinion = 1.0; + model.convergence_tol = 1e-6; + model.max_iterations = 100; + model.agents[0].data = 0.0; + model.agents[1].data = 1.0; do { model.iteration(); } while( !model.finished() ); - fmt::print( "N_iterations = {} (with convergence_tol {})\n", model.n_iterations, model.convergence_tol ); + INFO( fmt::format( "N_iterations = {} (with convergence_tol {})\n", model.n_iterations, model.convergence_tol ) ); for( size_t i = 0; i < n_agents; i++ ) { - fmt::print( "Opinion {} = {}\n", i, model.agents[i].opinion ); - REQUIRE_THAT( model.agents[i].opinion, WithinAbs( 0.5, model.convergence_tol * 10.0 ) ); + INFO( fmt::format( "Opinion {} = {}\n", i, model.agents[i].data ) ); + REQUIRE_THAT( model.agents[i].data, WithinAbs( 0.5, model.convergence_tol * 10.0 ) ); } } \ No newline at end of file diff --git a/test/test_network.cpp b/test/test_network.cpp new file mode 100644 index 0000000..36f9fc8 --- /dev/null +++ b/test/test_network.cpp @@ -0,0 +1,93 @@ +#include "network.hpp" +#include "network_generation.hpp" +#include +#include +#include +#include +#include + +TEST_CASE( "Testing the network class" ) +{ + using namespace Seldon; + + // Generate some network + const int n_agents = 20; + const int n_connections = 10; + std::mt19937 gen( 0 ); + auto network = generate_n_connections( n_agents, n_connections, gen ); + + // Does n_agents work? + REQUIRE( network->n_agents() == n_agents ); + + // Check that the function for setting neighbours and a single weight work + // Agent 3 + std::vector buffer_n_get{}; // buffer for getting neighbours + std::vector buffer_w_get{}; // buffer for getting the weights + std::vector neigh{ { 0, 10 } }; // new neighbours + std::vector weight{ 0.5, 0.5 }; // new weights (const) + network->set_neighbours_and_weights( 3, neigh, 0.5 ); + network->get_weights( 3, buffer_w_get ); + REQUIRE_THAT( buffer_w_get, Catch::Matchers::UnorderedRangeEquals( buffer_w_get ) ); + + // Change the connections for agent 3 + std::vector buffer_n{ { 0, 10, 15 } }; // new neighbours + std::vector buffer_w{ 0.1, 0.2, 0.3 }; // new weights + network->set_neighbours_and_weights( 3, buffer_n, buffer_w ); + + // Make sure the changes worked + network->get_neighbours( 3, buffer_n_get ); + network->get_weights( 3, buffer_w_get ); + + REQUIRE_THAT( buffer_n_get, Catch::Matchers::UnorderedRangeEquals( buffer_n ) ); + REQUIRE_THAT( buffer_w_get, Catch::Matchers::UnorderedRangeEquals( buffer_w ) ); + + // Check that the push_back function works for agent 3 + buffer_n.push_back( 5 ); // new neighbour + buffer_w.push_back( 1.0 ); // new weight for this new connection + network->push_back_neighbour_and_weight( 3, 5, 1.0 ); // new connection added with weight + // Check that the change worked for the push_back function + network->get_neighbours( 3, buffer_n_get ); + network->get_weights( 3, buffer_w_get ); + REQUIRE_THAT( buffer_n_get, Catch::Matchers::UnorderedRangeEquals( buffer_n ) ); + REQUIRE_THAT( buffer_w_get, Catch::Matchers::UnorderedRangeEquals( buffer_w ) ); + + // Now we test the transpose() function + + // First record all the old edges as tuples (i,j,w) where this edge goes from j -> i with weight w + std::set> old_edges; + for( size_t i_agent = 0; i_agent < network->n_agents(); i_agent++ ) + { + network->get_neighbours( i_agent, buffer_n ); + network->get_weights( i_agent, buffer_w ); + + for( size_t i_neighbour = 0; i_neighbour < buffer_n.size(); i_neighbour++ ) + { + auto neighbour = buffer_n[i_neighbour]; + auto weight = buffer_w[i_neighbour]; + std::tuple edge{ i_agent, neighbour, weight }; + old_edges.insert( edge ); + } + } + + network->transpose(); + + // Now we go over the transposed network and try to re-identify all edges + for( size_t i_agent = 0; i_agent < network->n_agents(); i_agent++ ) + { + network->get_neighbours( i_agent, buffer_n ); + network->get_weights( i_agent, buffer_w ); + + for( size_t i_neighbour = 0; i_neighbour < buffer_n.size(); i_neighbour++ ) + { + auto neighbour = buffer_n[i_neighbour]; + auto weight = buffer_w[i_neighbour]; + std::tuple edge{ + neighbour, i_agent, weight + }; // 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 + } + } + + REQUIRE( old_edges.empty() ); +} \ No newline at end of file diff --git a/test/test_sampling.cpp b/test/test_sampling.cpp new file mode 100644 index 0000000..3bc64d6 --- /dev/null +++ b/test/test_sampling.cpp @@ -0,0 +1,128 @@ +#include "util/math.hpp" +#include +#include +#include +#include +#include +#include +#include +#include + +double compute_p( size_t k, size_t n ) +{ + if( k == 0 ) + { + return 0.0; + } + else + { + double p = 1.0 / ( double( n ) - 1.0 ); + return p + ( 1.0 - p ) * compute_p( k - 1, n - 1 ); + } +} + +TEST_CASE( "Testing sampling functions" ) +{ + std::random_device rd; + std::mt19937 gen( rd() ); + + SECTION( "draw_unique_k_from_n", "Drawing k numbers out of n" ) + { + + const size_t N_RUNS = 10000; + + const size_t k = 6; + const size_t n = 100; + const size_t ignore_idx = 11; + + std::vector histogram( n, 0 ); // Count how often each element occurs amongst all samples + + std::vector buffer{}; + for( size_t i = 0; i < N_RUNS; i++ ) + { + Seldon::draw_unique_k_from_n( ignore_idx, k, n, buffer, gen ); + for( const auto & n : buffer ) + { + histogram[n]++; + } + } + + // In each run there is a probability of p for each element to be selected + // That means for each histogram bin we have a binomial distribution with p + double p = compute_p( k, n ); + + size_t mean = N_RUNS * p; + // The variance of a binomial distribution is var = n*p*(1-p) + size_t sigma = std::sqrt( N_RUNS * p * ( 1.0 - p ) ); + + INFO( "Binomial distribution parameters" ); + INFO( fmt::format( " p = {}", p ) ); + INFO( fmt::format( " mean = {}", mean ) ); + INFO( fmt::format( " sigma = {}", sigma ) ); + + REQUIRE( histogram[ignore_idx] == 0 ); // The ignore_idx should never be selected + + size_t number_outside_three_sigma = 0; + for( const auto & n : histogram ) + { + if( n == 0 ) + { + continue; + } + INFO( fmt::format( " n = {}", n ) ); + INFO( fmt::format( " mean = {}", mean ) ); + INFO( fmt::format( " sigma = {}", sigma ) ); + + if( std::abs( double( n ) - double( mean ) ) > 3.0 * sigma ) + { + number_outside_three_sigma++; + } + + REQUIRE_THAT( n, Catch::Matchers::WithinAbs( mean, 5 * sigma ) ); + } + + if( number_outside_three_sigma > 0.01 * N_RUNS ) + WARN( fmt::format( + "Many deviations beyond the 3 sigma range. {} out of {}", number_outside_three_sigma, N_RUNS ) ); + } + + SECTION( "weighted_reservior_sampling", "Testing weighted reservoir sampling with A_ExpJ algorithm" ) + { + + const size_t N_RUNS = 10000; + + const size_t k = 6; + const size_t n = 100; + const size_t ignore_idx = 11; + const size_t ignore_idx2 = 29; + + std::vector histogram( n, 0 ); // Count how often each element occurs amongst all samples + + auto weight_callback = []( size_t idx ) { + if( ( idx == ignore_idx ) | ( idx == ignore_idx2 ) ) + { + return 0.0; + } + else + { + std::abs( double( n / 2.0 ) - double( idx ) ); + } + }; + + std::vector buffer{}; + + for( size_t i = 0; i < N_RUNS; i++ ) + { + Seldon::reservoir_sampling_A_ExpJ( k, n, weight_callback, buffer, gen ); + for( const auto & n : buffer ) + { + histogram[n]++; + } + } + + REQUIRE( histogram[ignore_idx] == 0 ); // The ignore_idx should never be selected + REQUIRE( histogram[ignore_idx2] == 0 ); // The ignore_idx should never be selected + + // TODO: histogram and sigma test + } +} \ No newline at end of file diff --git a/test/test_tarjan.cpp b/test/test_tarjan.cpp index c7fb551..5938932 100644 --- a/test/test_tarjan.cpp +++ b/test/test_tarjan.cpp @@ -1,4 +1,5 @@ #include +#include #include "connectivity.hpp" #include @@ -49,7 +50,7 @@ TEST_CASE( "Test for Tarjan's algorithm for strongly connected networks", "[tarj // List of SCC // [[5, 4], [3], [2, 1, 0], [9, 8, 7, 6]] - fmt::print( "SCC = {}\n", tarjan_scc.scc_list ); + INFO( fmt::format( "SCC = {}\n", tarjan_scc.scc_list ) ); std::set> expected_scc{ { 5, 4 }, { 3 }, { 2, 1, 0 }, { 9, 8, 7, 6 } }; @@ -61,5 +62,5 @@ TEST_CASE( "Test for Tarjan's algorithm for strongly connected networks", "[tarj } // There should be 4 strongly connected components - REQUIRE( tarjan_scc.scc_list.size() == 4 ); + REQUIRE_THAT( tarjan_scc.scc_list, Catch::Matchers::SizeIs( 4 ) ); } \ No newline at end of file