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 powerful but painful to use in practice; passing parameters through the binaries by hand gets tedious fast when doing iterative work. This post shows how to wrap TISEAN binaries in R using the system() call, giving you 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\\RtmpIFNjFe"
## [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\\RtmpIFNjFe\\file4b285dd11260"# 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\\RtmpIFNjFe"
## [1] "C:/Users/admin/Desktop/tisean/TISEAN_3.0.0-windows/Tisean_3.0.0//bin/addnoise C:\\Users\\admin\\AppData\\Local\\Temp\\RtmpIFNjFe\\file4b2875d05a96 -show.output.on.console0 -v0.04 -o C:\\Users\\admin\\AppData\\Local\\Temp\\RtmpIFNjFe\\file4b287a8d422b"# 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\\RtmpIFNjFe"
## [1] "C:/Users/admin/Desktop/tisean/TISEAN_3.0.0-windows/Tisean_3.0.0//bin/delay C:\\Users\\admin\\AppData\\Local\\Temp\\RtmpIFNjFe\\file4b28199a4c39 -show.output.on.console1 -m2 -o C:\\Users\\admin\\AppData\\Local\\Temp\\RtmpIFNjFe\\file4b28f935f60"# 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\\RtmpIFNjFe"
## [1] "C:/Users/admin/Desktop/tisean/TISEAN_3.0.0-windows/Tisean_3.0.0//bin/ghkss C:\\Users\\admin\\AppData\\Local\\Temp\\RtmpIFNjFe\\file4b28a1334aa -show.output.on.consoleTRUE -m1,7 -q2 -r0.05 -k20 -i2 -o C:\\Users\\admin\\AppData\\Local\\Temp\\RtmpIFNjFe\\file4b2835ec43d1"# 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\\RtmpIFNjFe"
## [1] "C:/Users/admin/Desktop/tisean/TISEAN_3.0.0-windows/Tisean_3.0.0//bin/delay C:\\Users\\admin\\AppData\\Local\\Temp\\RtmpIFNjFe\\file4b288de6413 -show.output.on.console1 -m2 -o C:\\Users\\admin\\AppData\\Local\\Temp\\RtmpIFNjFe\\file4b28154390d"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())