From a0cd5ce1dcc4ba520a7e321b23e5a243ce2170ec Mon Sep 17 00:00:00 2001 From: Moritz Sallermann Date: Sat, 23 Mar 2024 15:16:20 +0000 Subject: [PATCH 1/5] 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 2/5] 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 3/5] 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 4/5] 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 5/5] 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 }