---
title: "hOUwie performance and stability improvements"
author: "James Boyko"
output:
  pdf_document:
    fig_caption: yes
vignette: >
   %\VignetteEngine{knitr::rmarkdown}
   %\VignetteIndexEntry{hOUwie performance and stability improvements}
   \usepackage[utf8]{inputenc}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
  echo = FALSE,
  warning = FALSE,
  message = FALSE,
  fig.align = "center",
  out.width = "90%"
)

performance <- data.frame(
  scenario = c(
    "Moderate tree",
    "More histories",
    "1,000-tip tree",
    "16 states / 74 parameters",
    "500 tips / 8 states"
  ),
  tips = c(96, 96, 1000, 250, 500),
  histories = c(25, 100, 25, 25, 50),
  states = c(4, 4, 4, 16, 8),
  parameters = c(8, 8, 8, 74, 22),
  previous_seconds = c(0.335, 1.283, 4.713, 1.908, 5.513),
  updated_seconds = c(0.131, 0.463, 1.372, 0.365, 1.267),
  stringsAsFactors = FALSE
)
performance$speedup <- performance$previous_seconds /
  performance$updated_seconds

performance_table <- data.frame(
  Scenario = performance$scenario,
  Tips = performance$tips,
  Histories = performance$histories,
  States = performance$states,
  Parameters = performance$parameters,
  `Previous version (s)` = sprintf("%.3f", performance$previous_seconds),
  `Updated version (s)` = sprintf("%.3f", performance$updated_seconds),
  Speedup = sprintf("%.2fx", performance$speedup),
  check.names = FALSE
)
```

hOUwie now runs between `r sprintf("%.2f", min(performance$speedup))` and
`r sprintf("%.2f", max(performance$speedup))` times faster across the benchmark
set. The largest improvements appear on large trees and models with many states
and parameters.

## How much faster is it?

```{r performance-table}
knitr::kable(
  performance_table,
  align = c("l", rep("r", 7)),
  caption = "Median elapsed time across repeated runs. Lower times are better."
)
```

```{r performance-plot, fig.width=8, fig.height=4.8, fig.cap="Speedup of the updated implementation over the previous version."}
bar_colors <- ifelse(performance$speedup >= 4, "#2b8cbe", "#7bccc4")
old_par <- par(mar = c(8, 4.2, 1, 0.5))
bars <- barplot(
  performance$speedup,
  names.arg = performance$scenario,
  las = 2,
  ylim = c(0, max(performance$speedup) * 1.18),
  ylab = "Speedup (times faster)",
  col = bar_colors,
  border = NA
)
abline(h = 1, lty = 2, col = "grey50")
text(
  bars,
  performance$speedup,
  labels = sprintf("%.2fx", performance$speedup),
  pos = 3,
  cex = 0.9
)
par(old_par)
```

The moderate 96-tip analysis is 2.6 times faster. Increasing the number of
histories from 25 to 100 raises that to 2.8 times. On the larger and
more parameter-rich cases, the updated implementation is 3.4 to 5.2 times
faster. A 1,000-tip run falls from 4.713 seconds to 1.372 seconds, while the
16-state, 74-parameter case falls from 1.908 seconds to 0.365 seconds.

Every run produced the requested number of stochastic histories. The median
size of the returned history objects changed by less than one percent in every
case.

## Where the speedup comes from

Most of the improvement comes from preparing the tree structure once and
reusing it. A phylogenetic tree has the same topology throughout a fit. The
updated code prepares its paths, descendants, edge positions, and state
relationships up front, then keeps compact lookup tables for the calculations
that follow.

- Tree traversals and tip-to-root paths are prepared once instead of rebuilt for
  every stochastic history.

- Repeated searches through edge and node vectors are replaced with direct
  indexing.

- The likelihood calculation works with a compact representation of each
  history. Full painted trees are still assembled for the returned result, but
  they no longer need to be the working format for every intermediate step.

- Matrix and cache bookkeeping is lighter, especially when the number of states
  and parameters grows.

Together, these changes remove a large amount of repeated work from the inner
likelihood loop. The 16-state analysis benefits most: the updated implementation
spends more of its time on the statistical calculation and less on repeatedly
organizing the same tree.

## A steadier optimization problem

hOUwie estimates a likelihood by averaging over sampled evolutionary histories.
That sampling introduces a little numerical movement from one evaluation to the
next. An optimizer can mistake that movement for a meaningful improvement and
spend time following noise.

The updated implementation compares candidate parameter values with a shared
underlying random stream within each optimizer start. The histories still
respond to the candidate parameters; they are not frozen. The shared stream
simply makes one candidate's likelihood more directly comparable with the next.
This behavior is enabled by default and can be controlled with the
`common_random_numbers` argument to `hOUwie()`.

On an easy example, both approaches found the same answer. On a more difficult
hidden-state model, the steadier comparison produced better and more consistent
fits on fresh evaluations:

| Measure | Ordinary sampling | Shared random stream |
|:--|--:|--:|
| Mean fresh-evaluation log likelihood | -27.9244 | -27.2323 |
| Median fresh-evaluation log likelihood | -28.0759 | -27.1515 |
| Between-run standard deviation | 0.6145 | 0.2669 |
| Mean likelihood evaluations | 524.0 | 438.5 |
| Mean elapsed time (seconds) | 75.24 | 61.48 |

The steadier search used 16.3% fewer likelihood evaluations and 18.3% less wall
time. It found the better fresh-evaluated solution in four of six paired starts.
Multiple starts remain important when hidden states create weakly identified
parameters or local optima.

## Reliability improvements around the likelihood

The fitting path also includes the following reliability improvements:

- Cached likelihoods use the full parameter vector, so distinct candidate
  models cannot accidentally share a result.

- Cache capacity is bounded more safely for long optimization runs.

- Failed optimizer results are filtered before selecting the best fit.

- Trait data are matched to tree tips by species name rather than assumed row
  position.

- History identifiers work with ten or more states.

- Negative trait values are shifted consistently when a model requires positive
  values.

- Input checks and parameter bounds reject several malformed or impossible
  configurations earlier and more clearly.

The expanded suite contains 235 expectations across 80 tests, covering the
faster execution path and these edge cases.