Add new estimate for the expected number of inversions

[pribiller]

This pull request adds a function to compute the expected number of inversions separating two genomes, based on the number of cycles in their breakpoint graph. It implements Theorem 3 from Berestycki and Durrett (2006).
It also updates the existing documentation for all functions in the 'Similarity indexes' and 'Rearrangement distances' families to include a reference to the new function.

[charles-plessy]

Quick comments before going to exercise.

[charles-plessy]

This is actually where I would use return():

if ((N < 2) | (k  < 0)) return(NA) # Invalid value (negative number of inversions or not enough markers).
if            (k == 0)  return(N)  # No inversions. Genomes are the same.

This way the main part of the code does not need intendation and we save a few line of screen space.

[pribiller]

Thank you! I updated it.

[charles-plessy]

If it only depends on n would it be useful to memoise the results at build time?

[pribiller]

Hi Charles! Thanks for the comment. Yes, the purpose of this function is to enable memoization. In the main function (inversionEstimate_BD), the chromosome pairs are sorted by the number of aligned blocks. If two or more chromosome pairs share the same number of aligned blocks, this function is called only once, and its output is reused for all of them. I am not sure if this is what you meant, but the intent was to avoid recomputing the same table of values multiple times.

To add some context, before integrating it into GenomicBreaks, this setup step was computed considering all chromosome pairs across all pairwise alignments. In that setting, this sort of memoization provided a clearer time benefit. Here, however, the computation is restricted to chromosome pairs within a single alignment, so I do not expect significant gains, since these chromosomes will likely have different numbers of aligned blocks.

Because of that, rather than running this setup function (which requires something on the order of N operations, where N is the number of aligned blocks), I implemented another function that does a sort of binary search to find the most likely number of inversions. This second function, called when mode="single" in inversionEstimate_BD, requires something on the order of log N operations per chromosome pair. For most cases in GenomicBreaks, I expect this second method to be faster, so I set it as the default.

[charles-plessy]

What I mean is that if you do something like this for a larger interval than 1:100

inversionEstimate_BD_setup_list <- lapply(1:100, inversionEstimate_BD_setup)

and save this as a data object in the package like the structural variant examples, then you can call inversionEstimate_BD_setup_list[[n]] instead of inversionEstimate_BD_setup(n) when n is in your interval and it will just make a lookup in a pre-computed table instead of computing the value.

This said, it looks like inversionEstimate_BD_setup is not that slow unless n is larger than 20,000, so maybe this optimisation is not needed…

(Also I do not know how big would be the object and how slow would be the lookup).

[charles-plessy]

Also, if in the output of inversionEstimate_BD_setup, "k" is just a monotonic index, maybe you can just return a numeric vector with the values of "nb_cyc_k"

> inversionEstimate_BD_setup_ <- \(n) {l <- inversionEstimate_BD_setup(n) ; lapply(l, \(out) out$nb_cyc_k)}
> lapply(1:4, inversionEstimate_BD_setup_) |> as("NumericList")
NumericList of length 4
[[1]] 2 0.984462234842717
[[2]] 3 1.99062437898994 1.03427740641114 0.232605910104653
[[3]] 4 2.99701081160137 1.96892446968543 1.14265787960311 0.371790526805333
[[4]] 5 3.99904914769039 2.97044701099091 1.99661738721024 1.27460728026869 0.518574633717348

[pribiller]

Thanks for the clarification! Yes, having the saved pre-computed object is a good idea for speeding up the estimation, but I agree it might not be needed for our dataset. Perhaps we could first run the implementation as it is on some clades and check whether the performance is acceptable. Then, if needed, we could add the pre-computed object. I also have a few other ideas for optimizing this method that I can implement if it turns out to be necessary.

Regarding your other suggestions, I will add them to the code and push later.

[pribiller]

Hi Charles! k is monotonic increasing but not consecutive: as the number of events k increases, the expected number of cycles requires more events to be reduced by one. To avoid computing many values of k that yield approximately the same expected number of cycles, the method skips some ks.
For example, for inversionEstimate_BD_setup(8), the computed ks are: 1, ..., 8, 9, 11 (10 is not computed). As the input n increases, more k values are skipped.

[pribiller]

Hi Charles! Thank you for the suggestions. I have pushed an updated version, if you have any additional comments I can amend the commit.