Add new function to choose origin time of the new mutation:

0. If the motation's node's time can be represented exactly as
   int64, the use that representation as the time.
1. If the mutation's node's parent is NULL, choose
   the closest int64 <= the node's time.
2. Otherwise, choose a uniform int64 <= the node time and >
   the parent node's time.

Fixes #798

fwdpy11/_types/diploid_population.py

@@ -389,7 +389,11 @@ def add_mutation(
           exactly `ndescendants` alive nodes.
         * From this list of candidate nodes, one is chosen randomly
           to be the node where we place the new mutation.
-        * This node's time is the origin time of the new mutation.
+        * If the node's parent's time is NULL, then the mutations' origin
+          time is the first 64 bit integer ancestral to the node time,
+          which could be the node time itself.
+          Otherwise, a 64 bit integer is chosen uniformly between the
+          node time and the node's parent's time.
         * If `deme` is None, any set of `ndescendants` will be considered.
         * If `deme` is :math:`\geq 0`, all `ndescendants` alive nodes must be
           from that deme.

fwdpy11/src/functions/add_mutation.cc

@@ -158,6 +158,46 @@ namespace
     }
 }
 
+std::int64_t
+generate_mutation_time(const fwdpy11::GSLrng_t& rng,
+                       const std::vector<fwdpp::ts::node>& node_table,
+                       const std::int32_t mutation_node,
+                       const std::int32_t mutation_node_parent)
+{
+    auto mutation_node_time = node_table[mutation_node].time;
+    // Check for easy case where we can assign mutation time == node time
+    if (mutation_node_time < 0.0)
+        {
+            if (mutation_node_time - std::floor(mutation_node_time) == 0.0)
+                {
+                    return static_cast<std::int64_t>(std::floor(mutation_node_time));
+                }
+        }
+    else if (mutation_node_time - std::ceil(mutation_node_time) == 0.0)
+        {
+            return static_cast<std::int64_t>(std::ceil(mutation_node_time));
+        }
+
+    if (mutation_node_parent == fwdpp::ts::NULL_INDEX)
+        {
+            // Time is the closest int64_t to the node time
+            if (mutation_node_time < 0.0)
+                {
+                    return static_cast<std::int64_t>(std::floor(mutation_node_time));
+                }
+            return static_cast<std::int64_t>(std::ceil(mutation_node_time));
+        }
+    // choose a random time
+    auto candidate_time = gsl_ran_flat(rng.get(), mutation_node_time,
+                                       node_table[mutation_node_parent].time);
+    if (candidate_time < 0.0)
+        {
+            return static_cast<std::int64_t>(std::floor(mutation_node_time));
+        }
+    return static_cast<std::int64_t>(std::ceil(mutation_node_time));
+}
+
+// Returns size_t max value if no candidates are found.
 std::size_t
 add_mutation(const fwdpy11::GSLrng_t& rng, const double left, const double right,
              const fwdpp::ts::table_index_t ndescendants,
@@ -226,18 +266,21 @@ add_mutation(const fwdpy11::GSLrng_t& rng, const double left, const double right
         = static_cast<std::size_t>(gsl_ran_flat(rng.get(), 0, candidates.size()));
 
     auto candidate_data = std::move(candidates[candidate]);
-    // TODO: what to do if candidates.empty() == true?
+    auto mutation_node_time = generate_mutation_time(
+        rng, pop.tables->nodes, candidate_data.node, candidate_data.parent);
+    if (mutation_node_time > pop.tables->nodes[candidate_data.node].time)
+        {
+            throw std::runtime_error(
+                "origin time for mutation is invalid: it must be <= node time");
+        }
 
     // NOTE: can we do all of what we need using existing
     // mutate/recombine functions?  Gotta check the back-end.
 
     // The new variant gets a derived state of 1
     auto empty = fwdpp::empty_mutation_queue();
     new_mutation_key = fwdpy11::infsites_Mutation(
-        empty, pop.mutations, pop.mut_lookup, false,
-        // The new mutation's origin time == the node time
-        static_cast<fwdpy11::mutation_origin_time>(
-            pop.tables->nodes[candidate_data.node].time),
+        empty, pop.mutations, pop.mut_lookup, false, mutation_node_time,
         // The lambdas will simply transfer input parameters to the new object.
         [&rng, &candidate_data]() {
             return gsl_ran_flat(rng.get(), candidate_data.left, candidate_data.right);