Wrapping TISEAN Binaries in R
Chris Dissanayake
2023-02-17
TISEAN is an established software package for nonlinear time series analysis, commonly used in dynamical systems work. It’s painful to use in practice, mostly because passing parameters to the binaries by typing gets tedious fast when doing iterative work. This post shows how to wrap TISEAN binaries in R using the system() call, giving us a cleaner interface that plays nicely. The example uses a Henon map, a standard toy model which makes it a good test case. The workflow simulates a Henon map, adds nonlinear noise, applies Grassberger noise reduction via ghkss, and compares the embedded attractors before and after filtering using Takens theorem.
tisean_path <- "C:/Users/admin/Desktop/tisean/TISEAN_3.0.0-windows/Tisean_3.0.0/"
get_binaries <- function(tisean_path) {
if (is.null(tisean_path) || tisean_path == "") {
stop("Please provide a valid path to TISEAN.")
}
binary_location <- paste0(tisean_path, "/bin/")
fortran_location <- paste0(tisean_path, "/source_f/")
c_location <- paste0(tisean_path, "/source_c/")
aa <- list(
list.files(binary_location),
list.files(fortran_location),
list.files(c_location)
)
aa <- data.frame(
lapply(aa, "length<-", max(lengths(aa))),
stringsAsFactors = FALSE
)
names(aa) <- c("bin", "fortran", "c")
return(aa)
}
get_binaries(tisean_path = tisean_path)## bin fortran c
## 1 addnoise.exe addnoise.f ar-model.c
## 2 ar-model.exe any_s.f arima-model.c
## 3 ar-run.exe ar-run.f av-d2.c
## 4 arima-model.exe arguments.f boxcount.c
## 5 autocor.exe autocor.f corr.c
## 6 av-d2.exe c1.f d2.c
## 7 boxcount.exe c2d.f delay.c
## 8 c1.exe c2g.f extrema.c
## 9 c2d.exe c2naive.f false_nearest.c
## 10 c2g.exe c2t.f fsle.c
## 11 c2naive.exe choose.f ghkss.c
## 12 c2t.exe cluster.f histogram.c
## 13 choose.exe commandline.f lfo-ar.c
## 14 cluster.exe compare.f lfo-run.c
## 15 compare.exe d1.f lfo-test.c
## 16 corr.exe endtoend.f low121.c
## 17 d2.exe events.f lyap_k.c
## 18 delay.exe gpl.txt lyap_r.c
## 19 endtoend.exe help.f lyap_spec.c
## 20 events.exe henon.f lzo-gm.c
## 21 extrema.exe ikeda.f lzo-run.c
## 22 false_nearest.exe intervals.f lzo-test.c
## 23 fsle.exe istdio_temp.f Makefile.in
## 24 ghkss.exe lazy.f makenoise.c
## 25 henon.exe lorenz.f mem_spec.c
## 26 histogram.exe Makefile.in mutual.c
## 27 ikeda.exe neigh.f new.tgz
## 28 intervals.exe nmore.f nrlazy.c
## 29 lazy.exe normal.f nstat_z.c
## 30 lfo-ar.exe notch.f pca.c
## 31 lfo-run.exe others poincare.c
## 32 lfo-test.exe pc.f polyback.c
## 33 lorenz.exe predict.f polynom.c
## 34 low121.exe project.f polynomp.c
## 35 lyap_k.exe randomize polypar.c
## 36 lyap_r.exe rank.f rbf.c
## 37 lyap_spec.exe readfile.f recurr.c
## 38 lzo-gm.exe rms.f resample.c
## 39 lzo-run.exe slatec rescale.c
## 40 lzo-test.exe spectrum.f routines
## 41 makenoise.exe spikeauto.f sav_gol.c
## 42 mem_spec.exe spikespec.f xcor.c
## 43 mutual.exe store_spec.f xzero.c
## 44 notch.exe stp.f <NA>
## 45 nrlazy.exe surrogates.f <NA>
## 46 nstat_z.exe timerev.f <NA>
## 47 pc.exe tospec.f <NA>
## 48 pca.exe totospec.f <NA>
## 49 poincare.exe upo.f <NA>
## 50 polyback.exe upoembed.f <NA>
## 51 polynom.exe verbose.f <NA>
## 52 polynomp.exe wiener1.f <NA>
## 53 polypar.exe wiener2.f <NA>
## 54 predict.exe xc2.f <NA>
## 55 project.exe xreadfile.f <NA>
## 56 randomize_auto_exp_random.exe xrecur.f <NA>
## 57 randomize_autop_exp_random.exe <NA> <NA>
## 58 randomize_spikeauto_exp_random.exe <NA> <NA>
## 59 randomize_spikespec_exp_event.exe <NA> <NA>
## 60 randomize_uneven_exp_random.exe <NA> <NA>
## 61 rbf.exe <NA> <NA>
## 62 recurr.exe <NA> <NA>
## 63 resample.exe <NA> <NA>
## 64 rescale.exe <NA> <NA>
## 65 rms.exe <NA> <NA>
## 66 sav_gol.exe <NA> <NA>
## 67 spectrum.exe <NA> <NA>
## 68 spikeauto.exe <NA> <NA>
## 69 spikespec.exe <NA> <NA>
## 70 stp.exe <NA> <NA>
## 71 surrogates.exe <NA> <NA>
## 72 timerev.exe <NA> <NA>
## 73 upo.exe <NA> <NA>
## 74 upoembed.exe <NA> <NA>
## 75 wiener1.exe <NA> <NA>
## 76 wiener2.exe <NA> <NA>
## 77 xc2.exe <NA> <NA>
## 78 xcor.exe <NA> <NA>
## 79 xrecur.exe <NA> <NA>
## 80 xzero.exe <NA> <NA>call_routines <- function(
tisean_path,
time_series = NULL,
routine = "henon", # leave this default for troubleshooting
remove_header = FALSE,
show_output_on_console = FALSE,
...
) {
if (is.null(tisean_path) || tisean_path == "") {
stop("Please provide a valid path to TISEAN.")
}
# check system; OS x not tested.
if (!.Platform$OS.type == "unix") {
paste0(routine, ".exe")
}
dots <- list(...)
param_string <- paste0("-", paste0(names(dots), unlist(dots)), collapse = " ")
binary_location <- paste0(tisean_path, "/bin/", routine)
print(tempdir())
output_file <- tempfile()
output_file_name <- gsub("(^.*\\\\)(.*$)", "\\2", output_file)
if (!is.null(time_series)) {
# create a temp file name
input_file <- tempfile()
# dump data to the temp file
write.table(time_series, file = input_file, row.names = F)
shell_command <- paste(
binary_location,
input_file,
param_string,
"-o",
output_file
)
} else {
shell_command <- paste(binary_location, param_string, "-o", output_file)
}
# sanity check
print(shell_command)
system(shell_command, show.output.on.console = show_output_on_console)
# select the latest temp files from the directory.
# multi-file output corner case?
output_file <- file.info(dir(tempdir(), full.names = TRUE)) %>%
arrange(desc(ctime)) %>%
filter(grepl(output_file_name, rownames(.))) %>%
slice(1)
# remove temp file(s)
if (!is.null(time_series)) {
unlink(input_file)
}
det <- read.table(rownames(output_file), header = remove_header)
unlink(output_file)
return(det)
}Test using an example - a Hénon map.
- Simulate a Henon map.
- Add nonlinear noise to (1).
- Apply Takens theorem. This will be our reference.
- Remove noise from (2) using
ghkss. - Apply Takens theorem to (4), and plot over (3).
# henon map
henon_m <- call_routines(
tisean_path = tisean_path,
routine = "henon",
remove_header = FALSE,
show.output.on.console = FALSE,
l = 10000
)## [1] "C:\\Users\\admin\\AppData\\Local\\Temp\\Rtmpo1pbeM"
## [1] "C:/Users/admin/Desktop/tisean/TISEAN_3.0.0-windows/Tisean_3.0.0//bin/henon -show.output.on.console0 -l10000 -o C:\\Users\\admin\\AppData\\Local\\Temp\\Rtmpo1pbeM\\file57586f6f6b17"# add nonlinear noise
henon_m_n <- call_routines(
tisean_path = tisean_path,
time_series = henon_m,
remove_header = FALSE,
show.output.on.console = FALSE,
routine = "addnoise",
v = 0.04
)## [1] "C:\\Users\\admin\\AppData\\Local\\Temp\\Rtmpo1pbeM"
## [1] "C:/Users/admin/Desktop/tisean/TISEAN_3.0.0-windows/Tisean_3.0.0//bin/addnoise C:\\Users\\admin\\AppData\\Local\\Temp\\Rtmpo1pbeM\\file57586e127ffb -show.output.on.console0 -v0.04 -o C:\\Users\\admin\\AppData\\Local\\Temp\\Rtmpo1pbeM\\file575827177a62"# delay vectors
henon_m_n_d <- call_routines(
tisean_path = tisean_path,
time_series = henon_m_n,
remove_header = FALSE,
show.output.on.console = TRUE,
routine = "delay",
m = 2
)## [1] "C:\\Users\\admin\\AppData\\Local\\Temp\\Rtmpo1pbeM"
## [1] "C:/Users/admin/Desktop/tisean/TISEAN_3.0.0-windows/Tisean_3.0.0//bin/delay C:\\Users\\admin\\AppData\\Local\\Temp\\Rtmpo1pbeM\\file57583e4a3588 -show.output.on.console1 -m2 -o C:\\Users\\admin\\AppData\\Local\\Temp\\Rtmpo1pbeM\\file5758223f16a6"# Grassberger noise reduction
henon_m_n_ghkss <- call_routines(
tisean_path = tisean_path,
time_series = henon_m_n,
remove_header = FALSE,
show.output.on.console = TRUE,
routine = "ghkss",
m = "1,7",
q = 2,
r = 0.05,
k = 20,
i = 2
)## [1] "C:\\Users\\admin\\AppData\\Local\\Temp\\Rtmpo1pbeM"
## [1] "C:/Users/admin/Desktop/tisean/TISEAN_3.0.0-windows/Tisean_3.0.0//bin/ghkss C:\\Users\\admin\\AppData\\Local\\Temp\\Rtmpo1pbeM\\file5758fc03c4e -show.output.on.consoleTRUE -m1,7 -q2 -r0.05 -k20 -i2 -o C:\\Users\\admin\\AppData\\Local\\Temp\\Rtmpo1pbeM\\file57584f3f65f3"# delay vectors
henon_m_n_ghkss_d <- call_routines(
tisean_path = tisean_path,
time_series = henon_m_n_ghkss,
remove_header = FALSE,
show.output.on.console = TRUE,
routine = "delay",
m = 2
)## [1] "C:\\Users\\admin\\AppData\\Local\\Temp\\Rtmpo1pbeM"
## [1] "C:/Users/admin/Desktop/tisean/TISEAN_3.0.0-windows/Tisean_3.0.0//bin/delay C:\\Users\\admin\\AppData\\Local\\Temp\\Rtmpo1pbeM\\file575826b26423 -show.output.on.console1 -m2 -o C:\\Users\\admin\\AppData\\Local\\Temp\\Rtmpo1pbeM\\file57581db5557a"henon_m_n_d |>
ggplot(aes(V1, V2)) +
geom_point(aes(col = "Noisy Data"), alpha = 0.5, size = 0.1) +
geom_point(
data = henon_m_n_ghkss_d,
aes(V1, V2, col = "Filtered Data"),
size = 0.05,
alpha = 0.5
) +
coord_fixed() +
labs(
title = "Embedded Noisy Henon Attractor Filtered With ghkss",
subtitle = "10,000 iterations"
) +
theme(axis.title = element_blank(), legend.title = element_blank())