From fd8407b5b05f6dbf29e3759252463cbdb375ad2d Mon Sep 17 00:00:00 2001
From: Saras Windecker <smwindecker@users.noreply.github.com>
Date: Sun, 5 Jul 2026 13:13:57 +0300
Subject: [PATCH 01/14] Fix proportion_infections date-alignment check for
 epiwave_greta_timeseries

The date-alignment check in prepare_observation_data() assumed any already-
classed epiwave_timeseries object stores its dates in a flat $date column,
but epiwave_greta_timeseries objects (e.g. the IHR-from-CHR pattern used for
hospitalisation proportions in tests/test_workflow/testing.R and test2.R,
built via create_epiwave_greta_timeseries()) are a list wrapping a greta
array alongside the date tibble, with dates nested at $timeseries$date
instead. The check was comparing as.Date(NULL) against target_infection_dates
and always failing with "`proportion_infections` dates must match
`target_infection_dates`".

Found by actually running tests/test_workflow/testing.R's basic
(cases + hospitalisations) fit against real data in ../data -- this is a
core supported pattern (proportion driven by a greta array), not an edge
case, and would have broken on the first real multi-stream fit using it.

Verified: the same fit now runs to completion against real data.
---
 R/prepare_observation_data.R | 16 +++++++++++++---
 1 file changed, 13 insertions(+), 3 deletions(-)

diff --git a/R/prepare_observation_data.R b/R/prepare_observation_data.R
index 760b7e4..b995623 100644
--- a/R/prepare_observation_data.R
+++ b/R/prepare_observation_data.R
@@ -36,12 +36,22 @@ prepare_observation_data <- function (observation_data,
   }
 
   prop <- observation_data$proportion_infections
-  if (!('epiwave_timeseries' %in% class(prop))) {
+  if (!inherits(prop, 'epiwave_timeseries')) {
     prop <- create_epiwave_timeseries(
       dates = target_infection_dates,
       value = prop)
-  } else if (!identical(as.Date(prop$date), as.Date(target_infection_dates))) {
-    stop('`proportion_infections` dates must match `target_infection_dates`')
+  } else {
+    # epiwave_greta_timeseries stores dates nested under $timeseries$date
+    # (it wraps a greta array alongside the date tibble), rather than a
+    # flat $date column like epiwave_timeseries/epiwave_fixed_timeseries
+    prop_dates <- if (inherits(prop, 'epiwave_greta_timeseries')) {
+      prop$timeseries$date
+    } else {
+      prop$date
+    }
+    if (!identical(as.Date(prop_dates), as.Date(target_infection_dates))) {
+      stop('`proportion_infections` dates must match `target_infection_dates`')
+    }
   }
 
   case_vec <- as_matrix(observation_data$timeseries_data, target_infection_dates)

From c804792ae647d7376d1898fa3cd0fe3fe4f3dc81 Mon Sep 17 00:00:00 2001
From: Saras Windecker <smwindecker@users.noreply.github.com>
Date: Sun, 5 Jul 2026 13:20:59 +0300
Subject: [PATCH 02/14] Add as_epiwave_timeseries(), auto-coerce raw
 observation data tables

Adds as_epiwave_timeseries(data): a new exported helper that takes a plain
data.frame/tibble with date/value columns (partial date coverage is fine)
and classes it as epiwave_fixed_timeseries -- replacing the fragile,
unvalidated class(x) <- c('epiwave_fixed_timeseries', 'epiwave_timeseries',
class(x)) pattern users previously had to hand-write for every observation
stream. Already-classed epiwave_timeseries objects and plain numeric
values/vectors pass through unchanged, so it's safe to call unconditionally.

prepare_observation_data() now calls this automatically on timeseries_data
and size_vec (mirroring the auto-coercion that delay_from_infection/
proportion_infections already had), so users can pass a plain data.frame
straight into define_observation_data()/define_sero_data() without ever
manually setting a class themselves. This directly targets the messy
"before define_observation_data()" data prep in tests/test_workflow/
testing.R, where every observation stream needed its own hand-rolled
class(x) <- ... line.

Verified: a raw unclassed data.frame produces identical output to the
equivalent pre-classed object; a data.frame missing the required columns
errors with a clear message instead of failing deep inside as_matrix().
---
 NAMESPACE                       |  1 +
 R/create_epiwave_timeseries.R   | 39 +++++++++++++++++++++++++++++++++
 R/define_observation_data.R     | 10 ++++++---
 R/prepare_observation_data.R    |  9 ++++++--
 man/as_epiwave_timeseries.Rd    | 30 +++++++++++++++++++++++++
 man/define_observation_data.Rd  |  4 +++-
 man/define_sero_data.Rd         |  6 +++--
 man/prepare_observation_data.Rd |  3 +++
 8 files changed, 94 insertions(+), 8 deletions(-)
 create mode 100644 man/as_epiwave_timeseries.Rd

diff --git a/NAMESPACE b/NAMESPACE
index 9590b46..2fec4de 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -3,6 +3,7 @@
 S3method(as_matrix,epiwave_greta_timeseries)
 S3method(as_matrix,epiwave_timeseries)
 S3method(as_matrix,numeric)
+export(as_epiwave_timeseries)
 export(as_matrix)
 export(compute_reff)
 export(create_dow_priors)
diff --git a/R/create_epiwave_timeseries.R b/R/create_epiwave_timeseries.R
index b2abba4..049aa30 100644
--- a/R/create_epiwave_timeseries.R
+++ b/R/create_epiwave_timeseries.R
@@ -86,3 +86,42 @@ create_epiwave_greta_timeseries <- function(dates,
                             class(long_combined))
   long_combined
 }
+
+#' Coerce a date/value table to an epiwave_fixed_timeseries object
+#'
+#' @description Observation data (case counts, hospitalisation counts,
+#'   seroprevalence survey counts/sample sizes) is naturally a table of
+#'   `date`/`value` pairs with real per-date observations, rather than a
+#'   single value replicated across dates -- unlike delay distributions or
+#'   proportions, which usually apply uniformly. This function takes such a
+#'   table (with genuinely partial date coverage, e.g. missing weekends, is
+#'   fine -- gaps are handled downstream by `as_matrix()`) and classes it as
+#'   an `epiwave_fixed_timeseries` object, so it doesn't need to be
+#'   hand-classed with `class(x) <- c(...)`.
+#'
+#'   Already-classed `epiwave_timeseries` objects and plain numeric
+#'   values/vectors are returned unchanged, so this is safe to call
+#'   unconditionally on anything accepted as observation data.
+#'
+#' @param data an `epiwave_timeseries` object, a numeric value/vector, or a
+#'   data.frame/tibble with `date` and `value` columns
+#'
+#' @return an `epiwave_timeseries` or numeric object, ready for `as_matrix()`
+#' @export
+as_epiwave_timeseries <- function(data) {
+
+  if (inherits(data, "epiwave_timeseries") || is.numeric(data)) {
+    return(data)
+  }
+
+  if (!all(c("date", "value") %in% names(data))) {
+    stop("`data` must have `date` and `value` columns to be used as ",
+         "observation data")
+  }
+
+  data <- tibble::as_tibble(data[c("date", "value")])
+  data$date <- as.Date(data$date)
+
+  class(data) <- c("epiwave_fixed_timeseries", "epiwave_timeseries", class(data))
+  data
+}
diff --git a/R/define_observation_data.R b/R/define_observation_data.R
index eb77e95..31d256d 100644
--- a/R/define_observation_data.R
+++ b/R/define_observation_data.R
@@ -5,7 +5,9 @@
 #'  multiple jurisdictions later via `stack_jurisdictions()`.
 #'
 #' @param timeseries_data timeseries data for data of interest, for one
-#'  jurisdiction
+#'  jurisdiction. A plain data.frame/tibble with `date` and `value` columns
+#'  is fine (gaps in date coverage are fine too) -- it doesn't need to be
+#'  pre-classed, see `as_epiwave_timeseries()`
 #' @param delay_from_infection a `discrete_pmf`/`discrete_weights` object
 #'  (replicated across dates), or an already time-varying
 #'  `discrete_pmf_series`/`discrete_weights_series` object
@@ -34,11 +36,13 @@ define_observation_data <- function (timeseries_data,
 #'  `stack_jurisdictions()`.
 #'
 #' @param timeseries_data seroprevalence survey timeseries data, for one
-#'  jurisdiction
+#'  jurisdiction. A plain data.frame/tibble with `date` and `value` columns
+#'  is fine, see `as_epiwave_timeseries()`
 #' @param total_pop population size of this jurisdiction (a single numeric
 #'  value)
 #' @param size_vec sample size of the seroprevalence survey, for one
-#'  jurisdiction
+#'  jurisdiction. A plain data.frame/tibble with `date` and `value` columns
+#'  is fine, see `as_epiwave_timeseries()`
 #' @param delay_from_infection typically a `discrete_weights`/
 #'  `discrete_weights_series` object (e.g. probability of testing
 #'  seropositive by day since infection, which need not sum to 1 -- unlike
diff --git a/R/prepare_observation_data.R b/R/prepare_observation_data.R
index b995623..ca0802b 100644
--- a/R/prepare_observation_data.R
+++ b/R/prepare_observation_data.R
@@ -6,6 +6,9 @@
 #'
 #' @param observation_data data for one jurisdiction, one stream, as
 #'  returned by `define_observation_data()`/`define_sero_data()`.
+#'  `timeseries_data` (and `size_vec`, for sero streams) may be a plain
+#'  data.frame/tibble with `date` and `value` columns -- it doesn't need to
+#'  be pre-classed as `epiwave_timeseries`, see `as_epiwave_timeseries()`.
 #'  `delay_from_infection` may be a single `discrete_pmf`/`discrete_weights`
 #'  object (replicated across `target_infection_dates`), or an already
 #'  time-varying `discrete_pmf_series`/`discrete_weights_series` (aligned to
@@ -54,7 +57,8 @@ prepare_observation_data <- function (observation_data,
     }
   }
 
-  case_vec <- as_matrix(observation_data$timeseries_data, target_infection_dates)
+  timeseries_data <- as_epiwave_timeseries(observation_data$timeseries_data)
+  case_vec <- as_matrix(timeseries_data, target_infection_dates)
   prop_vec <- as_matrix(prop, target_infection_dates)
 
   convolution_matrix <- new_convolution_matrix(delays)
@@ -81,7 +85,8 @@ prepare_observation_data <- function (observation_data,
     out$total_pop <- observation_data$total_pop
   }
   if ('size_vec' %in% names(observation_data)) {
-    out$size_vec <- as_matrix(observation_data$size_vec, target_infection_dates)
+    size_vec <- as_epiwave_timeseries(observation_data$size_vec)
+    out$size_vec <- as_matrix(size_vec, target_infection_dates)
   }
 
   out
