Wrapping TISEAN Binaries in R
Chris Dissanayake
2023-02-17
It is possible to wrap TISEAN binaries using Rcpp or
system, and the example below uses the latter.
Here are all the files in my TISEAN installation path:
tisean_path <- "/home/chrisdissanayake/Desktop/tisean/TISEAN_3.0.1_unix/Tisean_3.0.1"
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 addnoise.f ar-model.c
## 2 ar-model any_s.f arima-model.c
## 3 ar-run ar-run.f av-d2.c
## 4 arima-model arguments.f boxcount.c
## 5 autocor autocor.f corr.c
## 6 av-d2 c1.f d2.c
## 7 boxcount c2d.f delay.c
## 8 c1 c2g.f extrema.c
## 9 c2d c2naive.f false_nearest.c
## 10 c2g c2t.f fsle.c
## 11 c2naive choose.f ghkss.c
## 12 c2t cluster.f histogram.c
## 13 choose commandline.f lfo-ar.c
## 14 cluster compare.f lfo-run.c
## 15 compare d1.f lfo-test.c
## 16 corr endtoend.f low121.c
## 17 d2 events.f lyap_k.c
## 18 delay gpl.txt lyap_r.c
## 19 endtoend help.f lyap_spec.c
## 20 events henon.f lzo-gm.c
## 21 extrema ikeda.f lzo-run.c
## 22 false_nearest intervals.f lzo-test.c
## 23 fsle istdio_temp.f Makefile.in
## 24 ghkss lazy.f makenoise.c
## 25 henon lorenz.f mem_spec.c
## 26 histogram Makefile.in mutual.c
## 27 ikeda neigh.f new.tgz
## 28 intervals nmore.f nrlazy.c
## 29 lazy normal.f nstat_z.c
## 30 lfo-ar notch.f pca.c
## 31 lfo-run pc.f poincare.c
## 32 lfo-test predict.f polyback.c
## 33 lorenz project.f polynom.c
## 34 low121 randomize polynomp.c
## 35 lyap_k rank.f polypar.c
## 36 lyap_r readfile.f rbf.c
## 37 lyap_spec rms.f recurr.c
## 38 lzo-gm slatec resample.c
## 39 lzo-run spectrum.f rescale.c
## 40 lzo-test spikeauto.f routines
## 41 makenoise spikespec.f sav_gol.c
## 42 mem_spec store_spec.f xcor.c
## 43 mutual stp.f xzero.c
## 44 notch surrogates.f <NA>
## 45 nrlazy timerev.f <NA>
## 46 nstat_z tospec.f <NA>
## 47 pc totospec.f <NA>
## 48 pca upo.f <NA>
## 49 poincare upoembed.f <NA>
## 50 polyback verbose.f <NA>
## 51 polynom wiener1.f <NA>
## 52 polynomp wiener2.f <NA>
## 53 polypar xc2.f <NA>
## 54 predict xreadfile.f <NA>
## 55 project xrecur.f <NA>
## 56 randomize_auto_exp_random <NA> <NA>
## 57 randomize_autop_exp_random <NA> <NA>
## 58 randomize_spikeauto_exp_random <NA> <NA>
## 59 randomize_spikespec_exp_event <NA> <NA>
## 60 randomize_uneven_exp_random <NA> <NA>
## 61 rbf <NA> <NA>
## 62 recurr <NA> <NA>
## 63 resample <NA> <NA>
## 64 rescale <NA> <NA>
## 65 rms <NA> <NA>
## 66 sav_gol <NA> <NA>
## 67 spectrum <NA> <NA>
## 68 spikeauto <NA> <NA>
## 69 spikespec <NA> <NA>
## 70 stp <NA> <NA>
## 71 surrogates <NA> <NA>
## 72 timerev <NA> <NA>
## 73 upo <NA> <NA>
## 74 upoembed <NA> <NA>
## 75 wiener1 <NA> <NA>
## 76 wiener2 <NA> <NA>
## 77 xc2 <NA> <NA>
## 78 xcor <NA> <NA>
## 79 xrecur <NA> <NA>
## 80 xzero <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 Hénon 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] "/tmp/RtmpIgQOFD"
## [1] "/home/chrisdissanayake/Desktop/tisean/TISEAN_3.0.1_unix/Tisean_3.0.1/bin/henon -show.output.on.console0 -l10000 -o /tmp/RtmpIgQOFD/file23f2788fc6c7"# 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] "/tmp/RtmpIgQOFD"
## [1] "/home/chrisdissanayake/Desktop/tisean/TISEAN_3.0.1_unix/Tisean_3.0.1/bin/addnoise /tmp/RtmpIgQOFD/file23f24b146896 -show.output.on.console0 -v0.04 -o /tmp/RtmpIgQOFD/file23f25e927258"# 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] "/tmp/RtmpIgQOFD"
## [1] "/home/chrisdissanayake/Desktop/tisean/TISEAN_3.0.1_unix/Tisean_3.0.1/bin/delay /tmp/RtmpIgQOFD/file23f22d1adda2 -show.output.on.console1 -m2 -o /tmp/RtmpIgQOFD/file23f21ab17971"# 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] "/tmp/RtmpIgQOFD"
## [1] "/home/chrisdissanayake/Desktop/tisean/TISEAN_3.0.1_unix/Tisean_3.0.1/bin/ghkss /tmp/RtmpIgQOFD/file23f27e886d2c -show.output.on.consoleTRUE -m1,7 -q2 -r0.05 -k20 -i2 -o /tmp/RtmpIgQOFD/file23f226a04bf8"# 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] "/tmp/RtmpIgQOFD"
## [1] "/home/chrisdissanayake/Desktop/tisean/TISEAN_3.0.1_unix/Tisean_3.0.1/bin/delay /tmp/RtmpIgQOFD/file23f275117f11 -show.output.on.console1 -m2 -o /tmp/RtmpIgQOFD/file23f27f996c0a"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()
)