1 /* 2 * A sample test case which can be used as a template. 3 */ 4 #include <iostream> 5 #include <cmath> 6 #include <cppunit/TestCase.h> 7 #include <cppunit/extensions/HelperMacros.h> 8 #include <stdexcept> 9 10 #include "../../src/model.h" 11 #include "../../src/forest.h" 12 #include "../../src/summary_statistics/tmrca.h" 13 14 class TestModel : public CppUnit::TestCase { 15 16 CPPUNIT_TEST_SUITE( TestModel ); 17 18 CPPUNIT_TEST( testSetGetMutationRate ); 19 CPPUNIT_TEST( testAddChangePositions ); 20 CPPUNIT_TEST( testSetGetRecombinationRate ); 21 CPPUNIT_TEST( testGetPopulationSize ); 22 CPPUNIT_TEST( testAddChangeTime ); 23 CPPUNIT_TEST( testAddSampleSizes ); 24 CPPUNIT_TEST( testAddPopulationSizes ); 25 CPPUNIT_TEST( testAddRelativePopulationSizes ); 26 CPPUNIT_TEST( testAddGrowthRates ); 27 CPPUNIT_TEST( testAddMigRates ); 28 CPPUNIT_TEST( testAddMigRate ); 29 CPPUNIT_TEST( testDebugConstructor ); 30 CPPUNIT_TEST( testIncreaseTime ); 31 CPPUNIT_TEST( testGetNextTime ); 32 CPPUNIT_TEST( testGetters ); 33 CPPUNIT_TEST( testHasFixedTimeEvent ); 34 CPPUNIT_TEST( testCheck ); 35 CPPUNIT_TEST( testPopSizeAfterGrowth ); 36 CPPUNIT_TEST( testAddSummaryStatistic ); 37 CPPUNIT_TEST( testSetLocusLength ); 38 CPPUNIT_TEST( testAddPopToVectorList ); 39 CPPUNIT_TEST( testAddPopToMatrixList ); 40 41 CPPUNIT_TEST_SUITE_END(); 42 43 public: testAddChangeTime()44 void testAddChangeTime() { 45 Model model = Model(); 46 std::vector<double> v1 = std::vector<double>(1, 1), 47 v2 = std::vector<double>(1, 2), 48 v3 = std::vector<double>(1, 3); 49 50 // Check basic adding first time; 51 CPPUNIT_ASSERT( model.addChangeTime(0) == 0 ); 52 CPPUNIT_ASSERT( model.change_times_.size() == 1 ); 53 CPPUNIT_ASSERT( model.change_times_.at(0) == 0 ); 54 CPPUNIT_ASSERT( model.pop_sizes_list_.size() == 1 ); 55 model.pop_sizes_list_[0] = v1; 56 57 // Check adding a time at the end; 58 CPPUNIT_ASSERT_EQUAL( (size_t)1, model.addChangeTime(3) ); 59 CPPUNIT_ASSERT( model.change_times_.size() == 2 ); 60 CPPUNIT_ASSERT( model.change_times_.at(1) == 3 ); 61 CPPUNIT_ASSERT( model.pop_sizes_list_[0] == v1 ); 62 CPPUNIT_ASSERT( model.pop_sizes_list_[1].empty() ); 63 model.pop_sizes_list_[1] = v3; 64 65 // Check adding a time in the middle; 66 CPPUNIT_ASSERT_EQUAL( (size_t)1, model.addChangeTime(2) ); 67 CPPUNIT_ASSERT( model.change_times_.size() == 3 ); 68 CPPUNIT_ASSERT( model.change_times_.at(0) == 0 ); 69 CPPUNIT_ASSERT( model.change_times_.at(1) == 2 ); 70 CPPUNIT_ASSERT( model.change_times_.at(2) == 3 ); 71 CPPUNIT_ASSERT( model.pop_sizes_list_[0] == v1 ); 72 CPPUNIT_ASSERT( model.pop_sizes_list_[1].empty() ); 73 CPPUNIT_ASSERT( model.pop_sizes_list_[2] == v3 ); 74 model.pop_sizes_list_[1] = v2; 75 76 CPPUNIT_ASSERT( model.pop_sizes_list_[0] == v1 ); 77 CPPUNIT_ASSERT( model.pop_sizes_list_[1] == v2 ); 78 CPPUNIT_ASSERT( model.pop_sizes_list_[2] == v3 ); 79 80 // Check that we don't add a time twice 81 CPPUNIT_ASSERT_EQUAL( (size_t)0, model.addChangeTime(0) ); 82 CPPUNIT_ASSERT_EQUAL( (size_t)1, model.addChangeTime(1) ); 83 CPPUNIT_ASSERT_EQUAL( (size_t)2, model.addChangeTime(2) ); 84 CPPUNIT_ASSERT_EQUAL( (size_t)3, model.addChangeTime(3) ); 85 CPPUNIT_ASSERT( model.change_times_.size() == 4 ); 86 CPPUNIT_ASSERT( model.pop_sizes_list_.size() == 4 ); 87 } 88 testAddSampleSizes()89 void testAddSampleSizes() { 90 Model model = Model(); 91 model.addSampleSizes(0.0, std::vector<size_t>(1,3)); 92 CPPUNIT_ASSERT_EQUAL( model.sample_size(), (size_t)3 ); 93 CPPUNIT_ASSERT_EQUAL( model.sample_population(2), (size_t)0 ); 94 CPPUNIT_ASSERT_EQUAL( model.sample_time(1), (double)0.0 ); 95 96 model = Model(0); 97 model.set_population_number(3); 98 model.addSampleSizes(0.0, std::vector<size_t>(3,1)); 99 CPPUNIT_ASSERT_EQUAL( model.sample_size(), (size_t)3 ); 100 CPPUNIT_ASSERT_EQUAL( model.sample_population(0), (size_t)0 ); 101 CPPUNIT_ASSERT_EQUAL( model.sample_population(1), (size_t)1 ); 102 CPPUNIT_ASSERT_EQUAL( model.sample_population(2), (size_t)2 ); 103 CPPUNIT_ASSERT_EQUAL( model.sample_time(0), (double)0.0 ); 104 CPPUNIT_ASSERT_EQUAL( model.sample_time(1), (double)0.0 ); 105 CPPUNIT_ASSERT_EQUAL( model.sample_time(2), (double)0.0 ); 106 107 std::vector<size_t> sample_size; 108 sample_size.push_back(0); 109 sample_size.push_back(0); 110 sample_size.push_back(2); 111 model.addSampleSizes(1.0, sample_size); 112 CPPUNIT_ASSERT_EQUAL( (size_t)5, model.sample_size() ); 113 CPPUNIT_ASSERT_EQUAL( (size_t)2, model.sample_population(3) ); 114 CPPUNIT_ASSERT_EQUAL( (size_t)2, model.sample_population(4) ); 115 CPPUNIT_ASSERT_EQUAL( 1.0, model.sample_time(3) ); 116 CPPUNIT_ASSERT_EQUAL( 1.0, model.sample_time(4) ); 117 118 Model model2 = Model(); 119 std::vector<size_t> sample_sizes; 120 sample_sizes.push_back(5); 121 sample_sizes.push_back(2); 122 model2.addSampleSizes(0.0, sample_sizes); 123 CPPUNIT_ASSERT_EQUAL( model2.sample_size(), (size_t)7 ); 124 CPPUNIT_ASSERT_EQUAL( model2.sample_population(0), (size_t)0 ); 125 CPPUNIT_ASSERT_EQUAL( model2.sample_population(4), (size_t)0 ); 126 CPPUNIT_ASSERT_EQUAL( model2.sample_population(5), (size_t)1 ); 127 CPPUNIT_ASSERT_EQUAL( model2.sample_population(6), (size_t)1 ); 128 CPPUNIT_ASSERT_EQUAL( model2.sample_time(0), (double)0.0 ); 129 CPPUNIT_ASSERT_EQUAL( model2.sample_time(4), (double)0.0 ); 130 CPPUNIT_ASSERT_EQUAL( model2.sample_time(6), (double)0.0 ); 131 sample_sizes.clear(); 132 sample_sizes.push_back(2); 133 sample_sizes.push_back(1); 134 model2.addSampleSizes(7.4, sample_sizes); 135 CPPUNIT_ASSERT_EQUAL( model2.sample_size(), (size_t)10 ); 136 CPPUNIT_ASSERT_EQUAL( model2.sample_population(7), (size_t)0 ); 137 CPPUNIT_ASSERT_EQUAL( model2.sample_population(8), (size_t)0 ); 138 CPPUNIT_ASSERT_EQUAL( model2.sample_population(9), (size_t)1 ); 139 CPPUNIT_ASSERT_EQUAL( model2.sample_time(0), (double)0.0 ); 140 CPPUNIT_ASSERT_EQUAL( model2.sample_time(6), (double)0.0 ); 141 CPPUNIT_ASSERT_EQUAL( model2.sample_time(8), (double)7.4 ); 142 } 143 testAddPopulationSizes()144 void testAddPopulationSizes() { 145 Model model = Model(); 146 model.set_population_number(2); 147 model.addPopulationSizes(1, std::vector<double>(2, 4)); 148 model.resetTime(); 149 model.increaseTime(); 150 CPPUNIT_ASSERT_EQUAL( 1.0, model.getCurrentTime() ); 151 CPPUNIT_ASSERT( model.population_size(0) == 4 ); 152 153 model.addPopulationSizes(1, std::vector<double>(2, 5), true, false); 154 model.increaseTime(); 155 CPPUNIT_ASSERT_EQUAL( 1.0 * 4 * model.default_pop_size(), model.getCurrentTime() ); 156 CPPUNIT_ASSERT( areSame(model.population_size(0), 5) ); 157 158 model.addPopulationSizes(2, 10, true, false); 159 model.increaseTime(); 160 CPPUNIT_ASSERT_EQUAL( 2.0 * 4 * model.default_pop_size(), model.getCurrentTime() ); 161 CPPUNIT_ASSERT( areSame(model.population_size(0), 10) ); 162 163 auto pop_sizes2 = std::vector<double>(); 164 pop_sizes2.push_back(7); 165 pop_sizes2.push_back(6); 166 model.addPopulationSizes(0.0, pop_sizes2); 167 model.resetTime(); 168 CPPUNIT_ASSERT( areSame(model.population_size(0), 7) ); 169 CPPUNIT_ASSERT( areSame(model.population_size(1), 6) ); 170 171 CPPUNIT_ASSERT_THROW( model.addPopulationSizes(1, std::vector<double>(1, 5)), std::logic_error ); 172 CPPUNIT_ASSERT_THROW( model.addPopulationSizes(1, std::vector<double>(3, 5)), std::logic_error ); 173 } 174 testAddRelativePopulationSizes()175 void testAddRelativePopulationSizes() { 176 Model model = Model(); 177 model.set_population_number(2); 178 model.addPopulationSizes(1, std::vector<double>(2, .5), false, true); 179 model.resetTime(); 180 model.increaseTime(); 181 182 CPPUNIT_ASSERT_EQUAL( 1.0, model.getCurrentTime() ); 183 CPPUNIT_ASSERT_EQUAL( 0.5 * model.default_pop_size(), model.population_size(0) ); 184 185 model.addPopulationSizes(1, std::vector<double>(2, .4), true, true); 186 model.increaseTime(); 187 CPPUNIT_ASSERT_EQUAL( 1.0 * 4 * model.default_pop_size(), model.getCurrentTime() ); 188 CPPUNIT_ASSERT_EQUAL( 0.4 * model.default_pop_size(), model.population_size(0) ); 189 190 model.addPopulationSizes(2, 10, true, true); 191 model.increaseTime(); 192 CPPUNIT_ASSERT_EQUAL( 2.0 * 4 * model.default_pop_size(), model.getCurrentTime() ); 193 CPPUNIT_ASSERT( areSame(10.0 * model.default_pop_size(), model.population_size(0)) ); 194 195 CPPUNIT_ASSERT_THROW( model.addPopulationSizes(3, 0, true, true), std::invalid_argument); 196 } 197 testAddGrowthRates()198 void testAddGrowthRates() { 199 Model model = Model(); 200 model.set_population_number(2); 201 202 model.addGrowthRates(1, std::vector<double>(2, 1.5)); 203 CPPUNIT_ASSERT( model.growth_rates_list_.at(1).at(0) == 1.5 ); 204 205 std::vector<double> growth_rates2 = std::vector<double>(); 206 growth_rates2.push_back(2.5); 207 growth_rates2.push_back(3.5); 208 model.addGrowthRates(0, growth_rates2); 209 CPPUNIT_ASSERT( model.growth_rates_list_.at(0).at(0) == 2.5 ); 210 CPPUNIT_ASSERT( model.growth_rates_list_.at(0).at(1) == 3.5 ); 211 212 CPPUNIT_ASSERT_THROW( model.addGrowthRates(1, std::vector<double>(1, 5)), std::logic_error ); 213 CPPUNIT_ASSERT_THROW( model.addGrowthRates(1, std::vector<double>(3, 5)), std::logic_error ); 214 } 215 testGetPopulationSize()216 void testGetPopulationSize() { 217 Model model = Model(2); 218 model.set_population_number(2); 219 model.addSymmetricMigration(0.0, 1); 220 model.addGrowthRates(1.0, std::vector<double>(2, 1.5)); 221 model.finalize(); 222 model.resetTime(); 223 224 CPPUNIT_ASSERT( areSame( model.default_pop_size(), model.population_size() ) ); 225 CPPUNIT_ASSERT( areSame( model.default_pop_size(), model.population_size(0) ) ); 226 CPPUNIT_ASSERT( areSame( model.default_pop_size(), model.population_size(1) ) ); 227 model.increaseTime(); 228 CPPUNIT_ASSERT_EQUAL( 1.0, model.getCurrentTime() ); 229 CPPUNIT_ASSERT( areSame( model.default_pop_size(), model.population_size(0) ) ); 230 CPPUNIT_ASSERT( areSame( model.default_pop_size()* std::exp( -1.5 ), model.population_size(0, 2.0) ) ); 231 } 232 testAddMigRates()233 void testAddMigRates() { 234 Model model = Model(2); 235 model.set_population_number(3); 236 237 std::vector<double> rates; 238 for (size_t i = 1; i < 9; ++i) { 239 rates.push_back(i); 240 } 241 CPPUNIT_ASSERT_THROW( model.addMigrationRates(1, rates), std::logic_error ); 242 rates.push_back(9); 243 model.addMigrationRates(1, rates); 244 245 model.resetTime(); 246 model.increaseTime(); 247 CPPUNIT_ASSERT_EQUAL( 2.0, model.migration_rate(0,1) ); 248 CPPUNIT_ASSERT_EQUAL( 3.0, model.migration_rate(0,2) ); 249 CPPUNIT_ASSERT_EQUAL( 4.0, model.migration_rate(1,0) ); 250 CPPUNIT_ASSERT_EQUAL( 7.0, model.migration_rate(2,0) ); 251 CPPUNIT_ASSERT_EQUAL( 6.0, model.migration_rate(1,2) ); 252 } 253 testAddMigRate()254 void testAddMigRate() { 255 Model model = Model(2); 256 model.set_population_number(2); 257 model.addMigrationRate(0.0, 0, 1, 0.5); 258 model.resetTime(); 259 CPPUNIT_ASSERT_EQUAL( 0.5, model.migration_rate(0,1) ); 260 CPPUNIT_ASSERT( std::isnan(model.migration_rate(1,0)) ); 261 262 model.addMigrationRate(0.0, 0, 1, 0.7); 263 CPPUNIT_ASSERT_EQUAL( 0.7, model.migration_rate(0,1) ); 264 CPPUNIT_ASSERT( std::isnan(model.migration_rate(1,0)) ); 265 266 model.addMigrationRate(0.0, 1, 0, 0.9); 267 CPPUNIT_ASSERT_EQUAL( 0.7, model.migration_rate(0,1) ); 268 CPPUNIT_ASSERT_EQUAL( 0.9, model.migration_rate(1,0) ); 269 } 270 testDebugConstructor()271 void testDebugConstructor() { 272 Model model = Model(7); 273 CPPUNIT_ASSERT_EQUAL( (size_t)7, model.sample_size() ); 274 CPPUNIT_ASSERT( model.growth_rate(0) == 0 ); 275 CPPUNIT_ASSERT_EQUAL( (size_t)0, model.current_time_idx_ ); 276 CPPUNIT_ASSERT_EQUAL( (size_t)0, model.current_seq_idx_ ); 277 } 278 testAddChangePositions()279 void testAddChangePositions() { 280 Model model = Model(7); 281 CPPUNIT_ASSERT_EQUAL((size_t)1, model.change_position_.size()); 282 CPPUNIT_ASSERT_EQUAL(0.0, model.change_position_[0]); 283 CPPUNIT_ASSERT_EQUAL((size_t)1, model.mutation_rates_.size()); 284 CPPUNIT_ASSERT_EQUAL((size_t)1, model.recombination_rates_.size()); 285 286 // 0.0 is already in => no change 287 CPPUNIT_ASSERT_EQUAL((size_t)0, model.addChangePosition(0.0)); 288 CPPUNIT_ASSERT_EQUAL((size_t)1, model.change_position_.size()); 289 CPPUNIT_ASSERT_EQUAL(0.0, model.change_position_[0]); 290 CPPUNIT_ASSERT_EQUAL((size_t)1, model.mutation_rates_.size()); 291 CPPUNIT_ASSERT_EQUAL((size_t)1, model.recombination_rates_.size()); 292 293 // Really add 1.0 ad the end 294 CPPUNIT_ASSERT_EQUAL((size_t)1, model.addChangePosition(1.0)); 295 CPPUNIT_ASSERT_EQUAL((size_t)2, model.change_position_.size()); 296 CPPUNIT_ASSERT_EQUAL(0.0, model.change_position_[0]); 297 CPPUNIT_ASSERT_EQUAL(1.0, model.change_position_[1]); 298 CPPUNIT_ASSERT_EQUAL((size_t)2, model.mutation_rates_.size()); 299 CPPUNIT_ASSERT_EQUAL((size_t)2, model.recombination_rates_.size()); 300 301 // Really add 1.0 again => nothing changes 302 CPPUNIT_ASSERT_EQUAL((size_t)1, model.addChangePosition(1.0)); 303 CPPUNIT_ASSERT_EQUAL((size_t)2, model.change_position_.size()); 304 CPPUNIT_ASSERT_EQUAL(0.0, model.change_position_[0]); 305 CPPUNIT_ASSERT_EQUAL(1.0, model.change_position_[1]); 306 CPPUNIT_ASSERT_EQUAL((size_t)2, model.mutation_rates_.size()); 307 CPPUNIT_ASSERT_EQUAL((size_t)2, model.recombination_rates_.size()); 308 309 // Add 0.5 in the middle 310 CPPUNIT_ASSERT_EQUAL((size_t)1, model.addChangePosition(0.5)); 311 CPPUNIT_ASSERT_EQUAL((size_t)3, model.change_position_.size()); 312 CPPUNIT_ASSERT_EQUAL(0.0, model.change_position_[0]); 313 CPPUNIT_ASSERT_EQUAL(0.5, model.change_position_[1]); 314 CPPUNIT_ASSERT_EQUAL(1.0, model.change_position_[2]); 315 CPPUNIT_ASSERT_EQUAL((size_t)3, model.mutation_rates_.size()); 316 CPPUNIT_ASSERT_EQUAL((size_t)3, model.recombination_rates_.size()); 317 318 model.setLocusLength(10.0); 319 CPPUNIT_ASSERT_THROW(model.addChangePosition(-0.5), std::invalid_argument); 320 CPPUNIT_ASSERT_THROW(model.addChangePosition(10.5), std::invalid_argument); 321 } 322 testIncreaseTime()323 void testIncreaseTime() { 324 Model model = Model(7); 325 model.addGrowthRates(1.0, std::vector<double>(1, 1.5)); 326 model.addGrowthRates(2.0, std::vector<double>(1, 1)); 327 328 CPPUNIT_ASSERT_EQUAL( (size_t)0, model.current_time_idx_ ); 329 CPPUNIT_ASSERT_NO_THROW( model.increaseTime() ); 330 331 CPPUNIT_ASSERT_EQUAL( (size_t)1, model.current_time_idx_ ); 332 CPPUNIT_ASSERT_NO_THROW( model.increaseTime() ); 333 334 CPPUNIT_ASSERT_EQUAL( (size_t)2, model.current_time_idx_ ); 335 336 CPPUNIT_ASSERT_THROW( model.increaseTime(), std::out_of_range ); 337 } 338 testGetNextTime()339 void testGetNextTime() { 340 Model model = Model(7); 341 CPPUNIT_ASSERT( model.getNextTime() == DBL_MAX ); 342 343 model.addGrowthRates(1.0, std::vector<double>(1, 1.5)); 344 model.addGrowthRates(2.0, std::vector<double>(1, 1)); 345 346 CPPUNIT_ASSERT_EQUAL( 1.0, model.getNextTime() ); 347 model.increaseTime(); 348 CPPUNIT_ASSERT_EQUAL( 2.0, model.getNextTime() ); 349 model.increaseTime(); 350 CPPUNIT_ASSERT( model.getNextTime() == DBL_MAX ); 351 } 352 testGetters()353 void testGetters() { 354 Model model = Model(7); 355 model.addGrowthRates(1.0, std::vector<double>(1, 1.5)); 356 model.addPopulationSizes(2.0, std::vector<double>(1, 5000)); 357 model.addGrowthRates(3.0, std::vector<double>(1, 1)); 358 model.addPopulationSizes(4.0, std::vector<double>(1, 10000)); 359 360 CPPUNIT_ASSERT_EQUAL( (size_t)7, model.sample_size() ); 361 CPPUNIT_ASSERT_EQUAL( 0.0, model.growth_rate(0) ); 362 363 CPPUNIT_ASSERT_NO_THROW( model.increaseTime() ); // 1.0 364 CPPUNIT_ASSERT_EQUAL( 1.5, model.growth_rate(0) ); 365 366 CPPUNIT_ASSERT_NO_THROW( model.increaseTime() ); // 2.0 367 CPPUNIT_ASSERT_EQUAL( 1.5, model.growth_rate(0) ); 368 CPPUNIT_ASSERT_EQUAL( 5000.0, model.population_size(0) ); 369 370 CPPUNIT_ASSERT_NO_THROW( model.increaseTime() ); // 3.0 371 CPPUNIT_ASSERT_EQUAL( 1.0, model.growth_rate(0) ); 372 CPPUNIT_ASSERT_EQUAL( 5000.0, model.population_size(0) ); 373 374 CPPUNIT_ASSERT_NO_THROW( model.increaseTime() ); // 4.0 375 CPPUNIT_ASSERT_EQUAL( 1.0, model.growth_rate(0) ); 376 CPPUNIT_ASSERT_EQUAL( 10000.0, model.population_size(0) ); 377 } 378 testHasFixedTimeEvent()379 void testHasFixedTimeEvent() { 380 Model model = Model(); 381 model.set_population_number(3); 382 383 model.addSingleMigrationEvent(10, 1, 0, .5); 384 model.resetTime(); 385 CPPUNIT_ASSERT( ! model.hasFixedTimeEvent(1) ); 386 CPPUNIT_ASSERT( ! model.hasFixedTimeEvent(10) ); 387 model.increaseTime(); 388 CPPUNIT_ASSERT( model.hasFixedTimeEvent(10) ); 389 CPPUNIT_ASSERT( ! model.hasFixedTimeEvent(1) ); 390 CPPUNIT_ASSERT( ! model.hasFixedTimeEvent(20) ); 391 } 392 testSetGetMutationRate()393 void testSetGetMutationRate() { 394 Model model = Model(5); 395 model.setLocusLength(10); 396 model.setMutationRate(0.001); 397 CPPUNIT_ASSERT_EQUAL( 0.001, model.mutation_rate() ); 398 399 model.setMutationRate(10, false, false); 400 CPPUNIT_ASSERT_EQUAL( 10.0, model.mutation_rate() ); 401 402 model.setMutationRate(10, true, false); 403 CPPUNIT_ASSERT_EQUAL( 1.0, model.mutation_rate() ); 404 405 model.setMutationRate(10, false, true); 406 CPPUNIT_ASSERT_EQUAL( 10.0/(4*model.default_pop_size()), model.mutation_rate() ); 407 408 model.setMutationRate(10, true, true); 409 CPPUNIT_ASSERT_EQUAL(10.0/(4*model.default_pop_size()*10), 410 model.mutation_rate() ); 411 412 // Test setting multiple rates 413 model.setMutationRate(10, false, false, 0.0); 414 model.setMutationRate(5, false, false, 2.0); 415 CPPUNIT_ASSERT_EQUAL(10.0, model.mutation_rate()); 416 model.increaseSequencePosition(); 417 CPPUNIT_ASSERT_EQUAL(5.0, model.mutation_rate()); 418 419 model.setMutationRate(7.5, false, false, 1.0); 420 model.resetSequencePosition(); 421 CPPUNIT_ASSERT_EQUAL(10.0, model.mutation_rate()); 422 model.increaseSequencePosition(); 423 CPPUNIT_ASSERT_EQUAL(7.5, model.mutation_rate()); 424 model.increaseSequencePosition(); 425 CPPUNIT_ASSERT_EQUAL(5.0, model.mutation_rate()); 426 427 // Test if the vector are filled during finalization 428 model.resetSequencePosition(); 429 model.setRecombinationRate(2.5, false, false, 0.5); 430 model.finalize(); 431 CPPUNIT_ASSERT_EQUAL(10.0, model.mutation_rate()); 432 model.increaseSequencePosition(); 433 CPPUNIT_ASSERT_EQUAL(10.0, model.mutation_rate()); 434 model.increaseSequencePosition(); 435 CPPUNIT_ASSERT_EQUAL(7.5, model.mutation_rate()); 436 model.increaseSequencePosition(); 437 CPPUNIT_ASSERT_EQUAL(5.0, model.mutation_rate()); 438 } 439 testSetGetRecombinationRate()440 void testSetGetRecombinationRate() { 441 Model model = Model(); 442 CPPUNIT_ASSERT( !model.has_recombination() ); 443 444 model.setLocusLength(101); 445 model.setRecombinationRate(0.001); 446 CPPUNIT_ASSERT_EQUAL( 0.001, model.recombination_rate() ); 447 CPPUNIT_ASSERT_EQUAL( (size_t)101, model.loci_length() ); 448 CPPUNIT_ASSERT( model.has_recombination() ); 449 450 model.setRecombinationRate(0.001, false, false); 451 CPPUNIT_ASSERT_EQUAL( 0.001, model.recombination_rate() ); 452 CPPUNIT_ASSERT_EQUAL( (size_t)101, model.loci_length() ); 453 CPPUNIT_ASSERT( model.has_recombination() ); 454 455 model.setRecombinationRate(0.001, true, false); 456 CPPUNIT_ASSERT_EQUAL( 0.00001, model.recombination_rate() ); 457 CPPUNIT_ASSERT_EQUAL( (size_t)101, model.loci_length() ); 458 CPPUNIT_ASSERT( model.has_recombination() ); 459 460 model.setRecombinationRate(0.001, false, true); 461 CPPUNIT_ASSERT_EQUAL( 0.001/(4*model.default_pop_size()), model.recombination_rate() ); 462 CPPUNIT_ASSERT_EQUAL( (size_t)101, model.loci_length() ); 463 464 model.setRecombinationRate(0.001, true, true); 465 CPPUNIT_ASSERT( areSame(0.00001/(4*model.default_pop_size()), model.recombination_rate()) ); 466 CPPUNIT_ASSERT_EQUAL( (size_t)101, model.loci_length() ); 467 468 // Test setting multiple rates 469 model.setRecombinationRate(10, false, false, 0.0); 470 model.setRecombinationRate(5, false, false, 2.0); 471 CPPUNIT_ASSERT_EQUAL(10.0, model.recombination_rate()); 472 model.increaseSequencePosition(); 473 CPPUNIT_ASSERT_EQUAL(5.0, model.recombination_rate()); 474 475 model.setRecombinationRate(7.5, false, false, 1.0); 476 model.resetSequencePosition(); 477 CPPUNIT_ASSERT_EQUAL(10.0, model.recombination_rate()); 478 model.increaseSequencePosition(); 479 CPPUNIT_ASSERT_EQUAL(7.5, model.recombination_rate()); 480 model.increaseSequencePosition(); 481 CPPUNIT_ASSERT_EQUAL(5.0, model.recombination_rate()); 482 483 // Check for error on invalid change positions 484 CPPUNIT_ASSERT_THROW( model.setRecombinationRate(-0.001, 100), std::invalid_argument); 485 model.setLocusLength(10); 486 CPPUNIT_ASSERT_THROW(model.setRecombinationRate(7.5, false, false, -1.0), std::invalid_argument); 487 CPPUNIT_ASSERT_THROW(model.setRecombinationRate(7.5, false, false, 11.0), std::invalid_argument); 488 } 489 testCheck()490 void testCheck() { 491 Model model = Model(1); 492 CPPUNIT_ASSERT_THROW( model.check(), std::invalid_argument ); 493 494 model = Model(2); 495 model.set_population_number(2); 496 CPPUNIT_ASSERT_THROW( model.check(), std::invalid_argument ); 497 model.addMigrationRate(20, 0, 1, 0.0); 498 CPPUNIT_ASSERT_THROW( model.finalize(), std::invalid_argument ); 499 model.addMigrationRate(10, 0, 1, 5.0); 500 CPPUNIT_ASSERT_NO_THROW( model.finalize() ); 501 502 model = Model(3); 503 model.set_population_number(2); 504 model.addSingleMigrationEvent(1, 0, 1, 1); 505 CPPUNIT_ASSERT_NO_THROW( model.check() ); 506 } 507 testPopSizeAfterGrowth()508 void testPopSizeAfterGrowth() { 509 // Growth only 510 Model model = Model(5); 511 model.set_population_number(2); 512 model.addSymmetricMigration(0, 1.0); 513 model.addPopulationSizes(0, 1000); 514 std::vector<double> growth_rates; 515 growth_rates.push_back(-0.5); 516 growth_rates.push_back(0.5); 517 model.addGrowthRates(1.0, growth_rates); 518 model.addGrowthRates(2.5, 2); 519 model.addSingleMigrationEvent(3.5, 1, 0, 0.5); 520 model.finalize(); 521 522 model.increaseTime(); 523 model.increaseTime(); 524 double n_1 = 1000*std::exp(0.5*1.5), 525 n_2 = 1000*std::exp(-0.5*1.5); 526 CPPUNIT_ASSERT( areSame( n_1, model.population_size(0) ) ); 527 CPPUNIT_ASSERT( areSame( n_2, model.population_size(1) ) ); 528 model.increaseTime(); 529 CPPUNIT_ASSERT( areSame( n_1*std::exp(-2), model.population_size(0) ) ); 530 CPPUNIT_ASSERT( areSame( n_2*std::exp(-2), model.population_size(1) ) ); 531 CPPUNIT_ASSERT( areSame( n_1*std::exp(-2*2), model.population_size(0, 4.5) ) ); 532 CPPUNIT_ASSERT( areSame( n_2*std::exp(-2*2), model.population_size(1, 4.5) ) ); 533 534 // Growth with a pop size change in between 535 model = Model(5); 536 model.set_population_number(2); 537 model.addSymmetricMigration(0, 1.0); 538 growth_rates.clear(); 539 growth_rates.push_back(-0.5); 540 growth_rates.push_back(0.5); 541 model.addGrowthRates(1.0, growth_rates); 542 model.addSymmetricMigration(1.25, 2.0); 543 model.addPopulationSize(1.5, 0, 500); 544 model.addGrowthRate(2.5, 1, 2.0); 545 model.addPopulationSize(3.0, 0, 500); 546 model.addGrowthRate(3.5, 0, 2.0); 547 model.finalize(); 548 //std::cout << model << std::endl; 549 550 551 model.resetTime(); 552 double size0 = model.default_pop_size(); 553 double size1 = model.default_pop_size(); 554 CPPUNIT_ASSERT( areSame( size0, model.population_size(0) ) ); 555 CPPUNIT_ASSERT( areSame( size1, model.population_size(1) ) ); 556 557 model.increaseTime(); 558 CPPUNIT_ASSERT_EQUAL( 1.0, model.getCurrentTime() ); 559 CPPUNIT_ASSERT( areSame( size0, model.population_size(0) ) ); 560 CPPUNIT_ASSERT( areSame( size1, model.population_size(1) ) ); 561 562 model.increaseTime(); 563 CPPUNIT_ASSERT_EQUAL( 1.25, model.getCurrentTime() ); 564 CPPUNIT_ASSERT( areSame( size0*std::exp(0.5*0.25), model.population_size(0) ) ); 565 CPPUNIT_ASSERT( areSame( size0*std::exp(-0.5*0.25), model.population_size(1) ) ); 566 567 model.increaseTime(); 568 CPPUNIT_ASSERT_EQUAL( 1.5, model.getCurrentTime() ); 569 CPPUNIT_ASSERT( areSame( 500, model.population_size(0) ) ); 570 CPPUNIT_ASSERT( areSame( model.default_pop_size()*std::exp(-0.5*0.5), model.population_size(1) ) ); 571 572 model.increaseTime(); 573 CPPUNIT_ASSERT_EQUAL( 2.5, model.getCurrentTime() ); 574 size0 = 500*std::exp(0.5*1.0); 575 size1 = model.default_pop_size()*std::exp(-0.5*1.5); 576 CPPUNIT_ASSERT( areSame( size0, model.population_size(0) ) ); 577 CPPUNIT_ASSERT( areSame( size1, model.population_size(1) ) ); 578 579 model.increaseTime(); 580 CPPUNIT_ASSERT_EQUAL( 3.0, model.getCurrentTime() ); 581 size0 = 500; 582 size1 *= exp(-2*0.5); 583 CPPUNIT_ASSERT( areSame( size0, model.population_size(0) ) ); 584 CPPUNIT_ASSERT( areSame( size1, model.population_size(1) ) ); 585 586 model.increaseTime(); 587 CPPUNIT_ASSERT_EQUAL( 3.5, model.getCurrentTime() ); 588 size0 *= exp(0.5*0.5); 589 size1 *= exp(-2*0.5); 590 CPPUNIT_ASSERT( areSame( size0, model.population_size(0) ) ); 591 CPPUNIT_ASSERT( areSame( size1, model.population_size(1) ) ); 592 593 model = Model(5); 594 model.set_population_number(2); 595 model.addSymmetricMigration(0, 1.0); 596 model.addPopulationSize(0, 1, 500); 597 model.finalize(); 598 //std::cout << model << std::endl; 599 model.resetTime(); 600 CPPUNIT_ASSERT( areSame( model.default_pop_size(), model.population_size(0) ) ); 601 CPPUNIT_ASSERT( areSame( 500.0, model.population_size(1) ) ); 602 } 603 testAddSummaryStatistic()604 void testAddSummaryStatistic() { 605 Model model = Model(5); 606 CPPUNIT_ASSERT(model.countSummaryStatistics() == 0); 607 model.addSummaryStatistic(std::make_shared<TMRCA>()); 608 CPPUNIT_ASSERT(model.countSummaryStatistics() == 1); 609 } 610 testSetLocusLength()611 void testSetLocusLength() { 612 Model model = Model(5); 613 model.setLocusLength(10); 614 CPPUNIT_ASSERT_EQUAL((size_t)10, model.loci_length() ); 615 model.setMutationRate(5.0, true, false); 616 CPPUNIT_ASSERT_EQUAL(0.5, model.mutation_rate()); 617 618 model.setLocusLength(2000); 619 CPPUNIT_ASSERT_EQUAL((size_t)2000, model.loci_length() ); 620 CPPUNIT_ASSERT_EQUAL(0.0025, model.mutation_rate()); 621 622 model.setLocusLength(5000); 623 CPPUNIT_ASSERT_EQUAL((size_t)5000, model.loci_length() ); 624 CPPUNIT_ASSERT_EQUAL(0.001, model.mutation_rate()); 625 } 626 testAddPopToVectorList()627 void testAddPopToVectorList() { 628 Model model = Model(5); 629 std::vector<std::vector<double> > vector_list; 630 vector_list.push_back(std::vector<double>(2, 1.0)); 631 vector_list.push_back(std::vector<double>()); 632 vector_list.push_back(std::vector<double>(2, 0.5)); 633 634 model.addPopToVectorList(vector_list); 635 CPPUNIT_ASSERT_EQUAL( (size_t)3, vector_list.at(0).size() ); 636 CPPUNIT_ASSERT_EQUAL( 1.0, vector_list.at(0).at(0) ); 637 CPPUNIT_ASSERT_EQUAL( 1.0, vector_list.at(0).at(1) ); 638 CPPUNIT_ASSERT( std::isnan(vector_list.at(0).at(2))); 639 640 CPPUNIT_ASSERT_EQUAL( (size_t)3, vector_list.at(2).size() ); 641 CPPUNIT_ASSERT_EQUAL( 0.5, vector_list.at(2).at(0) ); 642 CPPUNIT_ASSERT_EQUAL( 0.5, vector_list.at(2).at(1) ); 643 CPPUNIT_ASSERT( std::isnan(vector_list.at(2).at(2))); 644 } 645 testAddPopToMatrixList()646 void testAddPopToMatrixList() { 647 Model model = Model(5); 648 std::vector<std::vector<double> > vector_list; 649 vector_list.push_back(std::vector<double>(2, 1.0)); 650 vector_list.push_back(std::vector<double>()); 651 vector_list.push_back(std::vector<double>(2, 0.5)); 652 653 model.set_population_number(3); 654 model.addPopToMatrixList(vector_list, 2); 655 656 CPPUNIT_ASSERT_EQUAL( (size_t)6, vector_list.at(0).size() ); 657 CPPUNIT_ASSERT_EQUAL( 1.0, vector_list.at(0).at(0) ); 658 CPPUNIT_ASSERT_EQUAL( 1.0, vector_list.at(0).at(2) ); 659 CPPUNIT_ASSERT( std::isnan(vector_list.at(0).at(1))); 660 CPPUNIT_ASSERT( std::isnan(vector_list.at(0).at(3))); 661 CPPUNIT_ASSERT( std::isnan(vector_list.at(0).at(4))); 662 CPPUNIT_ASSERT( std::isnan(vector_list.at(0).at(5))); 663 664 CPPUNIT_ASSERT_EQUAL( (size_t)6, vector_list.at(2).size() ); 665 CPPUNIT_ASSERT_EQUAL( 0.5, vector_list.at(2).at(0) ); 666 CPPUNIT_ASSERT_EQUAL( 0.5, vector_list.at(2).at(2) ); 667 CPPUNIT_ASSERT( std::isnan(vector_list.at(0).at(1))); 668 CPPUNIT_ASSERT( std::isnan(vector_list.at(0).at(3))); 669 CPPUNIT_ASSERT( std::isnan(vector_list.at(0).at(4))); 670 CPPUNIT_ASSERT( std::isnan(vector_list.at(0).at(5))); 671 } 672 }; 673 674 //Uncomment this to activate the test 675 CPPUNIT_TEST_SUITE_REGISTRATION( TestModel ); 676