warn.tools_verbose <- FALSE cat('Running test scripts:\n') library(levelSets) pt_in_slice <- levelSets:::pt_in_slice inRect <- levelSets:::inRect flipRays <- levelSets:::flipRays pt_on_line2 <- levelSets:::pt_on_line2 warnBadResp <- levelSets:::warnBadResp pairIdxs <- levelSets:::pairIdxs linePosn0 <- levelSets:::linePosn0 linePosn1 <- levelSets:::linePosn1 lineVec <- levelSets:::lineVec sliceMatch <- levelSets:::sliceMatch internalOpt <- levelSets:::internalOpt inpPts <- levelSets:::inpPts inpPairs <- levelSets:::inpPairs pairInfo <- levelSets:::pairInfo respPts <- levelSets:::respPts num_matrix <- levelSets:::num_matrix naF <- levelSets:::naF cat('axisRays.r', sep='\n') # User-level functions to create standard sets of rays. # # axisRays # randomRays # flipRays (not public) oldopt <- lsetPkgOpt(reset=TRUE) #----- cat(" axisRays ... ") #== Testing #----- 1D example. spec <- fnSpec(inpnames="X1") scl <- inpScale(matrix(5, nrow=1, ncol=1), spec=spec) rays <- axisRays(spec, scale=scl) stopifnot(dim(rayVec(rays)) == c(1, 1), rayVec(rays) == sqrt(5)) # unit vector transformed to length sqrt(5) #----- 3D example to illustrate 'slice'ing. #layout(matrix(c(1:4), nrow=2, byrow=TRUE)) spec <- fnSpec(inpnames=c("X1", "X2", "X3")) slice <- c(NA, NA, -2) rays <- axisRays(spec, degree=2, slice=slice, origin=c(0, 0, -2)) # 4 rays #print(names(rays)) #plot(rays, eqAxes=TRUE) stopifnot(all(pt_in_slice(rayPosn1(rays), slicemat=rbind(slice), sliceidx=1, inptol=inpTol(spec))), length(rays) == 4) # Specifying 'rect' and putting the ray origin on its boundary can lead # to ray reversal or dropping rays: rect <- inpRect(rbind(c(-1, 0, -2), c(1, 2, 3)), refpt=c(1, 0, -2), # upper bdry of X1, lower of X2, X3 spec=spec) tools::assertWarning( rays <- axisRays(spec, degree=2, slice=slice, rect=rect) # origin on bdry # of 'rect' causes 'X1', 'X1(-)X2' rays to be reversed # and 'X1(+)X2' ray to be dropped , verbose=warn.tools_verbose) #print(names(rays)) stopifnot(all(pt_in_slice(rayPosn1(rays), slicemat=rbind(slice), sliceidx=1, inptol=inpTol(spec)))) stopifnot(all(inRect(rayPosn1(rays), rect=rect, bdryOnly=TRUE)), length(rays) == 3) #plot(rays, eqAxes=TRUE, main="Ray origin on boundary of 'rect'") # Can also specify 'scale', which is applied before ray vectors are # extended to the boundary of 'rect'. scl <- inpScale(rbind(c(1, 0.5, 0.5), c(0.5, 2, 0.5), c(0.5, 0.5, 4)), spec=spec) # unequal scales, non-zero corr rays <- axisRays(spec, degree=2, origin=c(1, 0, -2), slice=slice, scale=scl) #plot(rays, eqAxes=TRUE, main="Scaled and rotated rays ...") #plot(rect, add=TRUE) # Rays 1 and 4 and their reverses are outside 'rect': (chk <- inRect(rayPosn1(c(rays, flipRays(rays))), rect=rect)) stopifnot(identical(unname(chk), c(FALSE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE))) # So forcing boundary of 'rect' creates only 2 rays: tools::assertWarning( rays <- axisRays(spec, degree=2, slice=slice, rect=rect, scale=scl) , verbose=warn.tools_verbose) stopifnot(inherits(rays, "inpRays"), length(rays) == 2, inRect(rayPosn1(rays), rect=rect, bdryOnly=TRUE)) #plot(rays, eqAxes=TRUE, main="... confined to 'rect'") # Shows how two rays are extended to bdry of 'rect' #layout(matrix(1)) cat("OK\n") #----- cat(" randomRays ... ") #== Testing spec <- fnSpec(inpnames=c("X")) # one-dimensional parameter space set.seed(1) rays <- randomRays(5, spec=spec) # random +/-1's stopifnot(inherits(rays, "inpRays"), length(rays) == 5, all(rayPosn1(rays) %in% c(-1, 1))) posn1 <- rayPosn1(randomRays(5, spec=spec, origin=c(5))) # random 4's and 6's stopifnot(is.matrix(posn1), dim(posn1) == c(5, 1), all(posn1 %in% c(4, 6))) pRect <- inpRect(matrix(c(-1, 4), nrow=2), spec=spec) posn1b <- rayPosn1(randomRays(5, rect=pRect)) # 'spec' taken from 'rect' # random -1's and 4's stopifnot(length(posn1b) == 5, all(posn1b %in% c(-1, 4))) #== 2-D examples. #layout(matrix(1:4, nrow=2, byrow=TRUE)) #----- Rays pass through random points on circle centered at (3,2) with # radius 1: spec <- fnSpec(inpnames=c("A", "B")) set.seed(1) rslt <- randomRays(50, spec=spec, origin=c(3, 2)) #plot(rslt, eqAxes=TRUE, main="Random rays, unit length") # Effect of 'scale': # Ellipsoid with wider range for A axis relative to B (semi-major and # semi-minor axis lengths given by sqrt(diag(mat))). scl <- inpScale(matrix(c(5, 0, 0, 1), nrow=2), spec=spec) set.seed(1) rslt <- randomRays(50, origin=c(3, 2), scale=scl) #plot(rslt, eqAxes=TRUE, main="Ray lengths scaled to an ellipse") # Ellipsoid with wider range for A axis relative to B, and rotated: scl <- inpScale(matrix(c(5, 2, 2, 1), nrow=2), spec=spec) set.seed(1) rslt <- randomRays(50, origin=c(3, 2), scale=scl) #plot(rslt, eqAxes=TRUE, main="Rays scaled and rotated") # ... adjust lengths to put position 1.0 on the boundary of a rectangle: rect <- inpRect(rbind(c(1, 1), c(5, 2.5)), spec=spec) set.seed(1) rslt <- randomRays(50, origin=c(3, 2), scale=scl, rect=rect) #plot(rslt, eqAxes=TRUE, main="Ray position 1 on a 'rect' boundary") #-- When the sides of 'rect' have very different lengths, 'scale' makes # a big difference. #layout(matrix(c(1:3, 0), nrow=2, byrow=TRUE)) # To highlight the effect, create a long, narrow rect, and use an off-center # 'origin'. corners <- rbind(c(0.5, 3), c(5.5, 3.5)) rect <- inpRect(corners, spec=spec) # Default for 'scale' adjusts for the side lengths, so rays fill 'rect': set.seed(1) rslt <- randomRays(50, origin=c(1, 3.25), rect=rect) #plot(rslt, eqAxes=FALSE, main="Default 'scale' adjusts for 'rect' sides") # If we intentionally ignore differences in side lengths, rays may be sparse # in the long directions: set.seed(1) rslt <- randomRays(50, origin=c(1, 3.25), rect=rect, scale=inpScale(spec=spec)) #plot(rslt, eqAxes=FALSE, main="No adjustment for 'rect' sides") #-- When 'origin' is itself on the boundary of 'rect', direction vectors # that point outside 'rect' and lead to degenerate rays are instead reversed. set.seed(1) rslt <- randomRays(50, origin=c(2, 3.5), rect=rect) #plot(rslt, eqAxes=FALSE, main="Ray origin on boundary of 'rect'") #layout(matrix(1)) cat("OK\n") #----- cat(" flipRays ... ") #-- 1D example # Set up an input space. spec <- fnSpec(inpnames=c("a")) rays <- inpRays(rbind(0, 3), origin=c(1), spec=spec) flipped <- flipRays(rays) stopifnot(inherits(flipped, "inpRays"), dim(rayPosn1(flipped)) == c(2, 1), as.vector(rayPosn1(flipped)) == c(2, -1), as.vector(rayOrigin(flipped)) == c(1)) # Origin on rectangle boundary: pRect <- inpRect(rbind(-1, 2), spec=spec) rays <- inpRays(c(0), origin=c(2), rect=pRect, spec=spec) tools::assertWarning( flipped <- flipRays(rays), verbose=warn.tools_verbose) # omit 1 ray stopifnot(inherits(flipped, "inpRays"), length(flipped) == 0) #-- 3D examples spec <- fnSpec(inpnames=c("A", "B", "C")) rays <- inpRays(rbind(c(1, 1, 1), c(1, 2, 3)), origin=c(0, -1, 1), spec=spec) rays2 <- flipRays(rays) # Without a reference rectangle, flipRays() is its own inverse: stopifnot(isTRUE(all.equal(flipRays(rays2), rays))) # Edge case: rays2 <- flipRays(rays[0]) stopifnot(inherits(rays2, "inpRays"), length(rays2) == 0) # Include a rect. pRect <- inpRect(rbind(c(0,0,0), c(1,1,1)), spec=spec) rays <- inpRays(vecs=rbind(c(-1,-1,0), c(0,0,1), c(0, -1, 1)), origin=c(0.5, 0.4, 0.3), rect=pRect, spec=spec) frays <- flipRays(rays) # Flipped rays must have same origin as original, lie on the same lines, # and have position 1 on the boundary of 'pRect'. stopifnot(identical(rayOrigin(rays), rayOrigin(frays))) stopifnot(all(pt_on_line2(rayPosn1(frays), lineidx=seq_len(length(rays)), x0=rayOrigin(rays), x1=rayPosn1(rays), inptol=inpTol(rays)))) stopifnot(all(inRect(rayPosn1(frays), rect=pRect, bdryOnly=TRUE))) # Include a rect, and put the ray origin on its boundary: pRect <- inpRect(rbind(c(0,0,0), c(1,1,1)), refpt=c(1,1,0.5), spec=spec) rays <- inpRays(vecs=rbind(c(-1,-1,0), c(0,0,1), c(0, -1, 1)), rect=pRect, spec=spec) # (First ray points inside the rectangle, second lies along the # edge at the intersection of two faces, third lies in one of the # faces containing the origin.) # First and last rays are dropped because after flipping their only # intersection with the rectangle is the origin: tools::assertWarning( frays <- flipRays(rays), verbose=warn.tools_verbose) stopifnot(length(frays) == 1, rayPosn1(frays) == c(1, 1, 0)) stopifnot(identical(rayOrigin(rays), rayOrigin(frays))) stopifnot(all(pt_on_line2(rayPosn1(frays), lineidx=2, x0=rayOrigin(rays), x1=rayPosn1(rays), inptol=inpTol(rays)))) cat("OK\n") #----- lsetPkgOpt(oldopt) cat('bdryFromRays.r', sep='\n') # 17 Feb 2025 / Richard Raubertas # 24 Feb 2025 : 'fnObj' no longer has 'spec' arg. # 12 Feb 2025 : Move 4D test to 'test_extras' directory. oldopt <- lsetPkgOpt(reset=TRUE) #----- cat(" bdryFromRays ... ") #== Testing #----- 1D example: respf <- function(x, inclGrad) { resp <- -1*(x^2 + 2*x + 1) if (inclGrad) structure(resp, gradient= -1*(2*x + 2)) else resp } respF <- fnObj(c("X1"), respfn=respf, hasgrad=TRUE) pRays <- inpRays(matrix(c(5), ncol=1), origin=-2, spec=respF) rslt <- bdryFromRays(respF, rays=pRays, lsetthresh=-4) stopifnot(inherits(rslt, "lsetSegs"), length(rslt) == 1, isTRUE(all.equal(as.vector(ptCoord(rslt)), c(-3, 1), tolerance=1e-6))) summ <- summary(rslt) stopifnot(summ$nlines == 1, colSums(summ$lineInfo)[1:4] == c(1, 0, 2, 0)) # Extract lines/rays from an 'lsetSegs' object. stopifnot(identical(inpLines(rslt), pRays)) # No contour crossing between default positions -1, 0, and 1: pRays <- inpRays(matrix(c(0), ncol=1), origin=-1, spec=respF) #-- Only a censored segment is found if extra test positions are turned off: rslt <- bdryFromRays(respF, rays=pRays, lsetthresh=-4, control=list(initPosns2=numeric(0))) grad <- attr(respInfo(rslt, inclGrad=TRUE), "gradient") stopifnot(as.vector(ptCoord(rslt)) == c(-2, 0), is.matrix(grad), dim(grad) == c(2, 1), identical(dimnames(grad), list(NULL, "X1")), as.vector(grad) == c(2, -2)) summ <- summary(rslt) stopifnot(summ$nlines == 1, colSums(summ$lineInfo)[1:4] == c(1, 0, 0, 2)) #-- Both contour crossings found with default extra test positions: rslt <- bdryFromRays(respF, rays=pRays, lsetthresh=-4) stopifnot(inherits(rslt, "lsetSegs"), length(rslt) == 1, isTRUE(all.equal(as.vector(ptCoord(rslt)), c(-3, 1), tolerance=1e-6))) summ <- summary(rslt) stopifnot(summ$nlines == 1, colSums(summ$lineInfo)[1:4] == c(1, 0, 2, 0)) grad <- attr(respInfo(rslt, inclGrad=TRUE), "gradient") stopifnot(is.matrix(grad), dim(grad) == c(2, 1), identical(dimnames(grad), list(NULL, "X1")), as.vector(grad) == c(4, -4)) # Edge case: no rays. rslt <- bdryFromRays(respF, rays=pRays[0], lsetthresh=-4) stopifnot(inherits(rslt, "lsetSegs"), length(rslt) == 0) summ <- summary(rslt) stopifnot(summ$nlines == 0, colSums(summ$lineInfo)[1:4] == c(0, 0, 0, 0)) # Contour crossing exactly at an 'initPosn': pRays <- inpRays(matrix(c(-3), ncol=1), origin=-2, spec=respF) rslt <- bdryFromRays(respF, rays=pRays, lsetthresh=-4) stopifnot(inherits(rslt, "lsetSegs"), length(rslt) == 1, isTRUE(all.equal(as.vector(ptCoord(rslt)), c(1, -3), tolerance=1e-6))) #----- 2D Multimodal example (two elliptical contours). Ex <- dbl_ellipse_2dEx respF <- Ex$fobj fspec <- fnSpec(respF) # Small search region. refregion <- inpRect(rbind(c(-8, -8), c(8, 8)), spec=fspec) set.seed(1) pRays <- randomRays(70, spec=fspec, origin=c(2, 1), rect=refregion) # Default 'initPosns2' is enough to find boundary points on both ellipses ... rslt <- bdryFromRays(respF, rays=pRays, lsetthresh=-40) summ <- summary(rslt) # 0 censored, 21 of 70 rays have more than one segment stopifnot(summ$nlines == 70, colSums(summ$lineInfo)[1:4] == c(91, 0, 182, 0)) # ... but 'initPosns' alone is not. rslt <- bdryFromRays(respF, rays=pRays, lsetthresh=-40, control=list(initPosns2=numeric(0))) summ <- summary(rslt) # 0 censored, 21 of 70 rays have more than one segment stopifnot(summ$nlines == 70, colSums(summ$lineInfo)[1:4] == c(70, 0, 100, 40)) # Test 'segBoundingRect': rng <- apply(ptCoord(rslt), 2, range) bdryrect <- segBoundingRect(rslt) # just includes points, so matches range: stopifnot(isTRUE(all.equal(unname(rng), unname(rectCorners(bdryrect))))) # Expand censored points: bdryrect2 <- segBoundingRect(rslt, expand=c(1.0, 2.0)) stopifnot(round(prod(rectSize(bdryrect))) == 229, round(prod(rectSize(bdryrect2))) == 817) #-- Suppose we make poor choices of 'rect'. # (a) Far from level set. refregion2a <- inpRect(rbind(c(14, -2), c(15, -1)), spec=fspec) set.seed(1) pRays2a <- randomRays(70, spec=fspec, rect=refregion2a) rslt2a <- bdryFromRays(respF, rays=pRays2a, lsetthresh=-40) # empty summ <- summary(rslt2a) stopifnot(summ$nlines == 70, colSums(summ$lineInfo)[1:4] == c(0, 0, 0, 0)) # (b) Small and contained entirely in the level set. refregion2b <- inpRect(rbind(c(-1, -1), c(0, 0)), spec=fspec) set.seed(1) pRays2b <- randomRays(70, spec=fspec, rect=refregion2b) rslt2b <- bdryFromRays(respF, rays=pRays2b, lsetthresh=-40) summ <- summary(rslt2b) # all censored stopifnot(summ$nlines == 70, colSums(summ$lineInfo)[1:4] == c(70, 0, 0, 140)) #plot(rslt2b, rect=refregion2b) #-- Test feasible point checking. Ex <- dbl_ellipse_X1bnd_2dEx # Feasible region: X1 >= 0. respF <- Ex$fobj rslt <- bdryFromRays(respF, rays=pRays, lsetthresh=-40) summ <- summary(rslt) # 0 of 182 censored, 50 on feas bdry; 21 rays w/ multiple segs stopifnot(summ$nlines == 70, colSums(summ$lineInfo)[1:4] == c(91, 50, 132, 0)) #-- Use a 'rayProfiles' object as 'lines' argument. profs <- respProfiles(respF, lines=pRays, lsetthresh=-40, posns=c(-2, -1.5, -1, -0.5, 0, 0.5, 1, 1.5, 2)) #-- Just extract existing bdry intersections from 'profs': rslt2 <- lsetSegs(profs) # '.respProfiles' method summ <- summary(rslt2) stopifnot(summ$nlines == 70, colSums(summ$lineInfo)[1:4] == c(102, 50, 148, 6)) # One segment crosses the gap between ellipses and so is invalid. Can be # fixed with 'lsetSegsCheck': chk <- lsetSegsCheck(rslt2, fnobj=respF, posns=(1:4)/5) stopifnot(length(chk$validIn) == length(rslt2), sum(chk$validIn) == length(rslt2) - 1, length(chk$validOut) == length(rslt2) + 1, all(chk$validOut)) rslt2b <- chk$segs summ <- summary(rslt2b) stopifnot(summ$nlines == 70, colSums(summ$lineInfo)[1:4] == c(103, 50, 150, 6)) #-- Try changing 'lsetthresh' value from what was used in profiling. rslt3b <- bdryFromRays(respF, rays=profs, lsetthresh=-30) stopifnot(lsetThresh(rslt3b) == -30) # Should be the same as if we started with a profile using the new threshold: profs_alt <- respProfiles(respF, lines=pRays, lsetthresh=-30, posns=c(-2, -1.5, -1, -0.5, 0, 0.5, 1, 1.5, 2)) rslt3b_alt <- bdryFromRays(respF, rays=profs_alt) stopifnot(isTRUE(all.equal(rslt3b, rslt3b_alt, tolerance=1e-6)), isTRUE(all.equal(summary(rslt3b), summary(rslt3b_alt), tolerance=1e-6))) #-- Ray profiles that never intersect the level set (e.g., any level set # defined by a positive contour value): profs <- respProfiles(respF, lines=pRays, posns=c(-1, 0, 1), lsetthresh=10) rslt4 <- bdryFromRays(respF, rays=profs) # should be empty stopifnot(inherits(rslt4, "lsetSegs"), length(rslt4) == 0) summ <- summary(rslt4) stopifnot(summ$nlines == 70, colSums(summ$lineInfo)[1:4] == c(0, 0, 0, 0)) # Test 'segBoundingRect' with no level set points: tools::assertWarning({ bdryrect <- segBoundingRect(rslt4) } ) stopifnot(is.null(bdryrect)) cat("OK\n") #----- lsetPkgOpt(oldopt) cat('bdrySearch.r', sep='\n') # 24 Feb 2025 / Richard Raubertas # 24 Jul 2025 : Update names of 'control' parameters. # 27 Aug 2025 : Return component 'origins' renamed to 'srchOrigins'. # 12 Feb 2026 : Remove 4D tests with t-distribution example (now in # 'tests_extra'). # 01 Mar 2026 : No longer reference 'lvlSets_env' for 2D example. oldopt <- lsetPkgOpt(reset=TRUE) lsetPkgOpt(warn.showBadRespPts=0) #----- cat(" bdrySearch ... \n") #== Testing #----- 1D example. cat(" 1D examples ... ") respf <- function(x, ...) { x*cos(x) } #xx <- seq(-3, 10, by=0.25) #plot(xx, respf(xx), type="l") #abline(h=0.2, lty=2) fobj <- fnObj(c("X1"), respfn=respf) rect1 <- inpRect(cbind(c(-3, 10)), spec=fobj) srch1 <- bdrySearch(fobj, lsetthresh=0.2, rect=rect1) bdry1 <- srch1$segs summ <- summary(bdry1) stopifnot(inherits(bdry1, "lsetSegs"), length(bdry1) == 2, # 2 segments nrow(srch1$srchOrigins) == 1, # matrix not 'respPts' summ$nlines == 1, colSums(summ$lineInfo)[1:4] == c(2, 0, 4, 0)) # Upper segment is correct. Lower seg should really be # many segs; multiple oscillations below X1=1 are missed: if (FALSE) { xx <- seq(-30, 10, by=0.25) plot(xx, respf(xx), type="l") abline(h=0.2, lty=2) abline(v=as.vector(ptCoord(bdry1)), col="red") } #-- Let search go wild chasing oscillations: allow multiple search steps # even in 1D. tools::assertWarning( srch1b <- bdrySearch(fobj, lsetthresh=0.2, rect=rect1, control=list(maxSteps=5)) ) # 12 seg pairs reported inconsistent # With 'initPosns2' values > 1 and 'rect' updating the search region expands # dramatically over steps, and futilely, since the number of segments is # infinite. bdry1b <- srch1b$segs summ <- summary(bdry1b) # 7 segs on one line, no censoring stopifnot(inherits(bdry1b, "lsetSegs"), length(bdry1b) == 7, # 7 segments nrow(srch1b$srchOrigins) == 5, # matrix not 'respPts' summ$nlines == 1, colSums(summ$lineInfo)[1:4] == c(7, 0, 14, 0)) #??????? Add option in lsetSegsCheck to only check validity or only consistency temp <- lsetSegsCheck(bdry1b, fnobj=fobj, posns=(1:4)/5) # default 1 fix cycle stopifnot(length(temp$validIn) == length(bdry1b), sum(temp$validIn) == 1, length(temp$validOut) == 16, sum(temp$validOut) == 5, is.matrix(temp$inconOut), nrow(temp$inconOut) == 20) # Goes from 6/7 invalid to 11/16, with 20 inconsistent pairs temp2 <- lsetSegsCheck(temp$segs, fnobj=fobj, posns=(1:4)/5, fixCycles=5) stopifnot(length(temp2$validIn) == length(temp$segs), sum(temp2$validIn) == 5, length(temp2$validOut) == 20, sum(temp2$validOut) == 20, is.matrix(temp2$inconOut), nrow(temp2$inconOut) == 3) # Goes from 11/16 invalid to 0/20, with 3 inconsistent pairs # However, ~60 segs are expected given the range of seg box #-- Limit rect expansion by restricting 'initPosns2' to values <= 1. tools::assertWarning( srch1d <- bdrySearch(fobj, lsetthresh=0.2, rect=rect1, control=list(maxSteps=5, initPosns2=c(0.25, 0.5))) ) # 3 seg pairs reported inconsistent bdry1d <- srch1d$segs summ <- summary(bdry1d) stopifnot(inherits(bdry1d, "lsetSegs"), length(bdry1d) == 3, # 3 segments nrow(srch1d$srchOrigins) == 5, # matrix not 'respPts' summ$nlines == 1, colSums(summ$lineInfo)[1:4] == c(3, 0, 6, 0)) # Three (inconsistent) segs in a box only slightly larger than 'rect1' (due # to multiple search origins). Correct answer is three segs, but two # intermediate crossing pts are missed. temp <- lsetSegsCheck(bdry1d, posn=(1:4)/5, fnobj=fobj, fixCycles=3) stopifnot(length(temp$validIn) == length(bdry1d), sum(temp$validIn) == 1, length(temp$validOut) == 3, sum(temp$validOut) == 3, is.matrix(temp$inconOut), nrow(temp$inconOut) == 0) #table(temp$validIn) # 2 of 3 orignal segs are invalid #table(temp$validOut) # now 3 valid segs #print(temp$inconOut) # none are inconsistent #print(segInfo(temp$segs)) if (FALSE) { xx <- seq(-5, 10, by=0.25) plot(xx, respf(xx), type="l") abline(h=0.2, lty=2) abline(v=as.vector(ptCoord(temp$segs)), col="red") } # All segs correct. cat("OK\n") #----- 2D multimodal example (level set has two separate elliptical parts), cat(" 2D examples ... ") Ex <- dbl_ellipse_X1bnd_2dEx # Feasible region: X1 >= 0. respF <- Ex$fobj fspec <- fnSpec(respF) # Initial search region (default reference point is center, (0, 0)): refregion1 <- inpRect(rbind(c(-8, -8), c(8, 8)), spec=fspec) # Boundary search with defaults: srch1 <- bdrySearch(respF, lsetthresh=-40, rect=refregion1) bdry1 <- srch1$segs summ <- summary(bdry1) stopifnot(inherits(bdry1, "lsetSegs"), length(bdry1) == 9, nrow(srch1$srchOrigins) == 4, # matrix not 'respPts' summ$nlines == 8, colSums(summ$lineInfo)[1:4] == c(9, 3, 15, 0)) # Use the 'control' argument to increase the number of segments. srch1b <- bdrySearch(respF, lsetthresh=-40, rect=refregion1, control=list(axisDegree=2, # increases rays per origin maxSteps=8)) bdry1b <- srch1b$segs srch1b$srchOrigins # 8 origins used, 4 rays per origin summ <- summary(bdry1b) # 34 segments on 32 lines; 0 censored endpts; 2 lines w/ > 1 seg stopifnot(inherits(bdry1b, "lsetSegs"), length(bdry1b) == 34, nrow(srch1b$srchOrigins) == 8, # matrix not 'respPts' summ$nlines == 32, colSums(summ$lineInfo)[1:4] == c(34, 15, 53, 0)) #plot(bdry1b) # many segs cross the gap between ellipses (i.e., invalid) # Check segment validity, and attempt to fix any invalid ones. chk <- lsetSegsCheck(bdry1b, fnobj=respF, posns=(1:4)/5) stopifnot(table(chk$validIn) == c(11, 23), # 11 invalid segments table(chk$validOut) == c(45), # Fixed by splitting each into two is.matrix(chk$inconOut), nrow(chk$inconOut) == 0) bdry1c <- chk$segs summ <- summary(bdry1c) stopifnot(inherits(bdry1c, "lsetSegs"), length(bdry1c) == 45, summ$nlines == 32, colSums(summ$lineInfo)[1:4] == c(45, 15, 75, 0)) # Fewer invalid segments if we use a finer and wider grid of starting # positions for search. srch3 <- bdrySearch(respF, lsetthresh=-40, rect=refregion1, currentBdry=srch1, control=list(axisDegree=2, initPosns=(0:8)/4)) stopifnot(nrow(srch3$srchOrigins) == 8) bdry3 <- srch3$segs chk <- lsetSegsCheck(bdry3, fnobj=respF, posns=(1:4)/5) stopifnot(table(chk$validIn) == c(4, 26), # Four invalid segments table(chk$validOut) == c(34)) # Fixed by splitting each into two bdry3b <- chk$segs summ <- summary(bdry3b) stopifnot(inherits(bdry3b, "lsetSegs"), length(bdry3b) == 34, summ$nlines == 24, colSums(summ$lineInfo)[1:4] == c(34, 9, 59, 0)) # Try skipping the random rotation of ray vectors. srch4 <- bdrySearch(respF, lsetthresh=-40, rect=refregion1, control=list(axisDegree=2, maxSteps=8, rotateRays=FALSE)) bdry4 <- srch4$segs summ <- summary(bdry4) stopifnot(inherits(bdry4, "lsetSegs"), length(bdry4) == 35, nrow(srch4$srchOrigins) == 8, # matrix not 'respPts' summ$nlines == 32, colSums(summ$lineInfo)[1:4] == c(35, 12, 58, 0)) #plot(bdry4, rect=refregion1) # Not very good. chk <- lsetSegsCheck(bdry4, fnobj=respF, posns=(1:4)/5) stopifnot(table(chk$validIn) == c(7, 28), table(chk$validOut) == 42, nrow(chk$inconOut) == 0) #plot(chk$segs, rect=refregion1) # not bad # Conclusion: # a) Important to have right search region ('rect'), updating after initial # guess(es) if necessary. # b) Need grid of 'initPosns' fine enough to match scale of interesting # features of level set shape (concavity or multiple parts). # c) Need to increase number of boundary pairs to see detailed shape of # level set. Can be done by increasing the number of origins # ('maxSteps') initially or over several calls, and/or by # increasing 'axisDegree'. # d) 'rotateRays' seems to be useful. cat("OK\n") #----- lsetPkgOpt(oldopt) cat('fnObj.r', sep='\n') # fnObj # hasGrad # warnBadResp* # fnObjCheck* (???? incomplete and disabled) # *= not public # 12 Feb 2025 / Richard Raubertas # 24 Feb 2025 : 'fnObj' no longer has a 'spec' arg. Test new 'update.fnObj'. oldopt <- lsetPkgOpt(reset=TRUE) lsetPkgOpt(warn.showBadRespPts=0) # avoid printing matrices, data frames #----- cat(" fnObj ... ") #== Testing # Function of 3 input variables, with gradient available, and an # extra argument 'powr'. f <- function(x, inclGrad=FALSE, powr) { y <- x[, "X1"]^powr[1] + x[, "X2"]^powr[2] + x[, "X3"]^powr[3] if (inclGrad) { grad <- array(NA_real_, dim=dim(x), dimnames=dimnames(x)) for (j in 1:3) { grad[, j] <- { if (powr[j] == 0) 0 else powr[j] * (x[, j]^(powr[j] - 1)) } } attr(y, "gradient") <- grad } y } # Associated feasibility function, passed the same extra argument. feasf <- function(x, powr) { feas <- rep(TRUE, nrow(x)) # Want to avoid applying non-integer powers to negative numbers. fracpow <- (powr != round(powr)) for (j in which(fracpow)) { feas <- feas & (x[, j] >= 0) # NOTE: '>=', not '>'. Feasible set # should be _closed_. } feas } fnobj <- fnObj(c("X1", "X2", "X3"), respfn=f, feasfn=feasf, hasgrad=TRUE, powr=c(-1, 0, 2.5)) stopifnot(hasGrad(fnobj)) fspec <- fnSpec(fnobj) # Tests: pts with NA coord, or generating Inf resp: pts <- rbind(c(0, 0, 0), c(2, 0, -1), c(1, 6, 1), c(NA, 1, 1)) colnames(pts) <- inpNames(fnobj) tools::assertError(respInfo(pts, fnobj=fnobj), # NA pt coord not allowed verbose=warn.tools_verbose) pts[4, ] <- c(-1, 1, 1) tools::assertWarning(respInfo(pts, fnobj=fnobj), # Inf response verbose=warn.tools_verbose) pts[1, ] <- c(3, -2, 0) resp1 <- respInfo(pts, fnobj=fnobj, inclGrad=TRUE) stopifnot(is.data.frame(resp1), nrow(resp1) == nrow(pts), resp1$feas == c(TRUE, FALSE, TRUE, TRUE), is.na(resp1$resp) == !resp1$feas) grad <- attr(resp1, "gradient") stopifnot(dim(grad) == dim(pts), apply(grad, 1, function(y) {any(is.na(y))}) == !resp1$feas) resp2 <- respInfo(pts, fnobj=fnobj, inclGrad=FALSE) grad <- attr(resp2, "gradient") stopifnot(is.null(grad)) # Test '.fnSpec' method: fspec2 <- fnSpec(inpnames=c("X1", "X2", "X3")) fnobj2 <- fnObj(fspec2, respfn=f, feasfn=feasf, hasgrad=TRUE, powr=c(-1, 0, 2.5)) resp2 <- respInfo(pts, fnobj=fnobj2, inclGrad=TRUE) stopifnot(identical(fnSpec(fnobj), fnSpec(fnobj2)), identical(resp1, resp2)) #-- Feasibility function based on 'feasbnds'. fnobj <- fnObj(c("X1", "X2", "X3"), respfn=f, feasbnds=rbind(rep(0, 3), rep(NA, 3)), hasgrad=FALSE, powr=c(-1, 0, 2.5), warnInfresp=FALSE) stopifnot(!hasGrad(fnobj)) pts <- rbind(c(0, 0, 0), c(2, 0, -1), c(1, 6, 1), c(-2, 1, 1)) colnames(pts) <- inpNames(fnobj) resp3 <- respInfo(pts, fnobj=fnobj) # no warning about Inf response stopifnot(is.data.frame(resp3), nrow(resp3) == nrow(pts), resp3$feas == c(TRUE, FALSE, TRUE, FALSE), is.na(resp3$resp) == !resp3$feas) #-- Combine explicit feasibility function with 'feasbnds'. fnobj <- fnObj(c("X1", "X2", "X3"), respfn=f, feasfn=feasf, feasbnds=list(X2=c(2, NA)), hasgrad=TRUE, powr=c(-1, 0, 2.5), warnInfresp=0) pts <- rbind(c(0, 0, 0), c(2, 3, -1), c(1, 6, 1), c(0, 4, 1), c(1, 1, 1)) colnames(pts) <- inpNames(fnobj) # Negative value makes 2nd pt infeasible, lower bound makes 1st and 5th # infeasible. resp4 <- respInfo(pts, fnobj=fnobj) # no warning about Inf response stopifnot(is.data.frame(resp4), nrow(resp4) == nrow(pts), resp4$feas == c(FALSE, FALSE, TRUE, TRUE, FALSE), is.na(resp4$resp) == !resp4$feas) grad <- attr(resp4, "gradient") stopifnot(is.null(grad)) # b/c default 'inclGrad' is FALSE even w/ 'hasGrad=TRUE' #----- Testing 'warnNAresp', and 'update' method. respf <- function(x, inclGrad=FALSE) { resp <- rowSums(x) ifelse(x[, 1] + x[, 2] >= 2, resp, NA) # implicit feasible region } feasf <- function(x) { # explicit feasible region (x[, 1] - x[, 2] >= 1) } pts <- rbind(c(3, 1, 4), # feasible by all tests c(5, -1, 4), # not feas wrt box constraint c(1, 1, 4), # not feas wrt explicit constraint c(1, 0, 4)) # not feas wrt implicit constraint # No feasfn, silently accept implicitly infeas pts: fnobj <- fnObj(c("X1", "X2", "X3"), respfn=respf, warnNAresp=FALSE, feasbnds=list(X2=c(0, NA))) # box constraint feasible region resp <- respInfo(pts, fnobj=fnobj) stopifnot(identical(resp[, "resp"], c(8, NA, 6, NA)), resp[, "feas"] == !is.na(resp$resp), identical(resp[, "lset"], c(NA, FALSE, NA, FALSE))) # No feasfn, warn about implicitly infeas pts: fnobj2 <- fnObj(c("X1", "X2", "X3"), respfn=respf, warnNAresp=TRUE, feasbnds=list(X2=c(0, NA))) tools::assertWarning(resp2 <- respInfo(pts, fnobj=fnobj2), verbose=warn.tools_verbose) stopifnot(identical(resp2[, "resp"], c(8, NA, 6, NA)), resp2[, "feas"] == !is.na(resp2$resp), identical(resp2[, "lset"], c(NA, FALSE, NA, FALSE))) # No feasfn, treat implicitly infeas pts as error: fnobj <- update(fnobj, warnNAresp=2) tools::assertError(resp <- respInfo(pts, fnobj=fnobj), verbose=warn.tools_verbose) # Same, relying on default 'warnNAresp': fnobj <- fnObj(c("X1", "X2", "X3"), respfn=respf) tools::assertError(resp <- respInfo(pts, fnobj=fnobj), verbose=warn.tools_verbose) # Specify feasfn, silently accept implicitly infeas pts: fnobj <- fnObj(c("X1", "X2", "X3"), respfn=respf, feasfn=feasf, feasbnds=list(X2=c(0, NA)), warnNAresp=FALSE) resp <- respInfo(pts, fnobj=fnobj) stopifnot(identical(resp[, "resp"], c(8, NA, NA, NA)), resp[, "feas"] == !is.na(resp$resp), identical(resp[, "lset"], c(NA, FALSE, FALSE, FALSE))) # Specify feasfn, warn about implicitly infeas pts: fnobj <- update(fnobj, warnNAresp=TRUE) tools::assertWarning(resp <- respInfo(pts, fnobj=fnobj), verbose=warn.tools_verbose) stopifnot(identical(resp[, "resp"], c(8, NA, NA, NA)), resp[, "feas"] == !is.na(resp$resp), identical(resp[, "lset"], c(NA, FALSE, FALSE, FALSE)), fnobj$warnNAresp == 1) # Specify feasfn, treat implicitly infeas pts as error: fnobj <- update(fnobj, warnNAresp=2) tools::assertError(resp <- respInfo(pts, fnobj=fnobj), verbose=warn.tools_verbose) # No explicit feasible region at all, just implicit one. respf <- function(x) { resp <- rowSums(x) ifelse(x[, 1] + x[, 2] >= 2, resp, NA) # implicit feasible region } # Default issues error: fnobj <- fnObj(c("X1", "X2", "X3"), respfn=respf) tools::assertError(resp <- respInfo(pts, fnobj=fnobj), verbose=warn.tools_verbose) # User must intentionally disable error to rely on implicit feasible # region: fnobj <- update(fnobj, warnNAresp=0) resp <- respInfo(pts, fnobj=fnobj) stopifnot(identical(resp[, "resp"], c(8, 8, 6, NA)), resp[, "feas"] == !is.na(resp$resp), "feasFn" %in% names(fnobj) && is.null(fnobj$feasFn)) stopifnot(!hasGrad(fnobj)) # default cat("OK\n") #----- cat(" warnBadResp ... ") #== Testing pts <- rbind(c(1, 2), c(2, 3), c(3, 4), c(4, 5), c(5, 6)) resp <- c(4, Inf, -Inf, 2, NA) tools::assertWarning(warnBadResp(resp, pts, type="Inf"), verbose=warn.tools_verbose) # 2 pts printed tools::assertWarning(warnBadResp(resp, pts, type="NA"), verbose=warn.tools_verbose) # 1 pt printed #lsetPkgOpt(warn.showBadRespPts=1) tools::assertWarning(warnBadResp(resp, pts, type="Inf"), verbose=warn.tools_verbose) # 1 pt printed #lsetPkgOpt(warn.showBadRespPts=0) tools::assertWarning(warnBadResp(resp, pts, type="Inf"), verbose=warn.tools_verbose) # no pts printed #lsetPkgOpt(warn.showBadRespPts=1) tools::assertError(warnBadResp(resp, pts, type="Inf", level=2), verbose=warn.tools_verbose) # 1 pt printed cat("OK\n") #----- cat(" update.fnobj ... ") fnobj <- fnObj(c("X1", "X2"), respfn=function(x) { colSums(x) }, inptol=c(1e-6, 1e-8)) inptol1 <- inpTol(fnobj) stopifnot(is.matrix(inptol1), dim(inptol1) == c(2, 2), unname(inptol1) == cbind(c(1e-6, 1e-8), c(1e-6, 1e-8))) fnobj2 <- update(fnobj, inptol=c(1e-4, 2e-8)) inptol2 <- inpTol(fnobj2) stopifnot(is.matrix(inptol2), dim(inptol2) == c(2, 2), unname(inptol2) == cbind(c(1e-4, 2e-8), c(1e-4, 2e-8))) cat("OK\n") #----- lsetPkgOpt(oldopt) # restore original options cat('fnSpec.r', sep='\n') # Define 'fnSpec' class and related accessor functions. # fnSpec # print.fnSpec* # fnspec_accessor* # inpDim # inpNames # feasBnds # inpTol # respTol # *= not public # 11 Feb 2025 / Richard Raubertas oldopt <- lsetPkgOpt(reset=TRUE) #----- cat(" fnSpec ... ") #== Testing fspec <- fnSpec(inpnames=c("X1", "X2", "X3"), feasbnds=list(X2=c(0, 10), X3=c(NA, 5))) stopifnot(identical(fspec, fnSpec(fspec))) stopifnot(dim(inpTol(fspec)) == c(2, 3)) stopifnot(!is.matrix(respTol(fspec)), length(respTol(fspec)) == 2) stopifnot(inpDim(fspec) == 3) stopifnot(inpNames(fspec) == c("X1", "X2", "X3")) stopifnot(feasBnds(fspec) == rbind(c(-Inf, 0, -Inf), c(Inf, 10, 5))) # Updating: fspec <- fnSpec(inpnames=c("X1", "X2", "X3")) fspec2 <- fnSpec(fspec, feasbnds=list(X2=c(0, 10), X1=c(NA, 5))) stopifnot(identical(inpTol(fspec2), inpTol(fspec))) fspec2b <- fnSpec(fspec, feasbnds=feasBnds(fspec2)) stopifnot(identical(fspec2, fspec2b)) fspec3 <- fnSpec(fnSpec(fspec2, resptol=c(1e-6, 1e-7)), resptol=respTol(fspec2)) stopifnot(identical(fspec3, fspec2)) cat("OK\n") #----- lsetPkgOpt(oldopt) cat('inpLines.r', sep='\n') # inpLines (S3 generic) # inpLines.default (handles matrices of pt coords, and 'inpPts' objects) # inpLines.inpPairs* (also handles objects inheriting from 'inpPairs') # inpLines.inpLines* (also handles objects inheriting from 'inpLines') # inpLines.fnSpec (allows creation of empty 'inpLines' object) # inpLines.respProfiles # inpLines.lsetSegs (needed, otherwise 'inpPairs' method takes precedence # over 'respProfiles' method) # *= not public oldopt <- lsetPkgOpt(reset=TRUE) #----- cat(" inpLines ... ") #== Testing fspec <- fnSpec(inpnames=c("X1", "X2", "X3")) pts <- rbind(c(0, 1, 3), c(2, 2, 4), c(5, -1, 2), c(1, 1, 1)) # The default method of defining positions 0 and 1 for each line is # to use consecutive pairs of rows of 'pts': ilines <- inpLines(pts, spec=fspec) stopifnot(length(ilines) == 2) # Alternative way of specifying positions 0 and 1 is by two matrices: ilines2 <- inpLines(pts[c(2, 4), ], x0=pts[c(1, 3), ], spec=fspec) stopifnot(identical(ptCoord(ilines), ptCoord(ilines2, which="paired"))) names(ilines) <- c("foo", "bar") stopifnot(colnames(pairIdxs(ilines)) == c("foo", "bar")) # Downshifting to 'inpPairs' should remove names. stopifnot(is.null(colnames(pairIdxs(inpPairs(ilines))))) # Empty 'inpLines' object. ilines0 <- inpLines(fspec) stopifnot(identical(class(ilines0), class(ilines)), length(ilines0) == 0, identical(ilines0, ilines[0])) ilines00 <- inpLines(pts[0, , drop=FALSE], spec=fspec) stopifnot(identical(ilines0, ilines00)) # Non-default pairing info. pidxs <- rbind(c(1, 1, 1), c(2, 3, 4)) colnames(pidxs) <- LETTERS[1:3] ilines2 <- inpLines(pts, pairidxs=pidxs, spec=fspec) # Accessor functions: stopifnot(length(ilines2) == 3) stopifnot(dim(pairIdxs(ilines2)) == c(2, 3), colnames(pairIdxs(ilines2)) == LETTERS[1:3]) pairInfo(ilines2) # data frame stopifnot(identical(unname(linePosn0(ilines2)), pts[c(1, 1, 1), ])) stopifnot(identical(unname(linePosn1(ilines2)), pts[2:4, ])) stopifnot(isTRUE(all.equal(lineVec(ilines2), linePosn1(ilines2) - linePosn0(ilines2)))) stopifnot(dim(ptCoord(ilines2)) == c(nrow(pts), inpDim(fspec))) # 4 rows (raw point coords, ignoring pairing) stopifnot(dim(ptCoord(ilines2, which="paired")) == c(2*length(ilines2), inpDim(fspec))) # 2*number of pairs rows names(ilines2) <- c("foo", "bar", "baz") stopifnot(identical(linePosn0(ilines2), ptCoord(ilines2, which="1st"))) stopifnot(identical(linePosn1(ilines2), ptCoord(ilines2, which="2nd"))) # Combine multiple 'inpLines' objects. stopifnot(identical(ilines2, c(ilines2))) # Internal representation will not be identical, so this may be FALSE ... #identical(ilines2, c(ilines2[1:2], ilines2[0], ilines2[3], NULL)) # ... but actual point pairs will be the same: ptprs <- ptCoord(ilines2, which="paired") cptprs <- ptCoord(c(ilines2[1:2], ilines2[0], ilines2[3], NULL), which="paired") stopifnot(identical(ptprs, cptprs)) infoprs <- pairInfo(ilines2) cinfoprs <- pairInfo(c(ilines2[1:2], ilines2[0], ilines2[3], NULL)) stopifnot(identical(infoprs$foo, cinfoprs$foo)) # When 'x' is already an 'inpLines' object. stopifnot(identical(ilines, inpLines(ilines)), identical(ilines0, inpLines(ilines0))) # 'strictClass' argument. rays <- inpRays(pts, origin=c(2, 2, 2), spec=fspec) lines_from_rays <- inpLines(rays) # default 'strictClass=TRUE' stopifnot(!inherits(lines_from_rays, "inpRays")) still_rays <- inpLines(rays, strictClass=FALSE) stopifnot(identical(rays, still_rays), identical(ilines, inpLines(ilines, strictClass=FALSE))) # Degenerate lines. pts <- rbind(c(0, 1, 3), c(2, 2, 4), c(5, -1, 2), c(0, 1, 3) + 1e-8) # One pair tied exactly, another to within default tolerance: pairs <- rbind(c(1, 1, 1, 1), c(1, 2, 3, 4)) # Default is to drop tied pairs and warn: tools::assertWarning(lines <- inpLines(pts, pairidxs=pairs, spec=fspec), verbose=warn.tools_verbose) # (Only the 1st warning message is printed, but both tied pairs are dropped:) stopifnot(inherits(lines, "inpLines"), length(lines) == 2) # Drop tied pairs silently: lines <- inpLines(pts, pairidxs=pairs, spec=fspec, warnDegenerate=0) stopifnot(inherits(lines, "inpLines"), length(lines) == 2) # Make tied pairs an error: tools::assertError(lines <- inpLines(pts, pairidxs=pairs, spec=fspec, warnDegenerate=2), verbose=warn.tools_verbose) cat("OK\n") #----- lsetPkgOpt(oldopt) cat('inpObjAccessors.r', sep='\n') # rectCorners # rectRefPt # rectSize # sliceMatch (not public) oldopt <- lsetPkgOpt(reset=TRUE) #----- cat(" sliceMatch ... ") #== Testing fspec <- fnSpec(inpnames=c("a", "b")) inptol <- inpTol(fspec) slices <- sliceMat(rbind(c(1, NA), c(NA, 4)), spec=fspec) y <- rbind(c(1, 3), c(2, 2), c(1, 4), c(3, 4)) stopifnot(identical(sliceMatch(slices, y, inptol=inptol), list(`a=1`=c(1L, 3L), `b=4`=c(3L, 4L)))) stopifnot(identical(sliceMatch(slices, inpPts(y, spec=fspec), inptol=inptol), list(`a=1`=c(1L, 3L), `b=4`=c(3L, 4L)))) # same stopifnot(identical(sliceMatch(slices, inpPairs(y, spec=fspec), inptol=inptol), list(`a=1`=integer(0), `b=4`=c(2L)))) stopifnot(identical(sliceMatch(slices, inpLines(y, spec=fspec), inptol=inptol), list(`a=1`=integer(0), `b=4`=c(2L)))) # same cat("OK\n") #----- lsetPkgOpt(oldopt) cat('inpObjClasses.r', sep='\n') # Classes of objects in an input space. ('inpLines' and 'inpRays' are part of # the 'ptsClass' set of classes and are in separate files.) # inpRect # inpScale # 27 Jul 2025 / Richard Raubertas oldopt <- lsetPkgOpt(reset=TRUE) #----- cat(" inpRect ... ") #== Testing fspec <- fnSpec(inpnames=c("X1", "X2", "X3")) # Directly create a hyper-rectangle: rect1 <- inpRect(rbind(c(0, -1.1, -4), c(2, -0.5, 1000)), spec=fspec) # Accessor functions. corners <- rectCorners(rect1) stopifnot(dim(corners) == c(2, inpDim(fspec)), inpNames(rect1) == inpNames(fspec), colnames(corners) == inpNames(fspec)) refpt <- rectRefPt(rect1) stopifnot(dim(refpt) == c(1, inpDim(fspec)), colnames(refpt) == inpNames(fspec)) size <- rectSize(rect1) stopifnot(length(size) == inpDim(fspec), names(size) == inpNames(fspec), isTRUE(all.equal(size, corners[2, ] - corners[1, ]))) tools::assertError( # hyper-rectangle must be finite inpRect(rbind(c(0, -1.1, -4), c(2, -0.5, Inf)), spec=fspec), verbose=warn.tools_verbose) # Hyper-rectangle must have non-zero d-dimensional volume. tools::assertError( inpRect(rbind(c(0, -1.1, 4), c(2, -0.5, 4 + 1e-8)), spec=fspec), verbose=warn.tools_verbose) # Non-default reference point: rect1 <- inpRect(rbind(c(0, 1, 1), c(0.5, 3, 6)), refpt=c(0, 1, 4), spec=fspec) # Update reference point: rect2 <- update(rect1, refpt=c(0.2, 3, 3)) stopifnot(isTRUE(all.equal(as.vector(rectRefPt(rect2)), c(0.2, 3, 3)))) tools::assertError( # 'refpt' not inside the rectangle rect2 <- update(rect1, refpt=c(0.2, 0, 3)), verbose=warn.tools_verbose) # If refpt is only slightly outside the rectangle, the rectangle will # be adjusted to accomodate it. fspec <- fnSpec(inpnames=c("X1", "X2", "X3"), inptol=c(0, 1e-3)) # defines "slightly" rect1 <- inpRect(rbind(c(0, 1, 1), c(0.5, 3, 6)), refpt=c(0, 1, 4), spec=fspec) rect1b <- update(rect1, refpt=c(0.2, 3 + 5e-4, 3)) # Strictly includes the new refpt: stopifnot(inRect(c(0.2, 3 + 5e-4, 3), rect=rect1b, bdryOnly=TRUE, inptol=inpTol(fspec)/1000), all(rectCorners(rect1b)[2, ] >= c(0.2, 3 + 5e-4, 3))) # Not true for original rect: stopifnot(!all(rectCorners(rect1)[2, ] >= c(0.2, 3 + 5e-4, 3))) cat("OK\n") #----- cat(" inpScale ... ") #== Testing dat <- data.matrix(iris[, 1:4]) fspec <- fnSpec(inpnames=colnames(dat)) covmat <- cov(dat) sds <- sqrt(diag(covmat)) cormat <- cor(dat) # Check equivalence of alternate ways of specifying V. scl1 <- inpScale(sds, spec=fspec) # w = 1/sds stopifnot(isTRUE(all.equal(scl1$w, 1/sds)), scl1$V_is_diag, isTRUE(all.equal(diag(scl1$W), unname(1/sds))), isTRUE(all.equal(diag(scl1$Winv), unname(sds)))) scl1a <- inpScale(diag(sds^2, nrow=4), spec=fspec) scl2 <- inpScale(covmat, spec=fnSpec(scl1)) # 'spec' copied from 'scl1' scl2a <- inpScale(covmat, spec=fspec) # consistency check stopifnot(isTRUE(all.equal(scl1$W, scl1a$W)), identical(scl2, scl2a)) # Identity checks. stopifnot(isTRUE(all.equal(scl2$Winv %*% scl2$W, diag(nrow=ncol(covmat)))), isTRUE(all.equal(t(scl2$W) %*% scl2$W, unname(solve(covmat)))), isTRUE(all.equal(scl2$Winv %*% t(scl2$Winv), unname(covmat)))) # Mahalanobis distance and Euclidean distance (from origin in this example). Vinv <- solve(covmat) mdist <- apply(dat, 1, function(y) { sqrt( t(y) %*% Vinv %*% y) }) scl2 <- inpScale(covmat, spec=fspec) tdat <- dat %*% t(scl2$W) # t() since we are transforming _row_ vectors edist <- sqrt(rowSums(tdat^2)) # Euclidean dist using transformed vectors stopifnot(isTRUE(all.equal(mdist, edist))) # Coordinate system rescaling/rotation. chk1 <- sweep(dat, 2, scl1$w, "*") # re-scale each dimension to unit SD stopifnot(isTRUE(all.equal(cov(chk1), cormat))) chk2 <- dat %*% t(scl2$W) # re-scale and orthogonalize (using t() since # we are transforming _row_ vectors) stopifnot(isTRUE(all.equal(cov(chk2), diag(4), check.attributes=FALSE))) # Reverse the transformation: z_to_y <- scl2$Winv chk3 <- chk2 %*% t(z_to_y) # t() since we are transforming _row_ vectors dimnames(chk3) <- dimnames(dat) stopifnot(isTRUE(all.equal(chk3, dat))) # Empty/null 'inpScale' ('w', a vector of 1's, 'V' the identity matrix): scl <- inpScale(spec=fspec) stopifnot(identical(unname(scl$w), rep(1, inpDim(fspec))), names(scl$w) == inpNames(fspec), identical(unname(scl$V), diag(nrow=inpDim(fspec))), identical(dimnames(scl$V), list(inpNames(fspec), inpNames(fspec))), scl$V_is_diag) #----- 1D example. fspec <- fnSpec(inpnames="X1") scl <- inpScale(matrix(5, nrow=1, ncol=1), spec=fspec) stopifnot(identical(fnSpec(scl), fspec)) cat("OK\n") #----- lsetPkgOpt(oldopt) cat('inpRays.r', sep='\n') # 15 Feb 2025 / Richard Raubertas # 25 Feb 2025 : Include tests of new 'update' method, including preservation # of ray directions. 'rayOrigin' now includes colnames. oldopt <- lsetPkgOpt(reset=TRUE) #----- cat(" inpRays ... ") #== Testing fspec <- fnSpec(inpnames=c("X1", "X2", "X3")) d <- inpDim(fspec) #-- Rays defined by position 1 pts pts <- rbind(c(0, 1, 3), c(2, 2, 4), c(5, -1, 2)) orig <- c(0.2, 1.1, 5.5) rays <- inpRays(pts, origin=orig, spec=fspec) # Or by specifying their direction vectors and origin: vecs <- sweep(pts, 2, orig, "-") # equivalent to 'rayVec(rays)' except dimnames stopifnot(isTRUE(all.equal(vecs, unname(rayVec(rays))))) rays2 <- inpRays(vecs=vecs, origin=orig, spec=fspec) stopifnot(identical(rays, rays2)) #-- Ray names: rownames(pts) <- paste0("ray", seq_len(nrow(pts))) rays1 <- inpRays(pts, spec=fspec) stopifnot(names(rays1) == rownames(pts)) vecs <- sweep(pts, 2, orig, "-") rays2 <- inpRays(vecs=pts, spec=fspec) stopifnot(identical(names(rays2), names(rays1))) names(rays1) <- LETTERS[1:3] stopifnot(names(rays1) == LETTERS[1:3]) # changed #-- Subscripting: stopifnot(identical(rays[], rays)) stopifnot(length(rays[2:3]) == 2) stopifnot(length(rays[0]) == 0) #-- Accessors and methods. stopifnot(dim(rayOrigin(rays)) == c(1, d), dim(rayPosn1(rays)) == c(length(rays), d), dim(rayVec(rays)) == c(length(rays), d)) rayrect <- rayRect(rays) if (!is.null(rayrect)) { stopifnot(inherits(rayrect, "inpRect"), identical(rectRefPt(rayrect), rayOrigin(rays)), all(inRect(rayPosn1(rays), rect=rayrect, bdryOnly=TRUE))) } # Can update the origin or lengths of ray vectors with the (internal) # 'scale' method; see the examples for 'scale.inpRays'. rays2 <- inpRays(rbind(c(1, 1, 1), c(1, 2, 3)), origin=c(0, 0, 0), spec=fspec) # 'rayVecs' and 'rayPosn1' are equivalent when the ray origin is 0: stopifnot(isTRUE(all.equal(rayVec(rays2), rayPosn1(rays2)))) #-- Empty 'inpRays': rays0 <- inpRays(spec=fspec) # no rays, 0 origin stopifnot(inherits(rays0, "inpRays"), length(rays0) == 0) rays0 <- inpRays(spec=fspec, origin=orig) # no rays stopifnot(length(rays0) == 0, rayOrigin(rays0) == orig) stopifnot(identical(rays0, rays[0])) #-- Combining objects: stopifnot(identical(rays, c(rays[0], NULL, rays[1:2], rays[3]))) rays2 <- inpRays(vecs=rayVec(rays), origin=rayOrigin(rays) + 2, spec=fspec) tools::assertError(rays12 <- c(rays, rays2), # different origins verbose=warn.tools_verbose) stopifnot(identical(rays, inpRays(rays))) #-- Ray position 1 points can be (re)defined by each ray's intersection # with the boundary of a hyper-rectangle: rect1 <- inpRect(rbind(c(0,0,0), c(1,1,1)), spec=fspec) # (default reference pt is center: c(0.5, 0.5, 0.5)) vecs1 <- rbind(c(-1,-1,0), c(0,0,1), c(0, -1, 1)) rays1 <- inpRays(vecs=vecs1, rect=rect1) # default 'origin' and 'spec' taken from 'rect' stopifnot(all(inRect(rayPosn1(rays1), rect=rect1, bdryOnly=TRUE)), identical(rayOrigin(rays1), rectRefPt(rect1)), identical(rayRect(rays1), rect1)) # Explicit 'origin' argument overrides use of reference point from # 'rect', and updates the latter: rays1 <- inpRays(vecs=vecs1, origin=c(1, 1, 0.5), # origin is on boundary of rect rect=rect1, spec=fspec) stopifnot(rayOrigin(rays1) == c(1, 1, 0.5)) stopifnot(inRect(rayPosn1(rays1), rect=rect1, bdryOnly=TRUE)) stopifnot(identical(rectRefPt(rayRect(rays1)), rayOrigin(rays1))) # An existing 'inpRays' object can be updated with a new origin and/or # hyper-rectangle (which must contain the existing or new ray origin). # Position 1 points are on the boundary of the new rectangle, the new # rectangle reference pt is set to the ray origin, which is unchanged. rect2 <- inpRect(rbind(c(-1, -2, -3), c(1, 2, 3)), spec=fspec) rays2 <- update(rays1, rect=rect2) stopifnot(identical(rectCorners(rayRect(rays2)), rectCorners(rect2)), identical(rayOrigin(rays2), rayOrigin(rays1)), inRect(rayPosn1(rays2), rect=rect2, bdryOnly=TRUE), # Ray directions are unchanged: isTRUE(all.equal(rowSums(rayVec(rays1)*rayVec(rays2)), sqrt(rowSums(rayVec(rays1)^2) * rowSums(rayVec(rays2)^2))))) # Note: internally, ref pt of 'rect' is set to ray origin: stopifnot(identical(rectRefPt(rayRect(rays2)), rayOrigin(rays2))) # Change only 'origin' (which should also reset ref pt of its rect). rays3 <- update(rays1, origin=c(0.2, 0.3, 0.4)) stopifnot(rectCorners(rayRect(rays3)) == rectCorners(rect1), as.vector(rayOrigin(rays3)) == c(0.2, 0.3, 0.4), as.vector(rectRefPt(rayRect(rays3))) == c(0.2, 0.3, 0.4)) # Change both 'origin' and 'rect'. rays3 <- update(rays1, origin=c(0, -1, 2), rect=rect2) stopifnot(identical(rectRefPt(rayRect(rays3)), rayOrigin(rays3)), inRect(rayPosn1(rays3), rect=rect2, bdryOnly=TRUE), as.vector(rayOrigin(rays3)) == c(0, -1, 2), # Ray directions are unchanged: isTRUE(all.equal(rowSums(rayVec(rays1)*rayVec(rays3)), sqrt(rowSums(rayVec(rays1)^2) * rowSums(rayVec(rays3)^2))))) # Rays with different rect's can be combined, as long as they have the # same origin. The combination uses the rect of the first. rays12 <- c(rays1, rays2) stopifnot(length(rays12) == length(rays1) + length(rays2), identical(rayRect(rays12), rayRect(rays1))) #-- Degenerate rays are dropped: vecs2 <- rbind(c(-1,-1,0), c(0,0,0), c(0, -1, 1)) # 2nd vector has 0 length tools::assertWarning( rays <- inpRays(vecs=vecs2, origin=c(0.5, 0.5, 0.5), spec=fspec), verbose=warn.tools_verbose) stopifnot(length(rays) == 2) # not 3 # One can also end up with a degenerate ray if a rect is specified and # its reference point is on the boundary. In that case a ray's only # intersection with the rect boundary may be the origin: rect <- inpRect(rbind(c(0,0,0), c(1,1,1)), spec=fspec) orig <- c(1,1,0.5) # on the boundary of 'rect' vecs <- rbind(c(-1,-1,0), c(0,1,1), c(0, -1, 1)) # The second direction vector points away from the interior of 'rect', # relative to 'orig'. tools::assertWarning( rays <- inpRays(vecs=vecs, origin=orig, rect=rect), # (Second ray's only intersection with 'rect' is the origin.) verbose=warn.tools_verbose) stopifnot(length(rays) == 2) # Can force degenerate rays to be treated as an error: tools::assertError( rays <- inpRays(vecs=vecs, origin=orig, rect=rect, warnDegenerate=2), verbose=warn.tools_verbose) #-- Other warnings and errors. orig2 <- c(2, 2, 6) tools::assertError( inpRays(pts, origin=orig2, rect=rect1), # ray origin does not belong to rectangle verbose=warn.tools_verbose) cat("OK\n") #----- lsetPkgOpt(oldopt) cat('inpRaysMethods.r', sep='\n') # Methods specifically for 'inpRays' objects (where corresponding 'inpPts' # and 'inpPairs' methods in 'ptsClassMethods.r' are not adequate). # # scale.inpRays* (16 Jan 2026: works ok, has tests, but not currently # used, so disabled) # plot.inpRays # update.inpRays # *= not public # 13 Apr 2026 / Richard Raubertas oldopt <- lsetPkgOpt(reset=TRUE) #----- if (exists("scale.inpRays")) { cat(" scale.inpRays ... ") #== Testing fspec <- fnSpec(inpnames=LETTERS[1:3]) pRect <- inpRect(rbind(c(2, 0, -2), c(4, 4, 4)), spec=fspec) # 2 x 4 x 6 box centered at (3,2,1) pRays <- inpRays(rbind(c(2, 0, -2), c(4, 0, 1), c(3, 2, 3), c(10, 10, 10)), origin=c(3, 3, 3), spec=fspec) scRays1 <- scale(pRays, center=TRUE, scale=FALSE) scRays2 <- scale(pRays, center=c(3, 2, 1), scale=FALSE) scRays3 <- scale(pRays, center=pRect, scale=FALSE) stopifnot(rayOrigin(scRays1) == t(rep(0, 3)), identical(scRays2, scRays3)) chk <- rayVec(scale(pRays, center=FALSE, scale=TRUE)) # unit length stopifnot(isTRUE(all.equal(sqrt(rowSums(chk^2)), rep(1, nrow(chk))))) # Translate and scale ray vectors to match reference point and boundary # of a hyper-rectangle: scRays4 <- scale(pRays, center=c(3, 2, 1), scale=pRect) scRays5 <- scale(pRays, center=pRect, scale=pRect) stopifnot(identical(scRays4, scRays5)) raypts5 <- rayPosn1(scRays5) stopifnot(all(inBox(raypts5, corners=rectCorners(pRect), bdryOnly=TRUE, inptol=inpTol(fspec)))) scRays6 <- scale(pRays, center=FALSE, scale=pRect) raypts6 <- rayPosn1(scRays6) stopifnot(all(inBox(raypts6, corners=rectCorners(pRect), bdryOnly=TRUE, inptol=inpTol(fspec)))) # Scale lengths of ray vectors without changing origin. scRays7 <- scale(pRays, center=FALSE, scale=2) # *divides* by 2 len <- sqrt(rowSums(rayVec(pRays)^2)) len7 <- sqrt(rowSums(rayVec(scRays7)^2)) stopifnot(isTRUE(all.equal(len7, len/2))) factors <- c(1, 3, 0.5, -2) # negative scale reverses ray vector direction scRays8 <- scale(pRays, center=FALSE, scale=factors) # *divides* by 'factors' len <- sqrt(rowSums(rayVec(pRays)^2)) len8 <- sqrt(rowSums(rayVec(scRays8)^2)) stopifnot(isTRUE(all.equal(len8, len/abs(factors)))) # Not allowed to use scale to create infinite or degenerate ray vectors: factors[c(1, 3)] <- c(0, Inf) tools::assertError(scale(pRays, center=FALSE, scale=factors), verbose=warn.tools_verbose) #----- 2D example for rotation/scaling of ray vectors by a matrix. # Create ray vectors distributed along the ellipse # X1^2 + 2*X1*X2 + 2*(X2^2) = 8, # then translate to have origin at c(2, 1): fspec <- fnSpec(inpnames=c("X1", "X2")) X1 <- c(-4:4) X2 <- c((-X1 - sqrt(16 - X1^2))/2, (-X1 + sqrt(16 - X1^2))/2) X1 <- c(X1, X1) pRays <- inpRays(vecs=cbind(X1, X2), origin=c(2, 1), spec=fspec) # Plot them: #plot(pRays, eqAxes=TRUE, rayPar=list(col="blue")) # The untranslated ellipse can also be expressed as points X with a # Mahalanobis distance of 8 from the origin, given by the quadratic # form t(X) %*% V^(-1) %*% X, where symmetric, positive-definite matrix # V is V <- matrix(c(2, -1, -1, 1), nrow=2) stopifnot(isTRUE(all.equal(diag(cbind(X1, X2) %*% solve(V) %*% rbind(X1, X2)), rep(8, length(pRays))))) # all 8 # An 'inpScale' object that transforms this ellipse to a unit circle is: pScale <- inpScale(8*V, spec=fspec) # ('8*' because ultimately we want to divide vectors by sqrt(8) to get # unit length--'mat' acts quadratically not linearly on coordinates.) # Alternative: inpScale(x=sqrt(c(8, 8)), mat=V, spec=fspec) scRays1 <- scale(pRays, center=FALSE, scale=pScale) stopifnot(isTRUE(all.equal(rowSums(rayVec(scRays1)^2), rep(1, length(pRays))))) # unit length #plot(scRays1, add=TRUE, rayPar=list(col="red")) # All have unit length, but are not in the same directions as the original # vectors (rotated). To keep the same directions, instead scale each # vector individually by its length: raylens <- sqrt(rowSums(rayVec(pRays)^2)) scRays2 <- scale(pRays, center=FALSE, scale=raylens) stopifnot(isTRUE(all.equal(rowSums(rayVec(scRays2)^2), rep(1, length(pRays))))) # unit length #plot(scRays2, add=TRUE, rayPar=list(col="green")) cat("OK\n") } #----- cat(" update.inpRays ... ") fspec <- fnSpec(inpnames=c("X1", "X2")) set.seed(1) rays <- randomRays(10, spec=fspec) # ray vectors have unit length, # origin at (0, 0) rect2 <- inpRect(rbind(c(0, 0), c(1, 1)), spec=fspec) # unit square rays2 <- update(rays, origin=c(0.3, 0.4), rect=rect2) stopifnot(identical(as.vector(rayOrigin(rays2)), c(0.3, 0.4)), identical(rectCorners(rayRect(rays2)), rectCorners(rect2)), # Ref pt has been set to new origin: identical(rectRefPt(rayRect(rays2)), rayOrigin(rays2)), inRect(rayPosn1(rays2), rect=rect2, bdryOnly=TRUE)) rays3 <- update(rays2, origin=c(0.5, 0.6)) stopifnot(identical(as.vector(rayOrigin(rays3)), c(0.5, 0.6)), identical(rectCorners(rayRect(rays3)), rectCorners(rect2)), # Ref pt has been set to new origin: identical(rectRefPt(rayRect(rays3)), rayOrigin(rays3)), inRect(rayPosn1(rays3), rect=rect2, bdryOnly=TRUE)) cat("OK\n") #----- lsetPkgOpt(oldopt) cat('lsetPkgOpt.r', sep='\n') # User-level and internal package options. # # lsetPkgOpt # internalOpt # < create environment 'lvlSets_env'> # .onLoad # 09 Feb 2025 / Richard Raubertas # 17 Feb 2026 : Change 'tol.resp', 'tol.ptCoord' to 'resptol', 'inptol', for # consistency with other functions. oldopt <- lsetPkgOpt(reset=TRUE) #----- cat(" lsetPkgOpt ... ") #== Testing oldopts0 <- lsetPkgOpt(reset=TRUE) # don't assume current opts are default # Default 'factory-fresh' options: defopt <- lsetPkgOpt() #str(defopt) tol1 <- lsetPkgOpt("resptol") # plain vector stopifnot(is.vector(tol1), length(tol1) == 2) oldtol1 <- lsetPkgOpt(resptol=1e-4) # list of length 1 stopifnot(identical(tol1, oldtol1[[1]])) opts <- lsetPkgOpt(c("warn.invalidProfile", "warn.showBadRespPts")) stopifnot(is.list(opts), length(opts) == 2) oldopts <- lsetPkgOpt(warn.invalidProfile=0, warn.showBadRespPts=10) stopifnot(identical(opts, oldopts)) lsetPkgOpt(reset=TRUE) stopifnot(identical(lsetPkgOpt(), defopt)) oldopts2 <- lsetPkgOpt(list(warn.invalidProfile=0, warn.showBadRespPts=10), foo="bar") # 'foo' ignored stopifnot(identical(oldopts, oldopts2)) # Errors from non-existant options. tools::assertError(lsetPkgOpt(c("inptol", "foo")), verbose=warn.tools_verbose) tools::assertError(lsetPkgOpt(inptol=c(2e-5, 3e-8), foo="bar"), verbose=warn.tools_verbose) # Error from setting an invalid value. tools::assertError(lsetPkgOpt(inptol=c(2e-5, 3e-8), resptol=c(1e-1)), verbose=warn.tools_verbose) # Return to options where we started. lsetPkgOpt(oldopts0) stopifnot(identical(lsetPkgOpt(), oldopts0)) cat("OK\n") #----- cat(" internalOpt ... ") #== Testing oldopts0 <- internalOpt(reset=TRUE) # don't assume current opts are default # Default 'factory-fresh' options: defopt <- internalOpt() #str(defopt) stopifnot(is.list(oldopts0), names(oldopts0) == names(defopt)) oldval <- internalOpt(tol.line_position=c(1e-6, 1.5e-8)) # list w/ one comp stopifnot(is.list(oldval), length(oldval) == 1) chk <- internalOpt("tol.line_position") # new value (vector) stopifnot(!is.list(chk), identical(chk, c(1e-6, 1.5e-8))) internalOpt(reset=TRUE) stopifnot(identical(defopt, internalOpt())) tools::assertError(internalOpt("foo"), verbose=warn.tools_verbose) # non-existant option tools::assertError(internalOpt(foo=42), verbose=warn.tools_verbose) # non-existant option tools::assertError(internalOpt(tol.line_position=c(1e-2, 1e-8)), verbose=warn.tools_verbose) # invalid value tools::assertError(internalOpt(tol.line_position=c(1e-2, 1e-8)), verbose=warn.tools_verbose) # invalid value # Restore starting opts. internalOpt(oldopts0) cat("OK\n") #----- lsetPkgOpt(oldopt) cat('lsetSegs.r', sep='\n') # lsetSegs (S3 generic) # lsetSegs.default* (get '$segs' component from a list) # lsetSegs.lsetSegs* (no-op, but needed to avoid wrong method dispatch) # lsetSegs.respProfiles* # lsetSegs.respPairs* (create from 'respPairs' by adding 'inpLines') # lsetSegs.fnSpec* (create an empty object) # *= not public oldopt <- lsetPkgOpt(reset=TRUE) #----- cat(" lsetSegs.default ... ") Ex <- dbl_ellipse_X1bnd_2dEx # see '?dbl_ellipse_2dEx' fobj <- Ex$fobj segs <- Ex$segs chk <- lsetSegsCheck(segs, fnobj=fobj, posns=(1:4)/5) # list stopifnot(identical(lsetSegs(chk), Ex$segs_checked)) cat("OK\n") #----- cat(" lsetSegs.lsetSegs ... ") # Should be a no-op. Ex <- dbl_ellipse_X1bnd_2dEx # see '?dbl_ellipse_2dEx' segs <- Ex$segs stopifnot(identical(segs, lsetSegs(segs))) cat("OK\n") #----- cat(" lsetSegs.fnSpec ... ") #== Testing Ex <- dbl_ellipse_2dEx respF <- Ex$fobj fspec <- fnSpec(respF) thresh <- Ex$thresh lsegs <- lsetSegs(fspec, lsetthresh=thresh, checkArg=TRUE) # empty stopifnot(inherits(lsegs, "lsetSegs"), length(lsegs) == 0) # Non-empty 'lines', but no segments. lines <- axisRays(fspec) lsegs <- lsetSegs(fspec, lsetthresh=thresh, lines=lines, checkArg=TRUE) # empty stopifnot(inherits(lsegs, "lsetSegs"), length(lsegs) == 0) cat("OK\n") #----- lsetPkgOpt(oldopt) cat('lsetSegsCheck.r', sep='\n') #----- cat(" lsetSegsCheck ... ") #== Testing Ex <- dbl_ellipse_X1bnd_2dEx segs <- Ex$segs # 67 segs on 50 rays, 8 invalid segs chk <- lsetSegsCheck(segs, fnobj=Ex$fobj, posns=(1:4)/5) stopifnot(inherits(chk$segs, "lsetSegs"), identical(as.vector(table(chk$validIn)), c(8L, 59L)), identical(as.vector(table(chk$validOut)), c(75L))) # Should be able to handle an empty 'segs'. chk <- lsetSegsCheck(segs[0], fnobj=Ex$fobj, posns=(1:4)/5) stopifnot(inherits(chk$segs, "lsetSegs"), length(chk$segs) == 0, identical(as.vector(table(chk$validIn)), integer(0)), identical(as.vector(table(chk$validOut)), integer(0))) cat("OK\n") #----- cat('lsetSegsMethods.r', sep='\n') # Methods for 'lsetSegs' objects. (Most standard generics use the '.inpPts' # or '.inpPairs' methods for these objects.) # # plot.lsetSegs # pairs.lsetSegs # summary.lsetSegs # print.summary.lsetSegs* # plotSegsList (not strictly a method, but designed for 'lsetSegs' objects) # *= not public oldopt <- lsetPkgOpt(reset=TRUE) #----- cat(" summary.lsetSegs ... ") #== Testing # Test that function works even with empty input. Ex <- dbl_ellipse_2dEx respF <- Ex$fobj fspec <- fnSpec(respF) thresh <- Ex$thresh lsegs <- lsetSegs(fspec, lsetthresh=thresh, checkArg=TRUE) # empty summ <- summary(lsegs) stopifnot(summ$nlines == 0, dim(summ$bBox) == c(2, 2), all(is.na(summ$bBox)), length(summ$respRange) == 2, all(is.na(summ$respRange)), is.data.frame(summ$lineInfo), nrow(summ$lineInfo) == 0) # Non-empty 'lines', but no segments. lines <- axisRays(fspec) lsegs <- lsetSegs(fspec, lsetthresh=thresh, lines=lines, checkArg=TRUE) # empty summ <- summary(lsegs) stopifnot(summ$nlines == length(lines), dim(summ$bBox) == c(2, 2), all(is.na(summ$bBox)), length(summ$respRange) == 2, all(is.na(summ$respRange)), is.data.frame(summ$lineInfo), nrow(summ$lineInfo) == length(lines)) cat("OK\n") #----- lsetPkgOpt(oldopt) cat('respAccessors.r', sep='\n') # Functions to extract response-related components or properties of objects # in the 'ptsClass' set of classes. # lsetThresh # respInfo (generic and 'respPts'*, 'respPairs'*, 'lsetSegs', 'default' # methods) # *= not public # 24 Feb 2025 / Richard Raubertas # 12 Feb 2026 : Update 'lsetThresh' tests. oldopt <- lsetPkgOpt(reset=TRUE) #----- cat(" lsetThresh ... ") stopifnot(is.na(lsetThresh(NULL))) Ex <- dbl_ellipse_2dEx segs <- Ex$segs stopifnot(identical(lsetThresh(segs), Ex$thresh)) cat("OK\n") #----- cat(" respInfo ... ") respF <- function(x, inclGrad=FALSE) { resp <- x[, "X1"] + x[, "X2"]^2 + x[, "X3"]^3 if (inclGrad) { grad <- cbind(1, 2*x[, "X2"], 3*x[, "X3"]^2) structure(resp, gradient=grad) } else resp } fnobj <- fnObj(c("X1", "X2", "X3"), respfn=respF, hasgrad=TRUE) # NA point coordinates are not allowed. tools::assertError(respInfo(c(NA, 1, 3), fnobj=fnobj, inclGrad=TRUE), verbose=warn.tools_verbose) # Evaluate the response function at a set of points. pts <- rbind(c(0, 1, 3), c(2, 2, 4), c(5, -1, 2)) # three input space points rpts <- respInfo(pts, fnobj=fnobj, lsetthresh=28) stopifnot(is.data.frame(rpts), nrow(rpts) == nrow(pts), is.null(attr(rpts, "gradient"))) rpts <- respInfo(pts, fnobj=fnobj, lsetthresh=28, inclGrad=TRUE) grad <- attr(rpts, "gradient") stopifnot(is.data.frame(rpts), nrow(rpts) == nrow(pts), is.matrix(grad), dim(grad) == dim(pts), identical(dimnames(grad), list(NULL, inpNames(fnobj)))) # If 'fnobj' supports gradients but they were not actually calculated, # return all NA's. rpts <- respPts(pts, fnobj=fnobj, lsetthresh=28, inclGrad=FALSE) rinfo <- respInfo(rpts, inclGrad=TRUE) grad <- attr(rinfo, "gradient") stopifnot(is.data.frame(rinfo), nrow(rinfo) == nrow(pts), is.matrix(grad), dim(grad) == dim(pts), identical(dimnames(grad), list(NULL, inpNames(fnobj))), all(is.na(grad))) # Empty input. resp <- respInfo(pts[0, , drop=FALSE], fnobj=fnobj, inclGrad=TRUE) stopifnot(is.data.frame(resp), nrow(resp) == 0, is.matrix(attr(resp, "gradient")), nrow(attr(resp, "gradient")) == 0) cat("OK\n") #----- lsetPkgOpt(oldopt) cat('respProfiles.r', sep='\n') # respProfiles (S3 generic) # (no default method) # respProfiles.fnObj (calculate from scratch) # respProfiles.lsetSegs* (extract from an 'lsetSegs' object) # *= not public # 15 Feb 2025 / Richard Raubertas # 24 Feb 2025 : 'fnObj' no longer has 'spec' arg. # 16 Jul 2025 : Take advantage of, and test, 'update' method for 'fnObj'. # 13 Feb 2026 : 3D and 5D examples moved to 'tests_extra'. oldopt <- lsetPkgOpt(reset=TRUE) #----- cat(" respProfiles.fnObj ... \n") #== Testing #----- 1D example. cat(" 1D example ... ") respf <- function(x, ...) { x*cos(x) } fobj <- fnObj(c("X1"), respfn=respf) # default tolerances ray1 <- inpRays(c(1), spec=fobj) # unit vector rslt <- respProfiles(fobj, rays=ray1, posns=seq(-3, 10), lsetthresh=0.2) # Methods for 'respProfiles' objects: summ <- summary(rslt) stopifnot(length(rslt) == 1, is.list(summ), identical(inpLines(rslt), ray1), identical(summ$lsetThresh, lsetThresh(rslt)), is.matrix(summ$ptCounts), dim(summ$ptCounts) == c(summ$nlines, 5), colSums(summ$ptCounts) == c(19, 19, 0, 11, 5), is.matrix(summ$respRange), dim(summ$respRange) == c(summ$nlines, 2)) # Extract and summarize level set segments. lssegs <- lsetSegs(rslt) # '.respProfiles' method summ <- summary(lssegs) stopifnot(length(lssegs) == 3, is.list(summ), identical(inpLines(lssegs), ray1), identical(summ$lsetThresh, lsetThresh(lssegs)), is.matrix(summ$bBox), dim(summ$bBox) == c(2, inpDim(lssegs)), is.data.frame(summ$lineInfo), nrow(summ$lineInfo) == length(ray1), colSums(summ$lineInfo[, 1:4]) == c(length(lssegs), 0, 5, 1), summ$nlines == nrow(summ$lineInfo), length(summ$respRange) == 2) rpts <- respPts(rslt) summ <- summary(rpts) stopifnot(length(rpts) == 19, is.list(summ), identical(summ$lsetThresh, lsetThresh(rpts)), is.matrix(summ$ptCounts), dim(summ$ptCounts) == c(1, 3), colSums(summ$ptCounts) == c(19, 19, 11), length(summ$respRange) == 2) # Edge case summ <- summary(rpts[0]) stopifnot(is.list(summ), identical(summ$lsetThresh, lsetThresh(rpts)), is.matrix(summ$ptCounts), dim(summ$ptCounts) == c(1, 3), colSums(summ$ptCounts) == c(0, 0, 0), length(summ$respRange) == 2) # Implicit definition of feasible region (response function returns NA # to indicate infeasible points). Note that feasible region is closed. respf2 <- function(x, ...) { x <- as.vector(x) ifelse(x <= 1 | x >= 2, x*cos(x), NA) } fobj2 <- fnObj(inpNames(fobj), respfn=respf2, warnNAresp=FALSE) rslt2 <- respProfiles(fobj2, rays=ray1, lsetthresh=0.2, posns=c(seq(-3, 9), 1.3)) # need a test position in # the infeasible region to detect it rinfo2 <- respInfo(rslt2) stopifnot(nrow(rinfo2) == 18, rinfo2$feas == c(rep(TRUE, 7), FALSE, rep(TRUE, 10)), identical(rinfo2$feasbdry, c(NA, rep(FALSE, 5), TRUE, FALSE, TRUE, rep(FALSE, 8), NA)), rinfo2$lset == c(rep(TRUE, 3), rep(FALSE, 2), rep(TRUE, 2), rep(FALSE, 4), rep(TRUE, 5), rep(FALSE, 2)), identical(rinfo2$lsetbdry, c(NA, FALSE, TRUE, FALSE, FALSE, TRUE, rep(FALSE, 5), TRUE, rep(FALSE, 3), TRUE, FALSE, FALSE))) cat("OK\n") #----- Multimodal example in 2D (surface with two elliptical contours). cat(" 2D multimodal example ... ") respf <- function(x) { c1 <- x[, 1]^2 + 2*x[, 1]*x[, 2] + 2*(x[, 2]^2) # Same ellipse, but shifted to be centered at (8, 8): c2 <- (x[, 1]-8)^2 + 2*(x[, 1]-8)*(x[, 2]-8) + 2*((x[, 2]-8)^2) resp <- -1 * pmin(c1, c2) resp } # Feasible region: X1 >= 0. feasf <- function(x) { x <- num_matrix(x); (x[, 1] >= 0) } fobj <- fnObj(c("X1", "X2"), respfn=respf, feasfn=feasf) # Calculate ray profiles in four standard directions. rays1 <- axisRays(fobj, degree=2) rslt1 <- respProfiles(fobj, rays=rays1, posns=seq(-8, 15), lsetthresh=-40) summ <- summary(rslt1) stopifnot(length(rslt1) == length(rays1), identical(inpLines(rslt1), rays1), is.list(summ), identical(summ$lsetThresh, lsetThresh(rslt1)), is.matrix(summ$ptCounts), dim(summ$ptCounts) == c(summ$nlines, 5), colSums(summ$ptCounts) == c(101, 77, 3, 48, 8), is.matrix(summ$respRange), dim(summ$respRange) == c(summ$nlines, 2)) # '[' and 'c' methods test. stopifnot(identical(rslt1, c(rslt1[1], rslt1[0], NULL, rslt1[2:4], NULL))) # Extract level set segments. lssegs <- lsetSegs(rslt1) # .respProfiles method summ <- summary(lssegs) stopifnot(length(lssegs) == 6, is.list(summ), identical(inpLines(lssegs), rays1), identical(summ$lsetThresh, lsetThresh(lssegs)), is.matrix(summ$bBox), dim(summ$bBox) == c(2, inpDim(lssegs)), is.data.frame(summ$lineInfo), nrow(summ$lineInfo) == length(rays1), colSums(summ$lineInfo[, 1:4]) == c(length(lssegs), 3, 8, 1), summ$nlines == nrow(summ$lineInfo), length(summ$respRange) == 2) # Test that 'lsetbdry', 'feasbdry', 'lsetcens' are consistent. (Regression test) rinfo <- respInfo(lssegs) stopifnot(is.data.frame(rinfo), nrow(rinfo) == 2*length(lssegs), !anyNA(rinfo[, "feas"]), !anyNA(rinfo[, "lset"]), !any(rinfo[!rinfo$feas, "feasbdry"]), all(rinfo$feas[naF(rinfo$feasbdry)]), !any(rinfo[!rinfo$lset, "lsetbdry"]), all(rinfo$lset[naF(rinfo$lsetbdry)]), !anyNA(rinfo[, "lsetcens"]), all(rinfo$lsetcens == !(naF(rinfo$feasbdry) | naF(rinfo$lsetbdry)))) # Just feas bdry search, no level set. rslt2 <- respProfiles(fobj, rays=rays1, posns=seq(-8, 15), feasbdrySrch=TRUE) summ <- summary(rslt2) # same 'nfeasbdry' as for 'rslt' stopifnot(length(rslt2) == length(rays1), identical(inpLines(rslt2), rays1), is.list(summ), identical(summ$lsetThresh, lsetThresh(rslt2)), is.matrix(summ$ptCounts), dim(summ$ptCounts) == c(summ$nlines, 5), identical(unname(colSums(summ$ptCounts)), c(96, 72, 3, NA, NA)), is.matrix(summ$respRange), dim(summ$respRange) == c(summ$nlines, 2)) # Profile the lines that contain the rays. (Mostly equivalent.) lines1 <- inpLines(rays1) rslt1a <- respProfiles(fobj, lines=lines1, posns=seq(-8, 15), lsetthresh=-40) stopifnot(length(rslt1a) == length(rslt1), identical(summary(rslt1)$ptCounts, summary(rslt1a)$ptCounts)) # Edge cases: no rays; no posns. rslt3 <- respProfiles(fobj, rays=rays1[0], posns=seq(-8, 15), lsetthresh=-40) # no rays, no points summ <- summary(rslt3) stopifnot(inherits(rslt3, "respProfiles"), length(rslt3) == 0, nrow(summ$ptCounts) == 0) rslt4 <- respProfiles(fobj, rays=rays1, posns=numeric(0), lsetthresh=-40) # 4 rays, no points summ <- summary(rslt4) stopifnot(inherits(rslt4, "respProfiles"), length(rslt4) == 4, nrow(summ$ptCounts) == 4, all(summ$ptCounts == 0)) stopifnot(identical(rslt4, c(rslt4[1], NULL, rslt4[2:4], rslt4[0]))) if (requireNamespace("ggplot2", quietly=TRUE)) { tools::assertError(print(plot(rslt3)), verbose=warn.tools_verbose) # nothing to plot tools::assertError(print(plot(rslt4)), verbose=warn.tools_verbose) # nothing to plot } cat("OK\n") #----- lsetPkgOpt(oldopt) cat('respPts.r', sep='\n') # respPts* (S3 generic) # respPts.default* (handles matrices and 'inpPts' objects) # respPts.respPts* (extract or update an existing object. 'strictInfoCols' arg) # respPts.data.frame* (assemble 'respPts' manually. not public.) # respPts.fnSpec* (create an empty object) # respPts.respPairs* (includes 'which' arg. Also handles 'lsetSegs') # *= not public # 12 Feb 2025 / Richard Raubertas # 24 Feb 2025 : 'fnObj' no longer has 'spec' arg. oldopt <- lsetPkgOpt(reset=TRUE) #----- cat(" respPts ... ") #== Testing respF <- function(x, inclGrad=FALSE) { resp <- x[, "X1"] + x[, "X2"]^2 + x[, "X3"]^3 if (inclGrad) { grad <- cbind(1, 2*x[, "X2"], 3*x[, "X3"]^2) structure(resp, gradient=grad) } else resp } fnobj <- fnObj(c("X1", "X2", "X3"), respfn=respF, hasgrad=TRUE) coord <- rbind(c(0, 1, 3), c(2, 2, 4), c(5, -1, 2)) # three input space points # 'respPts' object created directly from a response function: rpts <- respPts(coord, fnobj=fnobj, lsetthresh=28, inclGrad=TRUE) # Construct the same 'respPts' from an 'inpPts' object. ipts <- inpPts(coord, spec=fnobj) rpts1b <- respPts(ipts, fnobj=fnobj, lsetthresh=28, inclGrad=TRUE) stopifnot(identical(rpts1b, rpts)) # Construct the same 'respPts' object starting with a data frame of # response function values: info <- respInfo(rpts, inclGrad=TRUE) # data frame w/ attribute stopifnot(is.matrix(attr(info, "gradient")), dim(attr(info, "gradient")) == dim(coord)) rpts2 <- respPts(info, spec=fnSpec(fnobj), ptcoord=coord, lsetthresh=28, respgrad=attr(info, "gradient")) stopifnot(identical(rpts, rpts2)) rpts2 <- respPts(info, spec=fnSpec(fnobj), ptcoord=coord, lsetthresh=28) stopifnot(identical(rpts, rpts2)) # tests default for 'respgrad' # Test code branch that skips most arg checking: rpts2b <- respPts(info, spec=fnSpec(fnobj), ptcoord=coord, lsetthresh=28, respgrad=attr(info, "gradient"), checkArg=FALSE) stopifnot(identical(rpts2, rpts2b)) # If one only wants the response values, the data frame can # be created directly, rather than going through 'respPts': respinfo <- respInfo(coord, fnobj=fnobj, lsetthresh=28) stopifnot(identical(respinfo, respInfo(rpts))) # Creating an empty 'respPts'. rpts0 <- respPts(matrix(0, nrow=0, ncol=inpDim(fnobj)), fnobj=fnobj) stopifnot(inherits(rpts0, "respPts"), length(rpts0) == 0) rpts0b <- respPts(fnSpec(fnobj)) stopifnot(identical(rpts0, rpts0b)) # Indexing. stopifnot(length(rpts[]) == length(rpts), length(rpts[2]) == 1, length(rpts[0]) == 0) # Combine multiple 'respPts' objects. stopifnot(identical(rpts, c(rpts[0], rpts[1], NULL, rpts[2:3]))) # Accessor functions. stopifnot(inherits(fnSpec(rpts), "fnSpec"), identical(unname(ptCoord(rpts)), coord)) # Updating 'lsetThresh', 'respTol'. rpts2 <- respPts(rpts, lsetthresh=28 + 1e-6) stopifnot(identical(lsetThresh(rpts2), c(28 + 1e-6)), respInfo(rpts2)$lset == c(TRUE, TRUE, FALSE)) # response value 28 is within tolerance of lsetthresh rpts2b <- respPts(rpts2, resptol=c(1e-8, 1e-8)) stopifnot(respInfo(rpts2b)$lset == c(FALSE, TRUE, FALSE)) # with tighter tolerance, response value 28 is no longer # considered to be in the level set # Can't combine objects with different level set definitions: tools::assertError(c(rpts, rpts2), verbose=warn.tools_verbose) tools::assertError(c(rpts2, rpts2b), verbose=warn.tools_verbose) # Can reverse updates: rpts3 <- respPts(rpts2, lsetthresh=28) # reverse update stopifnot(identical(rpts, rpts3)) rpts3b <- respPts(rpts2b, resptol=respTol(rpts)) stopifnot(identical(rpts2, rpts3b)) stopifnot(class(rpts) == c("respPts", "inpPts"), class(inpPts(rpts)) == "inpPts") stopifnot(identical(rpts, respPts(rpts))) cat("OK\n") #----- lsetPkgOpt(oldopt)