When indexing/counting mutations at the end of a simulation, use

tables to count if neutral variants are present.

fwdpy11/src/evolve_population/index_and_count_mutations.cc

@@ -2,28 +2,58 @@
 #include <fwdpy11/types/DiploidPopulation.hpp>
 #include <fwdpp/ts/count_mutations.hpp>
 #include <fwdpp/internal/sample_diploid_helpers.hpp>
+#include <stdexcept>
 #include "index_and_count_mutations.hpp"
 
 void
 index_and_count_mutations(const bool suppress_edge_table_indexing,
                           fwdpy11::DiploidPopulation& pop)
 {
     pop.mcounts_from_preserved_nodes.resize(pop.mutations.size(), 0);
-    if (!suppress_edge_table_indexing)
+    bool neutral_variants_present = false;
+    for (auto& m : pop.tables->mutations)
+        {
+            if (m.neutral)
+                {
+                    neutral_variants_present = true;
+                    break;
+                }
+        }
+
+    // FIXME: there must be a way to remove this awkward check?
+    if (!suppress_edge_table_indexing && !neutral_variants_present)
         {
             return;
         }
-    if (pop.tables->preserved_nodes.empty())
+    if (pop.tables->preserved_nodes.empty() || neutral_variants_present)
         {
             pop.tables->build_indexes();
             pop.fill_alive_nodes();
-            fwdpp::ts::count_mutations(*pop.tables, pop.mutations,
-                                       pop.alive_nodes, pop.mcounts);
+            if (pop.tables->preserved_nodes.empty())
+                {
+                    fwdpp::ts::count_mutations(*pop.tables, pop.mutations,
+                                               pop.alive_nodes, pop.mcounts);
+                }
+            else
+                {
+                    fwdpp::ts::count_mutations(*pop.tables, pop.mutations,
+                                               pop.alive_nodes, pop.mcounts,
+                                               pop.mcounts_from_preserved_nodes);
+                }
+            for (auto& mr : pop.tables->mutations)
+                {
+                    if (pop.mcounts[mr.key] + pop.mcounts_from_preserved_nodes[mr.key]
+                        == 0)
+                        {
+                            throw std::runtime_error(
+                                "mutation in table has count of zero.");
+                        }
+                }
         }
     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);
         }
 }