####=====================  Flocking (R) — tidy-friendly + annotated  =====================####
suppressPackageStartupMessages({ library(tidyverse) })  # mostly for pipes/readability
##---------------------- WORLD / TIME ----------------------##
N           <- 200  # number of birds (100–1200). High = richer patterns, slower.
steps       <- 400   # total frames to simulate/render
L           <- 100    # world side length; plot shows [0,L]×[0,L]
dt          <- 1      # time step; usually keep at 1
speed       <- 1.5   # units per step (~0.01–0.02*L typical). Higher = stretchier shapes.
seed        <- 42     # RNG seed for reproducibility
##---------------------- NEIGHBORS / INTERACTION ----------------------##
k_neighbors <- 3     # 3–10. Low = local swirl/fragmentation; high = fast global alignment.
use_separation <- TRUE   # avoid collisions (short-range repulsion)
r_rep       <- .05    # repulsion radius (~0.02–0.05*L). Higher = puffier spacing.
sep_weight  <- 0.30   # 0..1. Higher = stronger push-away when too close.
##---------------------- COHESION + SWIRL + NOISE ----------------------##
use_cohesion <- TRUE  # pull toward mid-range neighbors (keeps group compact)
r_att        <- 20    # attraction radius (~0.15–0.30*L). Higher = tighter balls.
coh_weight   <- 0.10  # 0.05–0.20. Blend toward neighbors’ centroid.
swirl_weight <- 0.08  # 0.04–0.12. Tangential bias around center (prevents straight-line lock-in).
eta_start    <- 0.25  # early noise amplitude (radians). Lower = cleaner; higher = more texture.
eta_end      <- 0.40  # late noise (raise slightly vs start to avoid polarized drift).
##---------------------- COLOR DYNAMICS (NEIGHBOR-DRIVEN) ----------------------##
adopt_prob   <- 0.8   # 0.3–0.9. Higher = faster “voter model” coarsening to few colors.
self_bias    <- 0     # 0–4 extra votes for current color. Higher = stickier colors.
# Optional guard to keep more colors alive (set TRUE to enable guarded logic below)
use_guarded_colors <- TRUE  # TRUE→ require strong majority + tiny mutation
adopt_thresh       <- 0.55    # if enabled: require ≥55% neighbors for adoption
mutate_prob        <- 0.00    # if enabled: small chance to hop to random color
##---------------------- PALETTE ----------------------##
pal_user <- c(         # your custom colors
  "#30C5D2","#3698B8","#3C6B9E",
  "#423E84","#471069"
)
brighten_factor <- 1.15  # >1 brightens colors a touch (handy on black bg)
##---------------------- FADE IN/OUT ----------------------##
fade_in     <- 100   # frames to ramp in from black at the start (0 = off)
fade_out    <- 100  # frames to ramp out to black at the end   (0 = off)
fade_gamma  <- 1.2   # easing: 1=linear; >1 gentle; <1 snappy
##---------------------- BOUNDARY / RENDER ----------------------##
boundary_mode <- "wrap"   # "wrap"=torus; "reflect"=bounce off edges; "vanish" = birds flying straight
plot_every    <- 4       # render every k-th step (2 halves file count)
black_bg      <- FALSE     # TRUE=black background
arrow_len     <- 1        # heading arrow length (plot units)
pt_size       <- 1        # point size
##---------------------- OUTPUT (unique run dir + params) ----------------------##
runs_root <- "output"                                    # <— root changed per your request
run_tag   <- format(Sys.time(), "%Y%m%d_%H%M%S")         # e.g., 20250813_213455
out_dir   <- file.path(runs_root, paste0("run_", run_tag))
dir.create(out_dir, recursive = TRUE, showWarnings = FALSE)
##====================== HELPERS ======================##
set.seed(seed)
wrap <- function(x, L) x %% L                                        # periodic wrap
torus_delta <- function(d, L) d - L * round(d / L)                   # shortest diff on torus
angle_mean <- function(th) atan2(mean(sin(th)), mean(cos(th)))       # circular mean
normalize2 <- function(v){ m <- sqrt(sum(v*v)); if (m==0) c(1,0) else v/m }
blend_angles <- function(a, b, w){                                   # slerp on unit circle
  va <- c(cos(a), sin(a)); vb <- c(cos(b), sin(b))
  v  <- normalize2((1-w)*va + w*vb); atan2(v[2], v[1])
}
make_palette <- function(n, black_bg = TRUE){                        # default palette
  lum <- if (black_bg) 85 else 35
  grDevices::hcl(h = seq(10, 370, length.out = n + 1)[-1], c = 70, l = lum)
}
majority_int <- function(v, n){                                      # majority among ints
  tab <- tabulate(v, nbins = n); winners <- which(tab == max(tab))
  if (length(winners) == 1) winners else sample(winners, 1)
}
# Safe brighten (no pmin() dimension drop)
brighten <- function(cols, f = 1.15){
  M <- col2rgb(cols) / 255; M <- M * f; M[M > 1] <- 1; M[M < 0] <- 0
  rgb(M[1,], M[2,], M[3,])
}
fade_alpha <- function(t, steps, fade_in, fade_out, gamma = 1){      # fade curve
  a_in  <- if (fade_in  > 0 && t <= fade_in)  (t / fade_in)^gamma else 1
  out_p <- if (fade_out > 0) max(0, (t - (steps - fade_out)) / fade_out) else 0
  a_out <- if (fade_out > 0 && t > steps - fade_out) (1 - out_p)^gamma else 1
  min(a_in, a_out)
}
##====================== INIT ======================##
start_time <- Sys.time()
# palette hookup
pal <- if (length(pal_user) > 0) pal_user else make_palette(10, black_bg)
if (black_bg) pal <- brighten(pal, brighten_factor)
n_colors <- length(pal)
# state
x     <- runif(N, 0, L)
y     <- runif(N, 0, L)
theta <- runif(N, -pi, pi)
col_idx <- sample.int(n_colors, N, replace = TRUE)
# write params next to frames
params <- list(
  timestamp = run_tag, out_dir = out_dir, seed = seed,
  N = N, steps = steps, L = L, dt = dt, speed = speed,
  k_neighbors = k_neighbors, use_separation = use_separation,
  r_rep = r_rep, sep_weight = sep_weight,
  use_cohesion = use_cohesion, r_att = r_att, coh_weight = coh_weight,
  swirl_weight = swirl_weight, eta_start = eta_start, eta_end = eta_end,
  adopt_prob = adopt_prob, self_bias = self_bias,
  use_guarded_colors = use_guarded_colors, adopt_thresh = adopt_thresh, mutate_prob = mutate_prob,
  fade_in = fade_in, fade_out = fade_out, fade_gamma = fade_gamma,
  black_bg = black_bg, arrow_len = arrow_len, pt_size = pt_size,
  boundary_mode = boundary_mode, palette = pal, palette_user = pal_user
)
if (requireNamespace("jsonlite", quietly = TRUE)) {
  jsonlite::write_json(params, file.path(out_dir, "params.json"),
                       pretty = TRUE, auto_unbox = TRUE)
} else {
  dput(params, file = file.path(out_dir, "params.R"))
}
##====================== SIMULATION ======================##
frame_id <- 0L
NN <- vector("list", N)  # neighbor indices cache
for (t in 1:steps) {
  theta_new <- theta
  
  for (i in 1:N) {
    # neighbor deltas (torus for wrap; Euclidean for reflect)
    if (boundary_mode == "wrap") {
      dx <- torus_delta(x - x[i], L); dy <- torus_delta(y - y[i], L)
    } else {
      dx <- x - x[i]; dy <- y - y[i]
    }
    d2 <- dx^2 + dy^2; d2[i] <- Inf
    
    nn_idx  <- head(order(d2), k_neighbors); NN[[i]] <- nn_idx
    desired <- angle_mean(theta[nn_idx])                 # alignment
    
    if (use_separation) {                                # short-range repulsion
      near <- which(d2 < r_rep^2)
      if (length(near) > 0) {
        w <- 1 / (d2[near] + 1e-6)
        sep_vec   <- c(-sum(dx[near]*w), -sum(dy[near]*w))
        sep_angle <- atan2(sep_vec[2], sep_vec[1])
        desired   <- blend_angles(desired, sep_angle, sep_weight)
      }
    }
    
    if (use_cohesion) {                                  # mid-range attraction
      far <- which(d2 > r_rep^2 & d2 < r_att^2)
      if (length(far) > 0) {
        cx <- mean(dx[far]); cy <- mean(dy[far])
        coh_angle <- atan2(cy, cx)
        desired   <- blend_angles(desired, coh_angle, coh_weight)
      }
    }
    
    swirl_angle <- atan2(L/2 - y[i], L/2 - x[i]) + pi/2  # gentle orbiting
    desired     <- blend_angles(desired, swirl_angle, swirl_weight)
    
    # time-varying angular noise
    s_frac  <- if (steps > 1) (t - 1) / (steps - 1) else 0
    eta_t   <- eta_start + s_frac * (eta_end - eta_start)
    desired <- desired + runif(1, -eta_t, eta_t)
    
    theta_new[i] <- atan2(sin(desired), cos(desired))    # wrap to [-pi,pi]
  }
  
  # move + boundaries
  dx_step <- dt * speed * cos(theta_new)
  dy_step <- dt * speed * sin(theta_new)
  
  if (boundary_mode == "wrap") {
    x <- wrap(x + dx_step, L); y <- wrap(y + dy_step, L); theta <- theta_new
  } else if (boundary_mode == "vanish") {
    # Birds fly straight, disappear if out of bounds
    x_new <- x + dx_step
    y_new <- y + dy_step
    
    alive <- x_new >= 0 & x_new <= L & y_new >= 0 & y_new <= L
    
    # keep only those inside box
    x <- x_new[alive]
    y <- y_new[alive]
    theta <- theta_new[alive]
    col_idx <- col_idx[alive]
    N <- sum(alive)  # shrink flock size as they vanish
  } else {  # reflect (bounce)
    x_new <- x + dx_step; y_new <- y + dy_step; th <- theta_new
    hit_left  <- x_new < 0; hit_right <- x_new > L
    if (any(hit_left))  { x_new[hit_left]  <- -x_new[hit_left];       th[hit_left]  <- pi - th[hit_left] }
    if (any(hit_right)) { x_new[hit_right] <- 2*L - x_new[hit_right]; th[hit_right] <- pi - th[hit_right] }
    hit_bot <- y_new < 0; hit_top <- y_new > L
    if (any(hit_bot)) { y_new[hit_bot] <- -y_new[hit_bot]; th[hit_bot] <- -th[hit_bot] }
    if (any(hit_top)) { y_new[hit_top] <- 2*L - y_new[hit_top]; th[hit_top] <- -th[hit_top] }
    x <- x_new; y <- y_new; theta <- atan2(sin(th), cos(th))
  } 
  
  # COLOR DYNAMICS
  if (!use_guarded_colors) {
    # simple neighbor-majority with optional inertia
    col_idx_new <- col_idx
    for (i in 1:N) {
      neigh_cols <- col_idx[NN[[i]]]
      votes <- c(neigh_cols, rep(col_idx[i], self_bias))
      winner <- majority_int(votes, n_colors)
      if (runif(1) < adopt_prob) col_idx_new[i] <- winner
    }
    col_idx <- col_idx_new
  } else {
    # guarded adoption + tiny mutation (helps keep >3 colors alive)
    col_idx_new <- col_idx
    for (i in 1:N) {
      neigh <- NN[[i]]
      counts_raw <- tabulate(col_idx[neigh], nbins = n_colors)
      winner_raw <- which.max(counts_raw)
      prop       <- counts_raw[winner_raw] / length(neigh)
      counts_bias <- counts_raw; counts_bias[col_idx[i]] <- counts_bias[col_idx[i]] + self_bias
      winner_biased <- which.max(counts_bias)
      if (runif(1) < adopt_prob && prop >= adopt_thresh) {
        col_idx_new[i] <- winner_biased
      } else if (runif(1) < mutate_prob) {
        choices <- setdiff(seq_len(n_colors), col_idx[i])
        if (length(choices)) col_idx_new[i] <- sample(choices, 1)
      }
    }
    col_idx <- col_idx_new
  }
  
  # DRAW (PNG frames)
  frame_id <- frame_id + 1L
  png(file.path(out_dir, sprintf("frame_%04d.png", frame_id)),
      width = 800, height = 800, bg = if (black_bg) "black" else "white")
  
  par(mar = c(2,2,2,2))
  plot(NA, xlim = c(0, L), ylim = c(0, L), xaxs = "i", yaxs = "i",
       xlab = "", ylab = "", axes = FALSE)
  
  opac <- fade_alpha(t, steps, fade_in, fade_out, fade_gamma)
  cols <- adjustcolor(pal[col_idx], alpha.f = opac)
  
  # points(x, y, pch = 16, cex = pt_size, col = cols)
  # segments(x, y, x - arrow_len*cos(theta), y - arrow_len*sin(theta), lwd = 1, col = cols)
  arrows(
    x0 = x, 
    y0 = y, 
    x1 = x + arrow_len * cos(theta), 
    y1 = y + arrow_len * sin(theta),
    length = 0.15,    # size of the arrowhead
    angle = 25,       # opening angle of the head
    col = cols,
    lwd = 1
  )
  
  dev.off()
}
##====================== EXPORT (optional) ======================##
# Make a GIF from the frames (writes next to PNGs)
if (requireNamespace("gifski", quietly = TRUE)) {
  pngs <- Sys.glob(file.path(out_dir, "frame_*.png")) |> sort()
  gifski::gifski(pngs, gif_file = file.path(out_dir, "flock.gif"),
                 width = 800, height = 800, delay = 1/30)
}
# MP4 alternative (smaller file) — uncomment if you use {av}
# if (requireNamespace("av", quietly = TRUE)) {
#   pngs <- Sys.glob(file.path(out_dir, "frame_*.png")) |> sort()
#   av::av_encode_video(pngs, output = file.path(out_dir, "flock.mp4"), framerate = 30)
# }
end_time <- Sys.time()
elapsed <- end_time - start_time
print(elapsed)
Back to Top