diff --git a/R/findZeros.R b/R/findZeros.R index 61864ee9..1370d500 100644 --- a/R/findZeros.R +++ b/R/findZeros.R @@ -1,13 +1,13 @@ #' Find zeros of functions -#' +#' #' Compute numerically zeros of a function or simultaneous zeros of multiple functions. -#' @param expr A formula. The right side names the variable with respect to which the zeros should be found. -#' The left side is an expression, e.g. `sin(x) ~ x`. -#' All free variables (all but the variable on the right side) named in the expression must be assigned +#' @param expr A formula. The right side names the variable with respect to which the zeros should be found. +#' The left side is an expression, e.g. `sin(x) ~ x`. +#' All free variables (all but the variable on the right side) named in the expression must be assigned #' a value via `\ldots` #' @param \dots Formulas corresponding to additional functions to use in simultaneous zero finding #' and/or specific numerical values for the free variables in the expression. -#' @param xlim The range of the dependent variable to search for zeros. `Inf` is a legitimate value, +#' @param xlim The range of the dependent variable to search for zeros. `Inf` is a legitimate value, #' but is interpreted in the numerical sense as the non-Inf largest floating point number. This can also #' be specified replacing `x` with the name of the variable. See the examples. #' @param near a value near which zeros are desired @@ -15,23 +15,23 @@ #' alternative to using `xlim` to specify the search space. #' @param nearest the number of nearest zeros to return. Fewer are returned if fewer are found. #' @param iterate maximum number of times to iterate the search. Subsequent searches take place with the range -#' of previously found zeros. Choosing a large number here is likely to kill performance without +#' of previously found zeros. Choosing a large number here is likely to kill performance without #' improving results, but a value of 1 (the default) or 2 works well when searching in `c(-Inf,Inf)` for #' a modest number of zeros near `near`. -#' @param npts How many sub-intervals to divide the `xlim` into when looking for candidates for zeros. +#' @param npts How many sub-intervals to divide the `xlim` into when looking for candidates for zeros. #' The default is usually good enough. -#' If `Inf` is involved, the intervals are logarithmically spaced up to the largest finite floating point number. +#' If `Inf` is involved, the intervals are logarithmically spaced up to the largest finite floating point number. #' There is no guarantee that all the roots will be found. #' @param sortBy specifies how the zeros found will be sorted. Options are 'byx', 'byy', or 'radial'. -#' +#' #' @details #' Searches numerically using `uniroot`. -#' +#' #' @return A dataframe of zero or more numerical values. Plugging these into the #' expression on the left side of the formula should result in values near zero. #' -#' @author Daniel Kaplan (\email{kaplan@@macalester.edu}) -#' +#' @author Daniel Kaplan (\email{kaplan@@macalester.edu}) +#' #' @examples #' findZeros( sin(t) ~ t, xlim=c(-10,10) ) #' # Can use tlim or t.lim instead of xlim if we prefer @@ -51,10 +51,10 @@ #' # Zeros in multiple dimensions (not run: these take a long time) #' # findZeros(x^2+y^2+z^2-5~x&y&z, nearest=3000, within = 5) #' # findZeros(x*y+z^2~z&y&z, z+y~x&y&z, npts=10) -#' @keywords calculus +#' @keywords calculus #' @export -findZeros <- function(expr, ..., xlim=c(near-within, near+within), near=0, within=Inf, +findZeros <- function(expr, ..., xlim=c(near-within, near+within), near=0, within=Inf, nearest=10, npts=1000, iterate=1, sortBy=c('byx', 'byy', 'radial')) { dots <- list(...) sortBy <- match.arg(sortBy) @@ -65,35 +65,35 @@ findZeros <- function(expr, ..., xlim=c(near-within, near+within), near=0, withi } else { # this is the original call ignore.limits <- FALSE } - + if( length(rhsVars) != 1 ){ if(any(within==Inf)) within[within==Inf]=100 return(unique(signif(findZerosMult(expr,..., npts=nearest, rad=within, near = near, sortBy = sortBy), 7))) } - + pfun <- function(x){ # removed . from name, was .x mydots <- dots mydots[[rhsVars]] <- x - eval( lhs(expr), envir=mydots, enclos=parent.frame() ) + eval( lhs(expr), envir=mydots, enclos=environment(expr) ) # was enclos=parent.env() } - + if (! ignore.limits ) { xlim <- inferArgs( dots=dots, vars=rhsVars, defaults=list(xlim=xlim) )[['xlim']] } tryCatch( xlim <- range(xlim), error = function(e) stop(paste('Improper limits value --', e))) - - if( xlim[1] >= xlim[2] ) + + if( xlim[1] >= xlim[2] ) stop(paste("Left limit (", xlim[1], ") must be less than right limit (", xlim[2], ").")) internal.near <- near - if ( internal.near < xlim[1] || internal.near > xlim[2]) { - internal.near <- mean(xlim[1],xlim[2]) + if ( internal.near < xlim[1] || internal.near > xlim[2]) { + internal.near <- mean(xlim[1],xlim[2]) } mx <- max(xlim - internal.near ) # max amount to add to internal.near when searching mn <- min(xlim - internal.near ) # min amount to add to internal.near when searching (will be negative) if (mx < 0) stop('Bug alert: near outside search interval.') if (mn > 0) stop('Bug alert: near outside search interval.') - + # Deal with very large numbers for the interval, e.g. Inf verybig <- .Machine$double.xmax plainbig <- npts^(.75) # linear spacing below this. @@ -102,16 +102,16 @@ findZeros <- function(expr, ..., xlim=c(near-within, near+within), near=0, withi rightseq <- NULL leftseq <- NULL middleseq <- NULL - if( mx > plainbig ) { + if( mx > plainbig ) { rightseq <- exp( seq( log(max(plainbig,mn)),log(mx),length=npts) ) } middleseq <- seq( max(-plainbig,mn), min(plainbig,mx), length=npts) if( mn < -plainbig ){ leftseq <- -exp( seq( log(-mn), log(-min(-plainbig,mx)), length=npts)) } - + searchx <- sort(unique(internal.near + c(0, leftseq, middleseq, rightseq))) - + y <- sapply( searchx, pfun ) testinds <- which( diff(sign(y)) != 0 ) if (length(testinds) < 1 ) { @@ -124,22 +124,22 @@ findZeros <- function(expr, ..., xlim=c(near-within, near+within), near=0, withi for (k in 1:N) { # look in subinterval k, i.e., between testinds[k] and testinds[k+1] if ( searchx[testinds[k]] < searchx[testinds[k]+1] ) { ur <- uniroot(function(qqzz){ sapply( qqzz, pfun) }, lower=searchx[testinds[k]], upper=searchx[testinds[k]+1]) - zeros[k] <- round( ur$root, digits=trunc(-log10(ur$estim.prec)) ) + zeros[k] <- round( ur$root, digits=2 + trunc(-log10(ur$estim.prec)) ) } else { warning("Potential bug alert: Attempting to search in region where signs of function at endpoints are equal. Skipping this interval.") } } } - + o <- order( abs(zeros - near) ) result <- sort(unique(zeros[o[1:min(nearest,length(zeros))]])) if ( iterate > 0 && length(result) > 1 ) { adjust <- min(diff(result)) # Note: negative value of iterate to indicate that we will be in a secondary iteration - return ( findZeros( expr, ..., xlim=range(c(result-adjust, result+adjust)), - near=near, within=within, - nearest=nearest, - npts=npts, + return ( findZeros( expr, ..., xlim=range(c(result-adjust, result+adjust)), + near=near, within=within, + nearest=nearest, + npts=npts, iterate= list(iterate=iterate - 1, ignore.limits = TRUE) ) ) } else { result <- data.frame(result) @@ -167,14 +167,14 @@ findZeros <- function(expr, ..., xlim=c(near-within, near+within), near=0, withi #'# cloud(z~x+y, data=sphere) #'@export -solve.formula <- function(form, ..., near=0, +solve.formula <- function(form, ..., near=0, within=Inf, nearest=10, npts=1000, iterate=1, sortBy=c('byx', 'byy', 'radial')){ dots = list(...) system = list(form) sortBy <- match.arg(sortBy) freeVars = list() - - + + #Separate formulae and free vars if(length(dots)>0){ for(i in (1:length(dots))){ @@ -187,7 +187,7 @@ solve.formula <- function(form, ..., near=0, } } } - + #change expression into formula. for(i in (1:length(system))){ formula = system[[i]] @@ -197,7 +197,7 @@ solve.formula <- function(form, ..., near=0, formula[[2]] <- exp system[[i]] <- formula } - - return(do.call(findZeros, c(system, freeVars, near=near, + + return(do.call(findZeros, c(system, freeVars, near=near, within=within, nearest=nearest, npts=npts, iterate=iterate, sortBy=sortBy))) } diff --git a/tests/testthat/test-findZeros.R b/tests/testthat/test-findZeros.R index cb6bcd8b..1bd5b8bd 100644 --- a/tests/testthat/test-findZeros.R +++ b/tests/testthat/test-findZeros.R @@ -1,8 +1,8 @@ context('Finding Zeros') test_that("zeros are found", { - expect_equivalent( round(findZeros(x^2 - 2 ~x)[,"x"],4), round(c(-sqrt(2), sqrt(2)),4) ) - expect_equivalent( round(findZeros(x^2 - 3 ~x)[,"x"],4), round(c(-sqrt(3), sqrt(3)),4) ) + expect_true( sum((findZeros(x^2 - 2 ~x)- c(-sqrt(2), sqrt(2)))^2) < 0.0001) + expect_true( sum((findZeros(x^2 - 3 ~x)- c(-sqrt(3), sqrt(3)))^2) < 0.0001) }) test_that("zeros are within search interval", { @@ -12,7 +12,7 @@ test_that("zeros are within search interval", { test_that("Can find zeros in two variables",{ Z = findZerosMult(a*x^2-v~a&x, v=8) expect_that(Z[1,]$a*Z[1,]$x^2-8, equals(0,tol=.001) ) - + Z = findZerosMult(a^2+x^2-8~a&x, npts = 1000, sortBy='radial') expect_that(Z[46,]$a^2+Z[46,]$x^2-8, equals(0, tol=0.001)) }) @@ -24,7 +24,7 @@ test_that("Works with Broyden",{ z = Z[,"z"] expect_that(x*y+z^2, equals(rep(0, 10), tol=0.001)) expect_that(z+y, equals(rep(0, 10), tol=0.001)) - + Z = findZeros(x*y+z^2+z^2*w~w&x&y&z, w*z+y~w&x&y&z, npts=10) w = Z[,"w"] x = Z[,"x"] @@ -32,7 +32,7 @@ test_that("Works with Broyden",{ z = Z[,"z"] expect_that(x*y+z^2+z^2*w, equals(rep(0,10), tol=0.001)) expect_that(w*z+y, equals(rep(0,10), tol=0.001)) - + Z = findZeros(x*y^2-9+z~x&y&z, x*y*z+z*w~x&y&z, w=10) x = Z[,"x"] y = Z[,"y"] @@ -40,7 +40,7 @@ test_that("Works with Broyden",{ w = 10 expect_that(x*y^2-9+z, equals(rep(0, 10), tol=0.001)) expect_that(x*y*z+z*w, equals(rep(0,10), tol=10)) - + }) # These were taken out by DTK on Aug 10, since they are somewhat stochastic. @@ -51,7 +51,7 @@ test_that("Works with Broyden",{ # z = Z[, "z"] # expect_that((x-500)^2+ (y-200)^2 + (z-300)^2 - 25, equals(rep(0, 1000), tol=0.01)) # }) -# +# # test_that("Solve function works properly", { # sphere = solve(x^2+y^2+z^2==5~x&y&z, within=5, nearest=1000) # x = sphere[,"x"]