Skip to content

Changing environment referred to in pfun() #782

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
dtkaplan wants to merge 2 commits into
base: master
Choose a base branch
from
Open

Changing environment referred to in pfun() #782

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
c0d28f0
Changing environment referred to in pfun()
dtkaplan Feb 17, 2021
fb0e955
Precision problem in findZeros.
dtkaplan Feb 17, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
76 changes: 38 additions & 38 deletions R/findZeros.R
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,37 +1,37 @@
#' 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
#' @param within only look for zeros at least this close to near. `near` and `within` provide an
#' 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
Expand All @@ -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)
Expand All @@ -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.
Expand All @@ -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 ) {
Expand All @@ -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)
Expand Down Expand Up @@ -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))){
Expand All @@ -187,7 +187,7 @@ solve.formula <- function(form, ..., near=0,
}
}
}

#change expression into formula.
for(i in (1:length(system))){
formula = system[[i]]
Expand All @@ -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)))
}
14 changes: 7 additions & 7 deletions tests/testthat/test-findZeros.R
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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", {
Expand All @@ -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))
})
Expand All @@ -24,23 +24,23 @@ 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"]
y = Z[,"y"]
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"]
z = Z[,"z"]
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.
Expand All @@ -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"]
Expand Down