WIP: iteration of activity driven model

include/models/ActivityDrivenModel.hpp

@@ -44,15 +44,15 @@ class ActivityAgentModel : public Model<Agent<ActivityAgentData>>
     // 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
+    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 K           = 3.0; // Social interaction strength; K>0
 
     double convergence_tol = 1e-12; // TODO: ??
 

src/models/ActivityDrivenModel.cpp

@@ -1,4 +1,6 @@
 #include "models/ActivityDrivenModel.hpp"
+#include "util/math.hpp"
+#include <cstddef>
 #include <random>
 #include <stdexcept>
 
@@ -32,5 +34,59 @@ void Seldon::ActivityAgentModel::iteration()
 {
     Model<AgentT>::iteration();
 
-    // throw std::runtime_error( "ActivityAgentModel::iteration():: Not implemented!" );
+    std::uniform_real_distribution<> dis_activation( 0.0, 1.0 ); // Opinion initial values
+    std::vector<size_t> contacted_agents{};
+
+    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 )
+        {
+            // Contact m_agents according to the homophily distribution
+            // TODO: use homophily stuff instead of uniform
+            Seldon::draw_unique_k_from_n(
+                idx_agent, m, network.n_agents(), contacted_agents,
+                gen ); // now contacted_agents contains the indices of the contacted agents
+            // Set the *outgoing* edges
+            network.set_neighbours_and_weights( idx_agent, contacted_agents, 1.0 );
+        }
+        else
+        {
+            network.set_neighbours_and_weights( idx_agent, {}, {} );
+        }
+
+        // TODO: implement reciprocity
+        // ...
+    }
+
+    network.transpose(); // transpose the network, so that we have incoming edges
+
+    // Integrate the ODE
+    auto neighbour_buffer = std::vector<size_t>();
+    size_t j_index        = 0;
+
+    agents_current_copy = agents; // Set the copy to the current state of agents. TODO: this is somewhat wasteful since
+                                  // activities are not changing
+
+    for( size_t i = 0; i < network.n_agents(); i++ )
+    {
+        network.get_neighbours( i, neighbour_buffer ); // Get the incoming neighbours
+        for( size_t j = 0; j < neighbour_buffer.size(); j++ )
+        {
+            // TODO: currently this uses an euler integration -> use RK4
+            j_index = neighbour_buffer[j];
+            agents_current_copy[i].data.opinion
+                += dt * ( -agents[i].data.opinion + K * ( std::tanh( alpha * agents[j_index].data.opinion ) ) );
+        }
+    }
+
+    // Update the agents from the copy
+    for( std::size_t i = 0; i < agents.size(); i++ )
+    {
+        max_opinion_diff
+            = std::max( max_opinion_diff, std::abs( agents[i].data.opinion - agents_current_copy[i].data.opinion ) );
+        agents[i] = agents_current_copy[i];
+    }
 }