diff --git a/man/as_epiwave_timeseries.Rd b/man/as_epiwave_timeseries.Rd
new file mode 100644
index 0000000..6f43b15
--- /dev/null
+++ b/man/as_epiwave_timeseries.Rd
@@ -0,0 +1,30 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/create_epiwave_timeseries.R
+\name{as_epiwave_timeseries}
+\alias{as_epiwave_timeseries}
+\title{Coerce a date/value table to an epiwave_fixed_timeseries object}
+\usage{
+as_epiwave_timeseries(data)
+}
+\arguments{
+\item{data}{an \code{epiwave_timeseries} object, a numeric value/vector, or a
+data.frame/tibble with \code{date} and \code{value} columns}
+}
+\value{
+an \code{epiwave_timeseries} or numeric object, ready for \code{as_matrix()}
+}
+\description{
+Observation data (case counts, hospitalisation counts,
+seroprevalence survey counts/sample sizes) is naturally a table of
+\code{date}/\code{value} pairs with real per-date observations, rather than a
+single value replicated across dates -- unlike delay distributions or
+proportions, which usually apply uniformly. This function takes such a
+table (with genuinely partial date coverage, e.g. missing weekends, is
+fine -- gaps are handled downstream by \code{as_matrix()}) and classes it as
+an \code{epiwave_fixed_timeseries} object, so it doesn't need to be
+hand-classed with \code{class(x) <- c(...)}.
+
+Already-classed \code{epiwave_timeseries} objects and plain numeric
+values/vectors are returned unchanged, so this is safe to call
+unconditionally on anything accepted as observation data.
+}
diff --git a/man/define_observation_data.Rd b/man/define_observation_data.Rd
index 5239028..e35d4e6 100644
--- a/man/define_observation_data.Rd
+++ b/man/define_observation_data.Rd
@@ -13,7 +13,9 @@ define_observation_data(
 }
 \arguments{
 \item{timeseries_data}{timeseries data for data of interest, for one
-jurisdiction}
+jurisdiction. A plain data.frame/tibble with \code{date} and \code{value} columns
+is fine (gaps in date coverage are fine too) -- it doesn't need to be
+pre-classed, see \code{as_epiwave_timeseries()}}
 
 \item{delay_from_infection}{a \code{discrete_pmf}/\code{discrete_weights} object
 (replicated across dates), or an already time-varying
diff --git a/man/define_sero_data.Rd b/man/define_sero_data.Rd
index 3113a75..363fd05 100644
--- a/man/define_sero_data.Rd
+++ b/man/define_sero_data.Rd
@@ -15,13 +15,15 @@ define_sero_data(
 }
 \arguments{
 \item{timeseries_data}{seroprevalence survey timeseries data, for one
-jurisdiction}
+jurisdiction. A plain data.frame/tibble with \code{date} and \code{value} columns
+is fine, see \code{as_epiwave_timeseries()}}
 
 \item{total_pop}{population size of this jurisdiction (a single numeric
 value)}
 
 \item{size_vec}{sample size of the seroprevalence survey, for one
-jurisdiction}
+jurisdiction. A plain data.frame/tibble with \code{date} and \code{value} columns
+is fine, see \code{as_epiwave_timeseries()}}
 
 \item{delay_from_infection}{typically a \code{discrete_weights}/
 \code{discrete_weights_series} object (e.g. probability of testing
diff --git a/man/prepare_observation_data.Rd b/man/prepare_observation_data.Rd
index 2d5579c..5840ec2 100644
--- a/man/prepare_observation_data.Rd
+++ b/man/prepare_observation_data.Rd
@@ -9,6 +9,9 @@ prepare_observation_data(observation_data, target_infection_dates)
 \arguments{
 \item{observation_data}{data for one jurisdiction, one stream, as
 returned by \code{define_observation_data()}/\code{define_sero_data()}.
+\code{timeseries_data} (and \code{size_vec}, for sero streams) may be a plain
+data.frame/tibble with \code{date} and \code{value} columns -- it doesn't need to
+be pre-classed as \code{epiwave_timeseries}, see \code{as_epiwave_timeseries()}.
 \code{delay_from_infection} may be a single \code{discrete_pmf}/\code{discrete_weights}
 object (replicated across \code{target_infection_dates}), or an already
 time-varying \code{discrete_pmf_series}/\code{discrete_weights_series} (aligned to

From 8057dc9e989a61dfa134f6f09c7119053850295b Mon Sep 17 00:00:00 2001
From: Saras Windecker <smwindecker@users.noreply.github.com>
Date: Sun, 5 Jul 2026 13:32:59 +0300
Subject: [PATCH 03/14] Remove manual class(x) <- c(...) boilerplate from
 testing.R/test2.R

define_observation_data() now accepts a plain date/value data.frame directly
(via as_epiwave_timeseries(), added earlier on this branch), so the
hand-rolled class(notif_dat) <- c('epiwave_fixed_timeseries',
'epiwave_timeseries', class(notif_dat)) lines are no longer needed.

Verified: re-ran testing.R's basic (cases + hospitalisations) fit against
real data in ../data with the classing removed -- notif_dat/hosp_dat stay
plain tbl_df objects throughout, and the fit runs to completion identically.
---
 tests/test_workflow/test2.R   | 8 ++------
 tests/test_workflow/testing.R | 8 ++------
 2 files changed, 4 insertions(+), 12 deletions(-)

diff --git a/tests/test_workflow/test2.R b/tests/test_workflow/test2.R
index f751180..0b7001c 100644
--- a/tests/test_workflow/test2.R
+++ b/tests/test_workflow/test2.R
@@ -23,9 +23,8 @@ notif_dat$date <- study_seq
 colnames(notif_dat) <- c('study_area', 'date')
 notif_dat <- notif_dat |>
   tidyr::pivot_longer(!date, names_to = "jurisdiction", values_to = "value")
-class(notif_dat) <- c('epiwave_fixed_timeseries',
-                      'epiwave_timeseries',
-                      class(notif_dat))
+# no manual class(notif_dat) <- c(...) needed -- define_observation_data()
+# accepts a plain date/value data.frame directly (see as_epiwave_timeseries())
 
 # hospitalisation counts
 hosp_dat <- 'simdata/sim_study_hosp.csv' |>
@@ -36,9 +35,6 @@ hosp_dat$date <- study_seq
 colnames(hosp_dat) <- c('study_area', 'date')
 hosp_dat <- hosp_dat |>
   tidyr::pivot_longer(!date, names_to = "jurisdiction", values_to = "value")
-class(hosp_dat) <- c('epiwave_fixed_timeseries',
-                      'epiwave_timeseries',
-                      class(hosp_dat))
 
 # jurisdictions
 jurisdictions <- unique(notif_dat$jurisdiction)
diff --git a/tests/test_workflow/testing.R b/tests/test_workflow/testing.R
index 6d2631d..2ab731c 100644
--- a/tests/test_workflow/testing.R
+++ b/tests/test_workflow/testing.R
@@ -47,9 +47,8 @@ notif_dat <- notif_file |>
 if (exists('jurisdictions')) {
   notif_dat <- notif_dat[notif_dat$jurisdiction %in% jurisdictions,]
 }
-class(notif_dat) <- c('epiwave_fixed_timeseries',
-                      'epiwave_timeseries',
-                      class(notif_dat))
+# no manual class(notif_dat) <- c(...) needed -- define_observation_data()
+# accepts a plain date/value data.frame directly (see as_epiwave_timeseries())
 
 # # hospitalisation counts
 hosp_file <- paste0(not_synced_folder, '/COVID_live_cases_in_hospital.rds')
@@ -60,9 +59,6 @@ hosp_dat <- hosp_file |>
 if (exists('jurisdictions')) {
   hosp_dat <- hosp_dat[hosp_dat$jurisdiction %in% jurisdictions,]
 }
-class(hosp_dat) <- c('epiwave_fixed_timeseries',
-                     'epiwave_timeseries',
-                     class(hosp_dat))
 
 # jurisdictions
 if (!exists('jurisdictions')) {

From bbcaf269616dd7b20c9011107e10e59cdc924077 Mon Sep 17 00:00:00 2001
From: Saras Windecker <smwindecker@users.noreply.github.com>
Date: Sun, 5 Jul 2026 13:35:56 +0300
Subject: [PATCH 04/14] Add self-contained single/multi-jurisdiction workflow
 scripts

testing.R/test2.R depend on local, not-synced data (../data, simdata/) that
isn't in the repo -- useful as personal real-data scripts, but not runnable
by anyone else who clones it, and don't exercise the multi-jurisdiction path
at all (the package has never had a multi-jurisdiction reprex before this
branch).

Adds two fully self-contained scripts (fabricated data, no external
dependencies) that demonstrate the current API end to end:

- single_jurisdiction_workflow.R: raw data.frames passed straight to
  define_observation_data() (no manual class(x) <- c(...)), delay
  distributions built via epiwave.params (discrete_pmf, combined with `+`),
  a greta-array-derived proportion (IHR-from-CHR), and an "advanced" section
  showing a time-varying discrete_pmf_series delay.
- multi_jurisdiction_workflow.R: two jurisdictions with deliberately
  staggered/non-identical date coverage, dow_model = TRUE to exercise
  hierarchical DOW pooling, and an explicit alignment check confirming a
  date only one jurisdiction has data for stays correctly non-NA/NA in the
  right columns.

Both run to completion (small MCMC settings for speed) with their
stopifnot() checks passing.
---
 .../multi_jurisdiction_workflow.R             |  90 ++++++++++++
 .../single_jurisdiction_workflow.R            | 139 ++++++++++++++++++
 2 files changed, 229 insertions(+)
 create mode 100644 tests/test_workflow/multi_jurisdiction_workflow.R
 create mode 100644 tests/test_workflow/single_jurisdiction_workflow.R

diff --git a/tests/test_workflow/multi_jurisdiction_workflow.R b/tests/test_workflow/multi_jurisdiction_workflow.R
new file mode 100644
index 0000000..3a1fea8
--- /dev/null
+++ b/tests/test_workflow/multi_jurisdiction_workflow.R
@@ -0,0 +1,90 @@
+## Multi-jurisdiction (hierarchical/partially-pooled) workflow
+##
+## A fully self-contained example fitting two jurisdictions jointly, using
+## fabricated data. Demonstrates:
+##   - jurisdiction is only ever a dimension at the fit_waves() step: each
+##     jurisdiction's data is prepared independently via
+##     define_observation_data()/define_observation_model() (exactly like the
+##     single-jurisdiction workflow), then combined into a *named* list keyed
+##     by jurisdiction and passed to fit_waves()
+##   - jurisdictions can have different/staggered date coverage -- every
+##     per-jurisdiction vector is aligned to the same shared
+##     target_infection_dates axis, so this is safe
+##   - partial pooling: create_infection_timeseries() shares GP kernel
+##     hyperparameters across jurisdictions (independent draws, shared
+##     lengthscale/variance), and create_dow_priors() shares a hierarchical
+##     day-of-week prior (nation-wide effect + per-jurisdiction deviation)
+##     when dow_model = TRUE
+##
+## For a single-jurisdiction example, see single_jurisdiction_workflow.R.
+
+devtools::load_all()
+library(epiwave.params)
+library(distributional)
+
+set.seed(7)
+
+target_infection_dates <- seq(as.Date("2024-01-01"), by = "day", length.out = 150)
+cases_delay <- as_discrete_pmf(distributional::dist_gamma(shape = 3, rate = 0.5))
+
+make_jurisdiction_cases <- function(start, end, lambda, dow_model) {
+  dates <- target_infection_dates[start:end]
+  counts <- rpois(length(dates), lambda = lambda)
+  define_observation_data(
+    timeseries_data = data.frame(date = dates, value = counts),
+    delay_from_infection = cases_delay,
+    proportion_infections = 0.5,
+    dow_model = dow_model)
+}
+
+# Jurisdiction A: cases observed days 10-90.
+# Jurisdiction B: cases observed days 60-140.
+# Deliberately staggered/non-identical coverage, to exercise the
+# master-date-axis alignment that makes joint fitting safe.
+observation_model_A <- define_observation_model(
+  target_infection_dates = target_infection_dates,
+  cases = make_jurisdiction_cases(10, 90, lambda = 50, dow_model = TRUE)
+)
+observation_model_B <- define_observation_model(
+  target_infection_dates = target_infection_dates,
+  cases = make_jurisdiction_cases(60, 140, lambda = 30, dow_model = TRUE)
+)
+
+# named list, named by jurisdiction -- this is the only place jurisdiction
+# becomes a dimension
+observations_by_jurisdiction <- list(
+  jurisdiction_a = observation_model_A,
+  jurisdiction_b = observation_model_B
+)
+
+fit_object <- fit_waves(
+  observations_by_jurisdiction = observations_by_jurisdiction,
+  infection_model_type = "gp_growth_rate",
+  n_chains = 2,
+  warmup = 200,
+  n_samples = 200
+)
+
+stopifnot(identical(dim(fit_object$infection_model), c(150L, 2L)))
+stopifnot(identical(fit_object$jurisdictions, c("jurisdiction_a", "jurisdiction_b")))
+
+## -- sanity check: staggered coverage stays correctly aligned ---------------
+##
+## A date only jurisdiction A has data for should show up as non-NA for A and
+## NA for B in the underlying case matrix, and vice versa -- confirming row i
+## means the same calendar date for every jurisdiction's column.
+
+case_mat <- observation_model_A$observation_model_data$cases$case_vec
+row_only_a <- 20   # within A's 10-90 window, outside B's 60-140 window
+row_only_b <- 135  # within B's window, outside A's
+
+stopifnot(!is.na(observation_model_A$observation_model_data$cases$case_vec[row_only_a]))
+stopifnot(is.na(observation_model_B$observation_model_data$cases$case_vec[row_only_a]))
+stopifnot(is.na(observation_model_A$observation_model_data$cases$case_vec[row_only_b]))
+stopifnot(!is.na(observation_model_B$observation_model_data$cases$case_vec[row_only_b]))
+
+rhats <- coda::gelman.diag(fit_object$fit, autoburnin = FALSE, multivariate = FALSE)
+cat("max rhat:", max(rhats$psrf[, 1], na.rm = TRUE), "\n")
+
+cat("\nMulti-jurisdiction workflow complete -- jurisdictions:",
+    fit_object$jurisdictions, "\n")
diff --git a/tests/test_workflow/single_jurisdiction_workflow.R b/tests/test_workflow/single_jurisdiction_workflow.R
new file mode 100644
index 0000000..19f9894
--- /dev/null
+++ b/tests/test_workflow/single_jurisdiction_workflow.R
@@ -0,0 +1,139 @@
+## Single-jurisdiction workflow
+##
+## A minimal, fully self-contained example of the epiwave workflow for one
+## jurisdiction, using fabricated data (no external/not-synced files needed).
+## Demonstrates the current API end to end:
+##   - raw data.frames passed straight to define_observation_data() (no
+##     manual class(x) <- c(...) needed -- see as_epiwave_timeseries())
+##   - delay distributions built with epiwave.params (discrete_pmf, combined
+##     via `+`/add_discrete()), including a time-varying discrete_pmf_series
+##   - proportion_infections as both a fixed value and a greta-array-derived
+##     value (the IHR-from-CHR pattern)
+##   - define_observation_data() -> define_observation_model() -> fit_waves()
+##
+## For a multi-jurisdiction (hierarchical/partially-pooled) example, see
+## multi_jurisdiction_workflow.R.
+
+devtools::load_all()
+library(epiwave.params)
+library(distributional)
+
+set.seed(42)
+
+## -- date range -----------------------------------------------------------
+
+target_infection_dates <- seq(as.Date("2024-01-01"), by = "day", length.out = 150)
+
+## -- delay distributions ----------------------------------------------------
+
+# incubation period + onset-to-notification delay, combined by convolution
+# (the `+` operator on discrete_pmf objects convolves them -- equivalent to
+# add_discrete(incubation, onset_to_notification))
+incubation <- as_discrete_pmf(distributional::dist_weibull(shape = 1.83, scale = 4.93))
+onset_to_notification <- as_discrete_pmf(distributional::dist_gamma(shape = 3, rate = 0.5))
+cases_delay <- incubation + onset_to_notification
+
+# hospitalisation delay: a single distributional object works directly too
+hosp_delay <- as_discrete_pmf(distributional::dist_weibull(shape = 2.51, scale = 10.17))
+
+## -- fabricated observation data --------------------------------------------
+
+# case counts: a plain date/value data.frame -- no epiwave_timeseries class
+# needed, define_observation_data() coerces it automatically
+case_dates <- target_infection_dates[20:130]
+case_counts <- rpois(length(case_dates), lambda = 50 + 20 * sin(seq_along(case_dates) / 10))
+notif_dat <- data.frame(date = case_dates, value = case_counts)
+
+# hospitalisation counts, a subset of infections some days later
+hosp_dates <- target_infection_dates[30:120]
+hosp_counts <- rpois(length(hosp_dates), lambda = 5)
+hosp_dat <- data.frame(date = hosp_dates, value = hosp_counts)
+
+## -- proportions -------------------------------------------------------------
+
+# case ascertainment rate: a fixed proportion
+car <- 0.5
+
+# hospitalisation proportion: derived from CHR (case-hospitalisation rate)
+# and CAR via a greta array, so it can be estimated jointly with the rest of
+# the model rather than fixed
+chr <- greta::uniform(0, 1)
+ihr <- create_epiwave_greta_timeseries(
+  dates = target_infection_dates,
+  car = car,
+  chr_prior = chr)
+
+## -- observation model, one jurisdiction -------------------------------------
+
+observation_model <- define_observation_model(
+  target_infection_dates = target_infection_dates,
+
+  cases = define_observation_data(
+    timeseries_data = notif_dat,
+    delay_from_infection = cases_delay,
+    proportion_infections = car,
+    dow_model = TRUE),
+
+  hospitalisations = define_observation_data(
+    timeseries_data = hosp_dat,
+    delay_from_infection = hosp_delay,
+    proportion_infections = ihr)
+)
+
+# a single-jurisdiction fit is a length-1 named list, named by jurisdiction
+observations_by_jurisdiction <- list(example_jurisdiction = observation_model)
+
+## -- fit ----------------------------------------------------------------
+
+# flat_prior is the fastest infection model, good for a quick check like this
+# one; for a real analysis, 'gp_infections'/'gp_growth_rate'/
+# 'gp_growth_rate_deriv' model infections as a Gaussian Process instead (see
+# create_infection_timeseries()'s docs for how the three GP formulations
+# differ). Warmup/n_samples are kept small here for speed, not convergence.
+fit_object <- fit_waves(
+  observations_by_jurisdiction = observations_by_jurisdiction,
+  infection_model_type = "flat_prior",
+  n_chains = 2,
+  warmup = 200,
+  n_samples = 200
+)
+
+stopifnot(identical(dim(fit_object$infection_model), c(150L, 1L)))
+
+rhats <- coda::gelman.diag(fit_object$fit, autoburnin = FALSE, multivariate = FALSE)
+cat("max rhat:", max(rhats$psrf[, 1], na.rm = TRUE), "\n")
+
+infection_draws <- greta::calculate(
+  fit_object$infection_model,
+  values = fit_object$fit,
+  nsim = 200)[[1]]
+cat("posterior median infections range:",
+    round(range(apply(infection_draws, 2, median))), "\n")
+
+## -- advanced: a time-varying delay distribution -----------------------------
+##
+## delay_from_infection can also be an already time-varying
+## discrete_pmf_series/discrete_weights_series, e.g. representing a change in
+## testing/reporting delay partway through the study period.
+
+notification_delay_early <- as_discrete_pmf(distributional::dist_gamma(shape = 3, rate = 0.5))
+notification_delay_late <- as_discrete_pmf(distributional::dist_gamma(shape = 2, rate = 0.6))
+
+time_varying_cases_delay <- new_discrete_series(
+  values = c(
+    rep(list(notification_delay_early), 75),
+    rep(list(notification_delay_late), 75)
+  ),
+  index = target_infection_dates
+)
+
+observation_model_time_varying <- define_observation_model(
+  target_infection_dates = target_infection_dates,
+  cases = define_observation_data(
+    timeseries_data = notif_dat,
+    delay_from_infection = time_varying_cases_delay,
+    proportion_infections = car,
+    dow_model = TRUE)
+)
+
+cat("\nSingle-jurisdiction workflow complete.\n")

From d67e8ee6bfb5deee34b0783fefa49158b48be691 Mon Sep 17 00:00:00 2001
From: Saras Windecker <smwindecker@users.noreply.github.com>
Date: Sun, 5 Jul 2026 23:14:34 +0300
Subject: [PATCH 05/14] Add dev/sero-integration-notes.md, exclude dev/ from
 package build

Captures why the seroprevalence pathway (define_sero_data(),
create_small_sero_model(), the discrete_weights widening in
new_convolution_matrix()/evaluate()) is being pulled back out of this MVP,
the recommended architecture for re-integrating it (merge into
define_observation_data()/create_observation_model() rather than parallel
functions with if-based dispatch), the real bugs found in the original sero
code during this work, and what validation is still needed before it ships.

dev/ is excluded from R CMD build via .Rbuildignore, matching the existing
pattern for .claude/.positai.
---
 .Rbuildignore                 |   1 +
 dev/sero-integration-notes.md | 103 ++++++++++++++++++++++++++++++++++
 2 files changed, 104 insertions(+)
 create mode 100644 dev/sero-integration-notes.md

diff --git a/.Rbuildignore b/.Rbuildignore
index 855fc34..73057a5 100644
--- a/.Rbuildignore
+++ b/.Rbuildignore
@@ -6,3 +6,4 @@
 ^README\.Rmd$
 ^\.positai$
 ^\.claude$
+^dev$
diff --git a/dev/sero-integration-notes.md b/dev/sero-integration-notes.md
new file mode 100644
index 0000000..2878389
--- /dev/null
+++ b/dev/sero-integration-notes.md
@@ -0,0 +1,103 @@
+# Seroprevalence pathway: integration notes
+
+Status as of this note: seroprevalence support was built and verified against
+synthetic data during the jurisdiction-dimension refactor, then deliberately
+pulled back out of the MVP. This file exists so that work isn't lost, and so
+whoever picks it up next doesn't need to re-derive the design from scratch.
+Delete this file once sero is actually re-integrated.
+
+## Why it was pulled
+
+1. No real-data validation exists. Everything sero-related was only ever run
+   against fabricated data in a scratch reprex. `tests/test_workflow/testing.R`
+   has *always* had a non-functional, commented-out sero block (referencing
+   `sero_dat`/`sero_size_mat`/`sero_conversion`, none of which are ever
+   assigned) -- there has never been a working real-data example for sero in
+   this package.
+2. The `discrete_weights` choice for seroconversion (see below) is a real
+   modelling decision, not just a software one, and deserves a second pair of
+   eyes (e.g. Tianxiao/Nick) before it ships.
+3. We didn't want `if`/dispatch logic for a not-yet-validated pathway sitting
+   in the MVP's core functions (`fit_waves()`, `stack_jurisdictions()`,
+   `prepare_observation_data()`). Rather than merge it in behind conditionals,
+   we're waiting until it can go in as a first-class part of the existing
+   function workflow.
+
+## The modelling decision: discrete_weights, not discrete_pmf
+
+Case/hospitalisation notification is a one-time event, correctly modelled as
+a normalised `discrete_pmf` (delay from infection to that single event).
+Seroconversion is different: a person may test positive for many consecutive
+days, so the "delay from infection" isn't a probability distribution over a
+single event -- it's a persistence/detectability curve that doesn't need to
+sum to 1. `epiwave.params::discrete_weights`/`discrete_weights_series` (via
+`as_discrete_weights()`) is the right object for this.
+
+`new_convolution_matrix()`/`evaluate()` were widened during this work to
+accept `discrete_weights`/`discrete_weights_series` alongside
+`discrete_pmf`/`discrete_pmf_series` -- mechanically identical (a
+day-difference matrix looked up against a `step` column), just using
+`$weight` instead of `$prob`, and *not* forcing normalisation. This was
+reverted when sero was pulled (no other current consumer), but re-widening
+it is a small, mechanical change -- see the `remove-jurisdiction-dimension`
+PR (#51) history for the exact diff if useful as a reference, since the
+widening itself was never in question, just its inclusion in the MVP.
+
+## Recommended architecture when re-integrating
+
+Do **not** re-add `define_sero_data()`/`create_small_sero_model()` as
+parallel functions with an `if`-based dispatch in `fit_waves()` (that's the
+shape we just removed). Instead:
+
+1. Add `total_pop`/`size_vec` as **optional** arguments to
+   `define_observation_data()` (default `NULL`), and delete the separate
+   `define_sero_data()` wrapper. `prepare_observation_data()` already treats
+   `total_pop`/`size_vec` as generically-optional fields (checks
+   `'total_pop' %in% names(observation_data)`) -- it doesn't care which
+   constructor produced the list, so this needs no changes there once
+   `total_pop`/`size_vec` are threaded through.
+2. Merge `create_small_sero_model()`'s binomial-likelihood body into
+   `create_observation_model()` as an internal branch, dispatched on
+   `is.null(observations$total_pop)` -- i.e. move the dispatch that
+   currently lived in `fit_waves()` (`model_fn <- if (!is.null(stream$total_pop))
+   create_small_sero_model else create_observation_model`) into the body of
+   a single `create_observation_model()`, and delete
+   `create_small_sero_model()`. `fit_waves()` then calls one function per
+   stream uniformly, no dispatch needed there at all.
+3. `stack_jurisdictions()`'s `stack_stream()` helper needs its `total_pop`/
+   `size_mat` stacking logic re-added (this was mechanical, not
+   conceptually hard -- see PR #51 history) -- gated on the same
+   `total_pop`-presence check, no new dispatch pattern needed.
+
+This was discussed and agreed as the target shape *before* pulling sero back
+out, specifically so this note doesn't need to re-litigate the "one function
+vs two" question -- it's already decided, just not implemented yet.
+
+## Real bugs found in the sero pathway during this work (fix again if the old code is reused as a reference)
+
+If anyone copies from the pre-removal sero code (e.g. from PR #51's history)
+rather than writing fresh, watch for these -- all were real, silent bugs in
+the *original* (pre-refactor) sero code, fixed once during this work:
+
+1. `prepare_observation_data()` stored the seroprevalence sample size as
+   `size_vec`, but `create_small_sero_model()` read a field called
+   `size_mat` -- always `NULL`. Make sure whatever field name
+   `create_observation_model()`'s binomial branch reads matches what
+   `stack_jurisdictions()` actually produces.
+2. The old `create_small_sero_model()` called
+   `new_convolution_matrix(delays, x, n_dates)`, a stale 3-argument call that
+   doesn't match the current 2-argument `new_convolution_matrix(pmf, n)`
+   signature -- would have errored if ever actually invoked.
+3. `fit_waves()`'s call to `create_small_sero_model()` was commented out
+   entirely -- the sero pathway was unreachable in the original code, not
+   just unvalidated.
+
+## What's still needed before this ships
+
+- A real (or realistic synthetic, committed as an actual test) seroprevalence
+  reprex end to end.
+- Domain sanity-check on the `discrete_weights` persistence-curve modelling
+  choice.
+- A committed test, not just a manual script -- there's no automated test
+  suite in this package at all yet (`tests/testthat/` isn't populated), so
+  this would ideally land alongside seeding that.

From 350f92e060e9e2bf146647a3cd0c33b5599aa79f Mon Sep 17 00:00:00 2001
From: Saras Windecker <smwindecker@users.noreply.github.com>
Date: Sun, 5 Jul 2026 23:22:09 +0300
Subject: [PATCH 06/14] Remove the seroprevalence pathway from the MVP
 (dispatch-free for now)

Per discussion: don't ship if-statement dispatch for a not-yet-validated
feature. Removes define_sero_data(), create_small_sero_model(), the
total_pop/size_vec plumbing through prepare_observation_data() and
stack_jurisdictions()'s stack_stream(), the model_fn dispatch in
fit_waves(), and the discrete_weights/discrete_weights_series widening in
new_convolution_matrix()/evaluate() (which had no other consumer once sero
is removed). inits_by_jurisdiction()'s mean_step() helper simplifies back to
a plain mean(x) call accordingly.

The design for how to re-integrate this cleanly (merge into
define_observation_data()/create_observation_model() rather than parallel
functions with if-based dispatch, as agreed before removing it), the real
bugs found in the original sero code, and what's still needed before it
ships are captured in dev/sero-integration-notes.md (previous commit).

fit_waves()/stack_jurisdictions()/prepare_observation_data() now only
handle the validated cases/hospitalisations pathway.
---
 R/create_epiwave_timeseries.R |  6 +--
 R/create_small_sero_model.R   | 82 -----------------------------------
 R/define_observation_data.R   | 46 +-------------------
 R/define_observation_model.R  |  4 +-
 R/evaluate.R                  | 38 ++++------------
 R/fit_waves.R                 | 13 ++----
 R/inits_by_jurisdiction.R     | 16 ++-----
 R/new_convolution_matrix.R    | 51 +++++++++-------------
 R/prepare_observation_data.R  | 48 ++++++++------------
 R/stack_jurisdictions.R       | 23 ++--------
 10 files changed, 64 insertions(+), 263 deletions(-)
 delete mode 100644 R/create_small_sero_model.R

diff --git a/R/create_epiwave_timeseries.R b/R/create_epiwave_timeseries.R
index 049aa30..c90f74d 100644
--- a/R/create_epiwave_timeseries.R
+++ b/R/create_epiwave_timeseries.R
@@ -89,9 +89,9 @@ create_epiwave_greta_timeseries <- function(dates,
 
 #' Coerce a date/value table to an epiwave_fixed_timeseries object
 #'
-#' @description Observation data (case counts, hospitalisation counts,
-#'   seroprevalence survey counts/sample sizes) is naturally a table of
-#'   `date`/`value` pairs with real per-date observations, rather than a
+#' @description Observation data (case counts, hospitalisation counts) is
+#'   naturally a table of `date`/`value` pairs with real per-date
+#'   observations, rather than a
 #'   single value replicated across dates -- unlike delay distributions or
 #'   proportions, which usually apply uniformly. This function takes such a
 #'   table (with genuinely partial date coverage, e.g. missing weekends, is
diff --git a/R/create_small_sero_model.R b/R/create_small_sero_model.R
deleted file mode 100644
index 94df7d4..0000000
--- a/R/create_small_sero_model.R
+++ /dev/null
@@ -1,82 +0,0 @@
-#' Create small seroprevalence observation model
-#'
-#' @description This function constructs the observation model and defines
-#'  likelihood over seroprevalence survey data. It mirrors
-#'  `create_observation_model()`: a timeseries of expected observations is
-#'  computed from the infection timeseries, the infection-to-seroconversion
-#'  delay distribution(s), and the proportion of infections to be observed.
-#'  Unlike case/hospitalisation notification (a one-time event best modelled
-#'  as a `discrete_pmf`), seroconversion is typically persistent -- an
-#'  infected person may test positive for many consecutive days -- so
-#'  `delay_from_infection` for a sero stream is usually a `discrete_weights`/
-#'  `discrete_weights_series` object (e.g. probability of testing seropositive
-#'  by day since infection), rather than a normalised `discrete_pmf`. The
-#'  expected observation timeseries, divided by jurisdiction population,
-#'  is treated as the probability of a binomial distribution over the survey
-#'  sample size, from which the observed seropositive count is drawn.
-#'
-#' @param data_id name of observation model data in list
-#' @param observation_model_data list of observation model data lists, as
-#'  produced by `stack_jurisdictions()`
-#' @param infection_model greta arrays that define the infection model
-#'
-#' @importFrom greta %*% as_data binomial sweep
-#'
-#' @return greta arrays of observation model
-#' @export
-create_small_sero_model <- function (data_id = 'sero',
-                                      observation_model_data,
-                                      infection_model
-) {
-
-  observations <- observation_model_data[[data_id]]
-
-  case_mat <- observations$case_mat
-  size_mat <- observations$size_mat
-  prop_mat <- observations$prop_mat
-  convolution_matrices <- observations$convolution_matrices
-  total_pop <- observations$total_pop
-
-  n_jurisdictions <- length(convolution_matrices)
-
-  expected_cases_list <- lapply(
-    1:n_jurisdictions,
-    function(x) {
-      (convolution_matrices[[x]]) %*% infection_model[, x] * prop_mat[, x]
-    })
-
-  expected_cases <- do.call(
-    cbind,
-    expected_cases_list
-  )
-
-  prob <- sweep(expected_cases, 2, total_pop, "/")
-
-  valid_mat <- case_mat
-  valid_mat[is.na(case_mat)] <- FALSE
-  valid_mat[!is.na(case_mat)] <- TRUE
-  valid_idx <- as.logical(as.numeric(valid_mat))
-
-  case_mat_array <- greta::as_data(
-    as.numeric(case_mat)[valid_idx])
-
-  size_mat_array <- greta::as_data(
-    as.numeric(size_mat)[valid_idx])
-
-  greta::distribution(case_mat_array) <- greta::binomial(
-    size_mat_array,
-    prob[valid_idx])
-
-  greta_arrays <- list(
-    prob,
-    convolution_matrices
-  )
-
-  names(greta_arrays) <- c(
-    paste0(data_id, '_prob'),
-    paste0(data_id, '_convolution_matrices')
-  )
-
-  return(greta_arrays)
-
-}
diff --git a/R/define_observation_data.R b/R/define_observation_data.R
index 31d256d..e9e3b31 100644
--- a/R/define_observation_data.R
+++ b/R/define_observation_data.R
@@ -8,9 +8,8 @@
 #'  jurisdiction. A plain data.frame/tibble with `date` and `value` columns
 #'  is fine (gaps in date coverage are fine too) -- it doesn't need to be
 #'  pre-classed, see `as_epiwave_timeseries()`
-#' @param delay_from_infection a `discrete_pmf`/`discrete_weights` object
-#'  (replicated across dates), or an already time-varying
-#'  `discrete_pmf_series`/`discrete_weights_series` object
+#' @param delay_from_infection a `discrete_pmf` object (replicated across
+#'  dates), or an already time-varying `discrete_pmf_series` object
 #' @param proportion_infections proportion data
 #' @param dow_model logical indicating whether to apply a DOW
 #'
@@ -28,44 +27,3 @@ define_observation_data <- function (timeseries_data,
   out
 
 }
-
-#' Define seroprevalence observation data
-#'
-#' @description Bundle one jurisdiction's seroprevalence survey data. Call
-#'  this once per jurisdiction; combine multiple jurisdictions later via
-#'  `stack_jurisdictions()`.
-#'
-#' @param timeseries_data seroprevalence survey timeseries data, for one
-#'  jurisdiction. A plain data.frame/tibble with `date` and `value` columns
-#'  is fine, see `as_epiwave_timeseries()`
-#' @param total_pop population size of this jurisdiction (a single numeric
-#'  value)
-#' @param size_vec sample size of the seroprevalence survey, for one
-#'  jurisdiction. A plain data.frame/tibble with `date` and `value` columns
-#'  is fine, see `as_epiwave_timeseries()`
-#' @param delay_from_infection typically a `discrete_weights`/
-#'  `discrete_weights_series` object (e.g. probability of testing
-#'  seropositive by day since infection, which need not sum to 1 -- unlike
-#'  case/hospitalisation notification, seroconversion is usually persistent
-#'  rather than a one-time event). A `discrete_pmf`/`discrete_pmf_series` is
-#'  also accepted if a normalised delay is more appropriate.
-#' @param proportion_infections proportion data
-#' @param dow_model logical indicating whether to apply a DOW
-#'
-#' @return list of observation data for one data type
-define_sero_data <- function (timeseries_data,
-                              total_pop,
-                              size_vec,
-                              delay_from_infection,
-                              proportion_infections,
-                              dow_model = FALSE) {
-
-  out <- list(timeseries_data = timeseries_data,
-              total_pop = total_pop,
-              size_vec = size_vec,
-              delay_from_infection = delay_from_infection,
-              proportion_infections = proportion_infections,
-              dow_model = dow_model)
-  out
-
-}
diff --git a/R/define_observation_model.R b/R/define_observation_model.R
index 35ad656..043e0e5 100644
--- a/R/define_observation_model.R
+++ b/R/define_observation_model.R
@@ -2,12 +2,12 @@
 #'
 #' @description
 #' Bundle one jurisdiction's observation streams (e.g. cases,
-#' hospitalisations, sero) together. Call this once per jurisdiction;
+#' hospitalisations) together. Call this once per jurisdiction;
 #' combine multiple jurisdictions later via `stack_jurisdictions()`.
 #'
 #' @param target_infection_dates sequence of infection dates
 #' @param ... observation data sets for this jurisdiction (as returned by
-#'  `define_observation_data()`/`define_sero_data()`), named by stream
+#'  `define_observation_data()`), named by stream
 #'
 #' @return list describing one jurisdiction's observation model
 #' @export
diff --git a/R/evaluate.R b/R/evaluate.R
index 89fe315..3abbb21 100644
--- a/R/evaluate.R
+++ b/R/evaluate.R
@@ -1,15 +1,14 @@
-#' Evaluate probability mass or weight for a matrix of delay values
+#' Evaluate probability mass for a matrix of delay values
 #'
-#' Looks up the probability mass or weight values corresponding to a matrix
-#' of delay values from a `discrete_pmf`/`discrete_weights` object,
-#' returning a numeric matrix of the same dimensions. Delay values not
-#' present in the lookup return 0.
+#' Looks up the probability mass values corresponding to a matrix of delay
+#' values from a `discrete_pmf` object, returning a numeric matrix of the
+#' same dimensions. Delay values not present in the lookup return 0.
 #'
-#' @param pmf a `discrete_pmf` or `discrete_weights` object
+#' @param pmf a `discrete_pmf` object
 #' @param day_diff a matrix of integer delay values
 #' @param ... unused
 #'
-#' @return a numeric matrix of probability mass/weight values with the same
+#' @return a numeric matrix of probability mass values with the same
 #'   dimensions as `day_diff`
 #'
 #' @noRd
@@ -27,19 +26,9 @@ evaluate.discrete_pmf <- function(pmf, day_diff, ...) {
   day_diff
 }
 
-#' @noRd
-evaluate.discrete_weights <- function(pmf, day_diff, ...) {
-  pmf <- as.data.frame(pmf)
-  day_diff[] <- pmf$weight[
-    match(unlist(day_diff), pmf$step)
-  ]
-  day_diff[is.na(day_diff)] <- 0
-  day_diff
-}
-
-#' Evaluate probability mass/weight column-wise for a time-varying series
+#' Evaluate probability mass column-wise for a time-varying PMF series
 #'
-#' @param pmf a `discrete_pmf_series` or `discrete_weights_series` object
+#' @param pmf a `discrete_pmf_series` object
 #' @param day_diff a matrix of integer delay values
 #' @param ... unused
 #'
@@ -53,14 +42,3 @@ evaluate.discrete_pmf_series <- function(pmf, day_diff, ...) {
   )
   do.call(cbind, con_list)
 }
-
-#' @noRd
-evaluate.discrete_weights_series <- function(pmf, day_diff, ...) {
-  con_list <- lapply(
-    seq_len(ncol(day_diff)),
-    function(col_idx) {
-      evaluate(pmf$values[[col_idx]], day_diff[, col_idx])
-    }
-  )
-  do.call(cbind, con_list)
-}
diff --git a/R/fit_waves.R b/R/fit_waves.R
index 2c2cc56..ca4f5bb 100644
--- a/R/fit_waves.R
+++ b/R/fit_waves.R
@@ -67,15 +67,10 @@ fit_waves <- function (observations_by_jurisdiction,
 
   # observation model objects in observations
   observation_model_data <- observations$observation_model_data
-  observation_models <- lapply(names(observation_model_data), function(id) {
-    stream <- observation_model_data[[id]]
-    model_fn <- if (!is.null(stream$total_pop)) {
-      create_small_sero_model
-    } else {
-      create_observation_model
-    }
-    model_fn(id, observation_model_data, incidence)
-  })
+  observation_models <- lapply(names(observation_model_data),
+                               create_observation_model,
+                               observation_model_data,
+                               incidence)
   names(observation_models) <- names(observation_model_data)
 
   # greta model fit
diff --git a/R/inits_by_jurisdiction.R b/R/inits_by_jurisdiction.R
index a120b1a..851516e 100644
--- a/R/inits_by_jurisdiction.R
+++ b/R/inits_by_jurisdiction.R
@@ -2,8 +2,8 @@
 #'
 #' @param obs_data numeric vector of observed data, aligned to
 #'   `target_infection_dates` (`NA` for dates without observations)
-#' @param delays a `discrete_pmf_series` or `discrete_weights_series` delay
-#'   object, aligned to `target_infection_dates`
+#' @param delays a `discrete_pmf_series` delay object, aligned to
+#'   `target_infection_dates`
 #' @param obs_prop numeric vector of proportions, aligned to
 #'   `target_infection_dates`, or a greta array of the same
 #' @param target_infection_dates infection date sequence
@@ -39,18 +39,8 @@ inits_by_jurisdiction <- function (obs_data,
   # to shift observation data for calculation of inits
   delays_juris <- delays[case_dates]
 
-  # discrete_weights aren't a proper distribution (don't sum to 1), so
-  # normalise to a discrete_pmf first to get a meaningful mean step
-  mean_step <- function (x) {
-    if (inherits(x, 'discrete_pmf')) {
-      mean(x)
-    } else {
-      mean(epiwave.params::normalise(x))
-    }
-  }
-
   expected_delay_vals <- unlist(lapply(delays_juris$values, function (x)
-    round(mean_step(x))
+    round(mean(x))
   ))
 
   max_delay_vals <- unlist(lapply(delays_juris$values, function (x)
diff --git a/R/new_convolution_matrix.R b/R/new_convolution_matrix.R
index 8d81e92..1c33d3b 100644
--- a/R/new_convolution_matrix.R
+++ b/R/new_convolution_matrix.R
@@ -1,35 +1,29 @@
-#' Construct a forward convolution matrix from a discrete PMF or weights
+#' Construct a forward convolution matrix from a discrete PMF
 #'
 #' @description A common operation in discrete-time epidemic models is forward
 #' convolution through a delay distribution. This can be represented as a
 #' linear operator in matrix form. This function constructs a forward
-#' convolution matrix from a `discrete_pmf`/`discrete_pmf_series` object, or
-#' from a `discrete_weights`/`discrete_weights_series` object (e.g. an
-#' unnormalised persistence/detectability curve, such as the probability of
-#' testing seropositive by day since infection). The resulting matrix can be
-#' multiplied by a time series vector to convolve it forward through a delay
-#' distribution.
+#' convolution matrix from a `discrete_pmf` or `discrete_pmf_series` object.
+#' The resulting matrix can be multiplied by a time series vector to convolve
+#' it forward through a delay distribution.
 #'
-#' When `pmf` is a single `discrete_pmf`/`discrete_weights`, the same object
-#' is applied at every timepoint and `n` must be supplied. When `pmf` is a
-#' `discrete_pmf_series`/`discrete_weights_series`, `n` is derived from the
-#' series index and does not need to be supplied.
+#' When `pmf` is a single `discrete_pmf`, the same PMF is applied at every
+#' timepoint and `n` must be supplied. When `pmf` is a `discrete_pmf_series`,
+#' `n` is derived from the series index and does not need to be supplied.
 #'
 #' @srrstats {PD3.5} Discrete summation is used following Cori et al. (2013,
 #'   AJE) and Gostic et al. (2020, PLOS Comp Bio), who demonstrate its
-#'   suitability for discrete-time epidemiological models. The input PMF or
-#'   weight function is finite by construction, guaranteeing a finite sum
-#'   (PD3.5a).
+#'   suitability for discrete-time epidemiological models. Normalisation of
+#'   the input PMF guarantees a finite sum (PD3.5a).
 #'
-#' @param pmf a `discrete_pmf`, `discrete_weights`, `discrete_pmf_series`, or
-#'   `discrete_weights_series` object. A single `discrete_pmf`/
-#'   `discrete_weights` is applied uniformly across all timepoints. A
-#'   `discrete_pmf_series`/`discrete_weights_series` enables time-varying
-#'   delays, with `n` derived from the series index.
+#' @param pmf a `discrete_pmf` or `discrete_pmf_series` object. A single
+#'   `discrete_pmf` is applied uniformly across all timepoints. A
+#'   `discrete_pmf_series` enables time-varying delays, with `n` derived
+#'   from the series index.
 #' @param n a positive integer giving the number of timepoints, determining
-#'   the dimensions of the output matrix. Required when `pmf` is a single
-#'   `discrete_pmf`/`discrete_weights` object; inferred from the timeseries
-#'   index when `pmf` is a series object.
+#'   the dimensions of the output matrix. Required when `pmf` is a
+#'   `discrete_pmf`; inferred from timeseries index when `pmf` is a
+#'   `discrete_pmf_series`.
 #'
 #' @return an `n` by `n` numeric matrix. To apply the convolution, multiply
 #'   the matrix by a time series vector of length `n`, e.g.
@@ -43,10 +37,9 @@
 #' @export
 new_convolution_matrix <- function(pmf, n = NULL) {
 
-  if (!inherits(pmf, c("discrete_pmf", "discrete_pmf_series",
-                       "discrete_weights", "discrete_weights_series"))) {
+  if (!inherits(pmf, c("discrete_pmf", "discrete_pmf_series"))) {
     cli::cli_abort(
-      "`pmf` must be a {.cls discrete_pmf}, {.cls discrete_pmf_series}, {.cls discrete_weights}, or {.cls discrete_weights_series} object."
+      "`pmf` must be a {.cls discrete_pmf} or {.cls discrete_pmf_series} object."
     )
   }
 
@@ -60,10 +53,10 @@ new_convolution_matrix <- function(pmf, n = NULL) {
 
 #' @noRd
 resolve_n <- function(pmf, n) {
-  if (inherits(pmf, c("discrete_pmf_series", "discrete_weights_series"))) {
+  if (inherits(pmf, "discrete_pmf_series")) {
     if (!is.null(n)) {
       cli::cli_warn(
-        "`n` is ignored when `pmf` is a series object; derived from `pmf$index`."
+        "`n` is ignored when `pmf` is a {.cls discrete_pmf_series}; derived from `pmf$index`."
       )
     }
     # length(pmf$index) is guaranteed valid by new_discrete_series construction
@@ -75,9 +68,7 @@ resolve_n <- function(pmf, n) {
 #' @noRd
 validate_n <- function(n) {
   if (is.null(n)) {
-    cli::cli_abort(
-      "`n` must be supplied when `pmf` is a single {.cls discrete_pmf} or {.cls discrete_weights} object."
-    )
+    cli::cli_abort("`n` must be supplied when `pmf` is a {.cls discrete_pmf}.")
   }
   if (!is.numeric(n) || length(n) != 1) {
     cli::cli_abort("`n` must be a single numeric value.")
diff --git a/R/prepare_observation_data.R b/R/prepare_observation_data.R
index ca0802b..faa6536 100644
--- a/R/prepare_observation_data.R
+++ b/R/prepare_observation_data.R
@@ -5,37 +5,35 @@
 #'  jurisdictions are combined later, by `stack_jurisdictions()`.
 #'
 #' @param observation_data data for one jurisdiction, one stream, as
-#'  returned by `define_observation_data()`/`define_sero_data()`.
-#'  `timeseries_data` (and `size_vec`, for sero streams) may be a plain
-#'  data.frame/tibble with `date` and `value` columns -- it doesn't need to
-#'  be pre-classed as `epiwave_timeseries`, see `as_epiwave_timeseries()`.
-#'  `delay_from_infection` may be a single `discrete_pmf`/`discrete_weights`
-#'  object (replicated across `target_infection_dates`), or an already
-#'  time-varying `discrete_pmf_series`/`discrete_weights_series` (aligned to
-#'  `target_infection_dates`)
+#'  returned by `define_observation_data()`.
+#'  `timeseries_data` may be a plain data.frame/tibble with `date` and
+#'  `value` columns -- it doesn't need to be pre-classed as
+#'  `epiwave_timeseries`, see `as_epiwave_timeseries()`.
+#'  `delay_from_infection` may be a single `discrete_pmf` object (replicated
+#'  across `target_infection_dates`), or an already time-varying
+#'  `discrete_pmf_series` (aligned to `target_infection_dates`)
 #' @param target_infection_dates full date sequence of infection timeseries,
 #'  shared across all jurisdictions in a fit
 #'
 #' @return a list describing this jurisdiction's stream: `convolution_matrix`,
-#'  `case_vec`, `prop_vec`, `dow_model`, `inits_values`, `observable_idx`, and
-#'  (for seroprevalence streams) `total_pop`/`size_vec`
+#'  `case_vec`, `prop_vec`, `dow_model`, `inits_values`, `observable_idx`
 #' @export
 prepare_observation_data <- function (observation_data,
                                       target_infection_dates) {
 
   delays <- observation_data$delay_from_infection
-  if (inherits(delays, c('discrete_pmf_series', 'discrete_weights_series'))) {
+  if (inherits(delays, 'discrete_pmf_series')) {
     delays <- delays[as.Date(target_infection_dates)]
     if (!identical(as.Date(delays$index), as.Date(target_infection_dates))) {
       stop('`delay_from_infection` dates must match `target_infection_dates`')
     }
-  } else if (inherits(delays, c('discrete_pmf', 'discrete_weights'))) {
+  } else if (inherits(delays, 'discrete_pmf')) {
     delays <- epiwave.params::new_discrete_series(
       values = delays,
       index = target_infection_dates)
   } else {
-    stop('`delay_from_infection` must be a discrete_pmf, discrete_weights, ',
-         'discrete_pmf_series, or discrete_weights_series object')
+    stop('`delay_from_infection` must be a discrete_pmf or ',
+         'discrete_pmf_series object')
   }
 
   prop <- observation_data$proportion_infections
@@ -74,20 +72,10 @@ prepare_observation_data <- function (observation_data,
   inits_values[inits$observable_idx] <- inits$inits_values
   observable_idx[inits$observable_idx] <- TRUE
 
-  out <- list(convolution_matrix = convolution_matrix,
-              case_vec = case_vec,
-              prop_vec = prop_vec,
-              dow_model = observation_data$dow_model,
-              inits_values = inits_values,
-              observable_idx = observable_idx)
-
-  if ('total_pop' %in% names(observation_data)) {
-    out$total_pop <- observation_data$total_pop
-  }
-  if ('size_vec' %in% names(observation_data)) {
-    size_vec <- as_epiwave_timeseries(observation_data$size_vec)
-    out$size_vec <- as_matrix(size_vec, target_infection_dates)
-  }
-
-  out
+  list(convolution_matrix = convolution_matrix,
+       case_vec = case_vec,
+       prop_vec = prop_vec,
+       dow_model = observation_data$dow_model,
+       inits_values = inits_values,
+       observable_idx = observable_idx)
 }
diff --git a/R/stack_jurisdictions.R b/R/stack_jurisdictions.R
index b5ba4a8..a5cf252 100644
--- a/R/stack_jurisdictions.R
+++ b/R/stack_jurisdictions.R
@@ -30,7 +30,8 @@ stack_jurisdictions <- function (observations_by_jurisdiction) {
   target_infection_dates <- observations_by_jurisdiction[[1]]$target_infection_dates
   dates_match <- vapply(
     observations_by_jurisdiction,
-    function(x) identical(as.Date(x$target_infection_dates), as.Date(target_infection_dates)),
+    function(x) identical(as.Date(x$target_infection_dates),
+                          as.Date(target_infection_dates)),
     logical(1))
   if (!all(dates_match)) {
     stop('All jurisdictions must share the same target_infection_dates')
@@ -112,27 +113,9 @@ stack_stream <- function (stream_id,
 
   convolution_matrices <- lapply(per_jurisdiction, `[[`, "convolution_matrix")
 
-  sero_flags <- vapply(per_jurisdiction, function(x) !is.null(x$total_pop), logical(1))
-  if (length(unique(sero_flags)) > 1) {
-    stop(sprintf(
-      'Jurisdictions disagree on whether "%s" is a seroprevalence stream ',
-      '(define_sero_data() vs define_observation_data())',
-      stream_id))
-  }
-
-  out <- list(
+  list(
     convolution_matrices = convolution_matrices,
     case_mat = case_mat,
     prop_mat = prop_mat
   )
-
-  if (isTRUE(sero_flags[1])) {
-    out$total_pop <- vapply(per_jurisdiction, `[[`, numeric(1), "total_pop")
-    names(out$total_pop) <- jurisdictions
-    out$size_mat <- do.call(cbind, lapply(per_jurisdiction, `[[`, "size_vec"))
-    colnames(out$size_mat) <- jurisdictions
-    rownames(out$size_mat) <- as.character(target_infection_dates)
-  }
-
-  out
 }

From 78f6128f8645533a5f012b87bbe3ff0cee5c3ab9 Mon Sep 17 00:00:00 2001
From: Saras Windecker <smwindecker@users.noreply.github.com>
Date: Sun, 5 Jul 2026 23:28:26 +0300
Subject: [PATCH 07/14] Remove testing.R's dead sero block, regenerate docs

testing.R's sero fit block called define_sero_data(), which no longer
exists after the previous commit -- it was already non-functional (sero_dat/
sero_size_mat/sero_conversion were never assigned), so this just removes
dead code rather than breaking anything that worked. Left short pointer
comments to dev/sero-integration-notes.md at both spots (data prep, and
where a third define_observation_data() block would go) for whoever picks
this back up.

devtools::document() drops create_small_sero_model/define_sero_data exports
and man pages, and the now-unused importFrom(greta, binomial).

Verified after all sero-removal changes: devtools::check() is 0 errors, 0
warnings (same 3 pre-existing NOTEs); testing.R's basic fit still runs
against real data in ../data; both single_jurisdiction_workflow.R and
multi_jurisdiction_workflow.R still run to completion with their
stopifnot() checks passing.
---
 NAMESPACE                       |  2 --
 man/as_epiwave_timeseries.Rd    |  6 ++--
 man/create_small_sero_model.Rd  | 39 ---------------------
 man/define_observation_data.Rd  |  5 ++-
 man/define_observation_model.Rd |  4 +--
 man/define_sero_data.Rd         | 46 -------------------------
 man/inits_by_jurisdiction.Rd    |  4 +--
 man/new_convolution_matrix.Rd   | 33 ++++++++----------
 man/prepare_observation_data.Rd | 18 +++++-----
 tests/test_workflow/testing.R   | 61 +++------------------------------
 10 files changed, 35 insertions(+), 183 deletions(-)
 delete mode 100644 man/create_small_sero_model.Rd
 delete mode 100644 man/define_sero_data.Rd

diff --git a/NAMESPACE b/NAMESPACE
index 2fec4de..f322b04 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -12,7 +12,6 @@ export(create_epiwave_greta_timeseries)
 export(create_epiwave_timeseries)
 export(create_infection_timeseries)
 export(create_observation_model)
-export(create_small_sero_model)
 export(define_observation_data)
 export(define_observation_model)
 export(estimate_reff)
@@ -36,7 +35,6 @@ importFrom(ggplot2,theme)
 importFrom(greta,"%*%")
 importFrom(greta,apply)
 importFrom(greta,as_data)
-importFrom(greta,binomial)
 importFrom(greta,dirichlet)
 importFrom(greta,lognormal)
 importFrom(greta,negative_binomial)
diff --git a/man/as_epiwave_timeseries.Rd b/man/as_epiwave_timeseries.Rd
index 6f43b15..c99869a 100644
--- a/man/as_epiwave_timeseries.Rd
+++ b/man/as_epiwave_timeseries.Rd
@@ -14,9 +14,9 @@ data.frame/tibble with \code{date} and \code{value} columns}
 an \code{epiwave_timeseries} or numeric object, ready for \code{as_matrix()}
 }
 \description{
-Observation data (case counts, hospitalisation counts,
-seroprevalence survey counts/sample sizes) is naturally a table of
-\code{date}/\code{value} pairs with real per-date observations, rather than a
+Observation data (case counts, hospitalisation counts) is
+naturally a table of \code{date}/\code{value} pairs with real per-date
+observations, rather than a
 single value replicated across dates -- unlike delay distributions or
 proportions, which usually apply uniformly. This function takes such a
 table (with genuinely partial date coverage, e.g. missing weekends, is
diff --git a/man/create_small_sero_model.Rd b/man/create_small_sero_model.Rd
deleted file mode 100644
index f7dc57e..0000000
--- a/man/create_small_sero_model.Rd
+++ /dev/null
@@ -1,39 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/create_small_sero_model.R
-\name{create_small_sero_model}
-\alias{create_small_sero_model}
-\title{Create small seroprevalence observation model}
-\usage{
-create_small_sero_model(
-  data_id = "sero",
-  observation_model_data,
-  infection_model
-)
-}
-\arguments{
-\item{data_id}{name of observation model data in list}
-
-\item{observation_model_data}{list of observation model data lists, as
-produced by \code{stack_jurisdictions()}}
-
-\item{infection_model}{greta arrays that define the infection model}
-}
-\value{
-greta arrays of observation model
-}
-\description{
-This function constructs the observation model and defines
-likelihood over seroprevalence survey data. It mirrors
-\code{create_observation_model()}: a timeseries of expected observations is
-computed from the infection timeseries, the infection-to-seroconversion
-delay distribution(s), and the proportion of infections to be observed.
-Unlike case/hospitalisation notification (a one-time event best modelled
-as a \code{discrete_pmf}), seroconversion is typically persistent -- an
-infected person may test positive for many consecutive days -- so
-\code{delay_from_infection} for a sero stream is usually a \code{discrete_weights}/
-\code{discrete_weights_series} object (e.g. probability of testing seropositive
-by day since infection), rather than a normalised \code{discrete_pmf}. The
-expected observation timeseries, divided by jurisdiction population,
-is treated as the probability of a binomial distribution over the survey
-sample size, from which the observed seropositive count is drawn.
-}
diff --git a/man/define_observation_data.Rd b/man/define_observation_data.Rd
index e35d4e6..26dd540 100644
--- a/man/define_observation_data.Rd
+++ b/man/define_observation_data.Rd
@@ -17,9 +17,8 @@ jurisdiction. A plain data.frame/tibble with \code{date} and \code{value} column
 is fine (gaps in date coverage are fine too) -- it doesn't need to be
 pre-classed, see \code{as_epiwave_timeseries()}}
 
-\item{delay_from_infection}{a \code{discrete_pmf}/\code{discrete_weights} object
-(replicated across dates), or an already time-varying
-\code{discrete_pmf_series}/\code{discrete_weights_series} object}
+\item{delay_from_infection}{a \code{discrete_pmf} object (replicated across
+dates), or an already time-varying \code{discrete_pmf_series} object}
 
 \item{proportion_infections}{proportion data}
 
diff --git a/man/define_observation_model.Rd b/man/define_observation_model.Rd
index 7260faf..b58e03a 100644
--- a/man/define_observation_model.Rd
+++ b/man/define_observation_model.Rd
@@ -10,13 +10,13 @@ define_observation_model(target_infection_dates = NULL, ...)
 \item{target_infection_dates}{sequence of infection dates}
 
 \item{...}{observation data sets for this jurisdiction (as returned by
-\code{define_observation_data()}/\code{define_sero_data()}), named by stream}
+\code{define_observation_data()}), named by stream}
 }
 \value{
 list describing one jurisdiction's observation model
 }
 \description{
 Bundle one jurisdiction's observation streams (e.g. cases,
-hospitalisations, sero) together. Call this once per jurisdiction;
+hospitalisations) together. Call this once per jurisdiction;
 combine multiple jurisdictions later via \code{stack_jurisdictions()}.
 }
diff --git a/man/define_sero_data.Rd b/man/define_sero_data.Rd
deleted file mode 100644
index 363fd05..0000000
--- a/man/define_sero_data.Rd
+++ /dev/null
@@ -1,46 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/define_observation_data.R
-\name{define_sero_data}
-\alias{define_sero_data}
-\title{Define seroprevalence observation data}
-\usage{
-define_sero_data(
-  timeseries_data,
-  total_pop,
-  size_vec,
-  delay_from_infection,
-  proportion_infections,
-  dow_model = FALSE
-)
-}
-\arguments{
-\item{timeseries_data}{seroprevalence survey timeseries data, for one
-jurisdiction. A plain data.frame/tibble with \code{date} and \code{value} columns
-is fine, see \code{as_epiwave_timeseries()}}
-
-\item{total_pop}{population size of this jurisdiction (a single numeric
-value)}
-
-\item{size_vec}{sample size of the seroprevalence survey, for one
-jurisdiction. A plain data.frame/tibble with \code{date} and \code{value} columns
-is fine, see \code{as_epiwave_timeseries()}}
-
-\item{delay_from_infection}{typically a \code{discrete_weights}/
-\code{discrete_weights_series} object (e.g. probability of testing
-seropositive by day since infection, which need not sum to 1 -- unlike
-case/hospitalisation notification, seroconversion is usually persistent
-rather than a one-time event). A \code{discrete_pmf}/\code{discrete_pmf_series} is
-also accepted if a normalised delay is more appropriate.}
-
-\item{proportion_infections}{proportion data}
-
-\item{dow_model}{logical indicating whether to apply a DOW}
-}
-\value{
-list of observation data for one data type
-}
-\description{
-Bundle one jurisdiction's seroprevalence survey data. Call
-this once per jurisdiction; combine multiple jurisdictions later via
-\code{stack_jurisdictions()}.
-}
diff --git a/man/inits_by_jurisdiction.Rd b/man/inits_by_jurisdiction.Rd
index a7f5b37..eb5f957 100644
--- a/man/inits_by_jurisdiction.Rd
+++ b/man/inits_by_jurisdiction.Rd
@@ -10,8 +10,8 @@ inits_by_jurisdiction(obs_data, delays, obs_prop, target_infection_dates)
 \item{obs_data}{numeric vector of observed data, aligned to
 \code{target_infection_dates} (\code{NA} for dates without observations)}
 
-\item{delays}{a \code{discrete_pmf_series} or \code{discrete_weights_series} delay
-object, aligned to \code{target_infection_dates}}
+\item{delays}{a \code{discrete_pmf_series} delay object, aligned to
+\code{target_infection_dates}}
 
 \item{obs_prop}{numeric vector of proportions, aligned to
 \code{target_infection_dates}, or a greta array of the same}
diff --git a/man/new_convolution_matrix.Rd b/man/new_convolution_matrix.Rd
index 31dbaf4..bfd58c7 100644
--- a/man/new_convolution_matrix.Rd
+++ b/man/new_convolution_matrix.Rd
@@ -2,21 +2,20 @@
 % Please edit documentation in R/new_convolution_matrix.R
 \name{new_convolution_matrix}
 \alias{new_convolution_matrix}
-\title{Construct a forward convolution matrix from a discrete PMF or weights}
+\title{Construct a forward convolution matrix from a discrete PMF}
 \usage{
 new_convolution_matrix(pmf, n = NULL)
 }
 \arguments{
-\item{pmf}{a \code{discrete_pmf}, \code{discrete_weights}, \code{discrete_pmf_series}, or
-\code{discrete_weights_series} object. A single \code{discrete_pmf}/
-\code{discrete_weights} is applied uniformly across all timepoints. A
-\code{discrete_pmf_series}/\code{discrete_weights_series} enables time-varying
-delays, with \code{n} derived from the series index.}
+\item{pmf}{a \code{discrete_pmf} or \code{discrete_pmf_series} object. A single
+\code{discrete_pmf} is applied uniformly across all timepoints. A
+\code{discrete_pmf_series} enables time-varying delays, with \code{n} derived
+from the series index.}
 
 \item{n}{a positive integer giving the number of timepoints, determining
-the dimensions of the output matrix. Required when \code{pmf} is a single
-\code{discrete_pmf}/\code{discrete_weights} object; inferred from the timeseries
-index when \code{pmf} is a series object.}
+the dimensions of the output matrix. Required when \code{pmf} is a
+\code{discrete_pmf}; inferred from timeseries index when \code{pmf} is a
+\code{discrete_pmf_series}.}
 }
 \value{
 an \code{n} by \code{n} numeric matrix. To apply the convolution, multiply
@@ -27,17 +26,13 @@ the matrix by a time series vector of length \code{n}, e.g.
 A common operation in discrete-time epidemic models is forward
 convolution through a delay distribution. This can be represented as a
 linear operator in matrix form. This function constructs a forward
-convolution matrix from a \code{discrete_pmf}/\code{discrete_pmf_series} object, or
-from a \code{discrete_weights}/\code{discrete_weights_series} object (e.g. an
-unnormalised persistence/detectability curve, such as the probability of
-testing seropositive by day since infection). The resulting matrix can be
-multiplied by a time series vector to convolve it forward through a delay
-distribution.
+convolution matrix from a \code{discrete_pmf} or \code{discrete_pmf_series} object.
+The resulting matrix can be multiplied by a time series vector to convolve
+it forward through a delay distribution.
 
-When \code{pmf} is a single \code{discrete_pmf}/\code{discrete_weights}, the same object
-is applied at every timepoint and \code{n} must be supplied. When \code{pmf} is a
-\code{discrete_pmf_series}/\code{discrete_weights_series}, \code{n} is derived from the
-series index and does not need to be supplied.
+When \code{pmf} is a single \code{discrete_pmf}, the same PMF is applied at every
+timepoint and \code{n} must be supplied. When \code{pmf} is a \code{discrete_pmf_series},
+\code{n} is derived from the series index and does not need to be supplied.
 }
 \examples{
 pmf <- epiwave.params::as_discrete_pmf(
diff --git a/man/prepare_observation_data.Rd b/man/prepare_observation_data.Rd
index 5840ec2..8cbb0d7 100644
--- a/man/prepare_observation_data.Rd
+++ b/man/prepare_observation_data.Rd
@@ -8,22 +8,20 @@ prepare_observation_data(observation_data, target_infection_dates)
 }
 \arguments{
 \item{observation_data}{data for one jurisdiction, one stream, as
-returned by \code{define_observation_data()}/\code{define_sero_data()}.
-\code{timeseries_data} (and \code{size_vec}, for sero streams) may be a plain
-data.frame/tibble with \code{date} and \code{value} columns -- it doesn't need to
-be pre-classed as \code{epiwave_timeseries}, see \code{as_epiwave_timeseries()}.
-\code{delay_from_infection} may be a single \code{discrete_pmf}/\code{discrete_weights}
-object (replicated across \code{target_infection_dates}), or an already
-time-varying \code{discrete_pmf_series}/\code{discrete_weights_series} (aligned to
-\code{target_infection_dates})}
+returned by \code{define_observation_data()}.
+\code{timeseries_data} may be a plain data.frame/tibble with \code{date} and
+\code{value} columns -- it doesn't need to be pre-classed as
+\code{epiwave_timeseries}, see \code{as_epiwave_timeseries()}.
+\code{delay_from_infection} may be a single \code{discrete_pmf} object (replicated
+across \code{target_infection_dates}), or an already time-varying
+\code{discrete_pmf_series} (aligned to \code{target_infection_dates})}
 
 \item{target_infection_dates}{full date sequence of infection timeseries,
 shared across all jurisdictions in a fit}
 }
 \value{
 a list describing this jurisdiction's stream: \code{convolution_matrix},
-\code{case_vec}, \code{prop_vec}, \code{dow_model}, \code{inits_values}, \code{observable_idx}, and
-(for seroprevalence streams) \code{total_pop}/\code{size_vec}
+\code{case_vec}, \code{prop_vec}, \code{dow_model}, \code{inits_values}, \code{observable_idx}
 }
 \description{
 Prepare the data objects needed for one jurisdiction's
diff --git a/tests/test_workflow/testing.R b/tests/test_workflow/testing.R
index 2ab731c..c7870f1 100644
--- a/tests/test_workflow/testing.R
+++ b/tests/test_workflow/testing.R
@@ -1,10 +1,9 @@
 ## header
-#library(epiwave)
 devtools::load_all()
 library(epiwave.params)
 library(dplyr)
 library(greta)
-# devtools::load_all()
+
 ## user modified
 infection_days <- seq(from = as.Date('2021-04-01'),
                  to = as.Date('2022-01-01'), 'days')
@@ -13,30 +12,11 @@ study_seq <- seq(from = as.Date('2021-06-01'),
 jurisdictions <- 'VIC' #c('NSW', 'VIC')
 
 # specific folder with not synced data
-# not_synced_folder <- '~/not_synced/data/' # AH folder
 not_synced_folder <- '../data'
 
-
 ## data prep
 
-## sero simulated data
-# sero_file <- readRDS(file = paste0(not_synced_folder, "sero_data_test.RDS"))
-#
-# sero_dat <- sero_file |>
-#   select(jurisdiction,value,date)
-#
-# class(sero_dat) <- c('epiwave_fixed_timeseries',
-#                      'epiwave_timeseries',
-#                      class(sero_dat))
-#
-# size_dummy <- sero_file |>
-#   select(jurisdiction,size,date) |>
-#   rename(value = size)
-#
-# class(size_dummy) <- c('epiwave_fixed_timeseries',
-#                      'epiwave_timeseries',
-#                      class(size_dummy))
-# sero_size_mat <- as_matrix(size_dummy)
+# seroprevalence support is deferred for now -- see dev/sero-integration-notes.md
 
 # notification counts
 notif_file <- paste0(not_synced_folder, '/COVID_live_cases.rds')
@@ -47,8 +27,6 @@ notif_dat <- notif_file |>
 if (exists('jurisdictions')) {
   notif_dat <- notif_dat[notif_dat$jurisdiction %in% jurisdictions,]
 }
-# no manual class(notif_dat) <- c(...) needed -- define_observation_data()
-# accepts a plain date/value data.frame directly (see as_epiwave_timeseries())
 
 # # hospitalisation counts
 hosp_file <- paste0(not_synced_folder, '/COVID_live_cases_in_hospital.rds')
@@ -84,11 +62,6 @@ incubation <- readRDS('tests/test_distributions/incubation.rds')
 gi <- readRDS('tests/test_distributions/gi.rds')
 onset_to_notification <- readRDS('tests/test_distributions/onset_to_notification.rds')
 
-# sero_conversion <- readRDS(paste0(not_synced_folder, "sero_curve_test.RDS"))
-# sero_conversion <- as_discrete_weights(
-#   x = sero_conversion$delay,
-#   y = sero_conversion$mass)
-
 hosp_dist <- distributional::dist_weibull(shape = 2.51, scale = 10.17)
 hosp_delay_ecdf <- as_discrete_pmf(hosp_dist)
 
@@ -121,34 +94,8 @@ fit_object <- fit_waves(
   infection_model_type = 'gp_growth_rate'# 'flat_prior'#,define_infection_model() #
 )
 
-jurisdiction_with_sero_observation_models <- define_observation_model(
-
-  target_infection_dates = infection_days,
-
-  cases = define_observation_data(
-    timeseries_data = notif_dat,
-    delay_from_infection = add_discrete(
-      incubation,
-      onset_to_notification),
-    proportion_infections = car,
-    dow_model = TRUE),
-  sero = define_sero_data(
-    timeseries_data = sero_dat,
-    total_pop = 8e6,
-    size_vec = sero_size_mat,
-    delay_from_infection = sero_conversion,
-    proportion_infections = 1
-  )
-)
-
-with_sero_observations_by_jurisdiction <- setNames(
-  list(jurisdiction_with_sero_observation_models),
-  jurisdictions)
-
-fit_object <- fit_waves(
-  observations_by_jurisdiction = with_sero_observations_by_jurisdiction,
-  infection_model_type = 'gp_growth_rate'# 'flat_prior'#,define_infection_model() #
-)
+# a sero stream would go here as a third define_observation_data() block once
+# that pathway is re-integrated -- see dev/sero-integration-notes.md
 
 rhats <- coda::gelman.diag(fit_object$fit, autoburnin = FALSE, multivariate = FALSE)
 max(rhats$psrf[, 1], na.rm = TRUE)

From 303159be3a15ee86c8d99a3d5d4c0a7962296f22 Mon Sep 17 00:00:00 2001
From: Saras Windecker <smwindecker@users.noreply.github.com>
Date: Sun, 5 Jul 2026 23:44:47 +0300
Subject: [PATCH 08/14] Simplify jurisdiction stacking ergonomics, defer
 flat_prior-only inits

Two changes, implemented together since they touch the same functions:

1. stack_jurisdictions() becomes the explicit, user-facing combining step,
   taking named ... args (mirroring define_observation_model()'s own
   cases=/hospitalisations= pattern) instead of a pre-built named list --
   e.g. stack_jurisdictions(VIC = obs_vic, NSW = obs_nsw). A single
   jurisdiction needs no combining step at all: define_observation_model()'s
   output (now stamped with class epiwave_observation_model) goes straight
   to fit_waves(), which detects whether it already received a stacked
   object (class epiwave_stacked_observations) or a raw single-jurisdiction
   one and wraps the latter internally (auto-labelled, since there's no
   jurisdiction identity to preserve for n=1). This removes the
   setNames(list(...), jurisdictions) step that was previously required even
   for a single jurisdiction.

2. GAM-based initial values (inits_by_jurisdiction(), and the cross-stream
   union/mean that combines them) are only ever consumed by
   infection_model_type = 'flat_prior' -- every GP-based infection model
   ignores them entirely (confirmed in create_infection_timeseries():
   observable_infection is referenced only inside the flat_prior branch).
   Previously this ran unconditionally for every stream x jurisdiction
   during data prep, running an mgcv::gam() fit that was silently discarded
   whenever a GP model was used instead. prepare_observation_data() no
   longer computes inits at all (keeps delays in its output instead, needed
   later); define_observation_model()/stack_jurisdictions() no longer
   eagerly combine them either. A new compute_flat_prior_inits() does this
   work on the fully-stacked object, called lazily by fit_waves() only
   inside its flat_prior branch, before create_infection_timeseries() (which
   needs the result) and reused for the MCMC initial values.

Verified: a standalone check exercising both a single-jurisdiction fit
(flat_prior AND gp_growth_rate, confirming the GP path never triggers inits
computation) and a two-jurisdiction stack_jurisdictions() combination (plus
the unnamed-args error path) all pass. devtools::check() is 0 errors, 0
warnings. testing.R's real-data fit (../data) verified against both
flat_prior and gp_growth_rate. Both synthetic workflow scripts updated to
the new API and re-verified end to end.
---
 R/compute_flat_prior_inits.R                  | 66 +++++++++++++++++++
 R/define_observation_model.R                  | 25 +++----
 R/fit_waves.R                                 | 27 +++++---
 R/prepare_observation_data.R                  | 21 ++----
 R/stack_jurisdictions.R                       | 61 +++++++++--------
 man/define_observation_model.Rd               |  9 ++-
 man/fit_waves.Rd                              | 11 ++--
 man/prepare_observation_data.Rd               |  4 +-
 man/stack_jurisdictions.Rd                    | 25 +++----
 .../multi_jurisdiction_workflow.R             | 17 +++--
 .../single_jurisdiction_workflow.R            | 13 ++--
 tests/test_workflow/test2.R                   |  9 +--
 tests/test_workflow/testing.R                 |  9 +--
 13 files changed, 183 insertions(+), 114 deletions(-)
 create mode 100644 R/compute_flat_prior_inits.R

diff --git a/R/compute_flat_prior_inits.R b/R/compute_flat_prior_inits.R
new file mode 100644
index 0000000..a44bfd2
--- /dev/null
+++ b/R/compute_flat_prior_inits.R
@@ -0,0 +1,66 @@
+#' Compute initial values for the flat_prior infection model
+#'
+#' @description GAM-smoothed initial values and the observable-date window
+#'  are only used by `infection_model_type = 'flat_prior'` -- GP-based
+#'  infection models (`gp_infections`/`gp_growth_rate`/
+#'  `gp_growth_rate_deriv`) never reference `observable_idx`/`inits_values`.
+#'  This is deliberately not computed during data preparation
+#'  (`prepare_observation_data()`/`define_observation_model()`/
+#'  `stack_jurisdictions()`), since doing so unconditionally would run an
+#'  `mgcv::gam()` fit per stream per jurisdiction for no purpose whenever a
+#'  GP model is used instead. `fit_waves()` calls this lazily, only inside
+#'  its `flat_prior` branch, on the fully stacked multi-jurisdiction object.
+#'
+#' @param observations a stacked observations object, as produced by
+#'  `stack_jurisdictions()`
+#'
+#' @return list with `incidence_observable` (`date x jurisdiction` logical
+#'  matrix, `TRUE` where at least one stream has nearby information to infer
+#'  an infection) and `incidence_observable_inits` (`date x jurisdiction`
+#'  numeric matrix, mean GAM-smoothed initial value across streams)
+#' @noRd
+compute_flat_prior_inits <- function (observations) {
+
+  target_infection_dates <- observations$target_infection_dates
+  jurisdictions <- observations$target_jurisdictions
+  n_dates <- length(target_infection_dates)
+  stream_ids <- names(observations$observation_model_data)
+
+  per_jurisdiction <- lapply(seq_along(jurisdictions), function (j) {
+
+    stream_inits <- lapply(stream_ids, function (stream_id) {
+      stream <- observations$observation_model_data[[stream_id]]
+      inits_by_jurisdiction(
+        stream$case_mat[, j],
+        stream$delays[[j]],
+        stream$prop_mat[, j],
+        target_infection_dates)
+    })
+
+    observable_idx <- rep(FALSE, n_dates)
+    inits_stream_mat <- matrix(NA_real_, n_dates, length(stream_ids))
+
+    for (s in seq_along(stream_inits)) {
+      idx <- stream_inits[[s]]$observable_idx
+      observable_idx[idx] <- TRUE
+      inits_stream_mat[idx, s] <- stream_inits[[s]]$inits_values
+    }
+
+    list(
+      observable_idx = observable_idx,
+      inits_values = rowMeans(inits_stream_mat, na.rm = TRUE)
+    )
+  })
+
+  incidence_observable <- do.call(
+    cbind, lapply(per_jurisdiction, `[[`, "observable_idx"))
+  incidence_observable_inits <- do.call(
+    cbind, lapply(per_jurisdiction, `[[`, "inits_values"))
+  colnames(incidence_observable) <- jurisdictions
+  colnames(incidence_observable_inits) <- jurisdictions
+
+  list(
+    incidence_observable = incidence_observable,
+    incidence_observable_inits = incidence_observable_inits
+  )
+}
diff --git a/R/define_observation_model.R b/R/define_observation_model.R
index 043e0e5..c4d8299 100644
--- a/R/define_observation_model.R
+++ b/R/define_observation_model.R
@@ -2,14 +2,17 @@
 #'
 #' @description
 #' Bundle one jurisdiction's observation streams (e.g. cases,
-#' hospitalisations) together. Call this once per jurisdiction;
-#' combine multiple jurisdictions later via `stack_jurisdictions()`.
+#' hospitalisations) together. Call this once per jurisdiction. A
+#' single-jurisdiction fit can pass this straight to `fit_waves()`; to
+#' combine multiple jurisdictions, pass several of these to
+#' `stack_jurisdictions()` first.
 #'
 #' @param target_infection_dates sequence of infection dates
 #' @param ... observation data sets for this jurisdiction (as returned by
 #'  `define_observation_data()`), named by stream
 #'
-#' @return list describing one jurisdiction's observation model
+#' @return list describing one jurisdiction's observation model, with class
+#'  `epiwave_observation_model`
 #' @export
 #'
 define_observation_model <- function (target_infection_dates = NULL, ...) {
@@ -20,18 +23,10 @@ define_observation_model <- function (target_infection_dates = NULL, ...) {
                                             prepare_observation_data,
                                             target_infection_dates)
 
-  incidence_observable <- Reduce(
-    `|`,
-    lapply(prepared_observation_model_data, function(x) x$observable_idx))
+  out <- list(observation_model_data = prepared_observation_model_data,
+             target_infection_dates = target_infection_dates)
 
-  inits_stream_mat <- do.call(
-    cbind,
-    lapply(prepared_observation_model_data, function(x) x$inits_values))
-  incidence_observable_inits <- rowMeans(inits_stream_mat, na.rm = TRUE)
-
-  list(observation_model_data = prepared_observation_model_data,
-       target_infection_dates = target_infection_dates,
-       incidence_observable_inits = incidence_observable_inits,
-       incidence_observable = incidence_observable)
+  class(out) <- c("epiwave_observation_model", class(out))
+  out
 
 }
diff --git a/R/fit_waves.R b/R/fit_waves.R
index ca4f5bb..b33f48d 100644
--- a/R/fit_waves.R
+++ b/R/fit_waves.R
@@ -17,11 +17,10 @@
 #'  gp_growth_rate_derivative: in now-cast/forecast the infection trajectory
 #'  would follow the most recent growth rate trend.
 #'
-#' @param observations_by_jurisdiction a named list of per-jurisdiction
-#'   observation model bundles (output of `define_observation_model()`),
-#'   named by jurisdiction. A single-jurisdiction fit is a length-1 list.
-#'   Combined internally via `stack_jurisdictions()`; jurisdictions sharing
-#'   one fit are partially pooled (shared GP kernel hyperparameters, and
+#' @param observations either a single jurisdiction's observation model
+#'   (output of `define_observation_model()`), or several jurisdictions
+#'   already combined via `stack_jurisdictions()`. Jurisdictions sharing one
+#'   fit are partially pooled (shared GP kernel hyperparameters, and
 #'   shared/hierarchical DOW priors where requested).
 #' @param infection_model_type options include 'flat_prior', 'infections',
 #'   'growth_rate', 'growth_rate_derivative'. See description for more info.
@@ -34,7 +33,7 @@
 #' @return list of infection_model, fit, infection_days, and jurisdictions
 #' @export
 #'
-fit_waves <- function (observations_by_jurisdiction,
+fit_waves <- function (observations,
                        infection_model_type = c('flat_prior',
                                                 'gp_infections',
                                                 'gp_growth_rate',
@@ -45,7 +44,9 @@ fit_waves <- function (observations_by_jurisdiction,
                        n_samples = 2000,
                        n_extra_samples = 1000) {
 
-  observations <- stack_jurisdictions(observations_by_jurisdiction)
+  if (!inherits(observations, 'epiwave_stacked_observations')) {
+    observations <- stack_jurisdictions_list(list(observations))
+  }
 
   # prep the model objects
   target_infection_dates <- observations$target_infection_dates # Ihat
@@ -54,7 +55,15 @@ fit_waves <- function (observations_by_jurisdiction,
   jurisdictions <- observations$target_jurisdictions
   n_jurisdictions <- length(jurisdictions)
 
-  observable_infection <- observations$incidence_observable
+  # GAM-based inits/observable window are only meaningful for flat_prior --
+  # GP-based infection models never reference them -- so this is computed
+  # lazily, only when actually needed (see compute_flat_prior_inits()).
+  if (infection_model_type == 'flat_prior') {
+    flat_prior_inits <- compute_flat_prior_inits(observations)
+    observable_infection <- flat_prior_inits$incidence_observable
+  } else {
+    observable_infection <- NULL
+  }
 
   # set up the greta model
   # infection timeseries model
@@ -78,7 +87,7 @@ fit_waves <- function (observations_by_jurisdiction,
 
   if (infection_model_type == 'flat_prior') {
 
-    incidence_observable_inits <- observations$incidence_observable_inits
+    incidence_observable_inits <- flat_prior_inits$incidence_observable_inits
 
     indexed_incidence_observable_inits <-
       incidence_observable_inits[observable_infection]
diff --git a/R/prepare_observation_data.R b/R/prepare_observation_data.R
index faa6536..5ee7384 100644
--- a/R/prepare_observation_data.R
+++ b/R/prepare_observation_data.R
@@ -16,7 +16,10 @@
 #'  shared across all jurisdictions in a fit
 #'
 #' @return a list describing this jurisdiction's stream: `convolution_matrix`,
-#'  `case_vec`, `prop_vec`, `dow_model`, `inits_values`, `observable_idx`
+#'  `case_vec`, `prop_vec`, `delays`, `dow_model`. Does not include GAM-based
+#'  initial values -- those are only needed for `infection_model_type =
+#'  'flat_prior'`, so they're computed lazily by `compute_flat_prior_inits()`
+#'  when `fit_waves()` actually needs them, rather than unconditionally here.
 #' @export
 prepare_observation_data <- function (observation_data,
                                       target_infection_dates) {
@@ -61,21 +64,9 @@ prepare_observation_data <- function (observation_data,
 
   convolution_matrix <- new_convolution_matrix(delays)
 
-  inits <- inits_by_jurisdiction(
-    case_vec,
-    delays,
-    prop_vec,
-    target_infection_dates)
-
-  inits_values <- rep(NA_real_, length(target_infection_dates))
-  observable_idx <- rep(FALSE, length(target_infection_dates))
-  inits_values[inits$observable_idx] <- inits$inits_values
-  observable_idx[inits$observable_idx] <- TRUE
-
   list(convolution_matrix = convolution_matrix,
        case_vec = case_vec,
        prop_vec = prop_vec,
-       dow_model = observation_data$dow_model,
-       inits_values = inits_values,
-       observable_idx = observable_idx)
+       delays = delays,
+       dow_model = observation_data$dow_model)
 }
diff --git a/R/stack_jurisdictions.R b/R/stack_jurisdictions.R
index a5cf252..12102b9 100644
--- a/R/stack_jurisdictions.R
+++ b/R/stack_jurisdictions.R
@@ -1,30 +1,44 @@
 #' Stack per-jurisdiction observation models together
 #'
-#' @description Combine a named list of per-jurisdiction observation model
-#'  bundles (each produced by `define_observation_model()`) into the
-#'  `date x jurisdiction` matrices used by the model-fitting functions. A
-#'  single-jurisdiction fit is just a length-1 list. This is where DOW
-#'  correction is actually applied (it needs to know `n_jurisdictions`, which
-#'  earlier, per-jurisdiction steps don't).
+#' @description Combine several per-jurisdiction observation model bundles
+#'  (each produced by `define_observation_model()`) into the
+#'  `date x jurisdiction` matrices used by the model-fitting functions. This
+#'  is where DOW correction is actually applied (it needs to know
+#'  `n_jurisdictions`, which earlier, per-jurisdiction steps don't).
+#'
+#'  Only needed when combining more than one jurisdiction -- a
+#'  single-jurisdiction fit can pass `define_observation_model()`'s output
+#'  straight to `fit_waves()` without calling this at all.
 #'
 #'  Every jurisdiction must supply the same set of observation streams, and
 #'  `dow_model` must agree across jurisdictions within a given stream; both
 #'  are intentional limitations for now rather than permanent design
 #'  decisions.
 #'
-#' @param observations_by_jurisdiction a named list of per-jurisdiction
-#'  observation model bundles (output of `define_observation_model()`),
-#'  named by jurisdiction
+#' @param ... per-jurisdiction observation model bundles (output of
+#'  `define_observation_model()`), named by jurisdiction, e.g.
+#'  `stack_jurisdictions(VIC = observation_model_vic, NSW = observation_model_nsw)`
 #'
 #' @return list of stacked observation model data, in the shape expected by
-#'  `fit_waves()`
+#'  `fit_waves()`, with class `epiwave_stacked_observations`
 #' @export
-stack_jurisdictions <- function (observations_by_jurisdiction) {
+stack_jurisdictions <- function (...) {
+  stack_jurisdictions_list(list(...))
+}
+
+#' @noRd
+stack_jurisdictions_list <- function (observations_by_jurisdiction) {
 
   jurisdictions <- names(observations_by_jurisdiction)
-  if (is.null(jurisdictions) || any(jurisdictions == "")) {
-    stop('`observations_by_jurisdiction` must be a fully named list, ',
-         'named by jurisdiction')
+
+  if (length(observations_by_jurisdiction) == 1 &&
+      (is.null(jurisdictions) || identical(jurisdictions, ""))) {
+    # a single jurisdiction has no real identity to preserve -- fit_waves()
+    # routes a raw define_observation_model() object through here unnamed
+    jurisdictions <- "1"
+  } else if (is.null(jurisdictions) || any(jurisdictions == "")) {
+    stop('when combining more than one jurisdiction, all arguments to ',
+         'stack_jurisdictions() must be named, named by jurisdiction')
   }
 
   target_infection_dates <- observations_by_jurisdiction[[1]]$target_infection_dates
@@ -48,15 +62,6 @@ stack_jurisdictions <- function (observations_by_jurisdiction) {
     stop('Every jurisdiction must supply the same set of observation streams')
   }
 
-  incidence_observable <- do.call(
-    cbind,
-    lapply(observations_by_jurisdiction, `[[`, "incidence_observable"))
-  incidence_observable_inits <- do.call(
-    cbind,
-    lapply(observations_by_jurisdiction, `[[`, "incidence_observable_inits"))
-  colnames(incidence_observable) <- jurisdictions
-  colnames(incidence_observable_inits) <- jurisdictions
-
   observation_model_data <- lapply(
     stream_ids,
     stack_stream,
@@ -65,13 +70,13 @@ stack_jurisdictions <- function (observations_by_jurisdiction) {
     target_infection_dates)
   names(observation_model_data) <- stream_ids
 
-  list(
+  out <- list(
     observation_model_data = observation_model_data,
     target_infection_dates = target_infection_dates,
-    target_jurisdictions = jurisdictions,
-    incidence_observable = incidence_observable,
-    incidence_observable_inits = incidence_observable_inits
+    target_jurisdictions = jurisdictions
   )
+  class(out) <- c("epiwave_stacked_observations", class(out))
+  out
 }
 
 #' Stack one observation stream across jurisdictions
@@ -112,9 +117,11 @@ stack_stream <- function (stream_id,
   }
 
   convolution_matrices <- lapply(per_jurisdiction, `[[`, "convolution_matrix")
+  delays <- lapply(per_jurisdiction, `[[`, "delays")
 
   list(
     convolution_matrices = convolution_matrices,
+    delays = delays,
     case_mat = case_mat,
     prop_mat = prop_mat
   )
diff --git a/man/define_observation_model.Rd b/man/define_observation_model.Rd
index b58e03a..04c8274 100644
--- a/man/define_observation_model.Rd
+++ b/man/define_observation_model.Rd
@@ -13,10 +13,13 @@ define_observation_model(target_infection_dates = NULL, ...)
 \code{define_observation_data()}), named by stream}
 }
 \value{
-list describing one jurisdiction's observation model
+list describing one jurisdiction's observation model, with class
+\code{epiwave_observation_model}
 }
 \description{
 Bundle one jurisdiction's observation streams (e.g. cases,
-hospitalisations) together. Call this once per jurisdiction;
-combine multiple jurisdictions later via \code{stack_jurisdictions()}.
+hospitalisations) together. Call this once per jurisdiction. A
+single-jurisdiction fit can pass this straight to \code{fit_waves()}; to
+combine multiple jurisdictions, pass several of these to
+\code{stack_jurisdictions()} first.
 }
diff --git a/man/fit_waves.Rd b/man/fit_waves.Rd
index df004a8..7af61f0 100644
--- a/man/fit_waves.Rd
+++ b/man/fit_waves.Rd
@@ -5,7 +5,7 @@
 \title{Fit waves to observations}
 \usage{
 fit_waves(
-  observations_by_jurisdiction,
+  observations,
   infection_model_type = c("flat_prior", "gp_infections", "gp_growth_rate",
     "gp_growth_rate_deriv"),
   n_chains = 4,
@@ -16,11 +16,10 @@ fit_waves(
 )
 }
 \arguments{
-\item{observations_by_jurisdiction}{a named list of per-jurisdiction
-observation model bundles (output of \code{define_observation_model()}),
-named by jurisdiction. A single-jurisdiction fit is a length-1 list.
-Combined internally via \code{stack_jurisdictions()}; jurisdictions sharing
-one fit are partially pooled (shared GP kernel hyperparameters, and
+\item{observations}{either a single jurisdiction's observation model
+(output of \code{define_observation_model()}), or several jurisdictions
+already combined via \code{stack_jurisdictions()}. Jurisdictions sharing one
+fit are partially pooled (shared GP kernel hyperparameters, and
 shared/hierarchical DOW priors where requested).}
 
 \item{infection_model_type}{options include 'flat_prior', 'infections',
diff --git a/man/prepare_observation_data.Rd b/man/prepare_observation_data.Rd
index 8cbb0d7..545e996 100644
--- a/man/prepare_observation_data.Rd
+++ b/man/prepare_observation_data.Rd
@@ -21,7 +21,9 @@ shared across all jurisdictions in a fit}
 }
 \value{
 a list describing this jurisdiction's stream: \code{convolution_matrix},
-\code{case_vec}, \code{prop_vec}, \code{dow_model}, \code{inits_values}, \code{observable_idx}
+\code{case_vec}, \code{prop_vec}, \code{delays}, \code{dow_model}. Does not include GAM-based
+initial values -- those are only needed for \code{infection_model_type = 'flat_prior'}, so they're computed lazily by \code{compute_flat_prior_inits()}
+when \code{fit_waves()} actually needs them, rather than unconditionally here.
 }
 \description{
 Prepare the data objects needed for one jurisdiction's
diff --git a/man/stack_jurisdictions.Rd b/man/stack_jurisdictions.Rd
index ae728f6..664dfb6 100644
--- a/man/stack_jurisdictions.Rd
+++ b/man/stack_jurisdictions.Rd
@@ -4,24 +4,27 @@
 \alias{stack_jurisdictions}
 \title{Stack per-jurisdiction observation models together}
 \usage{
-stack_jurisdictions(observations_by_jurisdiction)
+stack_jurisdictions(...)
 }
 \arguments{
-\item{observations_by_jurisdiction}{a named list of per-jurisdiction
-observation model bundles (output of \code{define_observation_model()}),
-named by jurisdiction}
+\item{...}{per-jurisdiction observation model bundles (output of
+\code{define_observation_model()}), named by jurisdiction, e.g.
+\code{stack_jurisdictions(VIC = observation_model_vic, NSW = observation_model_nsw)}}
 }
 \value{
 list of stacked observation model data, in the shape expected by
-\code{fit_waves()}
+\code{fit_waves()}, with class \code{epiwave_stacked_observations}
 }
 \description{
-Combine a named list of per-jurisdiction observation model
-bundles (each produced by \code{define_observation_model()}) into the
-\verb{date x jurisdiction} matrices used by the model-fitting functions. A
-single-jurisdiction fit is just a length-1 list. This is where DOW
-correction is actually applied (it needs to know \code{n_jurisdictions}, which
-earlier, per-jurisdiction steps don't).
+Combine several per-jurisdiction observation model bundles
+(each produced by \code{define_observation_model()}) into the
+\verb{date x jurisdiction} matrices used by the model-fitting functions. This
+is where DOW correction is actually applied (it needs to know
+\code{n_jurisdictions}, which earlier, per-jurisdiction steps don't).
+
+Only needed when combining more than one jurisdiction -- a
+single-jurisdiction fit can pass \code{define_observation_model()}'s output
+straight to \code{fit_waves()} without calling this at all.
 
 Every jurisdiction must supply the same set of observation streams, and
 \code{dow_model} must agree across jurisdictions within a given stream; both
diff --git a/tests/test_workflow/multi_jurisdiction_workflow.R b/tests/test_workflow/multi_jurisdiction_workflow.R
index 3a1fea8..f8e8371 100644
--- a/tests/test_workflow/multi_jurisdiction_workflow.R
+++ b/tests/test_workflow/multi_jurisdiction_workflow.R
@@ -2,11 +2,11 @@
 ##
 ## A fully self-contained example fitting two jurisdictions jointly, using
 ## fabricated data. Demonstrates:
-##   - jurisdiction is only ever a dimension at the fit_waves() step: each
-##     jurisdiction's data is prepared independently via
+##   - jurisdiction is only ever a dimension at the stack_jurisdictions()
+##     step: each jurisdiction's data is prepared independently via
 ##     define_observation_data()/define_observation_model() (exactly like the
-##     single-jurisdiction workflow), then combined into a *named* list keyed
-##     by jurisdiction and passed to fit_waves()
+##     single-jurisdiction workflow), then combined via stack_jurisdictions(),
+##     named by jurisdiction, before being passed to fit_waves()
 ##   - jurisdictions can have different/staggered date coverage -- every
 ##     per-jurisdiction vector is aligned to the same shared
 ##     target_infection_dates axis, so this is safe
@@ -50,15 +50,15 @@ observation_model_B <- define_observation_model(
   cases = make_jurisdiction_cases(60, 140, lambda = 30, dow_model = TRUE)
 )
 
-# named list, named by jurisdiction -- this is the only place jurisdiction
-# becomes a dimension
-observations_by_jurisdiction <- list(
+# named arguments to stack_jurisdictions(), named by jurisdiction -- this is
+# the only place jurisdiction becomes a dimension
+stacked <- stack_jurisdictions(
   jurisdiction_a = observation_model_A,
   jurisdiction_b = observation_model_B
 )
 
 fit_object <- fit_waves(
-  observations_by_jurisdiction = observations_by_jurisdiction,
+  observations = stacked,
   infection_model_type = "gp_growth_rate",
   n_chains = 2,
   warmup = 200,
@@ -74,7 +74,6 @@ stopifnot(identical(fit_object$jurisdictions, c("jurisdiction_a", "jurisdiction_
 ## NA for B in the underlying case matrix, and vice versa -- confirming row i
 ## means the same calendar date for every jurisdiction's column.
 
-case_mat <- observation_model_A$observation_model_data$cases$case_vec
 row_only_a <- 20   # within A's 10-90 window, outside B's 60-140 window
 row_only_b <- 135  # within B's window, outside A's
 
diff --git a/tests/test_workflow/single_jurisdiction_workflow.R b/tests/test_workflow/single_jurisdiction_workflow.R
index 19f9894..4c3f005 100644
--- a/tests/test_workflow/single_jurisdiction_workflow.R
+++ b/tests/test_workflow/single_jurisdiction_workflow.R
@@ -80,18 +80,19 @@ observation_model <- define_observation_model(
     proportion_infections = ihr)
 )
 
-# a single-jurisdiction fit is a length-1 named list, named by jurisdiction
-observations_by_jurisdiction <- list(example_jurisdiction = observation_model)
-
 ## -- fit ----------------------------------------------------------------
 
+# a single jurisdiction's define_observation_model() output goes straight to
+# fit_waves() -- no list()/naming step needed. For multiple jurisdictions,
+# see multi_jurisdiction_workflow.R's explicit stack_jurisdictions() step.
+#
 # flat_prior is the fastest infection model, good for a quick check like this
 # one; for a real analysis, 'gp_infections'/'gp_growth_rate'/
 # 'gp_growth_rate_deriv' model infections as a Gaussian Process instead (see
 # create_infection_timeseries()'s docs for how the three GP formulations
 # differ). Warmup/n_samples are kept small here for speed, not convergence.
 fit_object <- fit_waves(
-  observations_by_jurisdiction = observations_by_jurisdiction,
+  observations = observation_model,
   infection_model_type = "flat_prior",
   n_chains = 2,
   warmup = 200,
@@ -113,8 +114,8 @@ cat("posterior median infections range:",
 ## -- advanced: a time-varying delay distribution -----------------------------
 ##
 ## delay_from_infection can also be an already time-varying
-## discrete_pmf_series/discrete_weights_series, e.g. representing a change in
-## testing/reporting delay partway through the study period.
+## discrete_pmf_series, e.g. representing a change in testing/reporting
+## delay partway through the study period.
 
 notification_delay_early <- as_discrete_pmf(distributional::dist_gamma(shape = 3, rate = 0.5))
 notification_delay_late <- as_discrete_pmf(distributional::dist_gamma(shape = 2, rate = 0.6))
diff --git a/tests/test_workflow/test2.R b/tests/test_workflow/test2.R
index 0b7001c..634c402 100644
--- a/tests/test_workflow/test2.R
+++ b/tests/test_workflow/test2.R
@@ -85,16 +85,13 @@ jurisdiction_observation_models <- define_observation_model(
     proportion_infections = ihr)
 )
 
-# a single-jurisdiction fit is a length-1 named list, named by jurisdiction
-observations_by_jurisdiction <- setNames(
-  list(jurisdiction_observation_models),
-  jurisdictions)
-
 # set seed for DSE so we can reproduce
 set.seed(20250512)
 
+# a single jurisdiction's define_observation_model() output goes straight to
+# fit_waves() -- no list()/naming step needed
 fit_object <- fit_waves(
-  observations_by_jurisdiction = observations_by_jurisdiction,
+  observations = jurisdiction_observation_models,
   infection_model_type = 'flat_prior',#,# define_infection_model()'gp_growth_rate' #
   n_chains = 10,
   max_convergence_tries = 2,
diff --git a/tests/test_workflow/testing.R b/tests/test_workflow/testing.R
index c7870f1..b274974 100644
--- a/tests/test_workflow/testing.R
+++ b/tests/test_workflow/testing.R
@@ -84,13 +84,10 @@ jurisdiction_basic_observation_models <- define_observation_model(
 
 )
 
-# a single-jurisdiction fit is a length-1 named list, named by jurisdiction
-basic_observations_by_jurisdiction <- setNames(
-  list(jurisdiction_basic_observation_models),
-  jurisdictions)
-
+# a single jurisdiction's define_observation_model() output goes straight to
+# fit_waves() -- no list()/naming step needed
 fit_object <- fit_waves(
-  observations_by_jurisdiction = basic_observations_by_jurisdiction,
+  observations = jurisdiction_basic_observation_models,
   infection_model_type = 'gp_growth_rate'# 'flat_prior'#,define_infection_model() #
 )
 

From dc3854c40f5a02d26449aa51acd244c32bca9a99 Mon Sep 17 00:00:00 2001
From: Saras Windecker <smwindecker@users.noreply.github.com>
Date: Mon, 6 Jul 2026 10:44:04 +0300
Subject: [PATCH 09/14] Fix inits computation to use pre-DOW-correction
 proportion

compute_flat_prior_inits() was reading prop_mat from the stacked object,
but stack_stream() applies DOW correction to prop_mat before returning it --
so for any stream with dow_model = TRUE, inits_by_jurisdiction() was being
handed a greta array (the DOW-corrected proportion) instead of a plain
numeric one. It didn't error (there's already a greta_array branch in
inits_by_jurisdiction()), but that branch draws from the *prior* predictive
distribution of the DOW effect via greta::calculate(nsim = 100) -- no MCMC
has run yet at this point -- making the "deterministic smoothed guess"
inits computation silently stochastic, and adding an unnecessary
100-draw simulation per DOW-modelled stream x jurisdiction whenever
flat_prior is used.

Checked against the original, pre-refactor code to confirm this wasn't just
a style difference: it explicitly computed inits from prop_mat *before*
applying DOW correction, deliberately. My first refactor (PR #51) preserved
this ordering (DOW correction lived in stack_jurisdictions(), which ran
after prepare_observation_data()'s inits computation); deferring inits to
run after stacking (previous commit) inverted it by accident.

Fix: stack_stream() now also returns prop_mat_raw (captured before DOW
correction is applied), and compute_flat_prior_inits() reads that instead
of prop_mat. create_observation_model() is unaffected -- it still reads the
DOW-corrected prop_mat for the actual likelihood, as it should.

Verified: prop_mat_raw is a plain numeric column (not a greta_array) even
when dow_model = TRUE, so inits_by_jurisdiction() no longer needs
greta::calculate() at all for this case. Re-ran both synthetic workflow
scripts and testing.R's real-data fit (flat_prior + dow_model = TRUE, the
exact previously-affected combination) end to end. devtools::check() is
still 0 errors, 0 warnings.
---
 R/compute_flat_prior_inits.R |  6 +++++-
 R/stack_jurisdictions.R      | 11 ++++++++++-
 2 files changed, 15 insertions(+), 2 deletions(-)

diff --git a/R/compute_flat_prior_inits.R b/R/compute_flat_prior_inits.R
index a44bfd2..6378c34 100644
--- a/R/compute_flat_prior_inits.R
+++ b/R/compute_flat_prior_inits.R
@@ -30,10 +30,14 @@ compute_flat_prior_inits <- function (observations) {
 
     stream_inits <- lapply(stream_ids, function (stream_id) {
       stream <- observations$observation_model_data[[stream_id]]
+      # prop_mat_raw, not prop_mat: inits must use the pre-DOW-correction
+      # proportion (see stack_stream()) -- using the corrected one would make
+      # this depend on a prior-predictive draw of the DOW effect instead of
+      # being a deterministic approximation
       inits_by_jurisdiction(
         stream$case_mat[, j],
         stream$delays[[j]],
-        stream$prop_mat[, j],
+        stream$prop_mat_raw[, j],
         target_infection_dates)
     })
 
diff --git a/R/stack_jurisdictions.R b/R/stack_jurisdictions.R
index 12102b9..ebf538f 100644
--- a/R/stack_jurisdictions.R
+++ b/R/stack_jurisdictions.R
@@ -104,6 +104,14 @@ stack_stream <- function (stream_id,
   colnames(case_mat) <- colnames(prop_mat) <- jurisdictions
   rownames(case_mat) <- rownames(prop_mat) <- as.character(target_infection_dates)
 
+  # kept alongside the (possibly DOW-corrected) prop_mat below specifically
+  # for compute_flat_prior_inits(): GAM-based inits must be computed from the
+  # raw proportion, not the DOW-corrected one (DOW correction wasn't even
+  # known/applicable at the point inits used to be computed, pre-refactor;
+  # using the corrected version here would make inits depend on a
+  # prior-predictive draw of the DOW effect instead of being deterministic)
+  prop_mat_raw <- prop_mat
+
   dow_flags <- vapply(per_jurisdiction, `[[`, logical(1), "dow_model")
   if (length(unique(dow_flags)) > 1) {
     stop(sprintf(
@@ -123,6 +131,7 @@ stack_stream <- function (stream_id,
     convolution_matrices = convolution_matrices,
     delays = delays,
     case_mat = case_mat,
-    prop_mat = prop_mat
+    prop_mat = prop_mat,
+    prop_mat_raw = prop_mat_raw
   )
 }

From de399945879fd72aa98539c23cb204ac61710868 Mon Sep 17 00:00:00 2001
From: Saras Windecker <smwindecker@users.noreply.github.com>
Date: Mon, 6 Jul 2026 10:53:35 +0300
Subject: [PATCH 10/14] Make prepare_observation_data() internal (drop export)

It's called from exactly one place (define_observation_model()) and never
appears in any workflow script -- it's implementation detail, the same
category as stack_stream()/compute_flat_prior_inits(), which are already
@noRd. Exporting it overstated it as public API and was contributing to a
"why are there three equally-important functions here" feeling, when really
only define_observation_data()/define_observation_model() (plus
stack_jurisdictions()/fit_waves()) are meant to be called directly.

No logic changed -- devtools::check() still 0 errors/0 warnings, both
synthetic workflow scripts re-verified end to end.
---
 NAMESPACE                       |  1 -
 R/prepare_observation_data.R    |  5 +++--
 man/prepare_observation_data.Rd | 32 --------------------------------
 3 files changed, 3 insertions(+), 35 deletions(-)
 delete mode 100644 man/prepare_observation_data.Rd

diff --git a/NAMESPACE b/NAMESPACE
index f322b04..c68d349 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -20,7 +20,6 @@ export(implement_day_of_week)
 export(inits_by_jurisdiction)
 export(new_convolution_matrix)
 export(plot_infection_traj)
-export(prepare_observation_data)
 export(stack_jurisdictions)
 importFrom(cowplot,panel_border)
 importFrom(cowplot,theme_cowplot)
diff --git a/R/prepare_observation_data.R b/R/prepare_observation_data.R
index 5ee7384..422b1e8 100644
--- a/R/prepare_observation_data.R
+++ b/R/prepare_observation_data.R
@@ -1,7 +1,8 @@
 #' Prepare observation data
 #'
 #' @description Prepare the data objects needed for one jurisdiction's
-#'  observation model for a single data stream (e.g. cases). Multiple
+#'  observation model for a single data stream (e.g. cases), called
+#'  internally by `define_observation_model()` once per stream. Multiple
 #'  jurisdictions are combined later, by `stack_jurisdictions()`.
 #'
 #' @param observation_data data for one jurisdiction, one stream, as
@@ -20,7 +21,7 @@
 #'  initial values -- those are only needed for `infection_model_type =
 #'  'flat_prior'`, so they're computed lazily by `compute_flat_prior_inits()`
 #'  when `fit_waves()` actually needs them, rather than unconditionally here.
-#' @export
+#' @noRd
 prepare_observation_data <- function (observation_data,
                                       target_infection_dates) {
 
diff --git a/man/prepare_observation_data.Rd b/man/prepare_observation_data.Rd
deleted file mode 100644
index 545e996..0000000
--- a/man/prepare_observation_data.Rd
+++ /dev/null
@@ -1,32 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/prepare_observation_data.R
-\name{prepare_observation_data}
-\alias{prepare_observation_data}
-\title{Prepare observation data}
-\usage{
-prepare_observation_data(observation_data, target_infection_dates)
-}
-\arguments{
-\item{observation_data}{data for one jurisdiction, one stream, as
-returned by \code{define_observation_data()}.
-\code{timeseries_data} may be a plain data.frame/tibble with \code{date} and
-\code{value} columns -- it doesn't need to be pre-classed as
-\code{epiwave_timeseries}, see \code{as_epiwave_timeseries()}.
-\code{delay_from_infection} may be a single \code{discrete_pmf} object (replicated
-across \code{target_infection_dates}), or an already time-varying
-\code{discrete_pmf_series} (aligned to \code{target_infection_dates})}
-
-\item{target_infection_dates}{full date sequence of infection timeseries,
-shared across all jurisdictions in a fit}
-}
-\value{
-a list describing this jurisdiction's stream: \code{convolution_matrix},
-\code{case_vec}, \code{prop_vec}, \code{delays}, \code{dow_model}. Does not include GAM-based
-initial values -- those are only needed for \code{infection_model_type = 'flat_prior'}, so they're computed lazily by \code{compute_flat_prior_inits()}
-when \code{fit_waves()} actually needs them, rather than unconditionally here.
-}
-\description{
-Prepare the data objects needed for one jurisdiction's
-observation model for a single data stream (e.g. cases). Multiple
-jurisdictions are combined later, by \code{stack_jurisdictions()}.
-}

From 32952fdc5a2512a46199acc12349c8a1dc27cd6f Mon Sep 17 00:00:00 2001
From: Saras Windecker <smwindecker@users.noreply.github.com>
Date: Mon, 6 Jul 2026 11:02:00 +0300
Subject: [PATCH 11/14] Make inits_by_jurisdiction() internal (drop export)

Same situation as prepare_observation_data(): called from exactly one place
(compute_flat_prior_inits()), never from a workflow script. Pure internal
machinery.

No logic changed -- devtools::check() still 0 errors/0 warnings, both
synthetic workflow scripts re-verified end to end.
---
 NAMESPACE                    |  1 -
 R/inits_by_jurisdiction.R    |  5 ++++-
 man/inits_by_jurisdiction.Rd | 26 --------------------------
 3 files changed, 4 insertions(+), 28 deletions(-)
 delete mode 100644 man/inits_by_jurisdiction.Rd

diff --git a/NAMESPACE b/NAMESPACE
index c68d349..1a97cb0 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -17,7 +17,6 @@ export(define_observation_model)
 export(estimate_reff)
 export(fit_waves)
 export(implement_day_of_week)
-export(inits_by_jurisdiction)
 export(new_convolution_matrix)
 export(plot_infection_traj)
 export(stack_jurisdictions)
diff --git a/R/inits_by_jurisdiction.R b/R/inits_by_jurisdiction.R
index 851516e..f5cf4ba 100644
--- a/R/inits_by_jurisdiction.R
+++ b/R/inits_by_jurisdiction.R
@@ -1,5 +1,8 @@
 #' Define initial values for the infection timeseries
 #'
+#' @description Called internally by `compute_flat_prior_inits()`, once per
+#'  stream per jurisdiction.
+#'
 #' @param obs_data numeric vector of observed data, aligned to
 #'   `target_infection_dates` (`NA` for dates without observations)
 #' @param delays a `discrete_pmf_series` delay object, aligned to
@@ -11,7 +14,7 @@
 #' @importFrom mgcv gam predict.gam
 #'
 #' @return initial values as a list
-#' @export
+#' @noRd
 inits_by_jurisdiction <- function (obs_data,
                                    delays,
                                    obs_prop,
diff --git a/man/inits_by_jurisdiction.Rd b/man/inits_by_jurisdiction.Rd
deleted file mode 100644
index eb5f957..0000000
--- a/man/inits_by_jurisdiction.Rd
+++ /dev/null
@@ -1,26 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/inits_by_jurisdiction.R
-\name{inits_by_jurisdiction}
-\alias{inits_by_jurisdiction}
-\title{Define initial values for the infection timeseries}
-\usage{
-inits_by_jurisdiction(obs_data, delays, obs_prop, target_infection_dates)
-}
-\arguments{
-\item{obs_data}{numeric vector of observed data, aligned to
-\code{target_infection_dates} (\code{NA} for dates without observations)}
-
-\item{delays}{a \code{discrete_pmf_series} delay object, aligned to
-\code{target_infection_dates}}
-
-\item{obs_prop}{numeric vector of proportions, aligned to
-\code{target_infection_dates}, or a greta array of the same}
-
-\item{target_infection_dates}{infection date sequence}
-}
-\value{
-initial values as a list
-}
-\description{
-Define initial values for the infection timeseries
-}

From ca51fa477d146b2a5e6619d62c00aba5b4680d21 Mon Sep 17 00:00:00 2001
From: Saras Windecker <smwindecker@users.noreply.github.com>
Date: Mon, 6 Jul 2026 11:23:55 +0300
Subject: [PATCH 12/14] Delete
 create_epiwave_timeseries()/create_epiwave_fixed_timeseries(), simplify prop
 coercion

create_epiwave_fixed_timeseries() had zero call sites anywhere in the
package or workflow scripts, and did nothing as_epiwave_timeseries() doesn't
already do better (it takes a real date/value table directly, rather than
requiring separate dates + value arguments).

create_epiwave_timeseries() had exactly one remaining caller, in
prepare_observation_data()'s proportion_infections coercion -- but that
coercion step was itself redundant: as_matrix() already has an
as_matrix.numeric method that recycles a scalar or validates a same-length
vector against target_infection_dates, identically to what going through
create_epiwave_timeseries() first produced (verified: identical() output
both ways, for both the scalar and vector cases). prepare_observation_data()
now only branches on already-classed epiwave_timeseries objects (for date
validation); a bare numeric prop passes straight to as_matrix(), which
dispatches correctly on its own.

create_epiwave_timeseries.R now holds only create_epiwave_greta_timeseries()
(kept in epiwave rather than epiwave.params specifically because it depends
on greta, which epiwave.params should not) and as_epiwave_timeseries().

Verified: both synthetic workflow scripts and testing.R's real-data fit
produce identical output to before these changes (e.g. single-jurisdiction
posterior median range unchanged at 0-7185). devtools::check() still 0
errors, 0 warnings.
---
 NAMESPACE                              |  2 --
 R/create_epiwave_timeseries.R          | 48 --------------------------
 R/prepare_observation_data.R           | 11 +++---
 man/create_epiwave_fixed_timeseries.Rd | 23 ------------
 man/create_epiwave_timeseries.Rd       | 23 ------------
 5 files changed, 6 insertions(+), 101 deletions(-)
 delete mode 100644 man/create_epiwave_fixed_timeseries.Rd
 delete mode 100644 man/create_epiwave_timeseries.Rd

diff --git a/NAMESPACE b/NAMESPACE
index 1a97cb0..5f366be 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -7,9 +7,7 @@ export(as_epiwave_timeseries)
 export(as_matrix)
 export(compute_reff)
 export(create_dow_priors)
-export(create_epiwave_fixed_timeseries)
 export(create_epiwave_greta_timeseries)
-export(create_epiwave_timeseries)
 export(create_infection_timeseries)
 export(create_observation_model)
 export(define_observation_data)
diff --git a/R/create_epiwave_timeseries.R b/R/create_epiwave_timeseries.R
index c90f74d..8cd1f65 100644
--- a/R/create_epiwave_timeseries.R
+++ b/R/create_epiwave_timeseries.R
@@ -1,51 +1,3 @@
-#' Expand constant value into a long tibble
-#'
-#' @description The epiwave model functions expect data in a long format,
-#'   structured to have a value for every date. This function creates a
-#'   tibble of this structure from a single value that is replicated in each
-#'   cell.
-#'
-#' @param dates infection dates sequence
-#' @param value value to be replicated in each cell
-#'
-#' @importFrom tibble tibble
-#'
-#' @return a long tibble with the constant value replicated across all
-#'   dates
-#' @export
-create_epiwave_timeseries <- function(dates,
-                                      value) {
-  long_combined <- tibble::tibble(date = dates,
-                                  value = value)
-
-  class(long_combined) <- c("epiwave_timeseries", class(long_combined))
-  long_combined
-}
-
-#' Expand a fixed value into a long tibble
-#'
-#' @description The epiwave model functions expect data in a long format,
-#'   structured to have a value for every date. This function creates a
-#'   tibble of this structure from a single fixed numeric value that is
-#'   replicated in each cell.
-#'
-#' @param dates infection dates sequence
-#' @param value a fixed numeric value to be replicated in each cell
-#'
-#' @return a long tibble with the fixed value replicated across all
-#'   dates, with class `epiwave_fixed_timeseries`
-#' @export
-create_epiwave_fixed_timeseries <- function(dates,
-                                            value) {
-  timeseries <- create_epiwave_timeseries(
-    dates = dates,
-    value = value
-  )
-
-  class(timeseries) <- c("epiwave_fixed_timeseries", class(timeseries))
-  timeseries
-}
-
 #' Create a greta-compatible timeseries object
 #'
 #' @description The epiwave model functions expect data in a long format,
diff --git a/R/prepare_observation_data.R b/R/prepare_observation_data.R
index 422b1e8..35a8370 100644
--- a/R/prepare_observation_data.R
+++ b/R/prepare_observation_data.R
@@ -40,12 +40,13 @@ prepare_observation_data <- function (observation_data,
          'discrete_pmf_series object')
   }
 
+  # a bare numeric prop (scalar or a vector already in date order) needs no
+  # coercion at all -- as_matrix.numeric() recycles/validates it directly.
+  # Only already-classed epiwave_timeseries objects (including
+  # epiwave_greta_timeseries) need their dates checked against
+  # target_infection_dates first.
   prop <- observation_data$proportion_infections
-  if (!inherits(prop, 'epiwave_timeseries')) {
-    prop <- create_epiwave_timeseries(
-      dates = target_infection_dates,
-      value = prop)
-  } else {
+  if (inherits(prop, 'epiwave_timeseries')) {
     # epiwave_greta_timeseries stores dates nested under $timeseries$date
     # (it wraps a greta array alongside the date tibble), rather than a
     # flat $date column like epiwave_timeseries/epiwave_fixed_timeseries
diff --git a/man/create_epiwave_fixed_timeseries.Rd b/man/create_epiwave_fixed_timeseries.Rd
deleted file mode 100644
index b3d03d5..0000000
--- a/man/create_epiwave_fixed_timeseries.Rd
+++ /dev/null
@@ -1,23 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/create_epiwave_timeseries.R
-\name{create_epiwave_fixed_timeseries}
-\alias{create_epiwave_fixed_timeseries}
-\title{Expand a fixed value into a long tibble}
-\usage{
-create_epiwave_fixed_timeseries(dates, value)
-}
-\arguments{
-\item{dates}{infection dates sequence}
-
-\item{value}{a fixed numeric value to be replicated in each cell}
-}
-\value{
-a long tibble with the fixed value replicated across all
-dates, with class \code{epiwave_fixed_timeseries}
-}
-\description{
-The epiwave model functions expect data in a long format,
-structured to have a value for every date. This function creates a
-tibble of this structure from a single fixed numeric value that is
-replicated in each cell.
-}
diff --git a/man/create_epiwave_timeseries.Rd b/man/create_epiwave_timeseries.Rd
deleted file mode 100644
index 797b37a..0000000
--- a/man/create_epiwave_timeseries.Rd
+++ /dev/null
@@ -1,23 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/create_epiwave_timeseries.R
-\name{create_epiwave_timeseries}
-\alias{create_epiwave_timeseries}
-\title{Expand constant value into a long tibble}
-\usage{
-create_epiwave_timeseries(dates, value)
-}
-\arguments{
-\item{dates}{infection dates sequence}
-
-\item{value}{value to be replicated in each cell}
-}
-\value{
-a long tibble with the constant value replicated across all
-dates
-}
-\description{
-The epiwave model functions expect data in a long format,
-structured to have a value for every date. This function creates a
-tibble of this structure from a single value that is replicated in each
-cell.
-}

From 2b1751e9c4e3194cf5c3416d3dae79abcc8e3976 Mon Sep 17 00:00:00 2001
From: Saras Windecker <smwindecker@users.noreply.github.com>
Date: Mon, 6 Jul 2026 11:35:00 +0300
Subject: [PATCH 13/14] Rename create_epiwave_greta_timeseries() ->
 as_greta_timeseries(), un-export as_epiwave_timeseries()

Renamed to pair naturally with as_epiwave_timeseries() now that the two
create_*/dead functions removed from this file are gone -- as_* matches the
coercion-style naming already used elsewhere (epiwave.params::as_discrete_pmf(),
as_matrix()). The class it produces is renamed to match:
epiwave_greta_timeseries -> greta_timeseries (still inherits from
epiwave_timeseries, so existing inherits(x, 'epiwave_timeseries') checks
that need to catch both plain and greta-backed timeseries objects continue
to work unchanged). as_matrix.epiwave_greta_timeseries() is renamed to
as_matrix.greta_timeseries() to match, and
prepare_observation_data()'s inherits() check on proportion_infections
updated accordingly.

as_epiwave_timeseries() becomes internal (@noRd): checked and confirmed no
workflow script actually calls it by name -- users get its behaviour
indirectly (passing a raw data.frame to define_observation_data(), which
coerces internally). It remains available as an unexported helper rather
than being deleted, since it's still the one general-purpose coercion path,
just not meant as public API for now.

Verified: class(ihr) is now `greta_timeseries epiwave_timeseries list`,
dispatch to as_matrix.greta_timeseries() works correctly, both synthetic
workflow scripts produce identical output to before the rename (e.g.
single-jurisdiction posterior median range unchanged at 0-7185), and
testing.R's real-data fit runs to completion. devtools::check() still 0
errors, 0 warnings.
---
 NAMESPACE                                     |  5 ++--
 R/as_matrix.R                                 |  4 +--
 R/create_epiwave_timeseries.R                 | 16 +++++-----
 R/define_observation_data.R                   |  4 +--
 R/prepare_observation_data.R                  | 12 ++++----
 man/as_epiwave_timeseries.Rd                  | 30 -------------------
 ...a_timeseries.Rd => as_greta_timeseries.Rd} | 12 ++++----
 man/define_observation_data.Rd                |  4 +--
 .../single_jurisdiction_workflow.R            |  4 +--
 tests/test_workflow/test2.R                   |  4 +--
 tests/test_workflow/testing.R                 |  2 +-
 11 files changed, 33 insertions(+), 64 deletions(-)
 delete mode 100644 man/as_epiwave_timeseries.Rd
 rename man/{create_epiwave_greta_timeseries.Rd => as_greta_timeseries.Rd} (68%)

diff --git a/NAMESPACE b/NAMESPACE
index 5f366be..ca425d0 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,13 +1,12 @@
 # Generated by roxygen2: do not edit by hand
 
-S3method(as_matrix,epiwave_greta_timeseries)
 S3method(as_matrix,epiwave_timeseries)
+S3method(as_matrix,greta_timeseries)
 S3method(as_matrix,numeric)
-export(as_epiwave_timeseries)
+export(as_greta_timeseries)
 export(as_matrix)
 export(compute_reff)
 export(create_dow_priors)
-export(create_epiwave_greta_timeseries)
 export(create_infection_timeseries)
 export(create_observation_model)
 export(define_observation_data)
diff --git a/R/as_matrix.R b/R/as_matrix.R
index 805124d..2785850 100644
--- a/R/as_matrix.R
+++ b/R/as_matrix.R
@@ -51,11 +51,11 @@ as_matrix.epiwave_timeseries <- function (data, target_infection_dates, ...) {
 }
 
 #' @export
-as_matrix.epiwave_greta_timeseries <- function (data, target_infection_dates, ...) {
+as_matrix.greta_timeseries <- function (data, target_infection_dates, ...) {
 
   if (!identical(as.Date(data$timeseries$date), as.Date(target_infection_dates))) {
     stop('`data$timeseries$date` must match `target_infection_dates` ',
-         'exactly for an `epiwave_greta_timeseries` object.')
+         'exactly for a `greta_timeseries` object.')
   }
 
   ihr <- data$ihr
diff --git a/R/create_epiwave_timeseries.R b/R/create_epiwave_timeseries.R
index 8cd1f65..7201ac3 100644
--- a/R/create_epiwave_timeseries.R
+++ b/R/create_epiwave_timeseries.R
@@ -13,13 +13,13 @@
 #'
 #' @importFrom tibble tibble
 #'
-#' @return a list with class `epiwave_greta_timeseries` containing
-#'   components `timeseries` (a tibble of dates) and `ihr` (a greta array of
-#'   implied hospitalisation rates)
+#' @return a list with class `greta_timeseries` containing components
+#'   `timeseries` (a tibble of dates) and `ihr` (a greta array of implied
+#'   hospitalisation rates)
 #' @export
-create_epiwave_greta_timeseries <- function(dates,
-                                            car,
-                                            chr_prior) {
+as_greta_timeseries <- function(dates,
+                                car,
+                                chr_prior) {
   long_unique <- tibble::tibble(date = dates)
 
   dim(chr_prior) <- length(dates)
@@ -33,7 +33,7 @@ create_epiwave_greta_timeseries <- function(dates,
   long_combined <- list(timeseries = long_unique,
                         ihr = ihr_greta)
 
-  class(long_combined) <- c("epiwave_greta_timeseries",
+  class(long_combined) <- c("greta_timeseries",
                             "epiwave_timeseries",
                             class(long_combined))
   long_combined
@@ -59,7 +59,7 @@ create_epiwave_greta_timeseries <- function(dates,
 #'   data.frame/tibble with `date` and `value` columns
 #'
 #' @return an `epiwave_timeseries` or numeric object, ready for `as_matrix()`
-#' @export
+#' @noRd
 as_epiwave_timeseries <- function(data) {
 
   if (inherits(data, "epiwave_timeseries") || is.numeric(data)) {
diff --git a/R/define_observation_data.R b/R/define_observation_data.R
index e9e3b31..d9f3e6f 100644
--- a/R/define_observation_data.R
+++ b/R/define_observation_data.R
@@ -6,8 +6,8 @@
 #'
 #' @param timeseries_data timeseries data for data of interest, for one
 #'  jurisdiction. A plain data.frame/tibble with `date` and `value` columns
-#'  is fine (gaps in date coverage are fine too) -- it doesn't need to be
-#'  pre-classed, see `as_epiwave_timeseries()`
+#'  is fine (gaps in date coverage are fine too) -- it's coerced
+#'  automatically, no need to pre-class it
 #' @param delay_from_infection a `discrete_pmf` object (replicated across
 #'  dates), or an already time-varying `discrete_pmf_series` object
 #' @param proportion_infections proportion data
diff --git a/R/prepare_observation_data.R b/R/prepare_observation_data.R
index 35a8370..bde5f13 100644
--- a/R/prepare_observation_data.R
+++ b/R/prepare_observation_data.R
@@ -43,14 +43,14 @@ prepare_observation_data <- function (observation_data,
   # a bare numeric prop (scalar or a vector already in date order) needs no
   # coercion at all -- as_matrix.numeric() recycles/validates it directly.
   # Only already-classed epiwave_timeseries objects (including
-  # epiwave_greta_timeseries) need their dates checked against
-  # target_infection_dates first.
+  # greta_timeseries) need their dates checked against target_infection_dates
+  # first.
   prop <- observation_data$proportion_infections
   if (inherits(prop, 'epiwave_timeseries')) {
-    # epiwave_greta_timeseries stores dates nested under $timeseries$date
-    # (it wraps a greta array alongside the date tibble), rather than a
-    # flat $date column like epiwave_timeseries/epiwave_fixed_timeseries
-    prop_dates <- if (inherits(prop, 'epiwave_greta_timeseries')) {
+    # greta_timeseries stores dates nested under $timeseries$date (it wraps
+    # a greta array alongside the date tibble), rather than a flat $date
+    # column like epiwave_timeseries/epiwave_fixed_timeseries
+    prop_dates <- if (inherits(prop, 'greta_timeseries')) {
       prop$timeseries$date
     } else {
       prop$date
diff --git a/man/as_epiwave_timeseries.Rd b/man/as_epiwave_timeseries.Rd
deleted file mode 100644
index c99869a..0000000
--- a/man/as_epiwave_timeseries.Rd
+++ /dev/null
@@ -1,30 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/create_epiwave_timeseries.R
-\name{as_epiwave_timeseries}
-\alias{as_epiwave_timeseries}
-\title{Coerce a date/value table to an epiwave_fixed_timeseries object}
-\usage{
-as_epiwave_timeseries(data)
-}
-\arguments{
-\item{data}{an \code{epiwave_timeseries} object, a numeric value/vector, or a
-data.frame/tibble with \code{date} and \code{value} columns}
-}
-\value{
-an \code{epiwave_timeseries} or numeric object, ready for \code{as_matrix()}
-}
-\description{
-Observation data (case counts, hospitalisation counts) is
-naturally a table of \code{date}/\code{value} pairs with real per-date
-observations, rather than a
-single value replicated across dates -- unlike delay distributions or
-proportions, which usually apply uniformly. This function takes such a
-table (with genuinely partial date coverage, e.g. missing weekends, is
-fine -- gaps are handled downstream by \code{as_matrix()}) and classes it as
-an \code{epiwave_fixed_timeseries} object, so it doesn't need to be
-hand-classed with \code{class(x) <- c(...)}.
-
-Already-classed \code{epiwave_timeseries} objects and plain numeric
-values/vectors are returned unchanged, so this is safe to call
-unconditionally on anything accepted as observation data.
-}
diff --git a/man/create_epiwave_greta_timeseries.Rd b/man/as_greta_timeseries.Rd
similarity index 68%
rename from man/create_epiwave_greta_timeseries.Rd
rename to man/as_greta_timeseries.Rd
index a91e732..6e31cf2 100644
--- a/man/create_epiwave_greta_timeseries.Rd
+++ b/man/as_greta_timeseries.Rd
@@ -1,10 +1,10 @@
 % Generated by roxygen2: do not edit by hand
 % Please edit documentation in R/create_epiwave_timeseries.R
-\name{create_epiwave_greta_timeseries}
-\alias{create_epiwave_greta_timeseries}
+\name{as_greta_timeseries}
+\alias{as_greta_timeseries}
 \title{Create a greta-compatible timeseries object}
 \usage{
-create_epiwave_greta_timeseries(dates, car, chr_prior)
+as_greta_timeseries(dates, car, chr_prior)
 }
 \arguments{
 \item{dates}{infection dates sequence}
@@ -16,9 +16,9 @@ create_epiwave_greta_timeseries(dates, car, chr_prior)
 the case hospitalisation rate}
 }
 \value{
-a list with class \code{epiwave_greta_timeseries} containing
-components \code{timeseries} (a tibble of dates) and \code{ihr} (a greta array of
-implied hospitalisation rates)
+a list with class \code{greta_timeseries} containing components
+\code{timeseries} (a tibble of dates) and \code{ihr} (a greta array of implied
+hospitalisation rates)
 }
 \description{
 The epiwave model functions expect data in a long format,
diff --git a/man/define_observation_data.Rd b/man/define_observation_data.Rd
index 26dd540..0a9eda9 100644
--- a/man/define_observation_data.Rd
+++ b/man/define_observation_data.Rd
@@ -14,8 +14,8 @@ define_observation_data(
 \arguments{
 \item{timeseries_data}{timeseries data for data of interest, for one
 jurisdiction. A plain data.frame/tibble with \code{date} and \code{value} columns
-is fine (gaps in date coverage are fine too) -- it doesn't need to be
-pre-classed, see \code{as_epiwave_timeseries()}}
+is fine (gaps in date coverage are fine too) -- it's coerced
+automatically, no need to pre-class it}
 
 \item{delay_from_infection}{a \code{discrete_pmf} object (replicated across
 dates), or an already time-varying \code{discrete_pmf_series} object}
diff --git a/tests/test_workflow/single_jurisdiction_workflow.R b/tests/test_workflow/single_jurisdiction_workflow.R
index 4c3f005..bebda45 100644
--- a/tests/test_workflow/single_jurisdiction_workflow.R
+++ b/tests/test_workflow/single_jurisdiction_workflow.R
@@ -4,7 +4,7 @@
 ## jurisdiction, using fabricated data (no external/not-synced files needed).
 ## Demonstrates the current API end to end:
 ##   - raw data.frames passed straight to define_observation_data() (no
-##     manual class(x) <- c(...) needed -- see as_epiwave_timeseries())
+##     manual class(x) <- c(...) needed -- no need to pre-class)
 ##   - delay distributions built with epiwave.params (discrete_pmf, combined
 ##     via `+`/add_discrete()), including a time-varying discrete_pmf_series
 ##   - proportion_infections as both a fixed value and a greta-array-derived
@@ -58,7 +58,7 @@ car <- 0.5
 # and CAR via a greta array, so it can be estimated jointly with the rest of
 # the model rather than fixed
 chr <- greta::uniform(0, 1)
-ihr <- create_epiwave_greta_timeseries(
+ihr <- as_greta_timeseries(
   dates = target_infection_dates,
   car = car,
   chr_prior = chr)
diff --git a/tests/test_workflow/test2.R b/tests/test_workflow/test2.R
index 634c402..64d0fea 100644
--- a/tests/test_workflow/test2.R
+++ b/tests/test_workflow/test2.R
@@ -24,7 +24,7 @@ colnames(notif_dat) <- c('study_area', 'date')
 notif_dat <- notif_dat |>
   tidyr::pivot_longer(!date, names_to = "jurisdiction", values_to = "value")
 # no manual class(notif_dat) <- c(...) needed -- define_observation_data()
-# accepts a plain date/value data.frame directly (see as_epiwave_timeseries())
+# accepts a plain date/value data.frame directly (no need to pre-class)
 
 # hospitalisation counts
 hosp_dat <- 'simdata/sim_study_hosp.csv' |>
@@ -47,7 +47,7 @@ car <- 0.42 # minimum viable product scenario
 chr <- greta::uniform(0, 1)
 # ihr <- chr * car
 # wrapper for ihr specific flow
-ihr <- create_epiwave_greta_timeseries(
+ihr <- as_greta_timeseries(
   dates = infection_days,
   car = car,
   chr_prior = chr)
diff --git a/tests/test_workflow/testing.R b/tests/test_workflow/testing.R
index b274974..c24263e 100644
--- a/tests/test_workflow/testing.R
+++ b/tests/test_workflow/testing.R
@@ -52,7 +52,7 @@ car <- 0.75
 
 # create IHR
 chr <- greta::uniform(0, 1)
-ihr <- create_epiwave_greta_timeseries(
+ihr <- as_greta_timeseries(
   dates = infection_days,
   car = car,
   chr_prior = chr)

From 63dc618046d5309b192ef30d0f3f1c646b9efa30 Mon Sep 17 00:00:00 2001
From: Saras Windecker <smwindecker@users.noreply.github.com>
Date: Mon, 6 Jul 2026 11:47:31 +0300
Subject: [PATCH 14/14] Add usage examples to README

Two minimal examples -- single jurisdiction, and multiple jurisdictions
combined via stack_jurisdictions() -- mirroring
tests/test_workflow/single_jurisdiction_workflow.R and
multi_jurisdiction_workflow.R, trimmed to the essential shape for a README.
Both verified to run to completion exactly as written before adding them;
chunks are eval = FALSE in README.Rmd since running them needs a full
greta/Python setup, which shouldn't be a precondition for a routine README
re-knit.

README.Rmd had been missing from the repo (only README.md existed) despite
a pre-commit hook expecting both to move together; it's restored here from
an old, partially-filled scaffold (its Installation section still suggested
pak::pak() and its Example section was a never-filled-in placeholder).
Updated Installation to match the remotes::install_github() already
decided in README.md, and replaced the Example placeholder with the new
Usage section. Left the empty Citation section and placeholder Support
text untouched rather than inventing content, and left the informal draft
paragraph after the srr-tags chunk untouched (out of scope for this
change).

Also gitignores README.html, a preview-render artifact that isn't source.
---
 .gitignore |   1 +
 README.Rmd | 129 +++++++++++++++++++++++++++++++++++++++++++++++++++++
 README.md  | 122 +++++++++++++++++++++++++++++++++++++++++++++-----
 3 files changed, 241 insertions(+), 11 deletions(-)
 create mode 100644 README.Rmd

diff --git a/.gitignore b/.gitignore
index ef1c072..e7225c7 100644
--- a/.gitignore
+++ b/.gitignore
@@ -9,3 +9,4 @@ inst/doc
 figures/
 simdata/
 .positai
+README.html
diff --git a/README.Rmd b/README.Rmd
new file mode 100644
index 0000000..3e7ebea
--- /dev/null
+++ b/README.Rmd
@@ -0,0 +1,129 @@
+---
+output: github_document
+---
+
+```{r, include = FALSE}
+knitr::opts_chunk$set(
+  collapse = TRUE,
+  comment = "#>",
+  fig.path = "man/figures/README-",
+  out.width = "100%"
+)
+```
+
+# epiwave
+
+```{r srr-tags, eval = FALSE, echo = FALSE}
+#' @srrstats {G1.0, BS4.0} The statistical software `greta` forms the underlying bayesian inference used in this package. The citation provided for that package also documents the sampling algorithms. 
+#' @srrstats {G1.1} This package is the first implementation of a novel algorithm to estimate infection incidence timeseries. It relies underneath on existing sampling algorithms from the `greta` software. 
+#' @srrstatsTODO {G1.5} *Software should include all code necessary to reproduce results which form the basis of performance claims made in associated publications.*
+#' @srrstatsTODO {BS1.2} *Description of how to specify prior distributions, both in textual form describing the general principles of specifying prior distributions, along with more applied descriptions and examples, within:*
+#' @srrstatsTODO {BS1.2a} *The main package `README`, either as textual description or example code*
+```
+
+This package is the first implementation of a novel algorithm to estimate infection incidence timeseries. It relies underneath on existing sampling algorithms from the `greta` software. we get infection incidence, no other package does this. estimating r effective. we could benchmark but pred accuracy not super important for users. 
+
+## Package to estimate epidemiological parameters using a Gaussian Process for infection timeseries
+
+[![source](https://img.shields.io/badge/source-GitHub-success?style=flat&labelColor=gray)](https://github.com/idem-lab/epiwave)
+[![wip](https://www.repostatus.org/badges/latest/wip.svg)](https://www.repostatus.org/#wip)
+[![v.0.1.0-alpha](https://zenodo.org/badge/720975798.svg)](https://zenodo.org/doi/10.5281/zenodo.10398079)
+
+`epiwave` estimates a time-series of the daily number of new infections, using a time-series of daily reported case numbers and an estimated distribution of the time delay between the onset of infection and case notification.
+This model is unique because it integrates multiple data sources to estimate the underlying infection dynamics.
+Besides case, hospitalisation, and death notification data, these data sources could include wastewater data and serological data.
+From the estimated time-series of daily new infections we can then estimate useful epidemic parameters such as the effective reproduction number over time.
+This package will be especially useful when case incidence data, the most commonly used source of information for estimating effective reproduction number, is unreliable or challenging to use, but other data sources are available.
+
+## Installation
+
+`epiwave` can be installed from [GitHub](https://github.com/) with:
+
+``` r
+remotes::install_github('idem-lab/epiwave')
+```
+
+## Usage
+
+### One jurisdiction
+
+```{r example-single, eval = FALSE}
+library(epiwave)
+library(epiwave.params)
+library(distributional)
+
+dates <- seq(as.Date("2024-01-01"), by = "day", length.out = 150)
+
+# delay from infection to case notification
+cases_delay <- as_discrete_pmf(distributional::dist_gamma(shape = 3, rate = 0.5))
+
+# case counts: a plain date/value data.frame, no need to pre-class it
+notif_dat <- data.frame(
+  date = dates[20:130],
+  value = rpois(111, lambda = 50)
+)
+
+observation_model <- define_observation_model(
+  target_infection_dates = dates,
+  cases = define_observation_data(
+    timeseries_data = notif_dat,
+    delay_from_infection = cases_delay,
+    proportion_infections = 0.5,
+    dow_model = TRUE)
+)
+
+# a single jurisdiction goes straight to fit_waves(), no combining step needed
+fit_object <- fit_waves(
+  observations = observation_model,
+  infection_model_type = "flat_prior"
+)
+```
+
+### Multiple jurisdictions
+
+Jurisdictions are prepared independently, then combined explicitly via
+`stack_jurisdictions()`, named by jurisdiction. Jurisdictions sharing one
+fit are partially pooled (shared GP kernel hyperparameters, and a shared
+hierarchical day-of-week prior where requested).
+
+```{r example-multi, eval = FALSE}
+observation_model_b <- define_observation_model(
+  target_infection_dates = dates,
+  cases = define_observation_data(
+    timeseries_data = data.frame(date = dates[60:140], value = rpois(81, lambda = 30)),
+    delay_from_infection = cases_delay,
+    proportion_infections = 0.5,
+    dow_model = TRUE)
+)
+
+stacked <- stack_jurisdictions(
+  jurisdiction_a = observation_model,
+  jurisdiction_b = observation_model_b
+)
+
+fit_object <- fit_waves(
+  observations = stacked,
+  infection_model_type = "gp_growth_rate"
+)
+```
+
+## Citation
+
+When using this package, please cite the underlying statistical software, `greta` as well as the package itself: 
+
+
+
+## Contribution
+
+This is a work in progress. 
+If you see any mistakes in the package (branch `main`), let us know by logging a bug on the [issues](https://github.com/idem-lab/epiwave/issues) page. 
+
+## Code of Conduct
+
+Please note that the `epiwave` project is released with a [Contributor Code of Conduct](https://idem-lab.github.io/epiwave/CODE_OF_CONDUCT.html). By contributing to this project, you agree to abide by its terms.
+
+## Support
+
+This project was supported by ???? BDSS??? [the Australia-Aotearoa Consortium for Epidemic Forecasting and Analytics.]
+
+<a href="https://acefa-hubs.github.io/"><img src="man/figures/ACEFA.png" align = "center" height="150" alt="EpiStrainDynamics website" /></a>
diff --git a/README.md b/README.md
index 4f9c984..36efcf8 100644
--- a/README.md
+++ b/README.md
@@ -1,25 +1,125 @@
----
-editor_options: 
-  markdown: 
-    wrap: sentence
----
 
 # epiwave
 
+This package is the first implementation of a novel algorithm to
+estimate infection incidence timeseries. It relies underneath on
+existing sampling algorithms from the `greta` software. we get infection
+incidence, no other package does this. estimating r effective. we could
+benchmark but pred accuracy not super important for users.
+
 ## Package to estimate epidemiological parameters using a Gaussian Process for infection timeseries
 
 [![source](https://img.shields.io/badge/source-GitHub-success?style=flat&labelColor=gray)](https://github.com/idem-lab/epiwave)
 [![wip](https://www.repostatus.org/badges/latest/wip.svg)](https://www.repostatus.org/#wip)
 [![v.0.1.0-alpha](https://zenodo.org/badge/720975798.svg)](https://zenodo.org/doi/10.5281/zenodo.10398079)
 
-`epiwave` estimates a time-series of the daily number of new infections, using a time-series of daily reported case numbers and an estimated distribution of the time delay between the onset of infection and case notification.
-This model is unique because it integrates multiple data sources to estimate the underlying infection dynamics.
-Besides case, hospitalisation, and death notification data, these data sources could include wastewater data and serological data.
-From the estimated time-series of daily new infections we can then estimate useful epidemic parameters such as the effective reproduction number over time.
-This package will be especially useful when case incidence data, the most commonly used source of information for estimating effective reproduction number, is unreliable or challenging to use, but other data sources are available.
+`epiwave` estimates a time-series of the daily number of new infections,
+using a time-series of daily reported case numbers and an estimated
+distribution of the time delay between the onset of infection and case
+notification. This model is unique because it integrates multiple data
+sources to estimate the underlying infection dynamics. Besides case,
+hospitalisation, and death notification data, these data sources could
+include wastewater data and serological data. From the estimated
+time-series of daily new infections we can then estimate useful epidemic
+parameters such as the effective reproduction number over time. This
+package will be especially useful when case incidence data, the most
+commonly used source of information for estimating effective
+reproduction number, is unreliable or challenging to use, but other data
+sources are available.
+
+## Installation
 
-`epiwave` can be installed with GitHub, as follows:
+`epiwave` can be installed from [GitHub](https://github.com/) with:
 
 ``` r
 remotes::install_github('idem-lab/epiwave')
 ```
+
+## Usage
+
+### One jurisdiction
+
+``` r
+library(epiwave)
+library(epiwave.params)
+library(distributional)
+
+dates <- seq(as.Date("2024-01-01"), by = "day", length.out = 150)
+
+# delay from infection to case notification
+cases_delay <- as_discrete_pmf(distributional::dist_gamma(shape = 3, rate = 0.5))
+
+# case counts: a plain date/value data.frame, no need to pre-class it
+notif_dat <- data.frame(
+  date = dates[20:130],
+  value = rpois(111, lambda = 50)
+)
+
+observation_model <- define_observation_model(
+  target_infection_dates = dates,
+  cases = define_observation_data(
+    timeseries_data = notif_dat,
+    delay_from_infection = cases_delay,
+    proportion_infections = 0.5,
+    dow_model = TRUE)
+)
+
+# a single jurisdiction goes straight to fit_waves(), no combining step needed
+fit_object <- fit_waves(
+  observations = observation_model,
+  infection_model_type = "flat_prior"
+)
+```
+
+### Multiple jurisdictions
+
+Jurisdictions are prepared independently, then combined explicitly via
+`stack_jurisdictions()`, named by jurisdiction. Jurisdictions sharing
+one fit are partially pooled (shared GP kernel hyperparameters, and a
+shared hierarchical day-of-week prior where requested).
+
+``` r
+observation_model_b <- define_observation_model(
+  target_infection_dates = dates,
+  cases = define_observation_data(
+    timeseries_data = data.frame(date = dates[60:140], value = rpois(81, lambda = 30)),
+    delay_from_infection = cases_delay,
+    proportion_infections = 0.5,
+    dow_model = TRUE)
+)
+
+stacked <- stack_jurisdictions(
+  jurisdiction_a = observation_model,
+  jurisdiction_b = observation_model_b
+)
+
+fit_object <- fit_waves(
+  observations = stacked,
+  infection_model_type = "gp_growth_rate"
+)
+```
+
+## Citation
+
+When using this package, please cite the underlying statistical
+software, `greta` as well as the package itself:
+
+## Contribution
+
+This is a work in progress. If you see any mistakes in the package
+(branch `main`), let us know by logging a bug on the
+[issues](https://github.com/idem-lab/epiwave/issues) page.
+
+## Code of Conduct
+
+Please note that the `epiwave` project is released with a [Contributor
+Code of
+Conduct](https://idem-lab.github.io/epiwave/CODE_OF_CONDUCT.html). By
+contributing to this project, you agree to abide by its terms.
+
+## Support
+
+This project was supported by ???? BDSS??? \[the Australia-Aotearoa
+Consortium for Epidemic Forecasting and Analytics.\]
+
+<a href="https://acefa-hubs.github.io/"><img src="man/figures/ACEFA.png" align = "center" height="150" alt="EpiStrainDynamics website" /></a>