Fixes #639

Fixes #646

Add extensive tests of ancient samples w/ and w/o neutral mutations.

Change test to trigger more failures.

Add run time checks of mutation table consistency.

For final mutation counting, always use tables when simulating neutral
mutations.

After simplification, build recycling queue from simplification results
if simulating neutral mutations.

index_and_count_mutations removes mutations from tables whose counts
are adjusted to zero if reset_treeseqs_to_alive_nodes_after_simplification is true.

Only adjust mcounts_from_preserved_nodes using
last_preserved_generation_counts if
reset_treeseqs_to_alive_nodes_after_simplification is false.

Only adjust mutation counts by last_preserved_generation_counts if
not simulating neutral mutations.

fwdpy11/CMakeLists.txt

@@ -82,7 +82,8 @@ set(EVOLVE_POPULATION_SOURCES src/evolve_population/init.cc
     src/evolve_population/no_stopping.cc
     src/evolve_population/remove_extinct_mutations.cc
     src/evolve_population/track_ancestral_counts.cc
-    src/evolve_population/remove_extinct_genomes.cc)
+    src/evolve_population/remove_extinct_genomes.cc
+    src/evolve_population/runtime_checks.cc)
 
 set(DISCRETE_DEMOGRAPHY_SOURCES src/discrete_demography/init.cc
     src/discrete_demography/MigrationMatrix.cc

fwdpy11/src/evolve_population/index_and_count_mutations.cc

@@ -1,29 +1,46 @@
-#include <numeric>
+#include <algorithm>
 #include <fwdpy11/types/DiploidPopulation.hpp>
+#include <fwdpp/ts/table_collection_functions.hpp>
 #include <fwdpp/ts/count_mutations.hpp>
 #include <fwdpp/internal/sample_diploid_helpers.hpp>
-#include "index_and_count_mutations.hpp"
 
 void
-index_and_count_mutations(const bool suppress_edge_table_indexing,
+index_and_count_mutations(bool suppress_edge_table_indexing,
+                          bool simulating_neutral_variants,
+                          bool reset_treeseqs_to_alive_nodes_after_simplification,
                           fwdpy11::DiploidPopulation& pop)
 {
     pop.mcounts_from_preserved_nodes.resize(pop.mutations.size(), 0);
-    if (!suppress_edge_table_indexing)
+    if (!suppress_edge_table_indexing && !simulating_neutral_variants)
         {
             return;
         }
-    if (pop.tables->preserved_nodes.empty())
+    if (pop.tables->preserved_nodes.empty() || simulating_neutral_variants)
         {
             pop.tables->build_indexes();
             pop.fill_alive_nodes();
-            fwdpp::ts::count_mutations(*pop.tables, pop.mutations,
-                                       pop.alive_nodes, pop.mcounts);
+            fwdpp::ts::count_mutations(*pop.tables, pop.mutations, pop.alive_nodes,
+                                       pop.mcounts, pop.mcounts_from_preserved_nodes);
         }
     else
         {
-            fwdpp::fwdpp_internal::process_haploid_genomes(
-                pop.haploid_genomes, pop.mutations, pop.mcounts);
+            fwdpp::fwdpp_internal::process_haploid_genomes(pop.haploid_genomes,
+                                                           pop.mutations, pop.mcounts);
+        }
+    if (reset_treeseqs_to_alive_nodes_after_simplification)
+        {
+            auto itr = std::remove_if(
+                begin(pop.tables->mutations), end(pop.tables->mutations),
+                [&pop](const auto& mr) {
+                    return pop.mcounts[mr.key] + pop.mcounts_from_preserved_nodes[mr.key]
+                           == 0;
+                });
+            auto d = std::distance(itr, end(pop.tables->mutations));
+            pop.tables->mutations.erase(itr, end(pop.tables->mutations));
+            if (d)
+                {
+                    fwdpp::ts::rebuild_site_table(*pop.tables);
+                }
         }
 }
 

fwdpy11/src/evolve_population/index_and_count_mutations.hpp

@@ -1,11 +1,11 @@
 #ifndef FWDPP_TSEVOLUTION_INDEX_AND_COUNT_HPP
 #define FWDPP_TSEVOLUTION_INDEX_AND_COUNT_HPP
 
-#include <vector>
-#include <cstdint>
 #include <fwdpy11/types/DiploidPopulation.hpp>
 
-void index_and_count_mutations(const bool suppress_edge_table_indexing,
+void index_and_count_mutations(bool suppress_edge_table_indexing,
+                               bool simulating_neutral_variants,
+                               bool reset_treeseqs_to_alive_nodes_after_simplification,
                                fwdpy11::DiploidPopulation& pop);
 
 #endif

fwdpy11/src/evolve_population/remove_extinct_mutations.cc

@@ -5,6 +5,7 @@
 #include <limits>
 #include <stdexcept>
 #include <fwdpy11/types/Population.hpp>
+#include "runtime_checks.hpp"
 
 namespace
 {
@@ -40,6 +41,7 @@ namespace
 void
 remove_extinct_mutations(fwdpy11::Population& pop)
 {
+    check_mutation_table_consistency_with_count_vectors(pop, __FILE__, __LINE__);
     std::vector<fwdpp::uint_t> summed_counts(pop.mcounts);
     std::transform(begin(pop.mcounts_from_preserved_nodes),
                    end(pop.mcounts_from_preserved_nodes), begin(summed_counts),

fwdpy11/src/evolve_population/runtime_checks.cc

@@ -0,0 +1,30 @@
+#include <sstream>
+#include <fwdpy11/types/DiploidPopulation.hpp>
+
+static std::string
+strip_unix_path(const std::string file)
+{
+    auto pos = file.find_last_of('/');
+    if (pos != std::string::npos)
+        {
+            return std::string(begin(file) + static_cast<std::ptrdiff_t>(pos) + 1,
+                               end(file));
+        }
+    return file;
+}
+
+void
+check_mutation_table_consistency_with_count_vectors(const fwdpy11::Population& pop,
+                                                    std::string file, int line)
+{
+    for (auto& mr : pop.tables->mutations)
+        {
+            if (pop.mcounts[mr.key] + pop.mcounts_from_preserved_nodes[mr.key] == 0)
+                {
+                    std::ostringstream msg;
+                    msg << "mutation table is inconsistent with count vectors: "
+                        << strip_unix_path(file) << ", line " << line;
+                    throw std::runtime_error(msg.str());
+                }
+        }
+}

fwdpy11/src/evolve_population/runtime_checks.hpp

@@ -0,0 +1,6 @@
+#include <string>
+#include <fwdpy11/types/Population.hpp>
+
+void check_mutation_table_consistency_with_count_vectors(const fwdpy11::Population& pop,
+                                                         std::string file, int line);
+

fwdpy11/src/evolve_population/with_tree_sequences.cc

@@ -49,6 +49,7 @@
 #include "remove_extinct_mutations.hpp"
 #include "track_ancestral_counts.hpp"
 #include "remove_extinct_genomes.hpp"
+#include "runtime_checks.hpp"
 
 namespace py = pybind11;
 namespace ddemog = fwdpy11::discrete_demography;
@@ -105,19 +106,26 @@ simplification(
 void
 final_population_cleanup(
     bool suppress_edge_table_indexing, bool preserve_selected_fixations,
-    bool remove_extinct_mutations_at_finish, std::uint32_t last_preserved_generation,
+    bool remove_extinct_mutations_at_finish, bool simulating_neutral_variants,
+    bool reset_treeseqs_to_alive_nodes_after_simplification,
+    std::uint32_t last_preserved_generation,
     const std::vector<std::uint32_t> &last_preserved_generation_counts,
     fwdpy11::DiploidPopulation &pop)
 {
-    index_and_count_mutations(suppress_edge_table_indexing, pop);
-    if (pop.generation == last_preserved_generation)
+    index_and_count_mutations(suppress_edge_table_indexing, simulating_neutral_variants,
+                              reset_treeseqs_to_alive_nodes_after_simplification, pop);
+    check_mutation_table_consistency_with_count_vectors(pop, __FILE__, __LINE__);
+    if (pop.generation == last_preserved_generation &&
+            !reset_treeseqs_to_alive_nodes_after_simplification &&
+            !simulating_neutral_variants)
         {
             std::transform(begin(pop.mcounts_from_preserved_nodes),
                            end(pop.mcounts_from_preserved_nodes),
                            begin(last_preserved_generation_counts),
                            begin(pop.mcounts_from_preserved_nodes),
                            std::minus<std::uint32_t>());
         }
+    check_mutation_table_consistency_with_count_vectors(pop, __FILE__, __LINE__);
     if (!preserve_selected_fixations)
         {
             auto itr = std::remove_if(
@@ -431,7 +439,9 @@ evolve_with_tree_sequences(
                     new_edge_buffer.reset(nullptr);
                     final_population_cleanup(
                         suppress_edge_table_indexing, preserve_selected_fixations,
-                        remove_extinct_mutations_at_finish, last_preserved_generation,
+                        remove_extinct_mutations_at_finish, simulating_neutral_variants,
+                        reset_treeseqs_to_alive_nodes_after_simplification,
+                        last_preserved_generation,
                         last_preserved_generation_counts, pop);
                     std::ostringstream o;
                     o << "extinction at time " << pop.generation;
@@ -502,8 +512,20 @@ evolve_with_tree_sequences(
                                             end(simplification_rv.second),
                                             std::numeric_limits<std::size_t>::max()),
                                 end(simplification_rv.second));
-                            genetics.mutation_recycling_bin = fwdpp::ts::make_mut_queue(
-                                pop.mcounts, pop.mcounts_from_preserved_nodes);
+                            if (simulating_neutral_variants)
+                                {
+                                    genetics.mutation_recycling_bin
+                                        = fwdpp::ts::make_mut_queue(
+                                            simplification_rv.second,
+                                            pop.mutations.size());
+                                }
+                            else
+                                {
+                                    genetics.mutation_recycling_bin
+                                        = fwdpp::ts::make_mut_queue(
+                                            pop.mcounts,
+                                            pop.mcounts_from_preserved_nodes);
+                                }
                         }
                     else
                         {
@@ -577,7 +599,8 @@ evolve_with_tree_sequences(
     simplifier_state.reset(nullptr);
     new_edge_buffer.reset(nullptr);
     final_population_cleanup(suppress_edge_table_indexing, preserve_selected_fixations,
-                             remove_extinct_mutations_at_finish,
+                             remove_extinct_mutations_at_finish, simulating_neutral_variants, 
+                             reset_treeseqs_to_alive_nodes_after_simplification,
                              last_preserved_generation, last_preserved_generation_counts,
                              pop);
     ddemog::save_model_state(std::move(current_demographic_state), demography);

tests/conftest.py

@@ -21,10 +21,26 @@
 import fwdpy11
 
 
-@pytest.fixture
+@pytest.fixture(scope="function")
 def rng(request):
     try:
         seed = request.param.get("seed", None)
         return fwdpy11.GSLrng(seed)
-    except: # NOQA
+    except AttributeError:  # NOQA
         return fwdpy11.GSLrng(42)
+
+
+@pytest.fixture(scope="function")
+def pop(request):
+    popsize = request.param["N"]
+    genome_length = request.param["genome_length"]
+    return fwdpy11.DiploidPopulation(popsize, genome_length)
+
+
+@pytest.fixture(scope="function")
+def mslike_pop(request):
+    try:
+        N = request.param.get("N", None)
+        return fwdpy11.DiploidPopulation(N, 1.0)
+    except AttributeError:  # NOQA
+        return fwdpy11.DiploidPopulation(1000, 1.0)