From 83db4ae2bf95032a2639052b0b324382e9b2b095 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Sat, 16 May 2026 14:49:53 +0000 Subject: [PATCH 1/3] fix: magmap rank bug, magbin nn.dists typo, plot.magbin speed, docs, dead code, add tests and vignette sections Agent-Logs-Url: https://github.com/asgr/magicaxis/sessions/74e9c01b-f592-4d54-8e81-353d57d8b91c Co-authored-by: asgr <5617132+asgr@users.noreply.github.com> --- R/magband.R | 12 -- R/magbin.R | 41 ++-- R/magmap.R | 2 +- man/magclip.Rd | 12 +- man/magmap.Rd | 2 +- man/magrun.Rd | 4 +- tests/testthat/test-functions.R | 326 ++++++++++++++++++++++++++++++++ vignettes/magicaxis.Rmd | 159 +++++++++++++++- 8 files changed, 510 insertions(+), 48 deletions(-) create mode 100644 tests/testthat/test-functions.R diff --git a/R/magband.R b/R/magband.R index c95796b..375b1cd 100644 --- a/R/magband.R +++ b/R/magband.R @@ -75,18 +75,6 @@ return = temp } -".Last.magroj"<- - local({ - val <- list( - projection = "", - parameters = NULL, - orientation = NULL, - centre = NULL, - longlim = NULL, - latlim = NULL - ) - function(new) if(!missing(new)) val <<- new else val - }) magband=function(crosseq = 0, peaklat = 0, width=10, res=1000, ...){ longlo=.Last.magproj()$longlim[1] diff --git a/R/magbin.R b/R/magbin.R index 8e7237b..67abbfa 100644 --- a/R/magbin.R +++ b/R/magbin.R @@ -223,7 +223,7 @@ tempambig = output$nn.idx[distorder] remdupe = duplicated(tempambig, incomparables=NA) output$nn.idx[distorder][remdupe] = NA - output$output$nn.dists[distorder][remdupe] = NA + output$nn.dists[distorder][remdupe] = NA # if(exactcount){ # if(shape=='hex' | shape=='hexagon'){ @@ -356,30 +356,21 @@ plot.magbin = function(x, colramp=hcl.colors(21), colstretch='lin', sizestretch= ParmOff(magplot, c(list(NA, NA, xlim=xlim, ylim=ylim), dots)) } } - #magplot(NA, NA, xlim=x$xlim, ylim=x$ylim, ...) - for(i in 1:dim(x$bins)[1]){ - if(!is.na(colmap$map[i])){ - if(is.na(x$dustlim)){ - if(x$shape=='hex' | x$shape=='hexagon'){ - .drawhex(x$bins[i,'x'], x$bins[i,'y'], unitcell=x$step*sizemap[i], col=colramp[colmap$map[i]], border=NA, direction=x$direction) - } - if(x$shape=='sq' | x$shape=='square'){ - .drawsquare(x$bins[i,'x'], x$bins[i,'y'], unitcell=x$step*sizemap[i], col=colramp[colmap$map[i]], border=NA) - } - if(x$shape=='tri' | x$shape=='triangle' | x$shape=='trihex'){ - .drawtriangle(x$bins[i,'x'], x$bins[i,'y'], unitcell=x$step*sizemap[i], col=colramp[colmap$map[i]], border=NA, type=x$bins[i,'type'], direction=x$direction) - } - }else if(x$bins[i,'count'] > x$dustlim){ - if(x$shape=='hex' | x$shape=='hexagon'){ - .drawhex(x$bins[i,'x'], x$bins[i,'y'], unitcell=x$step*sizemap[i], col=colramp[colmap$map[i]], border=NA, direction=x$direction) - } - if(x$shape=='sq' | x$shape=='square'){ - .drawsquare(x$bins[i,'x'], x$bins[i,'y'], unitcell=x$step*sizemap[i], col=colramp[colmap$map[i]], border=NA) - } - if(x$shape=='tri' | x$shape=='triangle' | x$shape=='trihex'){ - .drawtriangle(x$bins[i,'x'], x$bins[i,'y'], unitcell=x$step*sizemap[i], col=colramp[colmap$map[i]], border=NA, type=x$bins[i,'type'], direction=x$direction) - } - } + # Determine the drawing function once outside the loop (shape is constant). + # After the dustlim filter above, all remaining bins already satisfy count > dustlim, + # so we only need to skip bins where the colour map produced NA. + draw_bins <- which(!is.na(colmap$map)) + if(x$shape == 'hex' | x$shape == 'hexagon'){ + for(i in draw_bins){ + .drawhex(x$bins[i,'x'], x$bins[i,'y'], unitcell=x$step*sizemap[i], col=colramp[colmap$map[i]], border=NA, direction=x$direction) + } + }else if(x$shape == 'sq' | x$shape == 'square'){ + for(i in draw_bins){ + .drawsquare(x$bins[i,'x'], x$bins[i,'y'], unitcell=x$step*sizemap[i], col=colramp[colmap$map[i]], border=NA) + } + }else if(x$shape == 'tri' | x$shape == 'triangle' | x$shape == 'trihex'){ + for(i in draw_bins){ + .drawtriangle(x$bins[i,'x'], x$bins[i,'y'], unitcell=x$step*sizemap[i], col=colramp[colmap$map[i]], border=NA, type=x$bins[i,'type'], direction=x$direction) } } if(!is.null(x$dust)){ diff --git a/R/magmap.R b/R/magmap.R index d296b6e..ef987a0 100644 --- a/R/magmap.R +++ b/R/magmap.R @@ -50,7 +50,7 @@ magmap = function( data, } else if (type == 'rank') { locut = 1 hicut = length(data[good]) - data[good][order(data[good])] = locut:hicut + data[good] <- rank(data[good], ties.method = 'average') } else{ stop(type, ' is not a valid type option!') } diff --git a/man/magclip.Rd b/man/magclip.Rd index a32ae53..f18b2a1 100644 --- a/man/magclip.Rd +++ b/man/magclip.Rd @@ -5,7 +5,7 @@ Magical sigma clipping } \description{ -This function does intelligent autoamtic sigma-clipping of data. This is optionally used by \code{\link{magplot}} and \code{\link{maghist}}. +This function does intelligent automatic sigma-clipping of data. This is optionally used by \code{\link{magplot}} and \code{\link{maghist}}. } \usage{ magclip(x, sigma = 'auto', clipiters = 5, sigmasel = 1, estimate = 'both', extra = TRUE) @@ -19,7 +19,7 @@ Numeric; the values to be clipped. This can reasonably be a vector, a matrix or The level of sigma clipping to be done. If set to default of 'auto' it will dynamically choose a sigma level to cut at based on the length of x (or the clipped version once iterations have started), i.e.: sigma=qnorm(1-2/length(x)). This have the effect of removing unlikely values based on the chance of them occurring, i.e. there is roughly a 50\% chance of a 3.5 / 4.6 sigma Normal fluctuation occurring when you have 10,000 / 10,000,000 values, hence choosing this value dynamically is usually the best option. } \item{clipiters}{ -The maximum number of sigma clipping iterations to attempt. It will break out sooner than this if the iterations have converged. The default of 5 is usually plenty (up to the contamination being towards the 50\% level). The number of actual iterations is returns as \option{clipiters}. +The maximum number of sigma clipping iterations to attempt. It will break out sooner than this if the iterations have converged. The default of 5 is usually plenty (up to the contamination being towards the 50\% level). The number of actual iterations is returned as \option{clipiters}. } \item{sigmasel}{ The quantile to use when trying to estimate the true standard-deviation of the Normal distribution. if contamination is low then the default of 1 is about optimal in terms of S/N, but you might need to make the value lower when contamination is very high. @@ -28,17 +28,17 @@ The quantile to use when trying to estimate the true standard-deviation of the N Character; determines which side/s of the distribution are used to estimate Normal properties. The default is to use both sides (both) giving better S/N, but if you know that your contamination only comes from positive flux sources (e.g., astronomical data when trying to select sky pixels) then you should only use the lower side to determine Normal statistics (lo). Similarly if the contamination is on the low side then you should use the higher side to determine Normal statistics (hi). } \item{extra}{ -Logical; if TRUE then \option{clip} and \option{range} are computed and returns, else these are set to NA to reduce computation and memory. +Logical; if TRUE then \option{clip} and \option{range} are computed and returned, else these are set to NA to reduce computation and memory. } } \details{ -If you know more sepcific details about your data then you should probably carry out a thorough likelihood analysis, but the ad-hoc clipping done in \code{magclip} works pretty well in practice. +If you know more specific details about your data then you should probably carry out a thorough likelihood analysis, but the ad-hoc clipping done in \code{magclip} works pretty well in practice. } \value{ A list containing three items: -\item{x}{Numeric vector; the cliped \option{x} values.} -\item{clip}{Locial; logic of which values were clipped with the same type and shape attributes as the input \option{x} (i.e. if the original \option{x} was a matrix then \option{clip} would also be a matrix that matches element to element).} +\item{x}{Numeric vector; the clipped \option{x} values.} +\item{clip}{Logical; logic of which values were clipped with the same type and shape attributes as the input \option{x} (i.e. if the original \option{x} was a matrix then \option{clip} would also be a matrix that matches element to element).} \item{range}{The data range of clipped \option{x} values returned.} \item{clipiters}{The number of iterations made, which might be less than the input \option{clipiters} since it might converge faster.} } diff --git a/man/magmap.Rd b/man/magmap.Rd index a82a643..ad94a4d 100644 --- a/man/magmap.Rd +++ b/man/magmap.Rd @@ -5,7 +5,7 @@ Value remapper } \description{ -This function allows the use to remap a vector of values onto a different system. For instance you might have values stretching from -10 to 100 which you want mapped from 0 to 2/3 so you can then sue the output as an input for point colour or size. It allows clipping of values, rejection of bad values, and log stretching. +This function allows the use to remap a vector of values onto a different system. For instance you might have values stretching from -10 to 100 which you want mapped from 0 to 2/3 so you can then use the output as an input for point colour or size. It allows clipping of values, rejection of bad values, and log stretching. } \usage{ magmap(data, locut = 0, hicut = 1, flip = FALSE, range = c(0, 2/3), type = "quan", diff --git a/man/magrun.Rd b/man/magrun.Rd index cb778e0..7abb2c8 100644 --- a/man/magrun.Rd +++ b/man/magrun.Rd @@ -71,10 +71,10 @@ The standard deviations in the x bins. This is a two column data.frame if 'diff' \item{ysd}{ The standard deviations in the y bins. This is a two column data.frame if 'diff' is set to FALSE, giving the y-sd and y+sd values, or a single vector if 'diff' is set to TRUE. If Nscale is set to TRUE then this is also divided by sqrt the contributing objects in each bin. } - \item{bincen}{ + \item{bincens}{ The bin centres used in the chosen binning direction. } - \item{binlim}{ + \item{binlims}{ The bin limits used in the chosen binning direction. } \item{Nbins}{ diff --git a/tests/testthat/test-functions.R b/tests/testthat/test-functions.R new file mode 100644 index 0000000..45616a3 --- /dev/null +++ b/tests/testthat/test-functions.R @@ -0,0 +1,326 @@ +## Tests for magclip, maghist, magrun, magerr, magcurve, magcon, magbar, +## magcutout, and magmap type='rank'. +## +## Each group covers: +## 1. Correct return values / output structure. +## 2. That edge-case or previously-buggy paths work without error. + +library(magicaxis) + +options(rgl.useNULL = TRUE) + +with_null_dev <- function(expr) { + pdf(nullfile()) + on.exit(dev.off()) + force(expr) +} + +# ── magclip ─────────────────────────────────────────────────────────────────── + +test_that("magclip returns a list with x, clip, range, clipiters", { + set.seed(1) + result <- magclip(c(rnorm(200), runif(20, 5, 10))) + expect_true(is.list(result)) + expect_true(all(c("x", "clip", "range", "clipiters") %in% names(result))) +}) + +test_that("magclip with extra=FALSE returns NA for clip and range", { + result <- magclip(rnorm(100), extra = FALSE) + expect_true(is.na(result$clip)) + expect_true(is.na(result$range)) +}) + +test_that("magclip removes obvious outliers", { + set.seed(2) + x <- c(rnorm(500), 100, -100) + result <- magclip(x, sigma = 3) + expect_true(max(result$x) < 50) + expect_true(min(result$x) > -50) +}) + +test_that("magclip estimate='lo' and estimate='hi' work without error", { + set.seed(3) + x <- c(rnorm(300), runif(30, 4, 8)) + expect_no_error(magclip(x, estimate = 'lo')) + expect_no_error(magclip(x, estimate = 'hi')) +}) + +test_that("magclip with sigma='auto' and very small input does not error", { + expect_no_error(magclip(c(1, 2, 3))) +}) + +test_that("magclip clip vector has same length as input x", { + set.seed(4) + x <- rnorm(100) + result <- magclip(x) + expect_equal(length(result$clip), length(x)) +}) + +# ── magmap type='rank' (previously broken) ───────────────────────────────── + +test_that("magmap type='rank' actually reorders data (bug fix)", { + set.seed(5) + d <- c(3, 1, 4, 1, 5, 9, 2, 6) + result <- magmap(d, type = 'rank', range = c(1, 8)) + # After ranking, the minimum ranked value should map near range[1] + # and the maximum near range[2]. + expect_true(min(result$map, na.rm = TRUE) <= 1.5) + expect_true(max(result$map, na.rm = TRUE) >= 7.5) + # The mapping must actually differ across elements (was all identical before fix) + expect_true(length(unique(result$map)) > 1) +}) + +test_that("magmap type='rank' with flip=TRUE inverts the mapping", { + d <- 1:10 + normal <- magmap(d, type = 'rank', range = c(0, 1))$map + flipped <- magmap(d, type = 'rank', range = c(0, 1), flip = TRUE)$map + expect_true(cor(normal, flipped) < -0.9) +}) + +# ── maghist ─────────────────────────────────────────────────────────────────── + +test_that("maghist returns a histogram object", { + result <- maghist(rnorm(500), plot = FALSE, verbose = FALSE) + expect_true(inherits(result, "histogram")) +}) + +test_that("maghist appends 'summary' and 'ranges' to the histogram", { + result <- maghist(rnorm(200), plot = FALSE, verbose = FALSE) + expect_true("summary" %in% names(result)) + expect_true("ranges" %in% names(result)) +}) + +test_that("maghist with xlim scalar performs sigma clipping without error", { + result <- maghist(c(rnorm(500), 100), xlim = 3, plot = FALSE, verbose = FALSE) + expect_true(inherits(result, "histogram")) +}) + +test_that("maghist with log='x' works on positive data", { + result <- maghist(10^runif(300, 0, 3), log = 'x', plot = FALSE, verbose = FALSE) + expect_true(inherits(result, "histogram")) +}) + +test_that("maghist cumsum=TRUE gives non-decreasing counts", { + result <- maghist(rnorm(300), cumsum = TRUE, plot = FALSE, verbose = FALSE) + expect_true(all(diff(result$counts) >= 0)) +}) + +test_that("maghist renders without error", { + with_null_dev({ + expect_no_error(maghist(rnorm(300), verbose = FALSE)) + }) +}) + +# ── magrun ─────────────────────────────────────────────────────────────────── + +test_that("magrun returns a list with correct element names", { + set.seed(6) + xy <- cbind(seq(0, 2, len = 200), rnorm(200)) + result <- magrun(xy) + expect_true(is.list(result)) + # Check the names that the code actually returns + expect_true(all(c("x", "y", "xquan", "yquan", "xsd", "ysd", + "bincens", "binlims", "Nbins") %in% names(result))) +}) + +test_that("magrun Nbins sums match input count", { + set.seed(7) + n <- 1000 + xy <- cbind(runif(n), runif(n)) + result <- magrun(xy, bins = 5, equalN = TRUE) + expect_equal(sum(result$Nbins), n) +}) + +test_that("magrun with type='mean' works", { + set.seed(8) + xy <- cbind(seq(0, 1, len = 100), rnorm(100)) + expect_no_error(magrun(xy, type = 'mean')) +}) + +test_that("magrun with custom breaks works", { + set.seed(9) + xy <- cbind(seq(0, 2, len = 200), rnorm(200)) + result <- magrun(xy, bins = c(0, 0.5, 1.0, 1.5, 2.0)) + expect_equal(length(result$x), 4) +}) + +test_that("magrun Nscale=TRUE reduces spread relative to Nscale=FALSE", { + set.seed(10) + xy <- cbind(seq(0, 1, len = 500), rnorm(500)) + r_noscale <- magrun(xy, bins = 5, Nscale = FALSE) + r_scale <- magrun(xy, bins = 5, Nscale = TRUE) + # Nscale divides by sqrt(N), so quantile widths should be narrower + width_no <- mean(r_noscale$yquan[, 2] - r_noscale$yquan[, 1]) + width_sc <- mean(r_scale$yquan[, 2] - r_scale$yquan[, 1]) + expect_true(width_sc < width_no) +}) + +# ── magerr ─────────────────────────────────────────────────────────────────── + +test_that("magerr draws y error bars without error", { + with_null_dev({ + magplot(1:5, 1:5) + expect_no_error(magerr(x = 1:5, y = 1:5, ylo = 0.2, yhi = 0.3)) + }) +}) + +test_that("magerr draws x error bars without error", { + with_null_dev({ + magplot(1:5, 1:5) + expect_no_error(magerr(x = 1:5, y = 1:5, xlo = 0.1)) + }) +}) + +test_that("magerr draws both x and y error bars without error", { + with_null_dev({ + magplot(1:10, 1:10) + expect_no_error(magerr(x = 1:10, y = 1:10, xlo = 0.15, ylo = 0.25)) + }) +}) + +test_that("magerr with poly=TRUE draws error polygon without error", { + with_null_dev({ + magplot(1:5, 1:5) + expect_no_error( + magerr(x = 1:5, y = 1:5, ylo = rep(0.2, 5), yhi = rep(0.3, 5), + poly = TRUE, col = 'lightgrey', border = NA) + ) + }) +}) + +test_that("magerr on log-x axis works without error", { + with_null_dev({ + magplot(10^(1:5), 1:5, log = 'x') + expect_no_error( + magerr(x = 10^(1:5), y = 1:5, xlo = 10^(1:5) * 0.2) + ) + }) +}) + +test_that("magerr accepts x as a two-column matrix", { + with_null_dev({ + magplot(1:5, 1:5) + xy <- cbind(1:5, 1:5) + expect_no_error(magerr(x = xy, ylo = 0.2)) + }) +}) + +# ── magcurve ───────────────────────────────────────────────────────────────── + +test_that("magcurve returns a list with x and y", { + with_null_dev({ + result <- magcurve(sin(x), from = 0, to = 2 * pi) + expect_true(is.list(result)) + expect_equal(names(result), c("x", "y")) + expect_equal(length(result$x), 101) + }) +}) + +test_that("magcurve add=TRUE adds to an existing plot without error", { + with_null_dev({ + magplot(0, 0, xlim = c(0, 6), ylim = c(-1.5, 1.5)) + expect_no_error(magcurve(cos(x), from = 0, to = 2 * pi, add = TRUE)) + }) +}) + +test_that("magcurve with log x axis works on positive domain", { + with_null_dev({ + expect_no_error(magcurve(log10(x), from = 1, to = 1000, log = 'x')) + }) +}) + +# ── magcon ─────────────────────────────────────────────────────────────────── + +test_that("magcon returns a list with x, y, z components", { + with_null_dev({ + set.seed(11) + x <- rnorm(300); y <- x + rnorm(300, sd = 0.5) + result <- magcon(x, y, h = c(0.3, 0.3), dobar = FALSE) + expect_true(is.list(result)) + expect_true(all(c("x", "y", "z") %in% names(result))) + }) +}) + +test_that("magcon doim=FALSE, docon=TRUE works without error", { + with_null_dev({ + set.seed(12) + x <- rnorm(200); y <- rnorm(200) + expect_no_error(magcon(x, y, h = c(0.3, 0.3), doim = FALSE, dobar = FALSE)) + }) +}) + +test_that("magcon with add=TRUE adds to existing plot", { + with_null_dev({ + set.seed(13) + x <- rnorm(200); y <- rnorm(200) + magplot(x, y) + expect_no_error(magcon(x, y, h = c(0.3, 0.3), add = TRUE, dobar = FALSE)) + }) +}) + +# ── magbar ──────────────────────────────────────────────────────────────────── + +test_that("magbar renders a vertical colour bar without error", { + with_null_dev({ + magplot(1, 1, xlim = c(0, 10), ylim = c(0, 10)) + expect_no_error(magbar(range = c(0, 100), orient = 'v')) + }) +}) + +test_that("magbar renders a horizontal colour bar without error", { + with_null_dev({ + magplot(1, 1, xlim = c(0, 10), ylim = c(0, 10)) + expect_no_error(magbar(range = c(0, 100), orient = 'h')) + }) +}) + +test_that("magbar with log=TRUE works without error", { + with_null_dev({ + magplot(1, 1, xlim = c(0, 10), ylim = c(0, 10)) + expect_no_error(magbar(range = c(1, 1000), log = TRUE)) + }) +}) + +test_that("magbar position argument variants work without error", { + with_null_dev({ + for (pos in c('topright', 'topleft', 'bottomright', 'bottomleft')) { + magplot(1, 1, xlim = c(0, 10), ylim = c(0, 10)) + expect_no_error(magbar(position = pos, range = c(0, 1))) + } + }) +}) + +# ── magcutout ──────────────────────────────────────────────────────────────── + +test_that("magcutout extracts correct submatrix from interior", { + m <- matrix(1:100, 10, 10) + result <- magcutout(m, loc = c(5, 5), box = c(4, 4)) + expect_equal(dim(result$image), c(4, 4)) +}) + +test_that("magcutout returns loc, loc.orig, loc.diff", { + m <- matrix(1:100, 10, 10) + result <- magcutout(m, loc = c(5, 5), box = c(4, 4)) + expect_true(all(c("image", "loc", "loc.orig", "loc.diff", "xsel", "ysel") %in% + names(result))) +}) + +test_that("magcutout with paddim=TRUE pads near-edge cutouts to requested box size", { + m <- matrix(1:100, 10, 10) + result <- magcutout(m, loc = c(1, 1), box = c(6, 6), paddim = TRUE) + expect_equal(dim(result$image), c(6, 6)) +}) + +test_that("magcutout with paddim=FALSE clips to available data", { + m <- matrix(1:100, 10, 10) + result <- magcutout(m, loc = c(1, 1), box = c(6, 6), paddim = FALSE) + # Should not exceed matrix boundaries + expect_true(all(dim(result$image) <= c(6, 6))) +}) + +test_that("magcutout plots without error", { + with_null_dev({ + m <- matrix(rnorm(400), 20, 20) + expect_no_error(magcutout(m, loc = c(10, 10), box = c(8, 8), plot = TRUE)) + }) +}) diff --git a/vignettes/magicaxis.Rmd b/vignettes/magicaxis.Rmd index 04b104b..c5f08bf 100644 --- a/vignettes/magicaxis.Rmd +++ b/vignettes/magicaxis.Rmd @@ -132,7 +132,164 @@ magplot(h, log = "y", xlab = "x (log)", ylab = "Count (log)") --- -## magcon – 2-D density contours +## magclip – automatic sigma-clipping + +`magclip()` iteratively removes outliers that are unlikely given a Normal +distribution. It is used internally by `magplot()`, `maghist()`, and `magbin()` +when a scalar `xlim`/`ylim` is supplied, but it is also useful on its own: + +```{r magclip-basic} +set.seed(42) +x <- c(rnorm(500), runif(30, 5, 10)) # 500 Normal + 30 outliers +result <- magclip(x, sigma = 3) +cat("Before clipping:", length(x), "points\n") +cat("After clipping:", length(result$x), "points\n") +cat("Clipped range:", round(result$range, 2), "\n") +``` + +The `estimate` argument controls which tail is used to estimate the Normal +dispersion. Use `'lo'` when contamination is on the high side only: + +```{r magclip-lo} +set.seed(43) +x_hi <- c(rnorm(500), runif(40, 4, 8)) +r_both <- magclip(x_hi) +r_lo <- magclip(x_hi, estimate = 'lo') +cat("Both sides kept:", length(r_both$x), "| Low-side only kept:", length(r_lo$x), "\n") +``` + +--- + +## magrun – running statistics + +`magrun()` computes running medians (or means/modes) and scatter measures along +a binned axis, ideal for overlaying a trend on noisy scatter data: + +```{r magrun-basic} +set.seed(44) +n <- 1000 +xx <- seq(0, 2, len = n) +yy <- xx + rnorm(n) +magplot(xx, yy, col = "lightgrey", pch = ".", xlab = "x", ylab = "y") +run <- magrun(cbind(xx, yy), bins = 10) +lines(run$x, run$y, col = "red", lwd = 2) +lines(run$x, run$yquan[, 1], lty = 2, col = "red") +lines(run$x, run$yquan[, 2], lty = 2, col = "red") +``` + +The `Nscale = TRUE` flag divides the scatter by √N, giving error-in-the-mean +intervals (useful for assessing the significance of a trend): + +```{r magrun-Nscale} +set.seed(45) +n <- 500 +xx <- seq(0, 1, len = n) +yy <- xx + rnorm(n) +run_sd <- magrun(cbind(xx, yy), bins = 8, Nscale = FALSE, diff = TRUE) +run_sem <- magrun(cbind(xx, yy), bins = 8, Nscale = TRUE, diff = TRUE) +magplot(xx, yy, col = "lightgrey", pch = ".", xlab = "x", ylab = "y") +magerr(run_sd$x, run_sd$y, ylo = run_sd$yquan[, 1], yhi = run_sd$yquan[, 2], + col = "royalblue", length = 0, lty = 2) +magerr(run_sem$x, run_sem$y, ylo = run_sem$yquan[, 1], yhi = run_sem$yquan[, 2], + col = "red") +legend("topleft", legend = c("scatter (1σ)", "SEM (1σ)"), + col = c("royalblue", "red"), lty = c(2, 1), bty = "n") +``` + +--- + +## magerr – error bars + +`magerr()` adds symmetric or asymmetric error bars in x and/or y to an +existing plot. Lower errors are supplied via `xlo`/`ylo` and upper errors +via `xhi`/`yhi` (defaulting to the lower value if omitted): + +```{r magerr-basic} +set.seed(46) +n <- 12 +x <- sort(runif(n)) +y <- 2 * x + rnorm(n, sd = 0.15) +magplot(x, y, xlab = "x", ylab = "y") +magerr(x, y, xlo = 0.04, ylo = 0.1, yhi = 0.2) +``` + +Error ellipses are drawn when `corxy` is supplied: + +```{r magerr-ellipse} +set.seed(47) +n <- 8 +x <- runif(n, 0.1, 0.9) +y <- runif(n, 0.1, 0.9) +sx <- runif(n, 0.03, 0.08) +sy <- runif(n, 0.05, 0.12) +rho <- runif(n, -0.7, 0.7) +magplot(x, y, xlab = "x", ylab = "y", pch = 16) +magerr(x, y, xlo = sx, ylo = sy, corxy = rho) +``` + +A shaded error polygon is produced with `poly = TRUE`: + +```{r magerr-poly} +set.seed(48) +x <- seq(0, 2 * pi, len = 30) +y <- sin(x) + rnorm(30, sd = 0.1) +magplot(x, y, type = "l", xlab = "x", ylab = "sin(x)") +magerr(x, y, ylo = 0.15, yhi = 0.15, poly = TRUE, col = adjustcolor("steelblue", 0.3), + border = NA) +``` + +--- + +## magcurve – draw a function curve + +`magcurve()` evaluates and plots any R expression over a range. It is a +drop-in replacement for `curve()` that uses `magplot()` for its base plot, +giving pretty axes automatically: + +```{r magcurve-basic} +magcurve(sin(x), from = 0, to = 2 * pi, xlab = "x", ylab = "sin(x)") +magcurve(cos(x), add = TRUE, col = "steelblue") +legend("topright", legend = c("sin", "cos"), col = c("black", "steelblue"), + lty = 1, bty = "n") +``` + +Log axes are supported: + +```{r magcurve-log} +magcurve(x^2, from = 1, to = 100, log = "xy", + xlab = "x (log)", ylab = expression(x^2 ~ "(log)")) +``` + +--- + +## magbar – colour bar legend + +`magbar()` adds a standalone colour bar to any existing plot. The `position` +argument accepts compass-style strings (`'topright'`, `'bottomleft'`, etc.); +`orient` switches between vertical (`'v'`) and horizontal (`'h'`): + +```{r magbar-basic} +n <- 300 +set.seed(49) +x <- rnorm(n); y <- rnorm(n); z <- x^2 + y^2 +colramp <- hcl.colors(21, "YlOrRd") +colvals <- magmap(z, range = c(1, length(colramp)))$map +magplot(x, y, col = colramp[round(colvals)], pch = 16, + xlab = "x", ylab = "y") +magbar(position = "topright", range = range(z), col = colramp, + orient = "v", title = expression(x^2 + y^2)) +``` + +A horizontal bar at the bottom: + +```{r magbar-horizontal} +magplot(x, y, col = colramp[round(colvals)], pch = 16, + xlab = "x", ylab = "y") +magbar(position = "bottomleft", range = range(z), col = colramp, + orient = "h", title = "z") +``` + +--- `magcon()` computes a 2-D kernel density estimate and draws contours at chosen probability levels (default 50%, 68%, 95%) together with an optional image From d3da9a69bfeac107928bc250c546a8c9354c10a7 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Sat, 16 May 2026 14:50:49 +0000 Subject: [PATCH 2/3] fix: restore magcon section heading in vignette Agent-Logs-Url: https://github.com/asgr/magicaxis/sessions/74e9c01b-f592-4d54-8e81-353d57d8b91c Co-authored-by: asgr <5617132+asgr@users.noreply.github.com> --- vignettes/magicaxis.Rmd | 2 ++ 1 file changed, 2 insertions(+) diff --git a/vignettes/magicaxis.Rmd b/vignettes/magicaxis.Rmd index c5f08bf..24cd233 100644 --- a/vignettes/magicaxis.Rmd +++ b/vignettes/magicaxis.Rmd @@ -291,6 +291,8 @@ magbar(position = "bottomleft", range = range(z), col = colramp, --- +## magcon – 2-D density contours + `magcon()` computes a 2-D kernel density estimate and draws contours at chosen probability levels (default 50%, 68%, 95%) together with an optional image and colour bar. From fb1e01fbac0d6ae7a1e3df7daab7a4d5bc6afafe Mon Sep 17 00:00:00 2001 From: Aaron Robotham Date: Mon, 18 May 2026 16:12:41 +0800 Subject: [PATCH 3/3] Some tweaks to up date colour scheme, simplify code. Created a magcorner alias for magtri too (since lots of people call them corner plots). --- R/magcon.R | 14 ++++++++------ R/magtri.R | 31 +++++++++++++++++++++++-------- man/magcon.Rd | 2 +- man/magtri.Rd | 22 ++++++++++++++-------- 4 files changed, 46 insertions(+), 23 deletions(-) diff --git a/R/magcon.R b/R/magcon.R index 958681e..95bbd6a 100644 --- a/R/magcon.R +++ b/R/magcon.R @@ -1,5 +1,7 @@ -magcon <- -function(x,y,h,doim=TRUE,docon=TRUE,dobar=TRUE,ngrid=100,add=FALSE,xlab='',ylab='',imcol=c(NA,rev(rainbow(1000,start=0,end=2/3))),conlevels=c(0.5,pnorm(1)-pnorm(-1),0.95), barposition='topright', barorient='v',bartitle='Contained %',bartitleshift=0,xlim=NULL,ylim=NULL,weights=NULL,...){ +magcon = function(x,y,h,doim=TRUE,docon=TRUE,dobar=TRUE,ngrid=100,add=FALSE,xlab='',ylab='', + imcol=c(NA,hcl.colors(100)), conlevels=c(0.5,pnorm(1)-pnorm(-1),0.95), + barposition='topright', barorient='v',bartitle='Contained %',bartitleshift=0, + xlim=NULL,ylim=NULL,weights=NULL,...){ if(missing(y)){ if(!is.null(dim(x))){ if(dim(x)[2]>=2){y=x[,2];x=x[,1]} @@ -26,14 +28,14 @@ function(x,y,h,doim=TRUE,docon=TRUE,dobar=TRUE,ngrid=100,add=FALSE,xlab='',ylab= tempcon=sm.density(cbind(x,y),h=h,weights=weights,display='none',ngrid=ngrid,xlim=xlim+c(-diff(xlim),diff(xlim)),ylim=ylim+c(-diff(ylim),diff(ylim)),verbose=FALSE) tempcon$x=tempcon$eval.points[,1] tempcon$y=tempcon$eval.points[,2] - tempcon$z=tempcon$estimate + tempcon$z=tempcon$estimate / sum(tempcon$estimate, na.rm=TRUE) #} temp=sort(tempcon$z) tempsum=cumsum(temp) - convfunc=approxfun(tempsum,temp) - levelmap=approxfun(convfunc(seq(0,1,len=1000)*max(tempsum)),seq(0,1,len=1000)) - tempcon$z=matrix(levelmap(tempcon$z),nrow=ngrid) + convfunc=approxfun(temp, tempsum) + #levelmap=approxfun(convfunc(seq(0,1,len=1000)),seq(0,1,len=1000)) + tempcon$z=matrix(convfunc(tempcon$z),nrow=ngrid) tempcon$z[is.na(tempcon$z)]=min(tempcon$z,na.rm=TRUE) if(doim){ if(add==FALSE){ diff --git a/R/magtri.R b/R/magtri.R index 434634d..ab001c2 100644 --- a/R/magtri.R +++ b/R/magtri.R @@ -1,4 +1,4 @@ -magtri=function(chains, samples=1000, thin=1, samptype='end', grid=FALSE, do.tick=FALSE, refvals=NULL, lab=NULL, ...){ +magtri = function(chains, samples=1000, thin=1, samptype='end', grid=FALSE, do.tick=FALSE, refvals=NULL, lab=NULL, draw.sig=FALSE, ...){ chains=as.data.frame(chains) chaincolnames=colnames(chains) Nsamp=dim(chains)[1] @@ -71,17 +71,22 @@ magtri=function(chains, samples=1000, thin=1, samptype='end', grid=FALSE, do.tic if(i>j){ plot.new() plot.window(xlim=xrange,ylim=yrange) - xtemp=chains[usesamps,i] - ytemp=chains[usesamps,j] + xtemp = chains[usesamps,i] + ytemp = chains[usesamps,j] if(sd(xtemp)==0){xtemp=xtemp+rnorm(samples,sd=1e-3)} if(sd(ytemp)==0){ytemp=ytemp+rnorm(samples,sd=1e-3)} ParmOff(magaxis, dots, side=1:2, grid=grid, grid.col='lightgrey', labels=FALSE, do.tick=do.tick, .pass_dots=FALSE) - ParmOff(magcon, dots, x=xtemp, y=ytemp, dobar=FALSE, doim=FALSE, add=TRUE, lty=c(2,1,3), xlim=xrange, ylim=yrange, h=c(diff(xrange),diff(yrange))/10, .pass_dots=FALSE) - ParmOff(points, dots, x=meanvec[i], y=meanvec[j], col='red', pch=4, cex=2, .pass_dots=FALSE) + ParmOff(magcon, dots, x=xtemp, y=ytemp, dobar=FALSE, doim=FALSE, add=TRUE, lty=c(2,1,3), xlim=xrange, ylim=yrange, h=c(sdvec[i], sdvec[j])/5, .pass_dots=TRUE, .clash='last') + ParmOff(points.default, dots, x=meanvec[i], y=meanvec[j], col='red', pch=4, cex=2, .pass_dots=TRUE, .clash = 'last') box() - abline(v=meanvec[i],lty=1,col='red') - abline(v=meanvec[i]-sdvec[i],lty=3,col='red') - abline(v=meanvec[i]+sdvec[i],lty=3,col='red') + if(draw.sig){ + abline(v=meanvec[i],lty=1,col='red') + abline(v=meanvec[i]-sdvec[i],lty=3,col='red') + abline(v=meanvec[i]+sdvec[i],lty=3,col='red') + abline(h=meanvec[j],lty=1,col='darkgreen') + abline(h=meanvec[j]-sdvec[j],lty=3,col='darkgreen') + abline(h=meanvec[j]+sdvec[j],lty=3,col='darkgreen') + } if(!is.null(refvals)){ abline(v=refvals[i],lty=1, col='blue') } @@ -99,6 +104,14 @@ magtri=function(chains, samples=1000, thin=1, samptype='end', grid=FALSE, do.tic ParmOff(points.default, dots_carry, x=chains[usesamps,i], y=chains[usesamps,j], pch='.', col='darkgrey', .pass_dots=TRUE) ParmOff(points.default, dots_carry, x=meanvec[i], y=meanvec[j], col='red', pch=4, cex=2, .pass_dots=TRUE) box() + if(draw.sig){ + abline(v=meanvec[i],lty=1,col='red') + abline(v=meanvec[i]-sdvec[i],lty=3,col='red') + abline(v=meanvec[i]+sdvec[i],lty=3,col='red') + abline(h=meanvec[j],lty=1,col='darkgreen') + abline(h=meanvec[j]-sdvec[j],lty=3,col='darkgreen') + abline(h=meanvec[j]+sdvec[j],lty=3,col='darkgreen') + } if(i==1){ if(is.null(lab)){ ParmOff(magaxis, dots, side=2, ylab=chaincolnames[j], .pass_dots=FALSE) @@ -114,3 +127,5 @@ magtri=function(chains, samples=1000, thin=1, samptype='end', grid=FALSE, do.tic rownames(output)=chaincolnames return(invisible(output)) } + +magcorner = magtri \ No newline at end of file diff --git a/man/magcon.Rd b/man/magcon.Rd index 20465a4..5b3df30 100644 --- a/man/magcon.Rd +++ b/man/magcon.Rd @@ -9,7 +9,7 @@ This function generates pretty images and contours that reflect the 2D quantile } \usage{ magcon(x, y, h, doim = TRUE, docon = TRUE, dobar = TRUE, ngrid = 100, add = FALSE, - xlab = '', ylab = '', imcol = c(NA,rev(rainbow(1000, start = 0, end = 2/3))), + xlab = '', ylab = '', imcol = c(NA, hcl.colors(100)), conlevels = c(0.5, pnorm(1) - pnorm(-1), 0.95), barposition = "topright", barorient = "v",bartitle = "Contained \%", bartitleshift = 0, xlim = NULL, ylim = NULL, weights = NULL,...) diff --git a/man/magtri.Rd b/man/magtri.Rd index 0794ffa..9cf6beb 100644 --- a/man/magtri.Rd +++ b/man/magtri.Rd @@ -1,15 +1,18 @@ \name{magtri} \alias{magtri} +\alias{magcorner} %- Also NEED an '\alias' for EACH other topic documented here. \title{ -High level triangle plotting code for MCMC chains +High level triangle (corner) plotting code for MCMC chains } \description{ -A very high level (minimal options) MCMC chain triangle (AKA corner) plot function. The default is deliberately spartan in terms of options, but the result should be a clear set of covariance plots that should give quick insight into the stationary sampling quality of a set of MCMC posterior chains. +A very high level (minimal options) MCMC chain triangle (AKA corner) plot function. The default is deliberately spartan in terms of options, but the result should be a clear set of covariance plots that should give quick insight into the stationary sampling quality of a set of MCMC posterior chains. \code{magcorner} is just an alias to \code{magtri}, since these days such plots are more commonly called 'corner' plots, especially in Python-world. } \usage{ -magtri(chains, samples = 1000, thin = 1, samptype = 'end', grid = FALSE, do.tick = FALSE, -refvals = NULL, lab = NULL, ...) +magtri(chains, samples = 1000, thin = 1, samptype = 'end', grid = FALSE, + do.tick = FALSE, refvals = NULL, lab = NULL, draw.sig = FALSE, ...) +magcorner(chains, samples = 1000, thin = 1, samptype = 'end', grid = FALSE, + do.tick = FALSE, refvals = NULL, lab = NULL, draw.sig = FALSE, ...) } %- maybe also 'usage' for other objects documented here. \arguments{ @@ -36,19 +39,22 @@ Numeric vector; this gives reference values to overdraw. If provided it must be } \item{lab}{ Character vector; optional over-ride for column names when plotting the grid of scatter plots. +} + \item{draw.sig}{ +Logical; should vertical and horizontal sigma regions be overlaid on the scatter and contour panels. } \item{\dots}{ Extra arguments are passed to \option{\link{magcon}} for plotting. } } \details{ - This interface is deliberately very high level with few options. It is really designed to allow quick exploratory views of posterior samples from MCMC chains, and publication grade plots should be designed by the user. That said, in many situations the plots generated are of pleasant, clear and publishable quality. +This interface is deliberately very high level with few options. It is really designed to allow quick exploratory views of posterior samples from MCMC chains, and publication grade plots should be designed by the user. That said, in many situations the plots generated are of pleasant, clear and publishable quality. - Other types of data can be plotted using this function of course, but the default setup is tuned towards being useful for MCMC posterior chain samples. +Other types of data can be plotted using this function of course, but the default setup is tuned towards being useful for MCMC posterior chain samples. - The contour levels shown are the defaults for magcon, i.e. they contain 50\% (lty=2), 68\% (lty=1) and 95\% (lty=3) of the posterior chains. +The contour levels shown are the defaults for magcon, i.e. they contain 50\% (lty=2), 68\% (lty=1) and 95\% (lty=3) of the posterior chains. - The red cross shows the mean for the sampled posterior chain. The red vertical dashed line traces this over the contour plots. The red dotted line shows the +/- SD range of the sampled posterior chain. +The red cross shows the mean for the sampled posterior chain. The red vertical dashed line traces this over the contour plots. The red dotted line shows the +/- SD range of the sampled posterior chain. } \value{ Outputs a two column matrix containing the means and standard deviations fo the parameters. Generally run for the side effect of producing nice projected plots.