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 7818df2..99f0a83 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", @@ -12,8 +12,11 @@ Authors@R: c( 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. +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 +34,7 @@ Imports: dplyr, tidyr, stats, + parallel, purrr, rlang, ggplot2 @@ -41,7 +45,6 @@ RoxygenNote: 7.3.3 Suggests: knitr, rmarkdown, - markdown, forcats, testthat (>= 3.0.0), parallelly diff --git a/NEWS.md b/NEWS.md index 46a99ee..956dd05 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,4 +1,4 @@ -# moire (development version) +# moire 3.7.0 ## New Features @@ -8,6 +8,20 @@ `$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. +* 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 ## New Features diff --git a/R/data.R b/R/data.R index db01475..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 @@ -22,7 +30,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/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..be44733 --- /dev/null +++ b/cran-comments.md @@ -0,0 +1,20 @@ +## 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. 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/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 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) +} 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)) +})