Zero division error (with a putative solution)

[hanbin973]

A ZeroDivisionError happens in the following example (tsdate version 0.2.1)

import tskit
import msprime
import tsinfer
import tsdate

# simulate tree sequence
demography = msprime.Demography()
demography.add_population(name="Panmictic", initial_size=100_000, growth_rate=0.01)
num_individuals = 15000
seq_length = 1e6
ts = msprime.sim_ancestry(samples={"Panmictic": num_individuals},
                          sequence_length=seq_length,
                          recombination_rate=1e-8,
                          demography=demography,
                          ploidy=2,
                          random_seed=1)
ts = msprime.sim_mutations(ts,
                           rate=1e-8,
                           random_seed=1)

# obtain inferred tree
sample_data = tsinfer.SampleData.from_tree_sequence(ts)
inferred_ts = tsinfer.infer(sample_data)

# date nodes
simplified_ts = tsdate.preprocess_ts(inferred_ts)
dated_ts = tsdate.variational_gamma(simplified_ts, mutation_rate=1e-8, rescaling_intervals=1000) # 1000 is default

The following error occurs at the last line

File ~/.conda/envs/jax_main/lib/python3.12/site-packages/tsdate/variational.py:759, in ExpectationPropagation.rescale(self, rescale_intervals, rescale_segsites, rescale_iterations, quantile_width, progress)
    755 _, unique = np.unique(rescaled_nodes_time, return_index=True)
    756 original_breaks = piecewise_scale_point_estimate(
    757     rescaled_breaks, rescaled_nodes_time[unique], nodes_time[unique]
    758 )
--> 759 self.node_posterior[:] = piecewise_scale_posterior(
    760     self.node_posterior,
    761     original_breaks,
    762     rescaled_breaks,
    763     quantile_width,
    764 )
    765 self.mutation_posterior[:] = piecewise_scale_posterior(
    766     self.mutation_posterior,
    767     original_breaks,
    768     rescaled_breaks,
    769     quantile_width,
    770 )

ZeroDivisionError: division by zero

This is solved if rescaling_intervals is set to a smaller number such as 100

dated_ts = tsdate.variational_gamma(simplified_ts, mutation_rate=1e-8, rescaling_intervals=100)