## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
library(rareflow)

## -----------------------------------------------------------------------------
Qobs <- c(0.05, 0.90, 0.05)
px <- function(z) c(0.3, 0.4, 0.3)

## -----------------------------------------------------------------------------
flow <- makeflow("planar", list(u = 0.1, w = 0.2, b = 0))
fit <- fitflowvariational(Qobs, pxgivenz = px, nmc = 500)
fit$elbo

## -----------------------------------------------------------------------------
b <- function(x) -x
dt <- 0.01
T <- 1000
set.seed(1)
Winc <- rnorm(T, sd = sqrt(dt))

## -----------------------------------------------------------------------------
theta <- rep(0.5, T)
girsanov_logratio(theta, Winc, dt)


## -----------------------------------------------------------------------------
set.seed(123)
dt  <- 0.01
T   <- 1000
Winc <- rnorm(T, sd = sqrt(dt))
theta_path <- rep(0, length(Winc))
fit_gir <- fitflow_girsanov(
  observed      = Qobs,
  base_pxgivenz = px,
  theta_path    = theta_path,
  Winc          = Winc,
  dt            = dt
)


## -----------------------------------------------------------------------------
V <- function(x) x^4/4 - x^2/2
xs <- seq(-2, 2, length.out = 400)

plot(xs, V(xs), type = "l", lwd = 2,
     main = "Double-Well Potential",
     ylab = "V(x)", xlab = "x")
abline(v = c(-1, 1), lty = 3, col = "gray")


## -----------------------------------------------------------------------------
b<- function(x) {
  v<- as.numeric(x)
  return(c(v - v^3)) #double-well drift
}
qp<- FW_quasipotential(-1, 1, drift = b, T = 200, dt = 0.01)
qp$action

## -----------------------------------------------------------------------------
plot(qp$path, type = "l", main = "Minimum-Action Path (Freidlin–Wentzell)")


## -----------------------------------------------------------------------------
V2 <- function(x, y) 0.25 * (x^2 + y^2 - 1)^2

xs <- seq(-2, 2, length.out = 200)
ys <- seq(-2, 2, length.out = 200)
grid <- expand.grid(x = xs, y = ys)
Z <- matrix(V2(grid$x, grid$y), nrow = length(xs))

contour(xs, ys, Z,
        nlevels = 20,
        main = "2D Radial Double-Well Potential",
        xlab = "x", ylab = "y")


## -----------------------------------------------------------------------------
dt <- 0.01
T <- 5000
x <- matrix(0, nrow = T, ncol = 2)

b2 <- function(v) {
  r2 <- sum(v^2)
  -(r2 - 1) * v
}

for (t in 1:(T-1)) {
  drift <- b2(x[t, ])
  noise <- sqrt(dt) * rnorm(2)
  x[t+1, ] <- x[t, ] + drift * dt + noise
}

plot(x[,1], x[,2], type = "l",
     main = "2D Diffusion in a Radial Double-Well",
     xlab = "x", ylab = "y")


## -----------------------------------------------------------------------------
# Potential and drift
V2 <- function(x, y) 0.25 * (x^2 + y^2 - 1)^2
b2 <- function(v) {
  r2 <- sum(v^2)
  -(r2 - 1) * v
}

# String method parameters
N <- 80          # number of points in the string
steps <- 200     # number of iterations
dt_sm <- 0.01    # step size

# Initial straight-line path
path <- cbind(seq(-1, 1, length.out = N), rep(0, N))

# String method iterations
for (k in 1:steps) {
  # Update each interior point
  for (i in 2:(N-1)) {
    drift <- b2(path[i, ])
    path[i, ] <- path[i, ] + dt_sm * drift
  }
  
  # Reparametrize to keep points evenly spaced
  d <- sqrt(rowSums(diff(path)^2))
  s <- c(0, cumsum(d))
  s <- s / max(s)
  path <- cbind(
    approx(s, path[,1], xout = seq(0,1,length.out=N))$y,
    approx(s, path[,2], xout = seq(0,1,length.out=N))$y
  )
}

plot(path[,1], path[,2], type="l", lwd=2,
     main="2D Minimum Action Path (String Method)",
     xlab="x", ylab="y")
points(c(-1,1), c(0,0), pch=19, col="red")


## -----------------------------------------------------------------------------
library(ggplot2)
library(gganimate)


## ----eval = FALSE-------------------------------------------------------------
# # Simulate a 2D trajectory
# dt <- 0.01
# T <- 2000
# x <- matrix(0, nrow = T, ncol = 2)
# 
# for (t in 1:(T-1)) {
#   drift <- b2(x[t, ])
#   noise <- sqrt(dt) * rnorm(2)
#   x[t+1, ] <- x[t, ] + drift * dt + noise
# }
# 
# df <- data.frame(
#   t = 1:T,
#   x = x[,1],
#   y = x[,2]
# )
# 
# p <- ggplot(df, aes(x, y)) +
#   geom_path(alpha = 0.4) +
#   geom_point(aes(frame = t), color = "red", size = 2) +
#   coord_equal() +
#   labs(title = "2D Diffusion Trajectory", x = "x", y = "y")
# 
# animate(p, nframes = 200, fps = 20)
# 

## -----------------------------------------------------------------------------
b <- function(x) x - x^3
dt <- 0.01
T <- 2000
x <- numeric(T)
for (t in 1:(T-1)) {
  x[t+1] <- x[t] + b(x[t])*dt + sqrt(dt)*rnorm(1)
}


## -----------------------------------------------------------------------------
rare_event <- mean(x > 1.5)
rare_event


## -----------------------------------------------------------------------------
Qobs <- c(mean(x < -0.5), mean(abs(x) <= 0.5), mean(x > 0.5))
px <- function(z) c(0.3, 0.4, 0.3)


## -----------------------------------------------------------------------------
fit <- fitflowvariational(Qobs, pxgivenz = px)
fit$elbo