Nothing
context('Test that xmin estimations are correct')
if ( exists('EXTENDED_TESTS') && EXTENDED_TESTS ) {
# Change dir if running tests manually
if ( file.exists('./tests/testthat') ) {
library(testthat)
setwd('./tests/testthat')
}
# Setup pli from Clauzet et al's
for ( s in dir('./pli-R-v0.0.3-2007-07-25',
full.names = TRUE, pattern = '*.R') ) {
source(s)
}
# Compile auxiliary binaries
system("cd ./pli-R-v0.0.3-2007-07-25/zeta-function/ && make")
system("cd ./pli-R-v0.0.3-2007-07-25/exponential-integral/ && make")
system("cd ./pli-R-v0.0.3-2007-07-25/ && \
gcc discpowerexp.c -lm -o discpowerexp && \
chmod +x discpowerexp")
test_that("xmins estimation is correct", {
parms <- expand.grid(expo = c(1.5, 1.3, 1.2),
rate = c(0.001, 0.005, 0.01, 0.1, 0.15, 0.2, 0.25, 0.3,
0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1))
estim_xmin <- function(df) {
pldat <- round(rpowerexp(10000, 1, df[ ,'expo'], df[, 'rate']))
est_xmin <- suppressWarnings( spatialwarnings:::xmin_estim(pldat) )
# est_xmin <- spatialwarnings:::xmin_estim(pldat)
# Create pl object and estimate its xmin
pl_obj <- poweRlaw::displ$new(pldat)
est_xmin_plpkg <- poweRlaw::estimate_xmin(pl_obj)[["xmin"]]
if ( !is.na(est_xmin) &&
!is.na(est_xmin_plpkg) &&
length(unique(pldat)) >= 5 ) {
# Note: In some pathological cases (few unique patches), there can be
# a small difference in xmin, so we use an acceptable error here.
fit_is_ok <- ( abs(est_xmin - est_xmin_plpkg) <= 1 )
if ( est_xmin != est_xmin_plpkg ) {
cat(" Ours: ", est_xmin, ' -> poweRlaw\'s: ', est_xmin_plpkg, " [",
length(unique(pldat)), " unique patches]", "\n", sep = "")
}
# In this case, inspect the fit provided by the poweRlaw package
if ( exists('GRAPHICAL') && GRAPHICAL && est_xmin != est_xmin_plpkg ) {
dev.new()
par(mfrow = c(1, 2))
plot(log10(cumpsd(pldat[pldat >= est_xmin])))
our_fit <- pl_fit(pldat, xmin = est_xmin)
xs <- unique(round(seq(min(pldat), max(pldat), length.out = 100)))
lines(log10(xs),
log10(ippl(xs, our_fit$plexpo, xmin = est_xmin)))
title("OUR FIT")
plot(log10(cumpsd(pldat[pldat >= est_xmin_plpkg])))
plpkg_expo <- poweRlaw::estimate_xmin(pl_obj)$pars
xs <- unique(round(seq(min(pldat), max(pldat), length.out = 100)))
lines(log10(xs),
log10(ippl(xs, plpkg_expo, xmin = est_xmin_plpkg)))
title("PWL FIT")
browser()
}
}
data.frame(df, est_xmin = est_xmin, fit_is_ok = fit_is_ok)
}
xmin_ests <- ddply(parms, ~ expo + rate, estim_xmin)
# Check that fits are overwhelmingly good. Sometimes the poweRlaw code will
# produce different fits, so we need to be able to accomodate some errors here.
# We put five as a maximum tolerable number of mismatches (out of 45)
number_of_ok_fits <- sum(xmin_ests[ ,"fit_is_ok"])
cat(paste0("Est. xmins OK: ", number_of_ok_fits, "/", nrow(parms), "\n"))
expect_true({
( nrow(parms) - number_of_ok_fits ) <= 5
})
})
# Remove auxiliary binaries now that tests are done
system("cd ./pli-R-v0.0.3-2007-07-25/zeta-function/ && rm zeta_func zeta_func.o")
system("cd ./pli-R-v0.0.3-2007-07-25/exponential-integral/ && rm exp_int exp_int.o")
system("cd ./pli-R-v0.0.3-2007-07-25/ && rm discpowerexp")
} else {
message('Skipping xmin estimation testing')
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.