## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 6, fig.height = 4)

## ----setup--------------------------------------------------------------------
library(koopman.dmd)
set.seed(1)

## -----------------------------------------------------------------------------
t <- seq(0, 10, length.out = 200)
X <- rbind(sin(t), cos(t))
dim(X)   # 2 variables, 200 time steps

## -----------------------------------------------------------------------------
d <- dmd(X, rank = 2, dt = t[2] - t[1])
summary(d)

## -----------------------------------------------------------------------------
dmd_spectrum(d)

## -----------------------------------------------------------------------------
dmd_stability(d)

## -----------------------------------------------------------------------------
pred <- predict(d, n_ahead = 50)
dim(pred)

## ----fig.alt = "Observed signal with the DMD forecast appended"---------------
plot(t, X[1, ], type = "l", xlab = "time", ylab = "x1",
     xlim = c(0, 13), main = "DMD forecast")
t_future <- seq(max(t) + (t[2] - t[1]), by = t[2] - t[1], length.out = 50)
lines(t_future, pred[1, ], col = "red", lwd = 2)
legend("bottomleft", c("observed", "forecast"), col = c("black", "red"),
       lty = 1, bty = "n")

## -----------------------------------------------------------------------------
dmd_error(d)
dmd_residual(d)

## -----------------------------------------------------------------------------
A0 <- matrix(c(0.9, 0, 0.1, 0.8), 2, 2)
B0 <- matrix(c(0.5, 1), 2, 1)
m <- 120
X1 <- matrix(0, 2, m); X2 <- matrix(0, 2, m); U <- matrix(0, 1, m)
x <- c(1, -0.5)
for (i in seq_len(m)) {
  u_i <- sin(0.7 * (i - 1)) + 0.5 * cos(2.3 * (i - 1) + 1)
  X1[, i] <- x
  U[, i] <- u_i
  x <- as.numeric(A0 %*% x + B0 * u_i)
  X2[, i] <- x
}

fit <- dmdc(X1, X2, U, rank_input = 3)
round(fit$a, 6)   # recovers A0
round(fit$b, 6)   # recovers B0

## -----------------------------------------------------------------------------
fit2 <- dmdc(X1, X2, U, rank_input = 2, known_B = B0)
round(fit2$a, 6)

## -----------------------------------------------------------------------------
dmdc_stability(fit)
pred <- predict(fit, U = U)      # replaying the training input reproduces X2
max(abs(pred - X2))

## -----------------------------------------------------------------------------
t2 <- seq(0, 10, length.out = 200)
Xn <- rbind(sin(t2), sin(t2)^2)

plain  <- dmd(Xn, rank = 2)
lifted <- dmd(Xn, lifting = "polynomial", lifting_param = 2)

plain_err  <- dmd_error(plain)
lifted_err <- dmd_error(lifted)
c(plain = plain_err$rmse, lifted = lifted_err$rmse)

## -----------------------------------------------------------------------------
tt <- seq(0, 4 * pi, length.out = 200)
y  <- matrix(sin(tt), nrow = 1)   # a single row

h <- hankel_dmd(y, delays = 20)
h

## -----------------------------------------------------------------------------
pred_h <- predict(h, n_ahead = 20)
dim(pred_h)

## -----------------------------------------------------------------------------
g <- gla(X, n_eigenvalues = 2)
g

## -----------------------------------------------------------------------------
traj <- generate_trajectory("standard", c(0.1, 0.2), 2000, epsilon = 0.9)
dim(traj)

## ----fig.alt = "Orbit of the Chirikov standard map in phase space"------------
plot(traj[1, ], traj[2, ], pch = ".", xlab = "x", ylab = "p",
     main = "Chirikov standard map, epsilon = 0.9")

## -----------------------------------------------------------------------------
regular <- harmonic_time_average("standard", c(0.5, 0.0), "sin_pi", 0.5, 2000,
                                 epsilon = 0.9)
chaotic <- harmonic_time_average("standard", c(0.1, 0.2), "sin_pi", 0.5, 2000,
                                 epsilon = 0.9)
c(regular = regular$magnitude, chaotic = chaotic$magnitude)

## -----------------------------------------------------------------------------
conv <- hta_convergence("standard", c(0.5, 0.0), "sin_pi", 0.5, 2000,
                        epsilon = 0.9)
conv$dynamics_type

## ----fig.alt = "Mesochronic harmonic plot of the standard map"----------------
mhp <- mesochronic_compute("standard", c(0, 1), c(0, 1), 40, "sin_pi", 0.5, 300,
                           epsilon = 0.9)
image(mhp$x_coords, mhp$y_coords, mhp$hta_matrix,
      col = hcl.colors(64, "YlGnBu", rev = TRUE),
      xlab = "x", ylab = "p", main = "Mesochronic harmonic plot")

## -----------------------------------------------------------------------------
labels <- classify_phase_space(as.vector(mhp$hta_matrix))
table(factor(labels, levels = 1:3,
             labels = c("resonating", "chaotic", "non-resonating")))