From c2b8e6597d37791c57b9d1ff944e3dc1c2946ec8 Mon Sep 17 00:00:00 2001 From: Maxwell Murphy Date: Tue, 15 Sep 2026 08:44:12 -0700 Subject: [PATCH 1/3] Release 3.7.0: seedable MCMC runs Bump version to 3.7.0 and close out the NEWS entry for run_mcmc(seed =). Use \doi{} for the namibia_data source so CRAN incoming checks are clean. --- DESCRIPTION | 2 +- NEWS.md | 2 +- R/data.R | 2 +- man/namibia_data.Rd | 2 +- 4 files changed, 4 insertions(+), 4 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 7818df2..b9dad92 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: moire Title: Multiplicity of Infection and Allele Frequency Recovery from Noisy Polyallelic Genetics Data -Version: 3.6.0 +Version: 3.7.0 Authors@R: c( person(given = "Maxwell", family = "Murphy", diff --git a/NEWS.md b/NEWS.md index 46a99ee..173e690 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,4 +1,4 @@ -# moire (development version) +# moire 3.7.0 ## New Features diff --git a/R/data.R b/R/data.R index db01475..68377ba 100644 --- a/R/data.R +++ b/R/data.R @@ -22,7 +22,7 @@ #' \item{locus}{Genetic locus} #' \item{allele}{Allele observed} #' } -#' @source \url{https://doi.org/10.7554/eLife.43510.018} +#' @source \doi{10.7554/eLife.43510.018} "namibia_data" #' Allele frequencies for different regions diff --git a/man/namibia_data.Rd b/man/namibia_data.Rd index 7f4656a..1010c60 100644 --- a/man/namibia_data.Rd +++ b/man/namibia_data.Rd @@ -17,7 +17,7 @@ A data frame with 7 columns and 97214 rows: } } \source{ -\url{https://doi.org/10.7554/eLife.43510.018} +\doi{10.7554/eLife.43510.018} } \usage{ namibia_data From 35c99737ce1ff70073f5f9538436c1ec7c03b5ab Mon Sep 17 00:00:00 2001 From: Maxwell Murphy Date: Tue, 15 Sep 2026 09:03:59 -0700 Subject: [PATCH 2/3] Prepare 3.7.0 for CRAN submission Add runnable examples and @return documentation to every exported function, cite the Bioinformatics paper in Description, declare parallel in Imports, drop the unused markdown Suggest, and credit the bundled spline.h author as cph. Fix swapped epsilon_pos/epsilon_neg docs in simulate_observed_allele(). Add cran-comments.md. --- .Rbuildignore | 2 + DESCRIPTION | 15 +++- NEWS.md | 9 +++ R/data.R | 8 ++ R/mcmc.R | 35 +++++++++ R/simulator.R | 67 ++++++++++++++++- R/summary.R | 90 +++++++++++++++++++++++ R/utils.R | 37 +++++++++- cran-comments.md | 23 ++++++ man/calculate_he.Rd | 6 ++ man/calculate_med_allele_freqs.Rd | 8 ++ man/calculate_naive_allele_frequencies.Rd | 7 ++ man/calculate_naive_coi.Rd | 6 ++ man/calculate_naive_coi_offset.Rd | 6 ++ man/load_delimited_data.Rd | 13 ++++ man/load_long_form_data.Rd | 21 ++++++ man/mcmc_results.Rd | 5 +- man/plot_chain_swaps.Rd | 6 +- man/rdirichlet.Rd | 7 ++ man/run_mcmc.Rd | 36 +++++++++ man/simulate_allele_frequencies.Rd | 7 ++ man/simulate_data.Rd | 16 +++- man/simulate_observed_allele.Rd | 14 +++- man/simulate_observed_genotype.Rd | 11 +++ man/simulate_sample_coi.Rd | 7 ++ man/simulate_sample_genotype.Rd | 11 +++ man/simulated_data.Rd | 5 +- man/summarize_allele_freq_fn.Rd | 10 +++ man/summarize_allele_freqs.Rd | 9 +++ man/summarize_coi.Rd | 9 +++ man/summarize_effective_coi.Rd | 9 +++ man/summarize_epsilon_neg.Rd | 8 ++ man/summarize_epsilon_pos.Rd | 8 ++ man/summarize_he.Rd | 8 ++ man/summarize_relatedness.Rd | 9 +++ 35 files changed, 533 insertions(+), 15 deletions(-) create mode 100644 cran-comments.md diff --git a/.Rbuildignore b/.Rbuildignore index 7ae6c9e..4ea2ab6 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -24,3 +24,5 @@ scratch\.R ^debug-data$ ^cpp$ ^CONTEXT\.md$ +^cran-comments\.md$ +^CRAN-SUBMISSION$ diff --git a/DESCRIPTION b/DESCRIPTION index b9dad92..b624c44 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -11,9 +11,16 @@ Authors@R: c( family = "Greenhouse", role = c("aut", "ths"), email = "bryan.greenhouse@ucsf.edu", - comment = c(ORCID = "0000-0003-0287-9111"))) -Description: A Markov Chain Monte Carlo (MCMC) based approach to Bayesian estimation of individual level multiplicity of infection, within host relatedness, and - population allele frequencies from polyallelic genetic data. + comment = c(ORCID = "0000-0003-0287-9111")), + person(given = "Tino", + family = "Kluge", + role = "cph", + comment = "Author of the bundled cubic spline library (src/include/spline)")) +Description: A Markov Chain Monte Carlo (MCMC) based approach to Bayesian + estimation of individual level multiplicity of infection, within host + relatedness, and population allele frequencies from polyallelic genetic + data. Implements the model described in Murphy and Greenhouse (2024) + . License: GPL (>= 3) Encoding: UTF-8 LazyData: true @@ -31,6 +38,7 @@ Imports: dplyr, tidyr, stats, + parallel, purrr, rlang, ggplot2 @@ -41,7 +49,6 @@ RoxygenNote: 7.3.3 Suggests: knitr, rmarkdown, - markdown, forcats, testthat (>= 3.0.0), parallelly diff --git a/NEWS.md b/NEWS.md index 173e690..f2ebc1a 100644 --- a/NEWS.md +++ b/NEWS.md @@ -8,6 +8,15 @@ `$seed` (printed if the run errors), with offsets as `$chain_seeds`. Replay needs the same sampler settings; `num_cores` does not affect the draws. +## Other Changes + +* Every exported function now has runnable examples and documents its return + value, in preparation for CRAN submission. +* Corrected the swapped `epsilon_pos` / `epsilon_neg` descriptions in + `simulate_observed_allele()`. +* `parallel` is now declared in `Imports`; `markdown` was dropped from + `Suggests` as nothing used it. + # moire 3.6.0 ## New Features diff --git a/R/data.R b/R/data.R index 68377ba..451f39b 100644 --- a/R/data.R +++ b/R/data.R @@ -2,10 +2,18 @@ #' #' A simulated dataset created using `simulate_data()` #' +#' @format A list as returned by [simulate_data()], with 100 samples at 30 +#' loci (15 with 5 alleles, 15 with 10) and a mean COI of 5: observed `data`, +#' `sample_ids`, `loci`, `is_missing`, and the simulated truth +#' (`allele_freqs`, `sample_cois`, `sample_relatedness`, `true_genotypes`). "simulated_data" #' MCMC results from using the packaged simulated data and calling `run_mcmc()` #' +#' @format A list as returned by [run_mcmc()] on [simulated_data]: a single +#' chain with 1000 burnin and 1000 sampling iterations, using 80 parallel +#' tempering replicas with adaptive temperatures. The `mcmc_demo` vignette +#' shows the call. "mcmc_results" #' Genetic and epidemiological data from Namibia diff --git a/R/mcmc.R b/R/mcmc.R index dc74c9f..c021f49 100644 --- a/R/mcmc.R +++ b/R/mcmc.R @@ -79,6 +79,41 @@ #' and whose `$diagnostics` field has Failure-locus counts and classification #' (Consistent failure cause vs Hard starting set). Inspect with #' `tryCatch(run_mcmc(...), error = function(e) e$diagnostics)`. +#' +#' @return List with elements: +#' \describe{ +#' \item{chains}{List with one entry per chain. Each holds the sampled +#' `coi`, `allele_freqs`, `eps_pos`, `eps_neg`, and `relatedness`, along +#' with log-likelihood traces, acceptance rates, and parallel tempering +#' diagnostics. Pass the whole result to the `summarize_*()` functions.} +#' \item{args}{The arguments the sampler was called with, including `data`.} +#' \item{seed}{Integer seed used for the run.} +#' \item{chain_seeds}{Integer seed used for each chain.} +#' \item{runtime}{Wall-clock time of the run as a `difftime`.} +#' } +#' +#' @examples +#' sim <- simulate_data( +#' mean_coi = 2, +#' num_samples = 10, +#' epsilon_pos = 0.01, +#' epsilon_neg = 0.1, +#' locus_freq_alphas = list(rep(1, 4), rep(1, 4), rep(1, 4)) +#' ) +#' +#' # A short run for illustration. Real analyses should use far more burnin +#' # and samples (see the defaults) and parallel tempering (pt_chains) to +#' # improve mixing. +#' res <- run_mcmc( +#' sim, +#' is_missing = sim$is_missing, +#' burnin = 100, +#' samples_per_chain = 100, +#' verbose = FALSE, +#' seed = 1 +#' ) +#' res$seed +#' summarize_coi(res) run_mcmc <- function(data, is_missing = FALSE, diff --git a/R/simulator.R b/R/simulator.R index 834e0e1..0968d92 100644 --- a/R/simulator.R +++ b/R/simulator.R @@ -7,6 +7,12 @@ #' #' @param n total number of draws #' @param alpha vector controlling the concentration of simplex +#' +#' @return Numeric matrix with `n` rows and `length(alpha)` columns; each row +#' is a draw from the simplex. +#' +#' @examples +#' rdirichlet(3, alpha = c(1, 1, 1)) rdirichlet <- function(n, alpha) { len_alpha <- length(alpha) d <- matrix(rgamma(len_alpha * n, alpha), ncol = len_alpha, byrow = TRUE) @@ -24,6 +30,12 @@ rdirichlet <- function(n, alpha) { #' #' @param alpha vector parameter controlling the Dirichlet distribution #' @param num_loci total number of loci to draw +#' +#' @return Numeric matrix with `length(alpha)` rows and `num_loci` columns; +#' each column is an allele frequency vector for one locus. +#' +#' @examples +#' simulate_allele_frequencies(alpha = c(1, 1, 1, 1), num_loci = 3) simulate_allele_frequencies <- function(alpha, num_loci) { dists <- rdirichlet(num_loci, alpha) sapply(seq_len(num_loci), function(x) { @@ -41,6 +53,12 @@ simulate_allele_frequencies <- function(alpha, num_loci) { #' #' @param num_samples the total number of biological samples to simulate #' @param mean_coi mean multiplicity of infection +#' +#' @return Integer vector of length `num_samples` with a COI of at least 1 +#' for each sample. +#' +#' @examples +#' simulate_sample_coi(num_samples = 10, mean_coi = 2) simulate_sample_coi <- function(num_samples, mean_coi) { qpois(runif(num_samples, dpois(0, mean_coi), 1), mean_coi) } @@ -59,7 +77,16 @@ simulate_sample_coi <- function(num_samples, mean_coi) { #' multinomial distribution #' @param internal_relatedness numeric 0-1 indicating the probability for a #' strain's allele to come from an existing lineage within host -#' @export +#' +#' @return List with one element per sample, each a 1-row integer matrix +#' counting the number of distinct strains carrying each allele. +#' +#' @examples +#' simulate_sample_genotype( +#' sample_cois = c(1, 2, 3), +#' locus_allele_dist = c(0.5, 0.3, 0.2), +#' internal_relatedness = c(0, 0, 0.5) +#' ) simulate_sample_genotype <- function(sample_cois, locus_allele_dist, internal_relatedness) { purrr::map2(sample_cois, internal_relatedness, function(coi, r) { related_strains <- rbinom(1, coi - 1, r) @@ -79,9 +106,18 @@ simulate_sample_genotype <- function(sample_cois, locus_allele_dist, internal_re #' #' @param alleles A numeric vector representing the number of strains #' contributing each allele -#' @param epsilon_pos expected number of false negatives -#' @param epsilon_neg expected number of false positives +#' @param epsilon_pos expected number of false positives +#' @param epsilon_neg expected number of false negatives #' @param missingness probability that the data is missing +#' +#' @return Binary numeric vector the same length as `alleles` indicating +#' which alleles were observed. All zeros if the observation is missing. +#' +#' @examples +#' simulate_observed_allele( +#' alleles = c(2, 0, 1), +#' epsilon_pos = 0.01, epsilon_neg = 0.1, missingness = 0 +#' ) simulate_observed_allele <- function(alleles, epsilon_pos, epsilon_neg, missingness) { # scale eps to the number of alleles so that given a fixed COI, there is a fixed # number of expected false positives or negatives across loci of varying @@ -118,6 +154,16 @@ simulate_observed_allele <- function(alleles, epsilon_pos, epsilon_neg, missingn #' @param epsilon_pos expected number of false positives #' @param epsilon_neg expected number of false negatives #' @param missingness probability of data being missing +#' +#' @return List the same length as `true_genotypes`, each element a binary +#' vector of observed alleles as returned by [simulate_observed_allele()]. +#' +#' @examples +#' true_genotypes <- list(c(1, 0, 1), c(0, 2, 0)) +#' simulate_observed_genotype( +#' true_genotypes, +#' epsilon_pos = 0.01, epsilon_neg = 0.1, missingness = 0 +#' ) simulate_observed_genotype <- function(true_genotypes, epsilon_pos, epsilon_neg, @@ -149,8 +195,21 @@ simulate_observed_genotype <- function(true_genotypes, #' @param internal_relatedness List of internal relatedness values to be used #' instead of simulating #' @param missingness probability of data being missing -#' @return Simulated data that is structured to go into the MCMC sampler +#' @return List structured for [run_mcmc()], with elements `data`, +#' `sample_ids`, `loci`, and `is_missing`, along with the simulated truth +#' (`allele_freqs`, `sample_cois`, `sample_relatedness`, `true_genotypes`) +#' and the `input` arguments. #' +#' @examples +#' sim <- simulate_data( +#' mean_coi = 2, +#' num_samples = 10, +#' epsilon_pos = 0.01, +#' epsilon_neg = 0.1, +#' locus_freq_alphas = list(rep(1, 4), rep(1, 4), rep(1, 4)) +#' ) +#' sim$sample_cois +#' str(sim$data, max.level = 1) simulate_data <- function(mean_coi = NULL, num_samples, epsilon_pos, diff --git a/R/summary.R b/R/summary.R index 51e030f..5f0821b 100644 --- a/R/summary.R +++ b/R/summary.R @@ -12,6 +12,13 @@ #' @param mcmc_results Result of calling run_mcmc() #' #' @param merge_chains boolean indicating that all chain results should be merged +#' +#' @return If `merge_chains` is `TRUE`, a named list with one numeric vector of +#' allele frequencies per locus. Otherwise a list with one such list per chain. +#' +#' @examples +#' med_freqs <- calculate_med_allele_freqs(mcmc_results) +#' med_freqs[[1]] calculate_med_allele_freqs <- function(mcmc_results, merge_chains = TRUE) { if (merge_chains) { chains <- mcmc_results$chains @@ -51,6 +58,11 @@ calculate_med_allele_freqs <- function(mcmc_results, merge_chains = TRUE) { #' element is a collection of observations across samples at a #' single genetic locus. #' @param offset Numeric offset -- n'th highest number of observed alleles +#' +#' @return Numeric vector with one naive COI estimate per sample. +#' +#' @examples +#' calculate_naive_coi_offset(simulated_data$data, offset = 2) calculate_naive_coi_offset <- function(data, offset) { num_alleles_by_locus <- lapply(data, function(locus) { lapply(locus, sum) @@ -79,6 +91,11 @@ calculate_naive_coi_offset <- function(data, offset) { #' @param data List of lists of numeric vectors, where each list #' element is a collection of observations across samples at a #' single genetic locus. +#' +#' @return Numeric vector with one naive COI estimate per sample. +#' +#' @examples +#' calculate_naive_coi(simulated_data$data) calculate_naive_coi <- function(data) { num_alleles_by_locus <- lapply(data, function(locus) { lapply(locus, sum) @@ -108,6 +125,12 @@ calculate_naive_coi <- function(data) { #' #' @param data List of lists of numeric vectors, where each list element #' is a collection of observations across samples at a single genetic locus +#' +#' @return List with one numeric vector of allele frequencies per locus. +#' +#' @examples +#' naive_freqs <- calculate_naive_allele_frequencies(simulated_data$data) +#' naive_freqs[[1]] calculate_naive_allele_frequencies <- function(data) { allele_freqs <- lapply(data, function(locus) { allele_counts <- purrr::reduce(locus, function(f1, f2) f1 + f2) @@ -120,6 +143,11 @@ calculate_naive_allele_frequencies <- function(data) { #' @export #' #' @param allele_freqs Simplex of allele frequencies +#' +#' @return Numeric scalar, the expected heterozygosity `1 - sum(p^2)`. +#' +#' @examples +#' calculate_he(c(0.5, 0.3, 0.2)) calculate_he <- function(allele_freqs) { if (any(is.na(allele_freqs))) { allele_freqs <- replace(allele_freqs, which(is.na(allele_freqs)), 0) @@ -146,6 +174,14 @@ calculate_he <- function(allele_freqs) { #' distribution to return #' @param naive_offset Offset used in calculate_naive_coi_offset #' @param merge_chains boolean indicating that all chain results should be merged +#' +#' @return Data frame with one row per sample (per chain if `merge_chains` is +#' `FALSE`) giving posterior quantiles and mean of COI, naive estimates, and +#' the posterior probability that the sample is polyclonal. +#' +#' @examples +#' coi_summary <- summarize_coi(mcmc_results) +#' head(coi_summary) summarize_coi <- function(mcmc_results, lower_quantile = .025, upper_quantile = .975, naive_offset = 2, merge_chains = TRUE) { naive_coi <- calculate_naive_coi(mcmc_results$args$data$data) @@ -221,6 +257,13 @@ summarize_coi <- function(mcmc_results, lower_quantile = .025, #' @param upper_quantile The upper quantile of the posterior #' distribution to return #' @param merge_chains boolean indicating that all chain results should be merged +#' +#' @return Data frame with one row per sample (per chain if `merge_chains` is +#' `FALSE`) giving posterior quantiles and mean of the false negative rate. +#' +#' @examples +#' eps_neg_summary <- summarize_epsilon_neg(mcmc_results) +#' head(eps_neg_summary) summarize_epsilon_neg <- function(mcmc_results, lower_quantile = .025, upper_quantile = .975, merge_chains = TRUE) { if (merge_chains) { epsilon_neg <- lapply(1:length(mcmc_results$args$data$sample_ids), function(x) c()) @@ -285,6 +328,13 @@ summarize_epsilon_neg <- function(mcmc_results, lower_quantile = .025, upper_qua #' @param upper_quantile The upper quantile of the posterior #' distribution to return #' @param merge_chains boolean indicating that all chain results should be merged +#' +#' @return Data frame with one row per sample (per chain if `merge_chains` is +#' `FALSE`) giving posterior quantiles and mean of the false positive rate. +#' +#' @examples +#' eps_pos_summary <- summarize_epsilon_pos(mcmc_results) +#' head(eps_pos_summary) summarize_epsilon_pos <- function(mcmc_results, lower_quantile = .025, upper_quantile = .975, merge_chains = TRUE) { if (merge_chains) { epsilon_pos <- lapply(1:length(mcmc_results$args$data$sample_ids), function(x) c()) @@ -350,6 +400,15 @@ summarize_epsilon_pos <- function(mcmc_results, lower_quantile = .025, upper_qua #' @param upper_quantile The upper quantile of the posterior distribution #' to return #' @param merge_chains boolean indicating that all chain results should be merged +#' +#' @return Data frame with one row per locus (per chain if `merge_chains` is +#' `FALSE`) giving posterior quantiles and mean of `fn` applied to the sampled +#' allele frequencies. +#' +#' @examples +#' # Posterior summary of the major allele frequency at each locus +#' major_allele <- summarize_allele_freq_fn(mcmc_results, fn = max) +#' head(major_allele) summarize_allele_freq_fn <- function(mcmc_results, fn, lower_quantile = .025, upper_quantile = .975, merge_chains = TRUE) { @@ -430,6 +489,13 @@ summarize_allele_freq_fn <- function(mcmc_results, fn, #' to return #' @param merge_chains Merge the results of multiple chains into a single #' summary +#' +#' @return Data frame with one row per locus (per chain if `merge_chains` is +#' `FALSE`) giving posterior quantiles and mean of expected heterozygosity. +#' +#' @examples +#' he_summary <- summarize_he(mcmc_results) +#' head(he_summary) summarize_he <- function(mcmc_results, lower_quantile = .025, upper_quantile = .975, merge_chains = TRUE) { @@ -468,6 +534,14 @@ names_or_idxs <- function(vec) { #' @param upper_quantile The upper quantile of the posterior distribution #' to return #' @param merge_chains boolean indicating that all chain results should be merged +#' +#' @return Data frame with one row per allele (per chain if `merge_chains` is +#' `FALSE`) giving posterior quantiles and mean of its frequency, along with +#' `locus` and `allele` identifiers. +#' +#' @examples +#' allele_freq_summary <- summarize_allele_freqs(mcmc_results) +#' head(allele_freq_summary) summarize_allele_freqs <- function(mcmc_results, lower_quantile = .025, upper_quantile = .975, @@ -587,6 +661,14 @@ summarize_allele_freqs <- function(mcmc_results, #' @param upper_quantile The upper quantile of the posterior #' distribution to return #' @param merge_chains boolean indicating that all chain results should be merged +#' +#' @return Data frame with one row per sample (per chain if `merge_chains` is +#' `FALSE`) giving posterior quantiles and mean of within-host relatedness, +#' computed over draws where COI is greater than 1. +#' +#' @examples +#' relatedness_summary <- summarize_relatedness(mcmc_results) +#' head(relatedness_summary) summarize_relatedness <- function(mcmc_results, lower_quantile = .025, upper_quantile = .975, merge_chains = TRUE) { if (merge_chains) { relatedness <- lapply(1:length(mcmc_results$args$data$sample_ids), function(x) c()) @@ -660,6 +742,14 @@ summarize_relatedness <- function(mcmc_results, lower_quantile = .025, upper_qua #' @param upper_quantile The upper quantile of the posterior #' distribution to return #' @param merge_chains boolean indicating that all chain results should be merged +#' +#' @return Data frame with one row per sample (per chain if `merge_chains` is +#' `FALSE`) giving posterior quantiles and mean of effective COI, +#' `(1 - relatedness) * (COI - 1) + 1`. +#' +#' @examples +#' effective_coi_summary <- summarize_effective_coi(mcmc_results) +#' head(effective_coi_summary) summarize_effective_coi <- function(mcmc_results, lower_quantile = .025, upper_quantile = .975, merge_chains = TRUE) { if (merge_chains) { effective_coi <- lapply(1:length(mcmc_results$args$data$sample_ids), function(x) c()) diff --git a/R/utils.R b/R/utils.R index 9892223..edf76c0 100644 --- a/R/utils.R +++ b/R/utils.R @@ -13,7 +13,27 @@ #' @param warn_uninformative boolean whether or not to print message when #' removing uninformative loci #' +#' @return List structured for [run_mcmc()] with elements `sample_ids`, +#' `data` (list of per-locus lists of binary allele vectors), `loci`, +#' `is_missing` (loci by samples logical matrix), and `uninformative_loci` +#' (loci removed for having a single allele). +#' #' @importFrom rlang .data +#' +#' @examples +#' df <- data.frame( +#' sample_id = c("S1", "S1", "S1", "S2", "S2"), +#' locus = c("L1", "L1", "L2", "L1", "L2"), +#' allele = c("A", "B", "A", "A", "B") +#' ) +#' dat <- load_long_form_data(df) +#' dat$loci +#' dat$data[["L1"]] +#' +#' # A subset of the bundled Namibia data +#' ids <- unique(namibia_data$sample_id)[1:20] +#' dat <- load_long_form_data(namibia_data[namibia_data$sample_id %in% ids, ]) +#' str(dat, max.level = 1) load_long_form_data <- function(df, warn_uninformative = TRUE) { uninformative_loci <- df |> dplyr::ungroup() |> @@ -98,7 +118,19 @@ load_long_form_data <- function(df, warn_uninformative = TRUE) { #' @param warn_uninformative boolean whether or not to print message when #' removing uninformative loci #' +#' @return List structured for [run_mcmc()]; see [load_long_form_data()]. +#' #' @importFrom rlang .data +#' +#' @examples +#' df <- data.frame( +#' sample_id = c("S1", "S2", "S3"), +#' L1 = c("A;B", "A", "B"), +#' L2 = c("A", "A;C", "C") +#' ) +#' dat <- load_delimited_data(df, sep = ";") +#' dat$sample_ids +#' dat$data[["L2"]] load_delimited_data <- function(data, sep = ";", warn_uninformative = TRUE) { df <- data |> tidyr::pivot_longer(-"sample_id", @@ -124,8 +156,11 @@ load_delimited_data <- function(data, sep = ";", warn_uninformative = TRUE) { #' @importFrom ggplot2 aes coord_cartesian geom_point geom_vline ggplot #' @importFrom rlang .data #' -#' @return list of ggplot objects +#' @return list of ggplot objects, one per chain #' +#' @examples +#' plots <- plot_chain_swaps(mcmc_results) +#' plots[[1]] plot_chain_swaps <- function(mcmc_results) { plots <- lapply(mcmc_results$chains, function(chain) { # swaps for a chain happen every 2 samples diff --git a/cran-comments.md b/cran-comments.md new file mode 100644 index 0000000..f71d5a7 --- /dev/null +++ b/cran-comments.md @@ -0,0 +1,23 @@ +## R CMD check results + +0 errors | 0 warnings | 0 notes + +* This is a new submission. + +## Possible NOTE: installed package size + +The installed size is about 6.5Mb, of which `data` is 4.3Mb. The `data` +directory holds a precomputed `mcmc_results` object (already `xz` compressed) +so that the vignette and the `summarize_*()` and `plot_*()` examples can run +in well under a second each rather than the minutes a real MCMC run takes. + +## Compiled code + +The package uses OpenMP when available and falls back to serial execution +when it is not (for example on CRAN's macOS builders). All OpenMP flags come +from `$(SHLIB_OPENMP_CXXFLAGS)`. Threading is controlled by the +`pt_num_threads` argument to `run_mcmc()`, which defaults to 1; tests and +examples never use more than 2 cores. + +`src/include/spline/spline.h` is a bundled header-only cubic spline library +by Tino Kluge, licensed GPL-2 or later and listed as `cph` in `Authors@R`. diff --git a/man/calculate_he.Rd b/man/calculate_he.Rd index 5d18b62..d28e931 100644 --- a/man/calculate_he.Rd +++ b/man/calculate_he.Rd @@ -9,6 +9,12 @@ calculate_he(allele_freqs) \arguments{ \item{allele_freqs}{Simplex of allele frequencies} } +\value{ +Numeric scalar, the expected heterozygosity \code{1 - sum(p^2)}. +} \description{ Calculate the expected heterozygosity from allele frequencies } +\examples{ +calculate_he(c(0.5, 0.3, 0.2)) +} diff --git a/man/calculate_med_allele_freqs.Rd b/man/calculate_med_allele_freqs.Rd index 45b173d..80b2bfd 100644 --- a/man/calculate_med_allele_freqs.Rd +++ b/man/calculate_med_allele_freqs.Rd @@ -12,6 +12,10 @@ calculate_med_allele_freqs(mcmc_results, merge_chains = TRUE) \item{merge_chains}{boolean indicating that all chain results should be merged} } +\value{ +If \code{merge_chains} is \code{TRUE}, a named list with one numeric vector of +allele frequencies per locus. Otherwise a list with one such list per chain. +} \description{ Calculate the geometric median of the posterior distribution of allele frequencies @@ -20,3 +24,7 @@ frequencies Returns the geometric median of the posterior distribution, defined as the point minimizing the L2 distance from each sampled point. } +\examples{ +med_freqs <- calculate_med_allele_freqs(mcmc_results) +med_freqs[[1]] +} diff --git a/man/calculate_naive_allele_frequencies.Rd b/man/calculate_naive_allele_frequencies.Rd index 79db6ef..093d519 100644 --- a/man/calculate_naive_allele_frequencies.Rd +++ b/man/calculate_naive_allele_frequencies.Rd @@ -10,6 +10,9 @@ calculate_naive_allele_frequencies(data) \item{data}{List of lists of numeric vectors, where each list element is a collection of observations across samples at a single genetic locus} } +\value{ +List with one numeric vector of allele frequencies per locus. +} \description{ Calculate naive allele frequencies } @@ -17,3 +20,7 @@ Calculate naive allele frequencies Estimate naive allele frequencies from the empirical distribution of alleles } +\examples{ +naive_freqs <- calculate_naive_allele_frequencies(simulated_data$data) +naive_freqs[[1]] +} diff --git a/man/calculate_naive_coi.Rd b/man/calculate_naive_coi.Rd index 9fa7ef2..e89f9cd 100644 --- a/man/calculate_naive_coi.Rd +++ b/man/calculate_naive_coi.Rd @@ -11,6 +11,9 @@ calculate_naive_coi(data) element is a collection of observations across samples at a single genetic locus.} } +\value{ +Numeric vector with one naive COI estimate per sample. +} \description{ Calculate naive COI } @@ -18,3 +21,6 @@ Calculate naive COI Estimates the complexity of infection using a naive approach that chooses the highest number of observed alleles. } +\examples{ +calculate_naive_coi(simulated_data$data) +} diff --git a/man/calculate_naive_coi_offset.Rd b/man/calculate_naive_coi_offset.Rd index ab31a50..087810c 100644 --- a/man/calculate_naive_coi_offset.Rd +++ b/man/calculate_naive_coi_offset.Rd @@ -13,6 +13,9 @@ single genetic locus.} \item{offset}{Numeric offset -- n'th highest number of observed alleles} } +\value{ +Numeric vector with one naive COI estimate per sample. +} \description{ Calculate naive COI offset } @@ -20,3 +23,6 @@ Calculate naive COI offset Estimates the complexity of infection using a naive approach that chooses the n'th highest number of observed alleles. } +\examples{ +calculate_naive_coi_offset(simulated_data$data, offset = 2) +} diff --git a/man/load_delimited_data.Rd b/man/load_delimited_data.Rd index 2363a2f..d9e4bf2 100644 --- a/man/load_delimited_data.Rd +++ b/man/load_delimited_data.Rd @@ -14,6 +14,9 @@ load_delimited_data(data, sep = ";", warn_uninformative = TRUE) \item{warn_uninformative}{boolean whether or not to print message when removing uninformative loci} } +\value{ +List structured for \code{\link[=run_mcmc]{run_mcmc()}}; see \code{\link[=load_long_form_data]{load_long_form_data()}}. +} \description{ Load delimited data } @@ -24,3 +27,13 @@ representing the observed alleles at that locus for that sample. Returned data contains vectors \code{sample_ids} and \code{loci} that are ordered as the results will be ordered from running the MCMC algorithm. } +\examples{ +df <- data.frame( + sample_id = c("S1", "S2", "S3"), + L1 = c("A;B", "A", "B"), + L2 = c("A", "A;C", "C") +) +dat <- load_delimited_data(df, sep = ";") +dat$sample_ids +dat$data[["L2"]] +} diff --git a/man/load_long_form_data.Rd b/man/load_long_form_data.Rd index 3c2eedc..3d79545 100644 --- a/man/load_long_form_data.Rd +++ b/man/load_long_form_data.Rd @@ -14,6 +14,12 @@ locus for a given sample.} \item{warn_uninformative}{boolean whether or not to print message when removing uninformative loci} } +\value{ +List structured for \code{\link[=run_mcmc]{run_mcmc()}} with elements \code{sample_ids}, +\code{data} (list of per-locus lists of binary allele vectors), \code{loci}, +\code{is_missing} (loci by samples logical matrix), and \code{uninformative_loci} +(loci removed for having a single allele). +} \description{ Load long form data } @@ -23,3 +29,18 @@ Long form data is a data frame with vectors \code{sample_ids} and \code{loci} that are ordered as the results will be ordered from running the MCMC algorithm. } +\examples{ +df <- data.frame( + sample_id = c("S1", "S1", "S1", "S2", "S2"), + locus = c("L1", "L1", "L2", "L1", "L2"), + allele = c("A", "B", "A", "A", "B") +) +dat <- load_long_form_data(df) +dat$loci +dat$data[["L1"]] + +# A subset of the bundled Namibia data +ids <- unique(namibia_data$sample_id)[1:20] +dat <- load_long_form_data(namibia_data[namibia_data$sample_id \%in\% ids, ]) +str(dat, max.level = 1) +} diff --git a/man/mcmc_results.Rd b/man/mcmc_results.Rd index fb35e16..1a69c36 100644 --- a/man/mcmc_results.Rd +++ b/man/mcmc_results.Rd @@ -5,7 +5,10 @@ \alias{mcmc_results} \title{MCMC results from using the packaged simulated data and calling \code{run_mcmc()}} \format{ -An object of class \code{list} of length 3. +A list as returned by \code{\link[=run_mcmc]{run_mcmc()}} on \link{simulated_data}: a single +chain with 1000 burnin and 1000 sampling iterations, using 80 parallel +tempering replicas with adaptive temperatures. The \code{mcmc_demo} vignette +shows the call. } \usage{ mcmc_results diff --git a/man/plot_chain_swaps.Rd b/man/plot_chain_swaps.Rd index c416449..c56a5e1 100644 --- a/man/plot_chain_swaps.Rd +++ b/man/plot_chain_swaps.Rd @@ -10,7 +10,7 @@ plot_chain_swaps(mcmc_results) \item{mcmc_results}{list of results from \code{run_mcmc}} } \value{ -list of ggplot objects +list of ggplot objects, one per chain } \description{ Plot chain swap acceptance rates @@ -20,3 +20,7 @@ Plot the swap acceptance rates for each chain. The x-axis is the temperature, and the y-axis is the swap acceptance rate. The dashed lines indicate the temperatures used for parallel tempering. } +\examples{ +plots <- plot_chain_swaps(mcmc_results) +plots[[1]] +} diff --git a/man/rdirichlet.Rd b/man/rdirichlet.Rd index 11949a3..00cd4fb 100644 --- a/man/rdirichlet.Rd +++ b/man/rdirichlet.Rd @@ -11,9 +11,16 @@ rdirichlet(n, alpha) \item{alpha}{vector controlling the concentration of simplex} } +\value{ +Numeric matrix with \code{n} rows and \code{length(alpha)} columns; each row +is a draw from the simplex. +} \description{ Dirichlet distribution } \details{ Implementation of random sampling from a Dirichlet distribution } +\examples{ +rdirichlet(3, alpha = c(1, 1, 1)) +} diff --git a/man/run_mcmc.Rd b/man/run_mcmc.Rd index 8f4eb03..603ca24 100644 --- a/man/run_mcmc.Rd +++ b/man/run_mcmc.Rd @@ -139,6 +139,19 @@ and printed if the run errors; per-chain offsets are \verb{$chain_seeds}. Replay needs the same sampler settings; \code{num_cores} does not affect the draws.} } +\value{ +List with elements: +\describe{ +\item{chains}{List with one entry per chain. Each holds the sampled +\code{coi}, \code{allele_freqs}, \code{eps_pos}, \code{eps_neg}, and \code{relatedness}, along +with log-likelihood traces, acceptance rates, and parallel tempering +diagnostics. Pass the whole result to the \verb{summarize_*()} functions.} +\item{args}{The arguments the sampler was called with, including \code{data}.} +\item{seed}{Integer seed used for the run.} +\item{chain_seeds}{Integer seed used for each chain.} +\item{runtime}{Wall-clock time of the run as a \code{difftime}.} +} +} \description{ Sample from the target distribution using MCMC } @@ -152,3 +165,26 @@ and whose \verb{$diagnostics} field has Failure-locus counts and classification \code{tryCatch(run_mcmc(...), error = function(e) e$diagnostics)}. } +\examples{ +sim <- simulate_data( + mean_coi = 2, + num_samples = 10, + epsilon_pos = 0.01, + epsilon_neg = 0.1, + locus_freq_alphas = list(rep(1, 4), rep(1, 4), rep(1, 4)) +) + +# A short run for illustration. Real analyses should use far more burnin +# and samples (see the defaults) and parallel tempering (pt_chains) to +# improve mixing. +res <- run_mcmc( + sim, + is_missing = sim$is_missing, + burnin = 100, + samples_per_chain = 100, + verbose = FALSE, + seed = 1 +) +res$seed +summarize_coi(res) +} diff --git a/man/simulate_allele_frequencies.Rd b/man/simulate_allele_frequencies.Rd index 548a8fd..79fe54a 100644 --- a/man/simulate_allele_frequencies.Rd +++ b/man/simulate_allele_frequencies.Rd @@ -11,6 +11,10 @@ simulate_allele_frequencies(alpha, num_loci) \item{num_loci}{total number of loci to draw} } +\value{ +Numeric matrix with \code{length(alpha)} rows and \code{num_loci} columns; +each column is an allele frequency vector for one locus. +} \description{ Simulate allele frequencies } @@ -18,3 +22,6 @@ Simulate allele frequencies Simulate allele frequency vectors as a draw from a Dirichlet distribution } +\examples{ +simulate_allele_frequencies(alpha = c(1, 1, 1, 1), num_loci = 3) +} diff --git a/man/simulate_data.Rd b/man/simulate_data.Rd index 6ffa0a7..bf258ad 100644 --- a/man/simulate_data.Rd +++ b/man/simulate_data.Rd @@ -47,8 +47,22 @@ instead of simulating} \item{missingness}{probability of data being missing} } \value{ -Simulated data that is structured to go into the MCMC sampler +List structured for \code{\link[=run_mcmc]{run_mcmc()}}, with elements \code{data}, +\code{sample_ids}, \code{loci}, and \code{is_missing}, along with the simulated truth +(\code{allele_freqs}, \code{sample_cois}, \code{sample_relatedness}, \code{true_genotypes}) +and the \code{input} arguments. } \description{ Simulate data generated according to the assumed model } +\examples{ +sim <- simulate_data( + mean_coi = 2, + num_samples = 10, + epsilon_pos = 0.01, + epsilon_neg = 0.1, + locus_freq_alphas = list(rep(1, 4), rep(1, 4), rep(1, 4)) +) +sim$sample_cois +str(sim$data, max.level = 1) +} diff --git a/man/simulate_observed_allele.Rd b/man/simulate_observed_allele.Rd index 34618c7..1ee845c 100644 --- a/man/simulate_observed_allele.Rd +++ b/man/simulate_observed_allele.Rd @@ -10,12 +10,16 @@ simulate_observed_allele(alleles, epsilon_pos, epsilon_neg, missingness) \item{alleles}{A numeric vector representing the number of strains contributing each allele} -\item{epsilon_pos}{expected number of false negatives} +\item{epsilon_pos}{expected number of false positives} -\item{epsilon_neg}{expected number of false positives} +\item{epsilon_neg}{expected number of false negatives} \item{missingness}{probability that the data is missing} } +\value{ +Binary numeric vector the same length as \code{alleles} indicating +which alleles were observed. All zeros if the observation is missing. +} \description{ Simulates the observation process } @@ -24,3 +28,9 @@ Takes a numeric value representing the number of strains contributing an allele and returns a binary vector indicating the presence or absence of the allele. } +\examples{ +simulate_observed_allele( + alleles = c(2, 0, 1), + epsilon_pos = 0.01, epsilon_neg = 0.1, missingness = 0 +) +} diff --git a/man/simulate_observed_genotype.Rd b/man/simulate_observed_genotype.Rd index 01ba35f..d90511c 100644 --- a/man/simulate_observed_genotype.Rd +++ b/man/simulate_observed_genotype.Rd @@ -21,6 +21,10 @@ to sim_observed_allele} \item{missingness}{probability of data being missing} } +\value{ +List the same length as \code{true_genotypes}, each element a binary +vector of observed alleles as returned by \code{\link[=simulate_observed_allele]{simulate_observed_allele()}}. +} \description{ Simulate observed genotypes } @@ -28,3 +32,10 @@ Simulate observed genotypes Simulate the observation process across a list of observation vectors } +\examples{ +true_genotypes <- list(c(1, 0, 1), c(0, 2, 0)) +simulate_observed_genotype( + true_genotypes, + epsilon_pos = 0.01, epsilon_neg = 0.1, missingness = 0 +) +} diff --git a/man/simulate_sample_coi.Rd b/man/simulate_sample_coi.Rd index 63acc2d..efc6660 100644 --- a/man/simulate_sample_coi.Rd +++ b/man/simulate_sample_coi.Rd @@ -11,9 +11,16 @@ simulate_sample_coi(num_samples, mean_coi) \item{mean_coi}{mean multiplicity of infection} } +\value{ +Integer vector of length \code{num_samples} with a COI of at least 1 +for each sample. +} \description{ Simulate sample COI } \details{ Simulate sample COIs from a zero-truncated Poisson distribution } +\examples{ +simulate_sample_coi(num_samples = 10, mean_coi = 2) +} diff --git a/man/simulate_sample_genotype.Rd b/man/simulate_sample_genotype.Rd index 89face2..72b9c6a 100644 --- a/man/simulate_sample_genotype.Rd +++ b/man/simulate_sample_genotype.Rd @@ -16,6 +16,10 @@ multinomial distribution} \item{internal_relatedness}{numeric 0-1 indicating the probability for a strain's allele to come from an existing lineage within host} } +\value{ +List with one element per sample, each a 1-row integer matrix +counting the number of distinct strains carrying each allele. +} \description{ Simulate sample genotype } @@ -23,3 +27,10 @@ Simulate sample genotype Simulates sampling the genetics at a single locus given an allele frequency distribution and a vector of sample COIs } +\examples{ +simulate_sample_genotype( + sample_cois = c(1, 2, 3), + locus_allele_dist = c(0.5, 0.3, 0.2), + internal_relatedness = c(0, 0, 0.5) +) +} diff --git a/man/simulated_data.Rd b/man/simulated_data.Rd index 59e75d4..9006486 100644 --- a/man/simulated_data.Rd +++ b/man/simulated_data.Rd @@ -5,7 +5,10 @@ \alias{simulated_data} \title{Simulated genotyping data} \format{ -An object of class \code{list} of length 9. +A list as returned by \code{\link[=simulate_data]{simulate_data()}}, with 100 samples at 30 +loci (15 with 5 alleles, 15 with 10) and a mean COI of 5: observed \code{data}, +\code{sample_ids}, \code{loci}, \code{is_missing}, and the simulated truth +(\code{allele_freqs}, \code{sample_cois}, \code{sample_relatedness}, \code{true_genotypes}). } \usage{ simulated_data diff --git a/man/summarize_allele_freq_fn.Rd b/man/summarize_allele_freq_fn.Rd index e2550e9..750df10 100644 --- a/man/summarize_allele_freq_fn.Rd +++ b/man/summarize_allele_freq_fn.Rd @@ -26,6 +26,11 @@ to return} \item{merge_chains}{boolean indicating that all chain results should be merged} } +\value{ +Data frame with one row per locus (per chain if \code{merge_chains} is +\code{FALSE}) giving posterior quantiles and mean of \code{fn} applied to the sampled +allele frequencies. +} \description{ Summarize Function of Allele Frequencies } @@ -33,3 +38,8 @@ Summarize Function of Allele Frequencies General function to summarize the posterior distribution of functions of the sampled allele frequencies } +\examples{ +# Posterior summary of the major allele frequency at each locus +major_allele <- summarize_allele_freq_fn(mcmc_results, fn = max) +head(major_allele) +} diff --git a/man/summarize_allele_freqs.Rd b/man/summarize_allele_freqs.Rd index 678aa98..8d0db63 100644 --- a/man/summarize_allele_freqs.Rd +++ b/man/summarize_allele_freqs.Rd @@ -22,6 +22,11 @@ to return} \item{merge_chains}{boolean indicating that all chain results should be merged} } +\value{ +Data frame with one row per allele (per chain if \code{merge_chains} is +\code{FALSE}) giving posterior quantiles and mean of its frequency, along with +\code{locus} and \code{allele} identifiers. +} \description{ Summarize allele frequencies } @@ -29,3 +34,7 @@ Summarize allele frequencies Summarize individual allele frequencies from the posterior distribution of sampled allele frequencies } +\examples{ +allele_freq_summary <- summarize_allele_freqs(mcmc_results) +head(allele_freq_summary) +} diff --git a/man/summarize_coi.Rd b/man/summarize_coi.Rd index 49b0452..cfcddaf 100644 --- a/man/summarize_coi.Rd +++ b/man/summarize_coi.Rd @@ -25,6 +25,11 @@ distribution to return} \item{merge_chains}{boolean indicating that all chain results should be merged} } +\value{ +Data frame with one row per sample (per chain if \code{merge_chains} is +\code{FALSE}) giving posterior quantiles and mean of COI, naive estimates, and +the posterior probability that the sample is polyclonal. +} \description{ Summarize COI } @@ -34,3 +39,7 @@ a dataframe that contains summaries of the posterior distribution of COI for each biological sample, as well as naive estimates of COI. } +\examples{ +coi_summary <- summarize_coi(mcmc_results) +head(coi_summary) +} diff --git a/man/summarize_effective_coi.Rd b/man/summarize_effective_coi.Rd index 518afc0..0a5c28d 100644 --- a/man/summarize_effective_coi.Rd +++ b/man/summarize_effective_coi.Rd @@ -22,6 +22,11 @@ distribution to return} \item{merge_chains}{boolean indicating that all chain results should be merged} } +\value{ +Data frame with one row per sample (per chain if \code{merge_chains} is +\code{FALSE}) giving posterior quantiles and mean of effective COI, +\code{(1 - relatedness) * (COI - 1) + 1}. +} \description{ Summarize effective COI } @@ -30,3 +35,7 @@ Summarize effective COI from MCMC. Returns a dataframe that contains summaries of the posterior distribution of effective COI for each biological sample. } +\examples{ +effective_coi_summary <- summarize_effective_coi(mcmc_results) +head(effective_coi_summary) +} diff --git a/man/summarize_epsilon_neg.Rd b/man/summarize_epsilon_neg.Rd index 8314646..3f6ca48 100644 --- a/man/summarize_epsilon_neg.Rd +++ b/man/summarize_epsilon_neg.Rd @@ -22,6 +22,10 @@ distribution to return} \item{merge_chains}{boolean indicating that all chain results should be merged} } +\value{ +Data frame with one row per sample (per chain if \code{merge_chains} is +\code{FALSE}) giving posterior quantiles and mean of the false negative rate. +} \description{ Summarize epsilon_neg } @@ -30,3 +34,7 @@ Summarize epsilon negative results from MCMC. Returns a dataframe that contains summaries of the posterior distribution of epsilon negative for each biological sample. } +\examples{ +eps_neg_summary <- summarize_epsilon_neg(mcmc_results) +head(eps_neg_summary) +} diff --git a/man/summarize_epsilon_pos.Rd b/man/summarize_epsilon_pos.Rd index 16a69f7..40b9f20 100644 --- a/man/summarize_epsilon_pos.Rd +++ b/man/summarize_epsilon_pos.Rd @@ -22,6 +22,10 @@ distribution to return} \item{merge_chains}{boolean indicating that all chain results should be merged} } +\value{ +Data frame with one row per sample (per chain if \code{merge_chains} is +\code{FALSE}) giving posterior quantiles and mean of the false positive rate. +} \description{ Summarize epsilon_pos } @@ -30,3 +34,7 @@ Summarize epsilon positive results from MCMC. Returns a dataframe that contains summaries of the posterior distribution of epsilon positive for each biological sample. } +\examples{ +eps_pos_summary <- summarize_epsilon_pos(mcmc_results) +head(eps_pos_summary) +} diff --git a/man/summarize_he.Rd b/man/summarize_he.Rd index 6adda00..d3ee87a 100644 --- a/man/summarize_he.Rd +++ b/man/summarize_he.Rd @@ -23,6 +23,10 @@ to return} \item{merge_chains}{Merge the results of multiple chains into a single summary} } +\value{ +Data frame with one row per locus (per chain if \code{merge_chains} is +\code{FALSE}) giving posterior quantiles and mean of expected heterozygosity. +} \description{ Summarize locus heterozygosity } @@ -30,3 +34,7 @@ Summarize locus heterozygosity Summarize locus heterozygosity from the posterior distribution of sampled allele frequencies. } +\examples{ +he_summary <- summarize_he(mcmc_results) +head(he_summary) +} diff --git a/man/summarize_relatedness.Rd b/man/summarize_relatedness.Rd index 1dcbcfa..e59c0a7 100644 --- a/man/summarize_relatedness.Rd +++ b/man/summarize_relatedness.Rd @@ -22,6 +22,11 @@ distribution to return} \item{merge_chains}{boolean indicating that all chain results should be merged} } +\value{ +Data frame with one row per sample (per chain if \code{merge_chains} is +\code{FALSE}) giving posterior quantiles and mean of within-host relatedness, +computed over draws where COI is greater than 1. +} \description{ Summarize relatedness } @@ -30,3 +35,7 @@ Summarize relatedness results from MCMC. Returns a dataframe that contains summaries of the posterior distribution of relatedness for each biological sample. } +\examples{ +relatedness_summary <- summarize_relatedness(mcmc_results) +head(relatedness_summary) +} From 6688582710061bd5a8d215f0a60601de3d0d9702 Mon Sep 17 00:00:00 2001 From: Maxwell Murphy Date: Tue, 15 Sep 2026 10:12:04 -0700 Subject: [PATCH 3/3] Replace bundled spline.h with a monotone cubic Hermite interpolator adapt_temp() now interpolates temperature as a function of cumulative communication barrier (Fritsch-Carlson / PCHIP slopes) and reads the new ladder off an even barrier grid, so there is no root solving, no monotonicity post-fix, and no band-matrix LU. Adjacent chains with zero barrier collapse into a single knot; if no barrier has been observed the update is skipped. Removes src/include/spline (GPL-2 third-party, locally modified without a change notice) and the cph entry for it. Adds testthat coverage for the adapted ladder. --- DESCRIPTION | 6 +- NEWS.md | 5 + cran-comments.md | 3 - src/include/spline/LICENSE | 340 --------- src/include/spline/README.md | 82 --- src/include/spline/src/spline.h | 1110 ------------------------------ src/mcmc.cpp | 77 ++- src/monotone_interpolator.h | 102 +++ tests/testthat/test-adapt-temp.R | 56 ++ 9 files changed, 206 insertions(+), 1575 deletions(-) delete mode 100644 src/include/spline/LICENSE delete mode 100644 src/include/spline/README.md delete mode 100644 src/include/spline/src/spline.h create mode 100644 src/monotone_interpolator.h create mode 100644 tests/testthat/test-adapt-temp.R diff --git a/DESCRIPTION b/DESCRIPTION index b624c44..99f0a83 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -11,11 +11,7 @@ Authors@R: c( family = "Greenhouse", role = c("aut", "ths"), email = "bryan.greenhouse@ucsf.edu", - comment = c(ORCID = "0000-0003-0287-9111")), - person(given = "Tino", - family = "Kluge", - role = "cph", - comment = "Author of the bundled cubic spline library (src/include/spline)")) + comment = c(ORCID = "0000-0003-0287-9111"))) Description: A Markov Chain Monte Carlo (MCMC) based approach to Bayesian estimation of individual level multiplicity of infection, within host relatedness, and population allele frequencies from polyallelic genetic diff --git a/NEWS.md b/NEWS.md index f2ebc1a..956dd05 100644 --- a/NEWS.md +++ b/NEWS.md @@ -16,6 +16,11 @@ `simulate_observed_allele()`. * `parallel` is now declared in `Imports`; `markdown` was dropped from `Suggests` as nothing used it. +* Parallel tempering temperature adaptation now interpolates temperature as a + function of cumulative communication barrier with a monotone cubic Hermite + scheme (`src/monotone_interpolator.h`), replacing the bundled third-party + `spline.h` and its root solving. Adapted ladders differ slightly from + earlier versions; if no barrier has been observed yet the update is skipped. # moire 3.6.0 diff --git a/cran-comments.md b/cran-comments.md index f71d5a7..be44733 100644 --- a/cran-comments.md +++ b/cran-comments.md @@ -18,6 +18,3 @@ when it is not (for example on CRAN's macOS builders). All OpenMP flags come from `$(SHLIB_OPENMP_CXXFLAGS)`. Threading is controlled by the `pt_num_threads` argument to `run_mcmc()`, which defaults to 1; tests and examples never use more than 2 cores. - -`src/include/spline/spline.h` is a bundled header-only cubic spline library -by Tino Kluge, licensed GPL-2 or later and listed as `cph` in `Authors@R`. diff --git a/src/include/spline/LICENSE b/src/include/spline/LICENSE deleted file mode 100644 index 8cdb845..0000000 --- a/src/include/spline/LICENSE +++ /dev/null @@ -1,340 +0,0 @@ - GNU GENERAL PUBLIC LICENSE - Version 2, June 1991 - - Copyright (C) 1989, 1991 Free Software Foundation, Inc., - 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA - Everyone is permitted to copy and distribute verbatim copies - of this license document, but changing it is not allowed. - - Preamble - - The licenses for most software are designed to take away your -freedom to share and change it. By contrast, the GNU General Public -License is intended to guarantee your freedom to share and change free -software--to make sure the software is free for all its users. This -General Public License applies to most of the Free Software -Foundation's software and to any other program whose authors commit to -using it. (Some other Free Software Foundation software is covered by -the GNU Lesser General Public License instead.) You can apply it to -your programs, too. - - When we speak of free software, we are referring to freedom, not -price. Our General Public Licenses are designed to make sure that you -have the freedom to distribute copies of free software (and charge for -this service if you wish), that you receive source code or can get it -if you want it, that you can change the software or use pieces of it -in new free programs; and that you know you can do these things. - - To protect your rights, we need to make restrictions that forbid -anyone to deny you these rights or to ask you to surrender the rights. -These restrictions translate to certain responsibilities for you if you -distribute copies of the software, or if you modify it. - - For example, if you distribute copies of such a program, whether -gratis or for a fee, you must give the recipients all the rights that -you have. You must make sure that they, too, receive or can get the -source code. And you must show them these terms so they know their -rights. - - We protect your rights with two steps: (1) copyright the software, and -(2) offer you this license which gives you legal permission to copy, -distribute and/or modify the software. - - Also, for each author's protection and ours, we want to make certain -that everyone understands that there is no warranty for this free -software. If the software is modified by someone else and passed on, we -want its recipients to know that what they have is not the original, so -that any problems introduced by others will not reflect on the original -authors' reputations. - - Finally, any free program is threatened constantly by software -patents. We wish to avoid the danger that redistributors of a free -program will individually obtain patent licenses, in effect making the -program proprietary. To prevent this, we have made it clear that any -patent must be licensed for everyone's free use or not licensed at all. - - The precise terms and conditions for copying, distribution and -modification follow. - - GNU GENERAL PUBLIC LICENSE - TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION - - 0. This License applies to any program or other work which contains -a notice placed by the copyright holder saying it may be distributed -under the terms of this General Public License. The "Program", below, -refers to any such program or work, and a "work based on the Program" -means either the Program or any derivative work under copyright law: -that is to say, a work containing the Program or a portion of it, -either verbatim or with modifications and/or translated into another -language. (Hereinafter, translation is included without limitation in -the term "modification".) Each licensee is addressed as "you". - -Activities other than copying, distribution and modification are not -covered by this License; they are outside its scope. The act of -running the Program is not restricted, and the output from the Program -is covered only if its contents constitute a work based on the -Program (independent of having been made by running the Program). -Whether that is true depends on what the Program does. - - 1. You may copy and distribute verbatim copies of the Program's -source code as you receive it, in any medium, provided that you -conspicuously and appropriately publish on each copy an appropriate -copyright notice and disclaimer of warranty; keep intact all the -notices that refer to this License and to the absence of any warranty; -and give any other recipients of the Program a copy of this License -along with the Program. - -You may charge a fee for the physical act of transferring a copy, and -you may at your option offer warranty protection in exchange for a fee. - - 2. You may modify your copy or copies of the Program or any portion -of it, thus forming a work based on the Program, and copy and -distribute such modifications or work under the terms of Section 1 -above, provided that you also meet all of these conditions: - - a) You must cause the modified files to carry prominent notices - stating that you changed the files and the date of any change. - - b) You must cause any work that you distribute or publish, that in - whole or in part contains or is derived from the Program or any - part thereof, to be licensed as a whole at no charge to all third - parties under the terms of this License. - - c) If the modified program normally reads commands interactively - when run, you must cause it, when started running for such - interactive use in the most ordinary way, to print or display an - announcement including an appropriate copyright notice and a - notice that there is no warranty (or else, saying that you provide - a warranty) and that users may redistribute the program under - these conditions, and telling the user how to view a copy of this - License. (Exception: if the Program itself is interactive but - does not normally print such an announcement, your work based on - the Program is not required to print an announcement.) - -These requirements apply to the modified work as a whole. If -identifiable sections of that work are not derived from the Program, -and can be reasonably considered independent and separate works in -themselves, then this License, and its terms, do not apply to those -sections when you distribute them as separate works. But when you -distribute the same sections as part of a whole which is a work based -on the Program, the distribution of the whole must be on the terms of -this License, whose permissions for other licensees extend to the -entire whole, and thus to each and every part regardless of who wrote it. - -Thus, it is not the intent of this section to claim rights or contest -your rights to work written entirely by you; rather, the intent is to -exercise the right to control the distribution of derivative or -collective works based on the Program. - -In addition, mere aggregation of another work not based on the Program -with the Program (or with a work based on the Program) on a volume of -a storage or distribution medium does not bring the other work under -the scope of this License. - - 3. You may copy and distribute the Program (or a work based on it, -under Section 2) in object code or executable form under the terms of -Sections 1 and 2 above provided that you also do one of the following: - - a) Accompany it with the complete corresponding machine-readable - source code, which must be distributed under the terms of Sections - 1 and 2 above on a medium customarily used for software interchange; or, - - b) Accompany it with a written offer, valid for at least three - years, to give any third party, for a charge no more than your - cost of physically performing source distribution, a complete - machine-readable copy of the corresponding source code, to be - distributed under the terms of Sections 1 and 2 above on a medium - customarily used for software interchange; or, - - c) Accompany it with the information you received as to the offer - to distribute corresponding source code. (This alternative is - allowed only for noncommercial distribution and only if you - received the program in object code or executable form with such - an offer, in accord with Subsection b above.) - -The source code for a work means the preferred form of the work for -making modifications to it. For an executable work, complete source -code means all the source code for all modules it contains, plus any -associated interface definition files, plus the scripts used to -control compilation and installation of the executable. However, as a -special exception, the source code distributed need not include -anything that is normally distributed (in either source or binary -form) with the major components (compiler, kernel, and so on) of the -operating system on which the executable runs, unless that component -itself accompanies the executable. - -If distribution of executable or object code is made by offering -access to copy from a designated place, then offering equivalent -access to copy the source code from the same place counts as -distribution of the source code, even though third parties are not -compelled to copy the source along with the object code. - - 4. You may not copy, modify, sublicense, or distribute the Program -except as expressly provided under this License. Any attempt -otherwise to copy, modify, sublicense or distribute the Program is -void, and will automatically terminate your rights under this License. -However, parties who have received copies, or rights, from you under -this License will not have their licenses terminated so long as such -parties remain in full compliance. - - 5. You are not required to accept this License, since you have not -signed it. However, nothing else grants you permission to modify or -distribute the Program or its derivative works. These actions are -prohibited by law if you do not accept this License. Therefore, by -modifying or distributing the Program (or any work based on the -Program), you indicate your acceptance of this License to do so, and -all its terms and conditions for copying, distributing or modifying -the Program or works based on it. - - 6. Each time you redistribute the Program (or any work based on the -Program), the recipient automatically receives a license from the -original licensor to copy, distribute or modify the Program subject to -these terms and conditions. You may not impose any further -restrictions on the recipients' exercise of the rights granted herein. -You are not responsible for enforcing compliance by third parties to -this License. - - 7. If, as a consequence of a court judgment or allegation of patent -infringement or for any other reason (not limited to patent issues), -conditions are imposed on you (whether by court order, agreement or -otherwise) that contradict the conditions of this License, they do not -excuse you from the conditions of this License. If you cannot -distribute so as to satisfy simultaneously your obligations under this -License and any other pertinent obligations, then as a consequence you -may not distribute the Program at all. For example, if a patent -license would not permit royalty-free redistribution of the Program by -all those who receive copies directly or indirectly through you, then -the only way you could satisfy both it and this License would be to -refrain entirely from distribution of the Program. - -If any portion of this section is held invalid or unenforceable under -any particular circumstance, the balance of the section is intended to -apply and the section as a whole is intended to apply in other -circumstances. - -It is not the purpose of this section to induce you to infringe any -patents or other property right claims or to contest validity of any -such claims; this section has the sole purpose of protecting the -integrity of the free software distribution system, which is -implemented by public license practices. Many people have made -generous contributions to the wide range of software distributed -through that system in reliance on consistent application of that -system; it is up to the author/donor to decide if he or she is willing -to distribute software through any other system and a licensee cannot -impose that choice. - -This section is intended to make thoroughly clear what is believed to -be a consequence of the rest of this License. - - 8. If the distribution and/or use of the Program is restricted in -certain countries either by patents or by copyrighted interfaces, the -original copyright holder who places the Program under this License -may add an explicit geographical distribution limitation excluding -those countries, so that distribution is permitted only in or among -countries not thus excluded. In such case, this License incorporates -the limitation as if written in the body of this License. - - 9. The Free Software Foundation may publish revised and/or new versions -of the General Public License from time to time. Such new versions will -be similar in spirit to the present version, but may differ in detail to -address new problems or concerns. - -Each version is given a distinguishing version number. If the Program -specifies a version number of this License which applies to it and "any -later version", you have the option of following the terms and conditions -either of that version or of any later version published by the Free -Software Foundation. If the Program does not specify a version number of -this License, you may choose any version ever published by the Free Software -Foundation. - - 10. If you wish to incorporate parts of the Program into other free -programs whose distribution conditions are different, write to the author -to ask for permission. For software which is copyrighted by the Free -Software Foundation, write to the Free Software Foundation; we sometimes -make exceptions for this. Our decision will be guided by the two goals -of preserving the free status of all derivatives of our free software and -of promoting the sharing and reuse of software generally. - - NO WARRANTY - - 11. BECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, THERE IS NO WARRANTY -FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW. EXCEPT WHEN -OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES -PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED -OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF -MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE ENTIRE RISK AS -TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU. SHOULD THE -PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING, -REPAIR OR CORRECTION. - - 12. IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING -WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY AND/OR -REDISTRIBUTE THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, -INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING -OUT OF THE USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED -TO LOSS OF DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY -YOU OR THIRD PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER -PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE -POSSIBILITY OF SUCH DAMAGES. - - END OF TERMS AND CONDITIONS - - How to Apply These Terms to Your New Programs - - If you develop a new program, and you want it to be of the greatest -possible use to the public, the best way to achieve this is to make it -free software which everyone can redistribute and change under these terms. - - To do so, attach the following notices to the program. It is safest -to attach them to the start of each source file to most effectively -convey the exclusion of warranty; and each file should have at least -the "copyright" line and a pointer to where the full notice is found. - - {description} - Copyright (C) {year} {fullname} - - This program is free software; you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation; either version 2 of the License, or - (at your option) any later version. - - This program is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License along - with this program; if not, write to the Free Software Foundation, Inc., - 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. - -Also add information on how to contact you by electronic and paper mail. - -If the program is interactive, make it output a short notice like this -when it starts in an interactive mode: - - Gnomovision version 69, Copyright (C) year name of author - Gnomovision comes with ABSOLUTELY NO WARRANTY; for details type `show w'. - This is free software, and you are welcome to redistribute it - under certain conditions; type `show c' for details. - -The hypothetical commands `show w' and `show c' should show the appropriate -parts of the General Public License. Of course, the commands you use may -be called something other than `show w' and `show c'; they could even be -mouse-clicks or menu items--whatever suits your program. - -You should also get your employer (if you work as a programmer) or your -school, if any, to sign a "copyright disclaimer" for the program, if -necessary. Here is a sample; alter the names: - - Yoyodyne, Inc., hereby disclaims all copyright interest in the program - `Gnomovision' (which makes passes at compilers) written by James Hacker. - - {signature of Ty Coon}, 1 April 1989 - Ty Coon, President of Vice - -This General Public License does not permit incorporating your program into -proprietary programs. If your program is a subroutine library, you may -consider it more useful to permit linking proprietary applications with the -library. If this is what you want to do, use the GNU Lesser General -Public License instead of this License. - diff --git a/src/include/spline/README.md b/src/include/spline/README.md deleted file mode 100644 index 65c4e8d..0000000 --- a/src/include/spline/README.md +++ /dev/null @@ -1,82 +0,0 @@ -## C++ cubic spline interpolation - -![cubic C2 spline](https://kluge.in-chemnitz.de/opensource/spline/cubic_c2_spline_git.png) - -This is a lightweight implementation of **cubic splines** -to interpolate points f(xi) = yi with -the following features. - -* available spline types: - * **cubic C2 splines**: global, twice continuously differentiable - * **cubic Hermite splines**: local, continuously differentiable (C1) -* boundary conditions: **first** and **second order** derivatives can be specified, **not-a-knot** condition, periodic condition is not implemented -* extrapolation - * linear: if first order derivatives are specified or 2nd order = 0 - * quadratic: if 2nd order derivatives not equal to zero specified -* monotonicity can be enforced (when input is monotonic as well) - -### Usage -The library is a header-only file with no external dependencies and can -be used like this: - -```C++ -#include -#include "spline.h" -... - std::vector X, Y; - ... - // default cubic spline (C^2) with natural boundary conditions (f''=0) - tk::spline s(X,Y); // X needs to be strictly increasing - double value=s(1.3); // interpolated value at 1.3 - double deriv=s.deriv(1,1.3); // 1st order derivative at 1.3 - std::vector solutions = s.solve(0.0); // solves s(x)=0.0 - ... -``` - -The constructor can take more arguments to define the spline, e.g.: -```C++ - // cubic Hermite splines (C^1) with enforced monotonicity and - // left curvature equal to 0.0 and right slope equal 1.0 - tk::spline s(X,Y,tk::spline::cspline_hermite, true, - tk::spline::second_deriv, 0.0, - tk::spline::first_deriv, 1.0); -``` -This is identical to (must be called in that order): -```C++ - tk::spline s; - s.set_boundary(tk::spline::second_deriv, 0.0, - tk::spline::first_deriv, 1.0); - s.set_points(X,Y); - s.make_monotonic(); -``` - - -### Spline types -Splines are piecewise polynomial functions to interpolate points -(xi, yi). In particular, cubic splines can -be represented as -* f(x) = ai + bi (x-xi) + ci (x-xi)2 + di (x-xi)3, for all x in [xi, xi+1) -* f(xi)=yi - -The following splines are available. - -* `tk::spline::cspline`: cubic C2 spline - * twice continuously differentiable, e.g. f'(xi) and f''(xi) exist - * this, together with boundary conditions uniquely determines the spline - * requires solving a sparse equation system - * is a global spline in the sense that changing an input point will impact the spline everywhere - * setting first order derivatives at the boundary will break C2 at the boundary -* `tk::spline::cspline_hermite`: cubic Hermite spline - * once continuously differentiable (C1) - * first order derivatives are specified by finite differences, e.g. on a uniform x-grid: - * f'(xi) = (yi+1-yi-1)/(xi+1-xi-1) - * is a local spline in the sense that changing grid points will only affect the spline around that grid point in a few adjacent segments - -A function to enforce monotonicity is available as well: -* `tk::spline::make_monotonic()`: will make the spline monotonic if input grid points are monotonic - * this function can only be called after `set_points(...)` has been called - * it will break C2 if the original spline was C2 and not already monotonic - * it will break boundary conditions if it was not monotonic in the first or last segment - -### References -https://kluge.in-chemnitz.de/opensource/spline/ diff --git a/src/include/spline/src/spline.h b/src/include/spline/src/spline.h deleted file mode 100644 index 7855d75..0000000 --- a/src/include/spline/src/spline.h +++ /dev/null @@ -1,1110 +0,0 @@ -/* - * spline.h - * - * simple cubic spline interpolation library without external - * dependencies - * - * --------------------------------------------------------------------- - * Copyright (C) 2011, 2014, 2016, 2021 Tino Kluge (ttk448 at gmail.com) - * - * This program is free software; you can redistribute it and/or - * modify it under the terms of the GNU General Public License - * as published by the Free Software Foundation; either version 2 - * of the License, or (at your option) any later version. - * - * This program is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with this program. If not, see . - * --------------------------------------------------------------------- - * - */ - -#ifndef TK_SPLINE_H -#define TK_SPLINE_H - -#include -#include -#include -#include -#include -#ifdef HAVE_SSTREAM -#include -#include -#endif // HAVE_SSTREAM - -// unnamed namespace only because the implementation is in this -// header file and we don't want to export symbols to the obj files -namespace -{ - -namespace tk -{ - -// spline interpolation -class spline -{ - public: - // spline types - enum spline_type - { - linear = 10, // linear interpolation - cspline = 30, // cubic splines (classical C^2) - cspline_hermite = 31 // cubic hermite splines (local, only C^1) - }; - - // boundary condition type for the spline end-points - enum bd_type - { - first_deriv = 1, - second_deriv = 2, - not_a_knot = 3 - }; - - protected: - std::vector m_x, m_y; // x,y coordinates of points - // interpolation parameters - // f(x) = a_i + b_i*(x-x_i) + c_i*(x-x_i)^2 + d_i*(x-x_i)^3 - // where a_i = y_i, or else it won't go through grid points - std::vector m_b, m_c, m_d; // spline coefficients - float m_c0; // for left extrapolation - spline_type m_type; - bd_type m_left, m_right; - float m_left_value, m_right_value; - bool m_made_monotonic; - void set_coeffs_from_b(); // calculate c_i, d_i from b_i - size_t find_closest(float x) const; // closest idx so that m_x[idx]<=x - - public: - // default constructor: set boundary condition to be zero curvature - // at both ends, i.e. natural splines - spline() - : m_type(cspline), - m_left(second_deriv), - m_right(second_deriv), - m_left_value(0.0), - m_right_value(0.0), - m_made_monotonic(false) - { - ; - } - spline(const std::vector &X, const std::vector &Y, - spline_type type = cspline, bool make_monotonic = false, - bd_type left = second_deriv, float left_value = 0.0, - bd_type right = second_deriv, float right_value = 0.0) - : m_type(type), - m_left(left), - m_right(right), - m_left_value(left_value), - m_right_value(right_value), - m_made_monotonic( - false) // false correct here: make_monotonic() sets it - { - this->set_points(X, Y, m_type); - if (make_monotonic) - { - this->make_monotonic(); - } - } - - // modify boundary conditions: if called it must be before set_points() - void set_boundary(bd_type left, float left_value, bd_type right, - float right_value); - - // set all data points (cubic_spline=false means linear interpolation) - void set_points(const std::vector &x, const std::vector &y, - spline_type type = cspline); - - // adjust coefficients so that the spline becomes piecewise monotonic - // where possible - // this is done by adjusting slopes at grid points by a non-negative - // factor and this will break C^2 - // this can also break boundary conditions if adjustments need to - // be made at the boundary points - // returns false if no adjustments have been made, true otherwise - bool make_monotonic(); - - // evaluates the spline at point x - float operator()(float x) const; - float deriv(int order, float x) const; - - // solves for all x so that: spline(x) = y - std::vector solve(float y, bool ignore_extrapolation = true) const; - float solve_one(float y, bool ignore_extrapolation = true) const; - - // returns the input data points - std::vector get_x() const { return m_x; } - std::vector get_y() const { return m_y; } - float get_x_min() const - { - assert(!m_x.empty()); - return m_x.front(); - } - float get_x_max() const - { - assert(!m_x.empty()); - return m_x.back(); - } - -#ifdef HAVE_SSTREAM - // spline info string, i.e. spline type, boundary conditions etc. - std::string info() const; -#endif // HAVE_SSTREAM -}; - -namespace internal -{ - -// band matrix solver -class band_matrix -{ - private: - std::vector > m_upper; // upper band - std::vector > m_lower; // lower band - public: - band_matrix(){}; // constructor - band_matrix(int dim, int n_u, int n_l); // constructor - ~band_matrix(){}; // destructor - void resize(int dim, int n_u, int n_l); // init with dim,n_u,n_l - int dim() const; // matrix dimension - int num_upper() const { return (int)m_upper.size() - 1; } - int num_lower() const { return (int)m_lower.size() - 1; } - // access operator - float &operator()(int i, int j); // write - float operator()(int i, int j) const; // read - // we can store an additional diagonal (in m_lower) - float &saved_diag(int i); - float saved_diag(int i) const; - void lu_decompose(); - std::vector r_solve(const std::vector &b) const; - std::vector l_solve(const std::vector &b) const; - std::vector lu_solve(const std::vector &b, - bool is_lu_decomposed = false); -}; - -float get_eps(); - -std::vector solve_cubic(float a, float b, float c, float d, - int newton_iter = 0); - -} // namespace internal - -// --------------------------------------------------------------------- -// implementation part, which could be separated into a cpp file -// --------------------------------------------------------------------- - -// spline implementation -// ----------------------- - -void spline::set_boundary(spline::bd_type left, float left_value, - spline::bd_type right, float right_value) -{ - assert(m_x.size() == 0); // set_points() must not have happened yet - m_left = left; - m_right = right; - m_left_value = left_value; - m_right_value = right_value; -} - -void spline::set_coeffs_from_b() -{ - assert(m_x.size() == m_y.size()); - assert(m_x.size() == m_b.size()); - assert(m_x.size() > 2); - size_t n = m_b.size(); - if (m_c.size() != n) m_c.resize(n); - if (m_d.size() != n) m_d.resize(n); - - for (size_t i = 0; i < n - 1; i++) - { - const float h = m_x[i + 1] - m_x[i]; - // from continuity and differentiability condition - m_c[i] = - (3.0 * (m_y[i + 1] - m_y[i]) / h - (2.0 * m_b[i] + m_b[i + 1])) / h; - // from differentiability condition - m_d[i] = ((m_b[i + 1] - m_b[i]) / (3.0 * h) - 2.0 / 3.0 * m_c[i]) / h; - } - - // for left extrapolation coefficients - m_c0 = (m_left == first_deriv) ? 0.0 : m_c[0]; -} - -void spline::set_points(const std::vector &x, - const std::vector &y, spline_type type) -{ - assert(x.size() == y.size()); - assert(x.size() >= 3); - // not-a-knot with 3 points has many solutions - if (m_left == not_a_knot || m_right == not_a_knot) assert(x.size() >= 4); - m_type = type; - m_made_monotonic = false; - m_x = x; - m_y = y; - int n = (int)x.size(); - // check strict monotonicity of input vector x - for (int i = 0; i < n - 1; i++) - { - assert(m_x[i] < m_x[i + 1]); - } - - if (type == linear) - { - // linear interpolation - m_d.resize(n); - m_c.resize(n); - m_b.resize(n); - for (int i = 0; i < n - 1; i++) - { - m_d[i] = 0.0; - m_c[i] = 0.0; - m_b[i] = (m_y[i + 1] - m_y[i]) / (m_x[i + 1] - m_x[i]); - } - // ignore boundary conditions, set slope equal to the last segment - m_b[n - 1] = m_b[n - 2]; - m_c[n - 1] = 0.0; - m_d[n - 1] = 0.0; - } - else if (type == cspline) - { - // classical cubic splines which are C^2 (twice cont differentiable) - // this requires solving an equation system - - // setting up the matrix and right hand side of the equation system - // for the parameters b[] - int n_upper = (m_left == spline::not_a_knot) ? 2 : 1; - int n_lower = (m_right == spline::not_a_knot) ? 2 : 1; - internal::band_matrix A(n, n_upper, n_lower); - std::vector rhs(n); - for (int i = 1; i < n - 1; i++) - { - A(i, i - 1) = 1.0 / 3.0 * (x[i] - x[i - 1]); - A(i, i) = 2.0 / 3.0 * (x[i + 1] - x[i - 1]); - A(i, i + 1) = 1.0 / 3.0 * (x[i + 1] - x[i]); - rhs[i] = (y[i + 1] - y[i]) / (x[i + 1] - x[i]) - - (y[i] - y[i - 1]) / (x[i] - x[i - 1]); - } - // boundary conditions - if (m_left == spline::second_deriv) - { - // 2*c[0] = f'' - A(0, 0) = 2.0; - A(0, 1) = 0.0; - rhs[0] = m_left_value; - } - else if (m_left == spline::first_deriv) - { - // b[0] = f', needs to be re-expressed in terms of c: - // (2c[0]+c[1])(x[1]-x[0]) = 3 ((y[1]-y[0])/(x[1]-x[0]) - f') - A(0, 0) = 2.0 * (x[1] - x[0]); - A(0, 1) = 1.0 * (x[1] - x[0]); - rhs[0] = 3.0 * ((y[1] - y[0]) / (x[1] - x[0]) - m_left_value); - } - else if (m_left == spline::not_a_knot) - { - // f'''(x[1]) exists, i.e. d[0]=d[1], or re-expressed in c: - // -h1*c[0] + (h0+h1)*c[1] - h0*c[2] = 0 - A(0, 0) = -(x[2] - x[1]); - A(0, 1) = x[2] - x[0]; - A(0, 2) = -(x[1] - x[0]); - rhs[0] = 0.0; - } - else - { - assert(false); - } - if (m_right == spline::second_deriv) - { - // 2*c[n-1] = f'' - A(n - 1, n - 1) = 2.0; - A(n - 1, n - 2) = 0.0; - rhs[n - 1] = m_right_value; - } - else if (m_right == spline::first_deriv) - { - // b[n-1] = f', needs to be re-expressed in terms of c: - // (c[n-2]+2c[n-1])(x[n-1]-x[n-2]) - // = 3 (f' - (y[n-1]-y[n-2])/(x[n-1]-x[n-2])) - A(n - 1, n - 1) = 2.0 * (x[n - 1] - x[n - 2]); - A(n - 1, n - 2) = 1.0 * (x[n - 1] - x[n - 2]); - rhs[n - 1] = 3.0 * (m_right_value - - (y[n - 1] - y[n - 2]) / (x[n - 1] - x[n - 2])); - } - else if (m_right == spline::not_a_knot) - { - // f'''(x[n-2]) exists, i.e. d[n-3]=d[n-2], or re-expressed in c: - // -h_{n-2}*c[n-3] + (h_{n-3}+h_{n-2})*c[n-2] - h_{n-3}*c[n-1] = 0 - A(n - 1, n - 3) = -(x[n - 1] - x[n - 2]); - A(n - 1, n - 2) = x[n - 1] - x[n - 3]; - A(n - 1, n - 1) = -(x[n - 2] - x[n - 3]); - rhs[0] = 0.0; - } - else - { - assert(false); - } - - // solve the equation system to obtain the parameters c[] - m_c = A.lu_solve(rhs); - - // calculate parameters b[] and d[] based on c[] - m_d.resize(n); - m_b.resize(n); - for (int i = 0; i < n - 1; i++) - { - m_d[i] = 1.0 / 3.0 * (m_c[i + 1] - m_c[i]) / (x[i + 1] - x[i]); - m_b[i] = - (y[i + 1] - y[i]) / (x[i + 1] - x[i]) - - 1.0 / 3.0 * (2.0 * m_c[i] + m_c[i + 1]) * (x[i + 1] - x[i]); - } - // for the right extrapolation coefficients (zero cubic term) - // f_{n-1}(x) = y_{n-1} + b*(x-x_{n-1}) + c*(x-x_{n-1})^2 - float h = x[n - 1] - x[n - 2]; - // m_c[n-1] is determined by the boundary condition - m_d[n - 1] = 0.0; - m_b[n - 1] = 3.0 * m_d[n - 2] * h * h + 2.0 * m_c[n - 2] * h + - m_b[n - 2]; // = f'_{n-2}(x_{n-1}) - if (m_right == first_deriv) - m_c[n - 1] = 0.0; // force linear extrapolation - } - else if (type == cspline_hermite) - { - // hermite cubic splines which are C^1 (cont. differentiable) - // and derivatives are specified on each grid point - // (here we use 3-point finite differences) - m_b.resize(n); - m_c.resize(n); - m_d.resize(n); - // set b to match 1st order derivative finite difference - for (int i = 1; i < n - 1; i++) - { - const float h = m_x[i + 1] - m_x[i]; - const float hl = m_x[i] - m_x[i - 1]; - m_b[i] = -h / (hl * (hl + h)) * m_y[i - 1] + - (h - hl) / (hl * h) * m_y[i] + - hl / (h * (hl + h)) * m_y[i + 1]; - } - // boundary conditions determine b[0] and b[n-1] - if (m_left == first_deriv) - { - m_b[0] = m_left_value; - } - else if (m_left == second_deriv) - { - const float h = m_x[1] - m_x[0]; - m_b[0] = 0.5 * (-m_b[1] - 0.5 * m_left_value * h + - 3.0 * (m_y[1] - m_y[0]) / h); - } - else if (m_left == not_a_knot) - { - // f''' continuous at x[1] - const float h0 = m_x[1] - m_x[0]; - const float h1 = m_x[2] - m_x[1]; - m_b[0] = -m_b[1] + 2.0 * (m_y[1] - m_y[0]) / h0 + - h0 * h0 / (h1 * h1) * - (m_b[1] + m_b[2] - 2.0 * (m_y[2] - m_y[1]) / h1); - } - else - { - assert(false); - } - if (m_right == first_deriv) - { - m_b[n - 1] = m_right_value; - m_c[n - 1] = 0.0; - } - else if (m_right == second_deriv) - { - const float h = m_x[n - 1] - m_x[n - 2]; - m_b[n - 1] = 0.5 * (-m_b[n - 2] + 0.5 * m_right_value * h + - 3.0 * (m_y[n - 1] - m_y[n - 2]) / h); - m_c[n - 1] = 0.5 * m_right_value; - } - else if (m_right == not_a_knot) - { - // f''' continuous at x[n-2] - const float h0 = m_x[n - 2] - m_x[n - 3]; - const float h1 = m_x[n - 1] - m_x[n - 2]; - m_b[n - 1] = -m_b[n - 2] + 2.0 * (m_y[n - 1] - m_y[n - 2]) / h1 + - h1 * h1 / (h0 * h0) * - (m_b[n - 3] + m_b[n - 2] - - 2.0 * (m_y[n - 2] - m_y[n - 3]) / h0); - // f'' continuous at x[n-1]: c[n-1] = 3*d[n-2]*h[n-2] + c[n-1] - m_c[n - 1] = (m_b[n - 2] + 2.0 * m_b[n - 1]) / h1 - - 3.0 * (m_y[n - 1] - m_y[n - 2]) / (h1 * h1); - } - else - { - assert(false); - } - m_d[n - 1] = 0.0; - - // parameters c and d are determined by continuity and differentiability - set_coeffs_from_b(); - } - else - { - assert(false); - } - - // for left extrapolation coefficients - m_c0 = (m_left == first_deriv) ? 0.0 : m_c[0]; -} - -bool spline::make_monotonic() -{ - assert(m_x.size() == m_y.size()); - assert(m_x.size() == m_b.size()); - assert(m_x.size() > 2); - bool modified = false; - const int n = (int)m_x.size(); - // make sure: input data monotonic increasing --> b_i>=0 - // input data monotonic decreasing --> b_i<=0 - for (int i = 0; i < n; i++) - { - int im1 = std::max(i - 1, 0); - int ip1 = std::min(i + 1, n - 1); - if (((m_y[im1] <= m_y[i]) && (m_y[i] <= m_y[ip1]) && m_b[i] < 0.0) || - ((m_y[im1] >= m_y[i]) && (m_y[i] >= m_y[ip1]) && m_b[i] > 0.0)) - { - modified = true; - m_b[i] = 0.0; - } - } - // if input data is monotonic (b[i], b[i+1], avg have all the same sign) - // ensure a sufficient criteria for monotonicity is satisfied: - // sqrt(b[i]^2+b[i+1]^2) <= 3 |avg|, with avg=(y[i+1]-y[i])/h, - for (int i = 0; i < n - 1; i++) - { - float h = m_x[i + 1] - m_x[i]; - float avg = (m_y[i + 1] - m_y[i]) / h; - if (avg == 0.0 && (m_b[i] != 0.0 || m_b[i + 1] != 0.0)) - { - modified = true; - m_b[i] = 0.0; - m_b[i + 1] = 0.0; - } - else if ((m_b[i] >= 0.0 && m_b[i + 1] >= 0.0 && avg > 0.0) || - (m_b[i] <= 0.0 && m_b[i + 1] <= 0.0 && avg < 0.0)) - { - // input data is monotonic - float r = sqrt(m_b[i] * m_b[i] + m_b[i + 1] * m_b[i + 1]) / - std::fabs(avg); - if (r > 3.0) - { - // sufficient criteria for monotonicity: r<=3 - // adjust b[i] and b[i+1] - modified = true; - m_b[i] *= (3.0 / r); - m_b[i + 1] *= (3.0 / r); - } - } - } - - if (modified == true) - { - set_coeffs_from_b(); - m_made_monotonic = true; - } - - return modified; -} - -// return the closest idx so that m_x[idx] <= x (return 0 if x::const_iterator it; - it = std::upper_bound(m_x.begin(), m_x.end(), x); // *it > x - size_t idx = std::max(int(it - m_x.begin()) - 1, 0); // m_x[idx] <= x - return idx; -} - -float spline::operator()(float x) const -{ - // polynomial evaluation using Horner's scheme - // TODO: consider more numerically accurate algorithms, e.g.: - // - Clenshaw - // - Even-Odd method by A.C.R. Newbery - // - Compensated Horner Scheme - size_t n = m_x.size(); - size_t idx = find_closest(x); - - float h = x - m_x[idx]; - float interpol; - if (x < m_x[0]) - { - // extrapolation to the left - interpol = (m_c0 * h + m_b[0]) * h + m_y[0]; - } - else if (x > m_x[n - 1]) - { - // extrapolation to the right - interpol = (m_c[n - 1] * h + m_b[n - 1]) * h + m_y[n - 1]; - } - else - { - // interpolation - interpol = ((m_d[idx] * h + m_c[idx]) * h + m_b[idx]) * h + m_y[idx]; - } - return interpol; -} - -float spline::deriv(int order, float x) const -{ - assert(order > 0); - size_t n = m_x.size(); - size_t idx = find_closest(x); - - float h = x - m_x[idx]; - float interpol; - if (x < m_x[0]) - { - // extrapolation to the left - switch (order) - { - case 1: - interpol = 2.0 * m_c0 * h + m_b[0]; - break; - case 2: - interpol = 2.0 * m_c0; - break; - default: - interpol = 0.0; - break; - } - } - else if (x > m_x[n - 1]) - { - // extrapolation to the right - switch (order) - { - case 1: - interpol = 2.0 * m_c[n - 1] * h + m_b[n - 1]; - break; - case 2: - interpol = 2.0 * m_c[n - 1]; - break; - default: - interpol = 0.0; - break; - } - } - else - { - // interpolation - switch (order) - { - case 1: - interpol = (3.0 * m_d[idx] * h + 2.0 * m_c[idx]) * h + m_b[idx]; - break; - case 2: - interpol = 6.0 * m_d[idx] * h + 2.0 * m_c[idx]; - break; - case 3: - interpol = 6.0 * m_d[idx]; - break; - default: - interpol = 0.0; - break; - } - } - return interpol; -} - -std::vector spline::solve(float y, bool ignore_extrapolation) const -{ - std::vector x; // roots for the entire spline - std::vector root; // roots for each piecewise cubic - const size_t n = m_x.size(); - - // left extrapolation - if (ignore_extrapolation == false) - { - root = internal::solve_cubic(m_y[0] - y, m_b[0], m_c0, 0.0, 1); - for (size_t j = 0; j < root.size(); j++) - { - if (root[j] < 0.0) - { - x.push_back(m_x[0] + root[j]); - } - } - } - - // brute force check if piecewise cubic has roots in their resp. segment - // TODO: make more efficient - for (size_t i = 0; i < n - 1; i++) - { - root = internal::solve_cubic(m_y[i] - y, m_b[i], m_c[i], m_d[i], 1); - for (size_t j = 0; j < root.size(); j++) - { - float h = (i > 0) ? (m_x[i] - m_x[i - 1]) : 0.0; - float eps = internal::get_eps() * 512.0 * std::min(h, 1.0f); - if ((-eps <= root[j]) && (root[j] < m_x[i + 1] - m_x[i])) - { - float new_root = m_x[i] + root[j]; - if (x.size() > 0 && x.back() + eps > new_root) - { - x.back() = new_root; // avoid spurious duplicate roots - } - else - { - x.push_back(new_root); - } - } - } - } - - // right extrapolation - if (ignore_extrapolation == false) - { - root = internal::solve_cubic(m_y[n - 1] - y, m_b[n - 1], m_c[n - 1], - 0.0, 1); - for (size_t j = 0; j < root.size(); j++) - { - if (0.0 <= root[j]) - { - x.push_back(m_x[n - 1] + root[j]); - } - } - } - - return x; -}; - -float spline::solve_one(float y, bool ignore_extrapolation) const -{ - std::vector root; // roots for each piecewise cubic - const size_t n = m_x.size(); - - // left extrapolation - if (ignore_extrapolation == false) - { - root = internal::solve_cubic(m_y[0] - y, m_b[0], m_c0, 0.0, 1); - for (float j : root) - { - if (j < 0.0) return m_x[0] + j; - } - } - - // brute force check if piecewise cubic has roots in their resp. segment - // TODO: make more efficient - for (size_t i = 0; i < n - 1; i++) - { - root = internal::solve_cubic(m_y[i] - y, m_b[i], m_c[i], m_d[i], 1); - for (float j : root) - { - float h = (i > 0) ? (m_x[i] - m_x[i - 1]) : 0.0; - float eps = internal::get_eps() * 512.0 * std::min(h, 1.0f); - if ((-eps <= j) && (j < m_x[i + 1] - m_x[i])) - { - return m_x[i] + j; - } - } - } - - // right extrapolation - if (ignore_extrapolation == false) - { - root = internal::solve_cubic(m_y[n - 1] - y, m_b[n - 1], m_c[n - 1], - 0.0, 1); - for (size_t j = 0; j < root.size(); j++) - { - if (0.0 <= root[j]) return m_x[n - 1] + root[j]; - } - } - - return std::numeric_limits::quiet_NaN(); -} - -#ifdef HAVE_SSTREAM -std::string spline::info() const -{ - std::stringstream ss; - ss << "type " << m_type << ", left boundary deriv " << m_left << " = "; - ss << m_left_value << ", right boundary deriv " << m_right << " = "; - ss << m_right_value << std::endl; - if (m_made_monotonic) - { - ss << "(spline has been adjusted for piece-wise monotonicity)"; - } - return ss.str(); -} -#endif // HAVE_SSTREAM - -namespace internal -{ - -// band_matrix implementation -// ------------------------- - -band_matrix::band_matrix(int dim, int n_u, int n_l) { resize(dim, n_u, n_l); } -void band_matrix::resize(int dim, int n_u, int n_l) -{ - assert(dim > 0); - assert(n_u >= 0); - assert(n_l >= 0); - m_upper.resize(n_u + 1); - m_lower.resize(n_l + 1); - for (size_t i = 0; i < m_upper.size(); i++) - { - m_upper[i].resize(dim); - } - for (size_t i = 0; i < m_lower.size(); i++) - { - m_lower[i].resize(dim); - } -} -int band_matrix::dim() const -{ - if (m_upper.size() > 0) - { - return m_upper[0].size(); - } - else - { - return 0; - } -} - -// defines the new operator (), so that we can access the elements -// by A(i,j), index going from i=0,...,dim()-1 -float &band_matrix::operator()(int i, int j) -{ - int k = j - i; // what band is the entry - assert((i >= 0) && (i < dim()) && (j >= 0) && (j < dim())); - assert((-num_lower() <= k) && (k <= num_upper())); - // k=0 -> diagonal, k<0 lower left part, k>0 upper right part - if (k >= 0) - return m_upper[k][i]; - else - return m_lower[-k][i]; -} -float band_matrix::operator()(int i, int j) const -{ - int k = j - i; // what band is the entry - assert((i >= 0) && (i < dim()) && (j >= 0) && (j < dim())); - assert((-num_lower() <= k) && (k <= num_upper())); - // k=0 -> diagonal, k<0 lower left part, k>0 upper right part - if (k >= 0) - return m_upper[k][i]; - else - return m_lower[-k][i]; -} -// second diag (used in LU decomposition), saved in m_lower -float band_matrix::saved_diag(int i) const -{ - assert((i >= 0) && (i < dim())); - return m_lower[0][i]; -} -float &band_matrix::saved_diag(int i) -{ - assert((i >= 0) && (i < dim())); - return m_lower[0][i]; -} - -// LR-Decomposition of a band matrix -void band_matrix::lu_decompose() -{ - int i_max, j_max; - int j_min; - float x; - - // preconditioning - // normalize column i so that a_ii=1 - for (int i = 0; i < this->dim(); i++) - { - assert(this->operator()(i, i) != 0.0); - this->saved_diag(i) = 1.0 / this->operator()(i, i); - j_min = std::max(0, i - this->num_lower()); - j_max = std::min(this->dim() - 1, i + this->num_upper()); - for (int j = j_min; j <= j_max; j++) - { - this->operator()(i, j) *= this->saved_diag(i); - } - this->operator()(i, i) = 1.0; // prevents rounding errors - } - - // Gauss LR-Decomposition - for (int k = 0; k < this->dim(); k++) - { - i_max = std::min(this->dim() - 1, - k + this->num_lower()); // num_lower not a mistake! - for (int i = k + 1; i <= i_max; i++) - { - assert(this->operator()(k, k) != 0.0); - x = -this->operator()(i, k) / this->operator()(k, k); - this->operator()(i, k) = -x; // assembly part of L - j_max = std::min(this->dim() - 1, k + this->num_upper()); - for (int j = k + 1; j <= j_max; j++) - { - // assembly part of R - this->operator()(i, j) = - this->operator()(i, j) + x * this->operator()(k, j); - } - } - } -} -// solves Ly=b -std::vector band_matrix::l_solve(const std::vector &b) const -{ - assert(this->dim() == (int)b.size()); - std::vector x(this->dim()); - int j_start; - float sum; - for (int i = 0; i < this->dim(); i++) - { - sum = 0; - j_start = std::max(0, i - this->num_lower()); - for (int j = j_start; j < i; j++) sum += this->operator()(i, j) * x[j]; - x[i] = (b[i] * this->saved_diag(i)) - sum; - } - return x; -} -// solves Rx=y -std::vector band_matrix::r_solve(const std::vector &b) const -{ - assert(this->dim() == (int)b.size()); - std::vector x(this->dim()); - int j_stop; - float sum; - for (int i = this->dim() - 1; i >= 0; i--) - { - sum = 0; - j_stop = std::min(this->dim() - 1, i + this->num_upper()); - for (int j = i + 1; j <= j_stop; j++) - sum += this->operator()(i, j) * x[j]; - x[i] = (b[i] - sum) / this->operator()(i, i); - } - return x; -} - -std::vector band_matrix::lu_solve(const std::vector &b, - bool is_lu_decomposed) -{ - assert(this->dim() == (int)b.size()); - std::vector x, y; - if (is_lu_decomposed == false) - { - this->lu_decompose(); - } - y = this->l_solve(b); - x = this->r_solve(y); - return x; -} - -// machine precision of a float, i.e. the successor of 1 is 1+eps -float get_eps() -{ - // return std::numeric_limits::epsilon(); // __DBL_EPSILON__ - return 2.2204460492503131e-16; // 2^-52 -} - -// solutions for a + b*x = 0 -std::vector solve_linear(float a, float b) -{ - std::vector x; // roots - if (b == 0.0) - { - if (a == 0.0) - { - // 0*x = 0 - x.resize(1); - x[0] = 0.0; // any x solves it but we need to pick one - return x; - } - else - { - // 0*x + ... = 0, no solution - return x; - } - } - else - { - x.resize(1); - x[0] = -a / b; - return x; - } -} - -// solutions for a + b*x + c*x^2 = 0 -std::vector solve_quadratic(float a, float b, float c, - int newton_iter = 0) -{ - if (c == 0.0) - { - return solve_linear(a, b); - } - // rescale so that we solve x^2 + 2p x + q = (x+p)^2 + q - p^2 = 0 - float p = 0.5 * b / c; - float q = a / c; - float discr = p * p - q; - const float eps = 0.5 * internal::get_eps(); - float discr_err = (6.0 * (p * p) + 3.0 * fabs(q) + fabs(discr)) * eps; - - std::vector x; // roots - if (fabs(discr) <= discr_err) - { - // discriminant is zero --> one root - x.resize(1); - x[0] = -p; - } - else if (discr < 0) - { - // no root - } - else - { - // two roots - x.resize(2); - x[0] = -p - sqrt(discr); - x[1] = -p + sqrt(discr); - } - - // improve solution via newton steps - for (size_t i = 0; i < x.size(); i++) - { - for (int k = 0; k < newton_iter; k++) - { - float f = (c * x[i] + b) * x[i] + a; - float f1 = 2.0 * c * x[i] + b; - // only adjust if slope is large enough - if (fabs(f1) > 1e-8) - { - x[i] -= f / f1; - } - } - } - - return x; -} - -// solutions for the cubic equation: a + b*x +c*x^2 + d*x^3 = 0 -// this is a naive implementation of the analytic solution without -// optimisation for speed or numerical accuracy -// newton_iter: number of newton iterations to improve analytical solution -// see also -// gsl: gsl_poly_solve_cubic() in solve_cubic.c -// octave: roots.m - via eigenvalues of the Frobenius companion matrix -std::vector solve_cubic(float a, float b, float c, float d, - int newton_iter) -{ - if (d == 0.0) - { - return solve_quadratic(a, b, c, newton_iter); - } - - // convert to normalised form: a + bx + cx^2 + x^3 = 0 - if (d != 1.0) - { - a /= d; - b /= d; - c /= d; - } - - // convert to depressed cubic: z^3 - 3pz - 2q = 0 - // via substitution: z = x + c/3 - std::vector z; // roots of the depressed cubic - float p = -(1.0 / 3.0) * b + (1.0 / 9.0) * (c * c); - float r = 2.0 * (c * c) - 9.0 * b; - float q = -0.5 * a - (1.0 / 54.0) * (c * r); - float discr = p * p * p - q * q; // discriminant - // calculating numerical round-off errors with assumptions: - // - each operation is precise but each intermediate result x - // when stored has max error of x*eps - // - only multiplication with a power of 2 introduces no new error - // - a,b,c,d and some fractions (e.g. 1/3) have rounding errors eps - // - p_err << |p|, q_err << |q|, ... (this is violated in rare cases) - // would be more elegant to use boost::numeric::interval - const float eps = internal::get_eps(); - float p_err = - eps * ((3.0 / 3.0) * fabs(b) + (4.0 / 9.0) * (c * c) + fabs(p)); - float r_err = eps * (6.0 * (c * c) + 18.0 * fabs(b) + fabs(r)); - float q_err = 0.5 * fabs(a) * eps + - (1.0 / 54.0) * fabs(c) * (r_err + fabs(r) * 3.0 * eps) + - fabs(q) * eps; - float discr_err = (p * p) * (3.0 * p_err + fabs(p) * 2.0 * eps) + - fabs(q) * (2.0 * q_err + fabs(q) * eps) + - fabs(discr) * eps; - - // depending on the discriminant we get different solutions - if (fabs(discr) <= discr_err) - { - // discriminant zero: one or two real roots - if (fabs(p) <= p_err) - { - // p and q are zero: single root - z.resize(1); - z[0] = 0.0; // triple root - } - else - { - z.resize(2); - z[0] = 2.0 * q / p; // single root - z[1] = -0.5 * z[0]; // float root - } - } - else if (discr > 0) - { - // three real roots: via trigonometric solution - z.resize(3); - float ac = (1.0 / 3.0) * acos(q / (p * sqrt(p))); - float sq = 2.0 * sqrt(p); - z[0] = sq * cos(ac); - z[1] = sq * cos(ac - 2.0 * M_PI / 3.0); - z[2] = sq * cos(ac - 4.0 * M_PI / 3.0); - } - else if (discr < 0.0) - { - // single real root: via Cardano's fromula - z.resize(1); - float sgnq = (q >= 0 ? 1 : -1); - float basis = fabs(q) + sqrt(-discr); - float C = sgnq * pow(basis, 1.0 / 3.0); // c++11 has std::cbrt() - z[0] = C + p / C; - } - for (size_t i = 0; i < z.size(); i++) - { - // convert depressed cubic roots to original cubic: x = z - c/3 - z[i] -= (1.0 / 3.0) * c; - // improve solution via newton steps - for (int k = 0; k < newton_iter; k++) - { - float f = ((z[i] + c) * z[i] + b) * z[i] + a; - float f1 = (3.0 * z[i] + 2.0 * c) * z[i] + b; - // only adjust if slope is large enough - if (fabs(f1) > 1e-8) - { - z[i] -= f / f1; - } - } - } - // ensure if a=0 we get exactly x=0 as root - // TODO: remove this fudge - if (a == 0.0) - { - assert(z.size() > 0); // cubic should always have at least one root - float xmin = fabs(z[0]); - size_t imin = 0; - for (size_t i = 1; i < z.size(); i++) - { - if (xmin > fabs(z[i])) - { - xmin = fabs(z[i]); - imin = i; - } - } - z[imin] = 0.0; // replace the smallest absolute value with 0 - } - std::sort(z.begin(), z.end()); - return z; -} - -} // namespace internal - -} // namespace tk - -} // namespace - -#endif /* TK_SPLINE_H */ diff --git a/src/mcmc.cpp b/src/mcmc.cpp index a31d877..a3cd5cc 100644 --- a/src/mcmc.cpp +++ b/src/mcmc.cpp @@ -1,7 +1,7 @@ #include "mcmc.h" -#include "include/spline/src/spline.h" #include "mcmc_utils.h" +#include "monotone_interpolator.h" #include "seed.h" #include @@ -253,53 +253,60 @@ void MCMC::swap_chains(int step, bool burnin) int MCMC::get_hot_chain() { return swap_indices[0]; } +// Re-space the temperature ladder so that every adjacent pair of chains has +// the same accumulated communication barrier (expected swap rejections). +// Temperatures are interpolated as a function of cumulative barrier, so the +// new ladder is read off an even barrier grid directly. void MCMC::adapt_temp() { - // swap rate starts at t = 1 so we need to reverse the swap rate - const std::vector reversed_swap_barriers{swap_barriers.rbegin(), - swap_barriers.rend()}; + const size_t n = temp_gradient.size(); - // cumulative swap rate starts from t = 0 - std::vector cumulative_swap_rate(params.pt_chains.size(), 0.0); - for (size_t i = 1; i < cumulative_swap_rate.size(); i++) + // Cumulative barrier from the cold end (temp_gradient[0] == 1). The + // per-swap normalisation cancels out of the interpolation, so raw sums + // are used. + std::vector barrier(n, 0.0); + for (size_t i = 1; i < n; i++) { - cumulative_swap_rate[i] = cumulative_swap_rate[i - 1] + - reversed_swap_barriers[i - 1] / (static_cast(num_swaps) / 2.0f); + barrier[i] = barrier[i - 1] + static_cast(swap_barriers[i - 1]); } - - // gradient starts at t = 1 so we need to reverse the gradient - const std::vector reversed_gradient(temp_gradient.rbegin(), temp_gradient.rend()); - - const tk::spline s(reversed_gradient, cumulative_swap_rate, tk::spline::cspline, true); - - // target swap rates - std::vector cumulative_swap_grid = std::vector(cumulative_swap_rate.size(), 0.0f); - - cumulative_swap_grid.back() = cumulative_swap_rate.back(); - float step = cumulative_swap_rate.back() / (cumulative_swap_grid.size() - 1); - - for (size_t i = 1; i < cumulative_swap_grid.size() - 1; i++) + if (!std::isfinite(barrier.back()) || barrier.back() <= 0.0) { - cumulative_swap_grid[i] = cumulative_swap_grid[0] + i * step; + return; // nothing observed yet; keep accumulating } - std::vector new_temp_gradient(temp_gradient.size()); - new_temp_gradient[0] = temp_gradient.back(); - new_temp_gradient[temp_gradient.size() - 1] = 1.0; - for (size_t i = 1; i < temp_gradient.size() - 1; i++) + // Adjacent chains that never rejected a swap share a barrier value. + // Collapse each such run into one knot at its mean temperature so the + // knots are strictly increasing in barrier. + std::vector knot_x; + std::vector knot_y; + for (size_t i = 0; i < n;) { - new_temp_gradient[i] = s.solve_one(cumulative_swap_grid[i]); + size_t j = i; + double temp_sum = 0.0; + while (j < n && barrier[j] == barrier[i]) + { + temp_sum += temp_gradient[j]; + ++j; + } + knot_x.push_back(barrier[i]); + knot_y.push_back(temp_sum / static_cast(j - i)); + i = j; + } + if (knot_x.size() < 2) + { + return; } - std::reverse(new_temp_gradient.begin(), new_temp_gradient.end()); + const MonotoneCubicInterpolator temp_at_barrier(knot_x, knot_y); - // check if any new temperatures are NAN, bail on the update if so and keep updating the swap barriers - for (size_t i = 0; i < new_temp_gradient.size(); i++) + std::vector new_temp_gradient(n); + new_temp_gradient.front() = temp_gradient.front(); + new_temp_gradient.back() = temp_gradient.back(); + const double step = barrier.back() / static_cast(n - 1); + for (size_t i = 1; i + 1 < n; i++) { - if (std::isnan(new_temp_gradient[i])) - { - return; - } + new_temp_gradient[i] = + static_cast(temp_at_barrier(static_cast(i) * step)); } // update the temperatures using the new gradient diff --git a/src/monotone_interpolator.h b/src/monotone_interpolator.h new file mode 100644 index 0000000..730f87c --- /dev/null +++ b/src/monotone_interpolator.h @@ -0,0 +1,102 @@ +#pragma once + +#ifndef MONOTONE_INTERPOLATOR_H_ +#define MONOTONE_INTERPOLATOR_H_ + +#include +#include +#include +#include +#include + +// Shape-preserving cubic Hermite interpolation through knots with strictly +// increasing x and monotone y (Fritsch & Carlson 1980, with the +// Fritsch-Butland slope formula used by PCHIP). If the knot values are +// monotone the interpolant is monotone. Evaluation is clamped to the knot +// range. MCMC::adapt_temp() uses this to invert the cumulative communication +// barrier: temperature as a function of barrier. +class MonotoneCubicInterpolator +{ + public: + // Requires x.size() == y.size() >= 2 and x strictly increasing. + MonotoneCubicInterpolator(std::vector x, std::vector y) + : x_(std::move(x)), y_(std::move(y)), m_(x_.size(), 0.0) + { + const std::size_t n = x_.size(); + std::vector h(n - 1); + std::vector delta(n - 1); + for (std::size_t k = 0; k + 1 < n; ++k) + { + h[k] = x_[k + 1] - x_[k]; + delta[k] = (y_[k + 1] - y_[k]) / h[k]; + } + + if (n == 2) + { + m_[0] = m_[1] = delta[0]; + return; + } + + for (std::size_t k = 1; k + 1 < n; ++k) + { + if (delta[k - 1] * delta[k] <= 0.0) + { + m_[k] = 0.0; // local extremum or flat: keep it flat + } + else + { + const double w1 = 2.0 * h[k] + h[k - 1]; + const double w2 = h[k] + 2.0 * h[k - 1]; + m_[k] = (w1 + w2) / (w1 / delta[k - 1] + w2 / delta[k]); + } + } + m_[0] = edge_slope(h[0], h[1], delta[0], delta[1]); + m_[n - 1] = edge_slope(h[n - 2], h[n - 3], delta[n - 2], delta[n - 3]); + } + + double operator()(double x) const + { + if (x <= x_.front()) + { + return y_.front(); + } + if (x >= x_.back()) + { + return y_.back(); + } + const std::size_t k = + static_cast(std::upper_bound(x_.begin(), x_.end(), x) - x_.begin()) - 1; + const double h = x_[k + 1] - x_[k]; + const double t = (x - x_[k]) / h; + const double t2 = t * t; + const double t3 = t2 * t; + const double h00 = 2.0 * t3 - 3.0 * t2 + 1.0; + const double h10 = t3 - 2.0 * t2 + t; + const double h01 = -2.0 * t3 + 3.0 * t2; + const double h11 = t3 - t2; + return h00 * y_[k] + h10 * h * m_[k] + h01 * y_[k + 1] + h11 * h * m_[k + 1]; + } + + private: + // Three-point one-sided estimate for an end knot, clamped so the end + // segment stays shape preserving. h0/d0 belong to the end segment. + static double edge_slope(double h0, double h1, double d0, double d1) + { + double m = ((2.0 * h0 + h1) * d0 - h0 * d1) / (h0 + h1); + if (m * d0 <= 0.0) + { + return 0.0; + } + if (d0 * d1 <= 0.0 && std::fabs(m) > 3.0 * std::fabs(d0)) + { + return 3.0 * d0; + } + return m; + } + + std::vector x_; + std::vector y_; + std::vector m_; +}; + +#endif // MONOTONE_INTERPOLATOR_H_ diff --git a/tests/testthat/test-adapt-temp.R b/tests/testthat/test-adapt-temp.R new file mode 100644 index 0000000..e01601b --- /dev/null +++ b/tests/testthat/test-adapt-temp.R @@ -0,0 +1,56 @@ +pt_data <- function() { + set.seed(7) + moire::simulate_data( + mean_coi = 3, + num_samples = 12, + epsilon_pos = 0.01, + epsilon_neg = 0.1, + locus_freq_alphas = rep(list(rep(1, 5)), 6), + missingness = 0 + ) +} + +run_pt <- function(adapt_temp, pt_chains = 8) { + data <- pt_data() + moire::run_mcmc( + data, + is_missing = data$is_missing, + verbose = FALSE, + burnin = 200, + samples_per_chain = 20, + pt_chains = pt_chains, + adapt_temp = adapt_temp, + pre_adapt_steps = 25, + temp_adapt_steps = 25, + seed = 3L + ) +} + +test_that("adapted temperature ladder is finite, monotone, and keeps its ends", { + res <- run_pt(adapt_temp = TRUE) + ladder <- res$chains[[1]]$temp_gradient + expect_length(ladder, 8) + expect_true(all(is.finite(ladder))) + expect_equal(ladder[1], 1) + expect_equal(ladder[length(ladder)], 0) + expect_true(all(diff(ladder) < 0)) + expect_false(isTRUE(all.equal(ladder, seq(1, 0, length.out = 8)))) +}) + +test_that("temperature ladder is untouched when adapt_temp = FALSE", { + res <- run_pt(adapt_temp = FALSE) + expect_equal( + res$chains[[1]]$temp_gradient, seq(1, 0, length.out = 8), + tolerance = 1e-6 + ) +}) + +test_that("adaptation works with a user supplied temperature vector", { + temps <- c(1, 0.7, 0.4, 0.2, 0.05) + res <- run_pt(adapt_temp = TRUE, pt_chains = temps) + ladder <- res$chains[[1]]$temp_gradient + expect_length(ladder, length(temps)) + expect_equal(ladder[1], 1) + expect_equal(ladder[length(ladder)], 0.05) + expect_true(all(diff(ladder) < 0)) +})