Random Walks Dashboard
Interactive dashboard with walker paths, step distributions, and fractal visualizations
#| '!! shinylive warning !!': |
#| shinylive does not work in self-contained HTML documents.
#| Please set `embed-resources: false` in your metadata.
#| standalone: true
#| viewerHeight: 1800
# ============================================================================
# FAST Package Loading (Pre-built webR binaries from GitHub Pages)
# ============================================================================
# Packages built by GitHub Actions and deployed to GitHub Pages
# See: https://github.com/JohnGavin/randomwalk/wiki/Shinylive-Lessons-Learned
cat("\n=== DASHBOARD SETUP ===\n")
cat("Installing packages from GitHub Pages...\n\n")
# Install randomwalk from GitHub Pages (built by deploy-pages.yaml workflow)
tryCatch({
webr::install(
"randomwalk",
repos = c(
"https://johngavin.github.io/randomwalk", # GitHub Pages webR repo
"https://repo.r-wasm.org" # Standard webR packages
)
)
cat("✅ randomwalk installed\n")
}, error = function(e) {
cat("❌ ERROR installing randomwalk:", e$message, "\n")
})
# Load shiny (usually pre-installed)
library(shiny)
cat("✅ shiny\n")
# Install ggplot2 dependencies and ggplot2
tryCatch({
library(ggplot2)
cat("✅ ggplot2\n")
}, error = function(e) {
cat("⚠️ ggplot2 not found, installing dependencies...\n")
# Install munsell first (required ggplot2 dependency often missing)
webr::install("munsell", repos = "https://repo.r-wasm.org")
cat("✅ munsell installed\n")
webr::install("ggplot2", repos = "https://repo.r-wasm.org")
library(ggplot2)
cat("✅ ggplot2 installed and loaded\n")
})
# Load randomwalk
tryCatch({
library(randomwalk)
cat("✅ randomwalk (v", as.character(packageVersion("randomwalk")), ")\n", sep = "")
}, error = function(e) {
cat("❌ ERROR loading randomwalk:", e$message, "\n")
cat("Available packages:\n")
print(installed.packages()[, c("Package", "Version")])
})
cat("\n⚠️ WebR Mode: Synchronous only (workers=0)\n")
cat(" Reason: mirai/nanonext not available in WebAssembly\n")
cat(" Note: Async works in native R with crew backend\n")
cat("\n=== SETUP COMPLETE ===\n\n")
# ============================================================================
# Helper Functions
# ============================================================================
# WebR environment detection (from Issue #129)
is_webr <- function() {
exists(".webr_env", envir = .GlobalEnv) ||
identical(Sys.getenv("WEBR"), "1")
}
# Format duration in seconds to mm:ss or hh:mm:ss format
format_duration <- function(seconds) {
if (is.null(seconds) || is.na(seconds)) return("00:00")
seconds <- round(seconds)
if (seconds < 60) {
return(sprintf("00:%02d", seconds))
} else if (seconds < 3600) {
mins <- floor(seconds / 60)
secs <- seconds %% 60
return(sprintf("%02d:%02d", mins, secs))
} else {
hours <- floor(seconds / 3600)
mins <- floor((seconds %% 3600) / 60)
secs <- seconds %% 60
return(sprintf("%02d:%02d:%02d", hours, mins, secs))
}
}
# ============================================================================
# Shiny UI
# ============================================================================
ui <- fluidPage(
# Add CSS for collapsible panels and progress display
tags$head(
tags$style(HTML("
.timer-display {
font-size: 24px;
font-weight: bold;
color: #007bff;
text-align: center;
padding: 10px;
background: #f0f8ff;
border-radius: 5px;
margin: 10px 0;
}
#run_sim:disabled {
opacity: 0.6;
cursor: not-allowed;
}
details {
margin-bottom: 15px;
border: 1px solid #ddd;
border-radius: 4px;
padding: 10px;
background: #f8f9fa;
}
details[open] {
background: white;
}
summary {
font-weight: bold;
cursor: pointer;
padding: 5px;
color: #333;
user-select: none;
}
summary:hover {
color: #007bff;
}
details > *:not(summary) {
margin-top: 10px;
}
.progress-popup {
position: fixed;
top: 50%;
left: 50%;
transform: translate(-50%, -50%);
background: rgba(0, 0, 0, 0.8);
color: white;
padding: 20px;
border-radius: 10px;
z-index: 9999;
font-size: 18px;
min-width: 300px;
text-align: center;
}
/* Dark theme for plots */
.dark-plot-container {
background-color: #2b2b2b;
padding: 10px;
border-radius: 5px;
}
"))
),
# TOP SECTION: COLLAPSIBLE PARAMETER GROUPS
fluidRow(
column(12,
# Use Shiny's wellPanel with conditional rendering for collapsible sections
wellPanel(
actionButton("toggle_sim_params", "▶ Simulation Parameters",
style = "background: none; border: none; font-weight: bold; font-size: 16px; padding: 0; margin-bottom: 10px;"),
conditionalPanel(
condition = "output.show_sim_params",
fluidRow(
column(4,
h5("Grid Settings"),
sliderInput("grid_size", "Grid Size:",
min = 20, max = 400, value = 100, step = 20),
helpText("Size of the square grid (pixels per side)")
),
column(8,
h5("Walker Settings"),
sliderInput("n_walkers", "Number of Walkers:",
min = 1, max = 1000, value = 200, step = 50),
helpText(HTML("Number of walkers to simulate<br><em>Max: 70% of grid pixels (updates automatically)</em>")),
sliderInput("max_steps", "Max Steps per Walker:",
min = 100, max = 10000, value = 5000, step = 100),
helpText("Maximum steps before walker terminates")
)
)
)
),
wellPanel(
actionButton("toggle_move_params", "▶ Movement Settings",
style = "background: none; border: none; font-weight: bold; font-size: 16px; padding: 0; margin-bottom: 10px;"),
conditionalPanel(
condition = "output.show_move_params",
fluidRow(
column(6,
selectInput("neighborhood", "Neighborhood:",
choices = c("4-hood", "8-hood"),
selected = "4-hood"),
helpText("4-hood: up/down/left/right, 8-hood: includes diagonals")
),
column(6,
selectInput("boundary", "Boundary Behavior:",
choices = c("terminate", "wrap"),
selected = "terminate"),
helpText("terminate: stop at edges, wrap: toroidal topology")
)
)
)
),
wellPanel(
actionButton("toggle_validation_params", "▶ Advanced Settings",
style = "background: none; border: none; font-weight: bold; font-size: 16px; padding: 0; margin-bottom: 10px;"),
conditionalPanel(
condition = "output.show_validation_params",
fluidRow(
column(4,
selectInput("sync_mode", "Sync Mode:",
choices = c("static", "chunked"),
selected = "static"),
helpText("static: Frozen snapshots (~5% collision). chunked: Batched updates (~15% collision)")
),
column(4,
checkboxInput("validate_strict", "Strict Validation",
value = FALSE),
helpText("ON: Stop on isolated pixel. OFF: Warn and continue")
),
column(4,
sliderInput("validate_percent", "Validation %:",
min = 0, max = 20, value = 5, step = 1),
helpText("Check every N% of walkers (0 = off)")
)
)
)
),
# Runtime Estimate and Run Button (always visible)
wellPanel(
style = "background-color: #f0f8ff; margin-top: 15px;",
fluidRow(
column(6,
h5("Estimated Runtime:"),
uiOutput("runtime_estimate")
),
column(6,
uiOutput("run_button_ui")
)
)
)
)
),
hr(),
# MAIN SECTION: FULL-WIDTH TABSET
fluidRow(
column(12,
tabsetPanel(
id = "tabs",
# PAGE 1: Fractal Graph (Black Pixels)
tabPanel("Fractal Graph",
h4("Grid Visualization - Black Pixels"),
plotOutput("fractal_plot", height = "600px"),
tags$p(style = "margin-top: 10px; margin-bottom: 10px; color: #666; font-style: italic; text-align: center;",
textOutput("fractal_caption", inline = TRUE)),
downloadButton("download_fractal", "Download Plot", class = "btn-sm", style = "margin-top: 10px;"),
hr(),
verbatimTextOutput("fractal_stats"),
hr(),
h5("Simulation Status:"),
verbatimTextOutput("status", placeholder = TRUE),
tags$div(
style = "background: #f0f0f0; padding: 10px; margin-top: 10px; font-family: monospace; font-size: 14px; border: 1px solid #ccc;",
h5("Live Progress:"),
verbatimTextOutput("progress_timer")
)),
# PAGE 2: Walker Paths (First 25 and Last 25)
tabPanel("Walker Paths",
h4("Walker Trajectories"),
tags$p(style = "color: #555; margin-bottom: 10px;",
"Shows walker trajectories from start (circle) to termination (square). ",
"Colors are stable - same walker keeps same color when adjusting slider."
),
fluidRow(
column(6,
sliderInput("first_n_paths", "First N to terminate:",
min = 0, max = 20, value = 5, step = 1)),
column(6,
sliderInput("last_n_paths", "Last N to terminate:",
min = 0, max = 20, value = 5, step = 1))
),
helpText("Select from walkers that FOUND THE FRACTAL (black_neighbor). ",
"These start from random interior positions. ",
"Boundary walkers (short paths) are excluded."),
plotOutput("paths_plot_first25", height = "500px"),
tags$p(style = "margin-top: 10px; color: #666; font-style: italic; text-align: center;",
textOutput("paths_caption", inline = TRUE))),
# PAGE 3: Distributions (Path Lengths by Termination Reason)
tabPanel("Distributions",
h4("Path Length Distributions"),
plotOutput("dist_overall", height = "300px"),
h5("Overall Distribution"),
hr(),
plotOutput("dist_by_reason", height = "300px"),
h5("By Termination Reason")),
# PAGE 4: Statistics - Quartile Analysis
tabPanel("Statistics",
h4("Quartile-Based Statistics"),
tags$div(
style = "background-color: #f5f5f5; padding: 10px; margin-bottom: 15px; border-radius: 5px;",
tags$p(strong("Terminology:")),
tags$ul(
tags$li(strong("Launch sequence (ID):"), " Order walkers were created (1, 2, ..., N)"),
tags$li(strong("Termination sequence:"), " Order walkers died (first to hit black = #1)"),
tags$li(strong("Path length (steps):"), " How many steps a walker took before terminating")
),
tags$p("See ", tags$a(href = "#", onclick = "Shiny.setInputValue('tabs', 'Notes');", "Notes tab"),
" for details on memory-efficient quartile calculation.")
),
helpText("Path length distributions partitioned by termination order quartiles (Q1=earliest terminators, Q4=latest)"),
hr(),
h5("Path Length by Termination Quartile"),
plotOutput("stats_quartile_boxplot", height = "350px"),
hr(),
h5("Histogram Panels by Quartile"),
plotOutput("stats_quartile_histograms", height = "400px"),
hr(),
h5("Distribution Shape Comparison"),
fluidRow(
column(6, plotOutput("stats_violin", height = "300px")),
column(6, plotOutput("stats_trend", height = "300px"))
),
hr(),
h5("Summary Statistics by Quartile"),
verbatimTextOutput("stats_quartile_summary")),
# PAGE 5: Survival Curve
tabPanel("Survival Curve",
h4("Walker Survival Analysis"),
# Interactive slider for step threshold
sliderInput("survival_step", "Show survival at step:",
min = 0, max = 1000, value = 500, step = 50),
helpText("Drag to see what proportion of walkers survived to this step"),
# Main survival curve plot
plotOutput("survival_curve", height = "400px"),
hr(),
# Competing risks breakdown
h5("Termination by Reason"),
plotOutput("termination_breakdown", height = "250px"),
hr(),
# Summary statistics
verbatimTextOutput("survival_stats"),
hr(),
# Hypothesis Testing Section
h4("Hypothesis Test: Survival Decay"),
tags$p("Tests whether later walkers terminate faster as black pixels accumulate."),
# Hypothesis plot
plotOutput("hypothesis_plot", height = "350px"),
# Hypothesis statistics
verbatimTextOutput("hypothesis_stats")),
# PAGE 6: Debug Info (renumbered from PAGE 5)
tabPanel("Debug",
h4("Debug Information"),
h5("Package Versions"),
verbatimTextOutput("debug_versions"),
hr(),
h5("Current Inputs"),
verbatimTextOutput("debug_inputs"),
hr(),
h5("Backend Information"),
verbatimTextOutput("debug_backend"),
hr(),
h5("Periodic Updates"),
verbatimTextOutput("debug_periodic"),
hr(),
h5("Simulation Status:"),
verbatimTextOutput("status_debug", placeholder = TRUE),
tags$div(
style = "background: #f0f0f0; padding: 10px; margin-top: 10px; font-family: monospace; font-size: 14px; border: 1px solid #ccc;",
h5("Live Progress:"),
verbatimTextOutput("progress_timer_debug")
)),
# PAGE 7: Notes
tabPanel("Notes",
h4("About This Dashboard"),
tags$p("Interactive random walk simulation via WebR/Shinylive."),
hr(),
h5("Sequence Terminology"),
tags$div(
style = "background-color: #f0f8ff; padding: 12px; border-radius: 5px; margin-bottom: 15px;",
tags$p(strong("Launch Sequence (Walker ID):"), " Order walkers are created (1, 2, ..., N)"),
tags$p(strong("Termination Sequence:"), " Order walkers stop moving (first to die = order 1)"),
tags$p(style = "font-style: italic; color: #666;",
"These often differ! Walker ID 10 might terminate last (order 1000).")
),
tags$table(style = "width: 100%; border-collapse: collapse; margin-bottom: 15px;",
tags$thead(
tags$tr(style = "background-color: #e8e8e8;",
tags$th(style = "padding: 8px; border: 1px solid #ddd;", "Sequence"),
tags$th(style = "padding: 8px; border: 1px solid #ddd;", "Pros"),
tags$th(style = "padding: 8px; border: 1px solid #ddd;", "Cons"),
tags$th(style = "padding: 8px; border: 1px solid #ddd;", "Best For")
)
),
tags$tbody(
tags$tr(
tags$td(style = "padding: 8px; border: 1px solid #ddd;", "Launch (ID)"),
tags$td(style = "padding: 8px; border: 1px solid #ddd;", "Known at creation; memory-efficient"),
tags$td(style = "padding: 8px; border: 1px solid #ddd;", "Arbitrary ordering"),
tags$td(style = "padding: 8px; border: 1px solid #ddd;", "Memory-constrained sims")
),
tags$tr(
tags$td(style = "padding: 8px; border: 1px solid #ddd;", "Termination"),
tags$td(style = "padding: 8px; border: 1px solid #ddd;", "Meaningful; reveals patterns"),
tags$td(style = "padding: 8px; border: 1px solid #ddd;", "Unknown until end"),
tags$td(style = "padding: 8px; border: 1px solid #ddd;", "Statistical analysis")
)
)
),
hr(),
h5("Memory-Efficient Quartile Calculation"),
tags$div(
style = "background-color: #f5f5f5; padding: 12px; border-radius: 5px; margin-bottom: 15px;",
tags$p(strong("Key Insight:"), " Quartiles do NOT require full paths!"),
tags$p("Each walker stores only: ", tags$code("id, steps, termination_reason, termination_order")),
tags$pre(style = "background-color: #fff; padding: 8px; border-radius: 3px; font-size: 12px;",
"df$quartile <- cut(df$termination_order,\n breaks = quantile(df$termination_order, probs = c(0, 0.25, 0.5, 0.75, 1)),\n labels = c(\"Q1\", \"Q2\", \"Q3\", \"Q4\"),\n include.lowest = TRUE)"
),
tags$p(style = "font-size: 12px; color: #666;",
"Uses O(4 bytes) per walker, not O(thousands) for paths.")
),
hr(),
h5("Path Storage Strategy"),
tags$p(strong("Current Approach (Launch ID):"), " Store paths for first/last 2000 walkers BY LAUNCH ID."),
tags$p("This is memory-efficient but means path visualization shows walkers by creation order, not termination order."),
tags$p(style = "color: #666; font-style: italic;",
"Future: Option E (chunked disk storage) would stream paths to parquet files, enabling retrieval by termination order. ",
"Not available in WebR due to arrow/parquet limitations in WASM."),
hr(),
h5("Execution Modes"),
tags$p("This browser dashboard runs in ", strong("sync sequential mode"), " (workers=0)."),
tags$p("For async parallel execution (workers > 0), use native R with crew/nanonext."),
hr(),
h5("Sync Mode Options (Advanced Settings)"),
tags$ul(
tags$li(strong("static:"), " Walkers use frozen grid snapshots. Lower collision rate (~5%). Faster but less realistic."),
tags$li(strong("chunked:"), " Grid updated in batches. Higher collision rate (~15%). More realistic pattern.")
),
tags$p(tags$em("Tip: Use 'chunked' for more interesting fractal patterns.")),
hr(),
h5("Validation Options (Advanced Settings)"),
tags$ul(
tags$li(strong("Strict Validation:"), " Stop immediately if isolated pixel detected."),
tags$li(strong("Validation %:"), " Check for isolated pixels every N% of walkers.")
),
hr(),
h5("Browser Limitations"),
tags$ul(
tags$li("Single-threaded execution (~400 steps/sec)"),
tags$li("First load: 15-30s (WebR package downloads)"),
tags$li("Recommended: grid ≤100, walkers ≤500, steps ≤1000"),
tags$li("For larger simulations, use native R with crew")
),
hr(),
tags$a(href="https://johngavin.github.io/randomwalk/", "Documentation"),
tags$span(" | "),
tags$a(href="https://github.com/JohnGavin/randomwalk", "GitHub"))
) # End of tabs
) # End of column
) # End of fluidRow
)
# ============================================================================
# Shiny Server
# ============================================================================
server <- function(input, output, session) {
# Dynamic walker constraint: Limit to 70% of grid pixels
observe({
max_walkers <- floor(0.7 * input$grid_size^2)
updateSliderInput(
session,
"n_walkers",
max = max_walkers,
value = min(input$n_walkers, max_walkers)
)
})
# Reactive values for simulation state and history
sim_result <- reactiveVal(NULL)
sim_count <- reactiveVal(0) # Count total simulations run
sim_state <- reactiveVal("idle") # States: idle, running, complete, error
sim_start_time <- reactiveVal(NULL)
sim_end_time <- reactiveVal(NULL)
sim_current_step <- reactiveVal(0) # Track current simulation step
sim_completed_walkers <- reactiveVal(0) # Track completed walkers
sim_black_pixels <- reactiveVal(0) # Track black pixels
# Reactive values for collapsible panels
show_sim_params <- reactiveVal(FALSE)
show_move_params <- reactiveVal(FALSE)
show_validation_params <- reactiveVal(FALSE)
# Toggle simulation parameters panel
observeEvent(input$toggle_sim_params, {
current <- show_sim_params()
show_sim_params(!current)
# Update button text
updateActionButton(session, "toggle_sim_params",
label = ifelse(!current, "▼ Simulation Parameters", "▶ Simulation Parameters"))
})
# Toggle movement settings panel
observeEvent(input$toggle_move_params, {
current <- show_move_params()
show_move_params(!current)
# Update button text
updateActionButton(session, "toggle_move_params",
label = ifelse(!current, "▼ Movement Settings", "▶ Movement Settings"))
})
# Toggle validation settings panel
observeEvent(input$toggle_validation_params, {
current <- show_validation_params()
show_validation_params(!current)
# Update button text
updateActionButton(session, "toggle_validation_params",
label = ifelse(!current, "▼ Validation Settings", "▶ Validation Settings"))
})
# Output for conditional panels
output$show_sim_params <- reactive({ show_sim_params() })
output$show_move_params <- reactive({ show_move_params() })
output$show_validation_params <- reactive({ show_validation_params() })
outputOptions(output, "show_sim_params", suspendWhenHidden = FALSE)
outputOptions(output, "show_move_params", suspendWhenHidden = FALSE)
outputOptions(output, "show_validation_params", suspendWhenHidden = FALSE)
# Runtime estimate (WebR performance varies: ~50-500 steps/sec depending on grid size & walker count)
output$runtime_estimate <- renderUI({
# Runtime calculation considers BOTH factors:
# 1. Simulation steps = max_steps (all walkers move simultaneously each step)
# 2. Per-step cost increases with more walkers due to collision checks & neighbor lookups
#
# Total cost = max_steps × (base_per_step_cost + walker_scaling × n_walkers)
# Rearranged: effective_steps_per_sec = base_steps_per_sec / (1 + walker_scaling × n_walkers)
grid_area <- input$grid_size^2
# Empirical base performance model (1 walker) based on grid size:
# Small grids (≤100x100): ~300 steps/sec
# Medium grids (100-200): ~150 steps/sec
# Large grids (>200): ~75 steps/sec
if (grid_area <= 10000) {
base_steps_per_sec <- 300
} else if (grid_area <= 40000) {
base_steps_per_sec <- 150
} else {
base_steps_per_sec <- 75
}
# Walker scaling factor: each walker adds ~0.2% overhead per step
# (collision checks, neighbor lookups, state updates scale with active walkers)
walker_scaling <- 0.002
# Effective throughput accounting for both factors
effective_steps_per_sec <- base_steps_per_sec / (1 + walker_scaling * input$n_walkers)
# Total simulation time
total_steps <- input$max_steps
est_seconds <- total_steps / effective_steps_per_sec
if (est_seconds < 60) {
est_text <- paste0("~", format_duration(est_seconds))
color <- "green"
} else if (est_seconds < 300) {
est_text <- paste0("~", format_duration(est_seconds))
color <- "orange"
} else {
est_text <- paste0("~", format_duration(est_seconds), " (consider smaller params)")
color <- "red"
}
# Calculate total walker steps for reference
total_walker_steps <- input$n_walkers * input$max_steps
tags$div(
style = paste0("color: ", color, "; font-weight: bold;"),
est_text,
tags$br(),
tags$small(style = "color: gray;",
sprintf("Max steps: %s | Walkers: %s\n(%sK total walker-steps, ~%.0f steps/sec effective)",
format(input$max_steps, big.mark=","),
format(input$n_walkers, big.mark=","),
format(round(total_walker_steps/1000), big.mark=","),
effective_steps_per_sec))
)
})
# Reactive timer for progress updates (following shiny-examples/142-reactive-timer pattern)
# This timer continues running even during synchronous simulation
autoInvalidate <- reactiveTimer(100) # Update every 100ms for smooth progress
# Observe timer to update estimated progress during simulation
observe({
autoInvalidate() # This triggers every 100ms
if (sim_state() == "running" && !is.null(sim_start_time())) {
elapsed <- as.numeric(difftime(Sys.time(), sim_start_time(), units = "secs"))
# Estimate progress based on typical runtime
est_steps_per_sec <- 200 / (1 + 0.002 * input$n_walkers) # Adjust for walker count
est_runtime <- input$max_steps / est_steps_per_sec
progress <- min(elapsed / est_runtime, 0.99)
# Update estimated counts for display
sim_current_step(floor(input$max_steps * progress))
sim_completed_walkers(floor(input$n_walkers * progress * 0.8))
sim_black_pixels(floor(input$n_walkers * progress * 0.7))
# Update status display with current progress
output$status <- renderText({
paste0(
"SIMULATION RUNNING...\n",
"Run #", sim_count(), "\n",
"━━━━━━━━━━━━━━━━━━━━━━\n",
"Started: ", format(sim_start_time(), "%H:%M:%S"), "\n",
"Elapsed: ", format_duration(elapsed), "\n",
"Progress: ", sprintf("%d%%", round(progress * 100)), "\n",
"━━━━━━━━━━━━━━━━━━━━━━\n",
"Grid: ", input$grid_size, "×", input$grid_size, "\n",
"Walkers: ", sim_completed_walkers(), "/", input$n_walkers, "\n",
"Black Pixels: ", sim_black_pixels(), "\n",
"Workers: 0 (sync sequential)"
)
})
}
})
# Real-time progress timer (following https://shinylive.io/r/examples/#timer pattern)
# This WILL work in Shinylive/WebR - reactive timer continues during computation
output$progress_timer <- renderText({
autoInvalidate() # Use the reactive timer that updates every 100ms
if (sim_state() == "running") {
# Get current time
current_time <- Sys.time()
elapsed <- if (!is.null(sim_start_time())) {
round(as.numeric(difftime(current_time, sim_start_time(), units = "secs")), 1)
} else {
0
}
# Format with current counts from reactive values
sprintf("⏱️ Time: %s | Walkers: %d/%d | Black Pixels: %d | Step: %d/%d",
format_duration(elapsed),
sim_completed_walkers(),
input$n_walkers,
sim_black_pixels(),
sim_current_step(),
input$max_steps)
} else if (sim_state() == "complete") {
result <- sim_result()
if (!is.null(result)) {
sprintf("✓ Complete | Black Pixels: %d | Grid: %d%%",
result$statistics$black_pixels,
round(result$statistics$black_percentage))
} else {
""
}
} else {
# Show current time to prove timer is working
sprintf("Ready | Current time: %s", format(Sys.time(), "%H:%M:%S"))
}
})
# Duplicate outputs for Debug tab (Shiny requires unique output IDs)
output$status_debug <- renderText({
# Mirror the status output
if (sim_state() == "running") {
paste0(
"SIMULATION RUNNING...\n",
"Run #", sim_count(), "\n",
"Started: ", format(sim_start_time(), "%H:%M:%S")
)
} else if (sim_state() == "complete" && !is.null(sim_result())) {
result <- sim_result()
paste0(
"SIMULATION COMPLETE\n",
"Run #", sim_count(), "\n",
"Black pixels: ", result$statistics$black_pixels, "\n",
"Elapsed: ", format_duration(result$statistics$elapsed_time_secs)
)
} else if (sim_state() == "error") {
"SIMULATION ERROR"
} else {
"Ready to run simulation..."
}
})
output$progress_timer_debug <- renderText({
autoInvalidate() # Use the reactive timer
if (sim_state() == "running") {
elapsed <- if (!is.null(sim_start_time())) {
round(as.numeric(difftime(Sys.time(), sim_start_time(), units = "secs")), 1)
} else { 0 }
sprintf("⏱️ Time: %s | Walkers: %d/%d | Black: %d",
format_duration(elapsed),
sim_completed_walkers(), input$n_walkers,
sim_black_pixels())
} else if (sim_state() == "complete" && !is.null(sim_result())) {
result <- sim_result()
sprintf("✓ Complete | Black: %d | Grid: %d%%",
result$statistics$black_pixels,
round(result$statistics$black_percentage))
} else {
sprintf("Ready | %s", format(Sys.time(), "%H:%M:%S"))
}
})
# Render button with walker/pixel counts
output$run_button_ui <- renderUI({
if (sim_state() == "running") {
# During simulation, show real-time walker/pixel counts
invalidateLater(200) # Update every 200ms
elapsed <- if (!is.null(sim_start_time())) {
as.numeric(difftime(Sys.time(), sim_start_time(), units = "secs"))
} else {
0
}
actionButton("run_sim_disabled",
sprintf("Running... W: %d/%d | B: %d | %s",
sim_completed_walkers(), input$n_walkers,
sim_black_pixels(), format_duration(elapsed)),
class = "btn-warning btn-lg",
style = "width: 100%; font-size: 14px;",
disabled = TRUE)
} else if (sim_state() == "complete" && !is.null(sim_result())) {
# Show walker count and black pixel count when complete
result <- sim_result()
total_secs <- as.numeric(difftime(sim_end_time(), sim_start_time(), units = "secs"))
mins <- floor(total_secs / 60)
secs <- floor(total_secs %% 60)
# Get walker and black pixel counts
walker_count <- result$statistics$completed_walkers
black_count <- result$statistics$black_pixels
actionButton("run_sim", sprintf("Run Simulation (Walkers: %d. Black: %d. Time: %02d:%02d)",
walker_count, black_count, mins, secs),
class = "btn-success btn-lg",
style = "width: 100%;")
} else if (sim_state() == "error") {
actionButton("run_sim", "Run Simulation (Error - Try Again)",
class = "btn-danger btn-lg",
style = "width: 100%;")
} else {
actionButton("run_sim", "Run Simulation",
class = "btn-primary btn-lg",
style = "width: 100%;")
}
})
# Run simulation when button clicked
observeEvent(input$run_sim, {
# Initialize state
sim_state("running")
sim_start_time(Sys.time())
sim_count(sim_count() + 1)
output$status <- renderText({
paste0(
"SIMULATION RUNNING...\n",
"Run #", sim_count(), "\n",
"Grid: ", input$grid_size, "×", input$grid_size, "\n",
"Walkers: ", input$n_walkers, "\n",
"Workers: 0 (sync sequential)\n",
"Started: ", format(sim_start_time(), "%H:%M:%S"), "\n",
"Status: Processing..."
)
})
# Validate walker count
max_allowed <- floor(0.7 * input$grid_size^2)
if (input$n_walkers > max_allowed) {
sim_state("error")
output$status <- renderText({
sprintf("ERROR: Too many walkers (%d) for grid size %dx%d.\nMaximum allowed: %d (70%% of grid pixels)",
input$n_walkers, input$grid_size, input$grid_size, max_allowed)
})
return()
}
# Reset progress tracking
sim_current_step(0)
sim_completed_walkers(0)
sim_black_pixels(0)
# Run the simulation
tryCatch({
result <- run_simulation(
grid_size = input$grid_size,
n_walkers = input$n_walkers,
neighborhood = input$neighborhood,
boundary = input$boundary,
max_steps = input$max_steps,
workers = 0,
sync_mode = input$sync_mode,
verbose = FALSE,
validate_percent = input$validate_percent,
validate_strict = input$validate_strict
)
# Store result
sim_result(result)
sim_completed_walkers(result$statistics$completed_walkers)
sim_black_pixels(result$statistics$black_pixels)
sim_current_step(input$max_steps)
# Mark complete
sim_state("complete")
sim_end_time(Sys.time())
# Calculate elapsed time
elapsed <- as.numeric(difftime(sim_end_time(), sim_start_time(), units = "secs"))
# Re-enable the button
updateActionButton(session, "run_sim",
label = "Run Simulation",
disabled = FALSE)
# Browser console logging for WebR diagnostics
message(sprintf("\n=== SIMULATION COMPLETE ==="))
message(sprintf("Grid size: %dx%d", input$grid_size, input$grid_size))
message(sprintf("Walkers: %d", input$n_walkers))
message(sprintf("Black pixels: %d", sum(result$grid == 1)))
message(sprintf("Elapsed time: %s", format_duration(elapsed)))
message(sprintf("==========================\n"))
output$status <- renderText({
paste0(
"SIMULATION COMPLETE\n",
"Run #", sim_count(), "\n",
"━━━━━━━━━━━━━━━━━━━━━━\n",
"Started: ", format(sim_start_time(), "%H:%M:%S"), "\n",
"Finished: ", format(sim_end_time(), "%H:%M:%S"), "\n",
"Elapsed: ", format_duration(elapsed), "\n",
"━━━━━━━━━━━━━━━━━━━━━━\n",
"Backend: synchronous (WebR browser)\n",
"Black Pixels: ", result$statistics$black_pixels, "\n",
"Walkers: ", result$statistics$completed_walkers, " completed"
)
})
}, error = function(e) {
sim_state("error")
sim_end_time(Sys.time())
# Re-enable the button
updateActionButton(session, "run_sim",
label = "Run Simulation",
disabled = FALSE)
output$status <- renderText({
paste0(
"SIMULATION ERROR\n",
"━━━━━━━━━━━━━━━━━━━━━━\n",
"Error: ", conditionMessage(e), "\n",
"━━━━━━━━━━━━━━━━━━━━━━\n",
"Please try again with different parameters."
)
})
})
})
# PAGE 1: Fractal Graph with arrival time coloring
output$fractal_plot <- renderPlot({
# If simulation is running, show progress in plot
if (sim_state() == "running") {
invalidateLater(300) # Update every 300ms
# Create a simple progress visualization
par(bg = "gray70", mar = c(5, 4, 4, 2))
plot(1, type = "n", xlim = c(0, 1), ylim = c(0, 1),
xlab = "", ylab = "", axes = FALSE,
main = "Simulation Running...")
# Show progress as text in plot
text(0.5, 0.7, sprintf("Walkers: %d / %d",
sim_completed_walkers(), input$n_walkers),
cex = 2, font = 2)
text(0.5, 0.5, sprintf("Black pixels: %d",
sim_black_pixels()),
cex = 2, font = 2)
elapsed <- if (!is.null(sim_start_time())) {
as.numeric(difftime(Sys.time(), sim_start_time(), units = "secs"))
} else { 0 }
text(0.5, 0.3, sprintf("Time: %s", format_duration(elapsed)),
cex = 1.5)
return()
}
req(sim_result())
result <- sim_result()
# Defensive checks for WebR environment (Issue #157)
req(!is.null(result$grid))
req(nrow(result$grid) > 0)
req(sum(result$grid == 1) > 0) # Ensure has black pixels
# Use enhanced plot with color-coded arrival times
p <- plot_grid_enhanced(result,
quantiles = 5,
color_scheme = "viridis")
# Store the plot for caption extraction
# The plot_grid_enhanced function returns plot with caption as attribute
p
})
# Display the caption separately
output$fractal_caption <- renderText({
req(sim_result())
result <- sim_result()
# Generate the plot to get the caption
p <- plot_grid_enhanced(result,
quantiles = 5,
color_scheme = "viridis")
# Extract and return the caption
caption <- attr(p, "caption")
if (!is.null(caption) && caption != "") {
caption
} else {
""
}
})
output$fractal_stats <- renderText({
req(sim_result())
stats <- sim_result()$statistics
paste0(
"Black Pixels: ", stats$black_pixels,
" (", sprintf("%d%%", round(stats$black_percentage)), " of grid)\n",
"Grid Size: ", stats$grid_size, "×", stats$grid_size,
" (", stats$grid_size^2, " total cells)"
)
})
# Download handler for fractal plot
output$download_fractal <- downloadHandler(
filename = function() {
paste0("fractal_dla_", format(Sys.time(), "%Y%m%d_%H%M%S"), ".png")
},
content = function(file) {
req(sim_result())
result <- sim_result()
p <- plot_grid_enhanced(result, quantiles = 5, color_scheme = "viridis")
ggplot2::ggsave(file, plot = p, width = 8, height = 8, dpi = 150)
}
)
# PAGE 2: Walker Paths
output$paths_plot_first25 <- renderPlot({
req(sim_result())
result <- sim_result()
n_walkers <- length(result$walkers)
if (n_walkers >= 1) {
first_n <- min(input$first_n_paths, n_walkers)
last_n <- min(input$last_n_paths, n_walkers)
# CRITICAL FIX: Only select from walkers that HAVE paths stored
# Paths are only stored for BLACK_NEIGHBOR terminators (walkers that found the fractal)
# Boundary walkers terminate quickly but don't have paths - skip them
# Filter to walkers with paths first
walkers_with_paths <- Filter(function(w) {
!is.null(w$path) && length(w$path) > 0
}, result$walkers)
n_with_paths <- length(walkers_with_paths)
if (n_with_paths == 0) {
# No paths stored - show message
plot.new()
text(0.5, 0.5, "No walker paths available.\nPaths are stored for walkers that found the fractal.",
cex = 1.2, col = "red")
return()
}
# Sort by termination order (among walkers with paths)
walkers_ordered <- walkers_with_paths[order(sapply(walkers_with_paths, function(w) {
if (!is.null(w$termination_order)) w$termination_order else w$id
}))]
# Select first N and last N from walkers WITH PATHS
first_n <- min(first_n, n_with_paths)
last_n <- min(last_n, n_with_paths)
walker_indices <- unique(c(
if(first_n > 0) 1:first_n else NULL,
if(last_n > 0) (n_with_paths - last_n + 1):n_with_paths else NULL
))
if (length(walker_indices) > 0) {
# Plot only selected walkers' paths
filtered_walkers <- walkers_ordered[walker_indices]
# Create custom plot with better visibility
grid_size <- result$parameters$grid_size
# Create plot with gray background for contrast
par(bg = "gray70")
plot(1, type = "n", xlim = c(1, grid_size), ylim = c(1, grid_size),
xlab = "", ylab = "",
main = paste(length(walker_indices), "of", n_with_paths, "Fractal-Finding Walkers (by termination order)"),
asp = 1, col.main = "black", col.lab = "black", col.axis = "black")
# Add subtle grid for visibility
grid(nx = grid_size, ny = grid_size, col = "gray70", lty = 1, lwd = 0.3)
# Create STABLE color mapping based on termination order (not position)
distinct_colors <- c("red", "blue", "darkgreen", "purple", "orange",
"brown", "deeppink", "darkcyan", "darkmagenta", "darkblue",
"darkred", "forestgreen", "black", "goldenrod", "darkviolet")
# Map each walker to a consistent color based on their termination order
# This ensures walker #20 always gets the same color regardless of selection
color_map <- list()
for (i in seq_along(walkers_ordered)) {
term_order <- walkers_ordered[[i]]$termination_order
# Use modulo to cycle through colors if more walkers than colors
color_index <- ((term_order - 1) %% length(distinct_colors)) + 1
color_map[[as.character(term_order)]] <- distinct_colors[color_index]
}
# Collect all paths with STABLE colors
all_paths <- list()
all_starts <- list()
all_ends <- list()
for (i in seq_along(filtered_walkers)) {
walker <- filtered_walkers[[i]]
# Get stable color for this walker based on termination order
walker_color <- color_map[[as.character(walker$termination_order)]]
# Check if path exists and is not empty
if (!is.null(walker$path) && length(walker$path) > 0) {
path_matrix <- do.call(rbind, walker$path)
# Use walker termination order as key to ensure uniqueness
key <- paste0("walker_", walker$termination_order)
all_paths[[key]] <- list(
x = path_matrix[, 2],
y = path_matrix[, 1],
col = walker_color
)
all_starts[[key]] <- list(
x = path_matrix[1, 2],
y = path_matrix[1, 1],
col = walker_color
)
}
# Always mark end position
if (!is.null(walker$pos)) {
key <- paste0("walker_", walker$termination_order)
all_ends[[key]] <- list(
x = walker$pos[2],
y = walker$pos[1],
col = walker_color
)
}
}
# Draw all paths
for (path_data in all_paths) {
if (!is.null(path_data)) {
points(path_data$x, path_data$y,
pch = 16, col = path_data$col, cex = 0.8)
}
}
# Draw all start positions
for (start_data in all_starts) {
if (!is.null(start_data)) {
points(start_data$x, start_data$y,
pch = 21, bg = start_data$col, col = "black", cex = 2, lwd = 2)
}
}
# Draw all end positions
for (end_data in all_ends) {
if (!is.null(end_data)) {
points(end_data$x, end_data$y,
pch = 22, bg = end_data$col, col = "black", cex = 2.5, lwd = 2)
}
}
# Caption moved to separate textOutput below plot
} else {
plot.new()
text(0.5, 0.5, "No walkers selected")
}
} else {
plot.new()
text(0.5, 0.5, "No walkers available")
}
})
# Caption for Walker Paths plot
output$paths_caption <- renderText({
"Symbols: Circle = Start | Square = End | Dots = Path. Colors by termination order."
})
# PAGE 3: Distributions
output$dist_overall <- renderPlot({
req(sim_result())
result <- sim_result()
# Set gray background
par(bg = "gray70")
# Defensive check for walkers
if (is.null(result$walkers) || length(result$walkers) == 0) {
plot.new()
text(0.5, 0.5, "Run a simulation first", cex = 1.5, col = "gray")
return()
}
# Extract path lengths safely
# NOTE: walker$path is a LIST of positions, not a matrix
# Use walker$steps for accurate step count
path_lengths <- tryCatch({
sapply(result$walkers, function(w) {
if (!is.null(w$steps)) {
w$steps
} else {
NA
}
})
}, error = function(e) {
numeric(0)
})
# Remove NAs and ensure numeric type
path_lengths <- as.numeric(path_lengths[!is.na(path_lengths)])
if (length(path_lengths) == 0 || !is.numeric(path_lengths)) {
plot.new()
text(0.5, 0.5, "No valid path data available", cex = 1.5, col = "gray")
return()
}
# Need at least 2 values for meaningful distribution
if (length(path_lengths) < 2) {
plot.new()
text(0.5, 0.5, paste0("Only ", length(path_lengths), " path - need at least 2"),
cex = 1.2, col = "gray")
return()
}
# Plot histogram with error handling
tryCatch({
hist(path_lengths,
breaks = min(30, max(5, length(unique(path_lengths)))),
col = "steelblue",
border = "white",
main = "Overall Path Length Distribution",
xlab = "Path Length (steps)",
ylab = "Frequency")
abline(v = median(path_lengths), col = "red", lwd = 2, lty = 2)
legend("topright",
legend = c(paste("Median:", median(path_lengths)),
paste("Mean:", sprintf("%.1f", mean(path_lengths)))),
col = c("red", "black"),
lty = c(2, 1),
lwd = 2)
}, error = function(e) {
plot.new()
text(0.5, 0.5, paste0("Error creating plot:\n", e$message), cex = 1, col = "red")
})
})
output$dist_by_reason <- renderPlot({
req(sim_result())
result <- sim_result()
# Defensive check for walkers
if (is.null(result$walkers) || length(result$walkers) == 0) {
plot.new()
text(0.5, 0.5, "Run a simulation first", cex = 1.5, col = "gray")
return()
}
# Extract termination reasons and path lengths safely
# NOTE: walker$path is a LIST, not a matrix - use walker$steps
path_lengths <- tryCatch({
sapply(result$walkers, function(w) {
if (!is.null(w$steps)) w$steps else NA
})
}, error = function(e) numeric(0))
reasons <- tryCatch({
sapply(result$walkers, function(w) {
if (!is.null(w$termination_reason)) w$termination_reason else "unknown"
})
}, error = function(e) character(0))
# Remove NAs and ensure same length
valid_idx <- !is.na(path_lengths)
path_lengths <- as.numeric(path_lengths[valid_idx])
reasons <- as.character(reasons[valid_idx])
# Defensive check: ensure vectors are same length
min_len <- min(length(path_lengths), length(reasons))
if (min_len < length(path_lengths) || min_len < length(reasons)) {
path_lengths <- path_lengths[1:min_len]
reasons <- reasons[1:min_len]
}
# Check if we have valid data
if (length(path_lengths) == 0 || length(reasons) == 0 ||
!is.numeric(path_lengths) || !is.character(reasons)) {
plot.new()
text(0.5, 0.5, "No valid termination data available", cex = 1.5, col = "gray")
return()
}
# Need at least 2 data points
if (length(path_lengths) < 2) {
plot.new()
text(0.5, 0.5, "Need at least 2 walkers for distribution", cex = 1.2, col = "gray")
return()
}
# Create data frame with error handling
walkers_data <- tryCatch({
data.frame(
path_length = path_lengths,
reason = reasons,
stringsAsFactors = FALSE
)
}, error = function(e) {
NULL
})
if (is.null(walkers_data)) {
plot.new()
text(0.5, 0.5, "Error creating data frame", cex = 1.2, col = "red")
return()
}
# Plot boxplot by termination reason with error handling
tryCatch({
print(
ggplot(walkers_data, aes(x = reason, y = path_length, fill = reason)) +
geom_boxplot(alpha = 0.7) +
geom_jitter(width = 0.2, alpha = 0.3) +
labs(title = "Path Lengths by Termination Reason",
x = "Termination Reason",
y = "Path Length (steps)") +
theme_minimal() +
theme(legend.position = "none",
panel.background = element_rect(fill = "gray70"),
plot.background = element_rect(fill = "gray70"))
)
}, error = function(e) {
plot.new()
text(0.5, 0.5, paste0("Error creating plot:\n", e$message), cex = 1, col = "red")
})
})
# PAGE 4: Statistics - Quartile-based Analysis
# Helper function to prepare quartile data
get_quartile_data <- reactive({
req(sim_result())
result <- sim_result()
# Extract steps and termination order
steps <- sapply(result$walkers, function(w) w$steps)
term_order <- sapply(result$walkers, function(w) w$termination_order)
# Create data frame
df <- data.frame(
steps = steps,
termination_order = term_order,
stringsAsFactors = FALSE
)
# Assign quartiles based on termination order
df$quartile <- cut(df$termination_order,
breaks = quantile(df$termination_order, probs = c(0, 0.25, 0.5, 0.75, 1)),
labels = c("Q1 (First 25%)", "Q2 (25-50%)", "Q3 (50-75%)", "Q4 (Last 25%)"),
include.lowest = TRUE)
df
})
# Boxplot by quartile
output$stats_quartile_boxplot <- renderPlot({
req(sim_result())
df <- get_quartile_data()
if (is.null(df) || nrow(df) < 4) {
plot.new()
text(0.5, 0.5, "Need at least 4 walkers for quartile analysis", cex = 1.2, col = "gray")
return()
}
tryCatch({
print(
ggplot(df, aes(x = quartile, y = steps)) +
geom_boxplot(fill = "steelblue", alpha = 0.7, outlier.shape = 21) +
geom_jitter(width = 0.15, alpha = 0.3, size = 1) +
labs(
title = "Path Length Distribution by Termination Quartile",
subtitle = "Q1 = earliest terminators, Q4 = latest terminators",
x = "Termination Quartile",
y = "Path Length (steps)"
) +
theme_minimal(base_size = 12) +
theme(legend.position = "none",
axis.text.x = element_text(angle = 15, hjust = 1),
panel.background = element_rect(fill = "gray70"),
plot.background = element_rect(fill = "gray70"))
)
}, error = function(e) {
plot.new()
text(0.5, 0.5, paste0("Error: ", e$message), cex = 1, col = "red")
})
})
# Histogram panels by quartile
output$stats_quartile_histograms <- renderPlot({
req(sim_result())
df <- get_quartile_data()
if (is.null(df) || nrow(df) < 4) {
plot.new()
text(0.5, 0.5, "Need at least 4 walkers for quartile analysis", cex = 1.2, col = "gray")
return()
}
tryCatch({
print(
ggplot(df, aes(x = steps)) +
geom_histogram(bins = 20, fill = "steelblue", alpha = 0.8, color = "white") +
facet_wrap(~ quartile, scales = "free_y", ncol = 2) +
labs(
title = "Path Length Histograms by Termination Quartile",
x = "Path Length (steps)",
y = "Count"
) +
theme_minimal(base_size = 11) +
theme(legend.position = "none",
strip.text = element_text(face = "bold"),
panel.background = element_rect(fill = "gray70"),
plot.background = element_rect(fill = "gray70"))
)
}, error = function(e) {
plot.new()
text(0.5, 0.5, paste0("Error: ", e$message), cex = 1, col = "red")
})
})
# Violin plot for distribution shape comparison
output$stats_violin <- renderPlot({
req(sim_result())
df <- get_quartile_data()
if (is.null(df) || nrow(df) < 4) {
plot.new()
text(0.5, 0.5, "Need at least 4 walkers", cex = 1.2, col = "gray")
return()
}
tryCatch({
print(
ggplot(df, aes(x = quartile, y = steps)) +
geom_violin(fill = "steelblue", alpha = 0.6, trim = FALSE) +
geom_boxplot(width = 0.15, fill = "white", alpha = 0.8) +
labs(
title = "Distribution Shape by Quartile",
x = "Termination Quartile",
y = "Path Length (steps)"
) +
theme_minimal(base_size = 11) +
theme(axis.text.x = element_text(angle = 15, hjust = 1),
panel.background = element_rect(fill = "gray70"),
plot.background = element_rect(fill = "gray70"))
)
}, error = function(e) {
plot.new()
text(0.5, 0.5, paste0("Error: ", e$message), cex = 1, col = "red")
})
})
# Trend plot showing mean/median across quartiles
output$stats_trend <- renderPlot({
req(sim_result())
df <- get_quartile_data()
if (is.null(df) || nrow(df) < 4) {
plot.new()
text(0.5, 0.5, "Need at least 4 walkers", cex = 1.2, col = "gray")
return()
}
tryCatch({
# Calculate summary statistics by quartile
summary_df <- do.call(rbind, lapply(levels(df$quartile), function(q) {
q_data <- df[df$quartile == q, ]
data.frame(
quartile = q,
mean_steps = mean(q_data$steps),
median_steps = median(q_data$steps),
q_num = match(q, levels(df$quartile))
)
}))
print(
ggplot(summary_df, aes(x = q_num)) +
geom_line(aes(y = mean_steps, color = "Mean"), linewidth = 1.2) +
geom_point(aes(y = mean_steps, color = "Mean"), size = 3) +
geom_line(aes(y = median_steps, color = "Median"), linewidth = 1.2, linetype = "dashed") +
geom_point(aes(y = median_steps, color = "Median"), size = 3) +
scale_x_continuous(breaks = 1:4, labels = c("Q1", "Q2", "Q3", "Q4")) +
scale_color_manual(values = c("Mean" = "steelblue", "Median" = "darkred")) +
labs(
title = "Path Length Trend Across Quartiles",
x = "Termination Quartile",
y = "Path Length (steps)",
color = "Statistic"
) +
theme_minimal(base_size = 11) +
theme(legend.position = "bottom",
panel.background = element_rect(fill = "gray70"),
plot.background = element_rect(fill = "gray70"))
)
}, error = function(e) {
plot.new()
text(0.5, 0.5, paste0("Error: ", e$message), cex = 1, col = "red")
})
})
# Summary statistics by quartile
output$stats_quartile_summary <- renderPrint({
req(sim_result())
df <- get_quartile_data()
if (is.null(df) || nrow(df) < 4) {
cat("Need at least 4 walkers for quartile analysis\n")
return()
}
cat("=== QUARTILE SUMMARY STATISTICS ===\n\n")
for (q in levels(df$quartile)) {
q_data <- df[df$quartile == q, ]
cat(sprintf("%s (n=%d):\n", q, nrow(q_data)))
cat(sprintf(" Min: %d steps\n", min(q_data$steps)))
cat(sprintf(" Q1: %d steps\n", as.integer(quantile(q_data$steps, 0.25))))
cat(sprintf(" Median: %d steps\n", as.integer(median(q_data$steps))))
cat(sprintf(" Q3: %d steps\n", as.integer(quantile(q_data$steps, 0.75))))
cat(sprintf(" Max: %d steps\n", max(q_data$steps)))
cat(sprintf(" Mean: %.1f steps\n", mean(q_data$steps)))
cat(sprintf(" SD: %.1f steps\n", sd(q_data$steps)))
cat("\n")
}
# Overall comparison
cat("=== QUARTILE COMPARISON ===\n\n")
means <- tapply(df$steps, df$quartile, mean)
cat("Mean steps by quartile:\n")
for (i in seq_along(means)) {
cat(sprintf(" %s: %.1f\n", names(means)[i], means[i]))
}
# Trend indicator
if (means[4] < means[1]) {
cat("\n⬇️ DECREASING TREND: Later walkers survive fewer steps\n")
cat(sprintf(" (Q4 mean %.0f%% of Q1 mean)\n", 100 * means[4] / means[1]))
} else if (means[4] > means[1]) {
cat("\n⬆️ INCREASING TREND: Later walkers survive more steps\n")
cat(sprintf(" (Q4 mean %.0f%% of Q1 mean)\n", 100 * means[4] / means[1]))
} else {
cat("\n➡️ FLAT: No clear trend across quartiles\n")
}
})
# ============================================================================
# PAGE 5: Survival Curve - Server Logic
# ============================================================================
# Dynamic slider max based on actual max steps in data
observe({
req(sim_result())
result <- sim_result()
max_step <- max(sapply(result$walkers, function(w) w$steps))
updateSliderInput(session, "survival_step",
max = max_step,
value = min(input$survival_step, max_step))
})
# Survival curve plot (Kaplan-Meier style)
output$survival_curve <- renderPlot({
req(sim_result())
result <- sim_result()
# Build survival data
steps <- sapply(result$walkers, function(w) w$steps)
n_total <- length(steps)
# Calculate survival proportion at each unique step value
unique_steps <- sort(unique(steps))
survival_at_step <- sapply(unique_steps, function(s) {
sum(steps > s) / n_total
})
# Current slider position
current_step <- input$survival_step
current_survival <- sum(steps > current_step) / n_total
# Plot using ggplot2
df <- data.frame(step = unique_steps, survival = survival_at_step)
ggplot(df, aes(x = step, y = survival)) +
geom_step(color = "steelblue", linewidth = 1.2) +
geom_vline(xintercept = current_step, linetype = "dashed", color = "red") +
geom_point(aes(x = current_step, y = current_survival),
color = "red", size = 3) +
annotate("text", x = current_step, y = min(current_survival + 0.08, 1),
label = sprintf("%.1f%% alive", current_survival * 100),
color = "red", hjust = 0, size = 4) +
labs(
title = "Walker Survival Curve",
subtitle = sprintf("Red line: step %d (%.1f%% still walking)",
current_step, current_survival * 100),
x = "Steps",
y = "Proportion Surviving"
) +
scale_y_continuous(labels = scales::percent_format(), limits = c(0, 1)) +
theme_minimal(base_size = 14) +
theme(
plot.title = element_text(hjust = 0.5, face = "bold"),
plot.subtitle = element_text(hjust = 0.5, color = "gray40"),
panel.background = element_rect(fill = "gray70"),
plot.background = element_rect(fill = "gray70")
)
})
# Termination breakdown by reason
output$termination_breakdown <- renderPlot({
req(sim_result())
result <- sim_result()
reasons <- sapply(result$walkers, function(w) w$termination_reason)
reason_counts <- table(reasons)
df <- data.frame(
reason = names(reason_counts),
count = as.integer(reason_counts)
)
df$pct <- df$count / sum(df$count) * 100
# Color mapping for different termination reasons
color_map <- c(
"black_neighbor" = "steelblue",
"hit_boundary" = "orange",
"max_steps" = "gray50",
"touched_black" = "darkblue"
)
# Create fill colors vector matching data
df$fill_color <- ifelse(df$reason %in% names(color_map),
color_map[df$reason],
"gray70")
ggplot(df, aes(x = reorder(reason, -count), y = count, fill = reason)) +
geom_col() +
geom_text(aes(label = sprintf("%.1f%%", pct)), vjust = -0.5, size = 3.5) +
scale_fill_manual(values = color_map, na.value = "gray70") +
labs(
title = "How Walkers Terminated",
x = "Termination Reason",
y = "Count"
) +
theme_minimal(base_size = 12) +
theme(
legend.position = "none",
axis.text.x = element_text(angle = 30, hjust = 1),
plot.title = element_text(hjust = 0.5, face = "bold"),
panel.background = element_rect(fill = "gray70"),
plot.background = element_rect(fill = "gray70")
)
})
# Survival statistics text
output$survival_stats <- renderPrint({
req(sim_result())
result <- sim_result()
steps <- sapply(result$walkers, function(w) w$steps)
cat("=== SURVIVAL STATISTICS ===\n\n")
cat(sprintf("Total walkers: %d\n", length(steps)))
cat(sprintf("Median survival: %d steps\n", median(steps)))
cat(sprintf("Mean survival: %.1f steps\n", mean(steps)))
cat(sprintf("25th percentile: %d steps\n", as.integer(quantile(steps, 0.25))))
cat(sprintf("75th percentile: %d steps\n", as.integer(quantile(steps, 0.75))))
cat("\n")
cat("At current slider position:\n")
current_survival <- sum(steps > input$survival_step)
cat(sprintf(" %d walkers (%.1f%%) survived past step %d\n",
current_survival,
current_survival / length(steps) * 100,
input$survival_step))
})
# Hypothesis test: Geometric vs Exponential decay models
output$hypothesis_plot <- renderPlot({
req(sim_result())
result <- sim_result()
# Build survival data with termination order
survdat <- data.frame(
order = sapply(result$walkers, function(w) {
if (!is.null(w$termination_order)) w$termination_order else NA
}),
steps = sapply(result$walkers, function(w) w$steps)
)
# Remove NA orders (if termination_order not available)
survdat <- survdat[!is.na(survdat$order), ]
if (nrow(survdat) < 10) {
# Not enough data for hypothesis testing
plot.new()
text(0.5, 0.5, "Need at least 10 walkers with termination order for hypothesis test",
cex = 1.2, col = "gray")
return()
}
survdat <- survdat[order(survdat$order), ]
# H0: Geometric decay model
# Geometric: P(survive n steps) = (1-p)^n, E[steps] = (1-p)/p
# Fit p via method of moments: p_hat = 1 / (1 + mean(steps))
geom_p <- 1 / (1 + mean(survdat$steps))
# Predicted mean at each order position (linear decay in expected value)
survdat$geom_pred <- (1 / geom_p) * (1 - survdat$order / max(survdat$order) * 0.5)
# H1: Exponential decay model
# E[steps | order] = a * exp(-lambda * order)
survdat$log_steps <- log(survdat$steps + 1)
exp_model <- lm(log_steps ~ order, data = survdat)
survdat$exp_pred <- exp(fitted(exp_model))
# Non-parametric: Kolmogorov-Smirnov tests
geom_ks <- tryCatch({
ks.test(survdat$steps, function(x) pgeom(x, prob = geom_p))
}, error = function(e) list(p.value = 0, statistic = 1))
exp_rate <- 1 / mean(survdat$steps)
exp_ks <- tryCatch({
ks.test(survdat$steps, function(x) pexp(x, rate = exp_rate))
}, error = function(e) list(p.value = 0, statistic = 1))
# Determine which model fits better (higher p-value = better fit)
geom_better <- geom_ks$p.value > exp_ks$p.value
# Plot
ggplot(survdat, aes(x = order, y = steps)) +
geom_point(alpha = 0.4, size = 1.5, color = "gray40") +
geom_line(aes(y = geom_pred), color = "blue", linewidth = 1.2, linetype = "dashed") +
geom_line(aes(y = exp_pred), color = "red", linewidth = 1.2) +
annotate("text", x = max(survdat$order) * 0.6, y = max(survdat$steps) * 0.9,
label = sprintf("H0 Geometric: KS p=%.3f\nH1 Exponential: KS p=%.3f\nBetter fit: %s",
geom_ks$p.value, exp_ks$p.value,
ifelse(geom_better, "Geometric", "Exponential")),
hjust = 0, size = 3.5,
color = ifelse(geom_better, "blue", "red")) +
labs(
title = "Survival Decay: Geometric vs Exponential",
subtitle = sprintf("Best fit: %s model (higher KS p-value = better fit)",
ifelse(geom_better, "Geometric (H0)", "Exponential (H1)")),
x = "Walker Termination Order",
y = "Steps Survived",
caption = "Blue dashed: Geometric (H0) | Red solid: Exponential (H1)"
) +
theme_minimal(base_size = 12) +
theme(
panel.background = element_rect(fill = "gray70"),
plot.background = element_rect(fill = "gray70"),
plot.title = element_text(hjust = 0.5, face = "bold"),
plot.subtitle = element_text(hjust = 0.5,
color = ifelse(geom_better, "blue", "red"))
)
})
# Hypothesis test summary text
output$hypothesis_stats <- renderPrint({
req(sim_result())
result <- sim_result()
# Build survival data with termination order
survdat <- data.frame(
order = sapply(result$walkers, function(w) {
if (!is.null(w$termination_order)) w$termination_order else NA
}),
steps = sapply(result$walkers, function(w) w$steps)
)
survdat <- survdat[!is.na(survdat$order), ]
if (nrow(survdat) < 10) {
cat("Insufficient data for hypothesis testing.\n")
cat("Need at least 10 walkers with termination_order field.\n")
return()
}
cat("=== HYPOTHESIS TEST: GEOMETRIC vs EXPONENTIAL DECAY ===\n\n")
cat("H0 (Geometric): Step counts follow a geometric distribution\n")
cat("H1 (Exponential): Step counts follow an exponential distribution\n\n")
# Geometric model
geom_p <- 1 / (1 + mean(survdat$steps))
cat("--- Geometric Model (H0) ---\n")
cat(sprintf("Estimated p: %.4f\n", geom_p))
cat(sprintf("Expected mean: %.1f steps\n", (1 - geom_p) / geom_p))
# Exponential model
exp_rate <- 1 / mean(survdat$steps)
survdat$log_steps <- log(survdat$steps + 1)
exp_model <- lm(log_steps ~ order, data = survdat)
cat("\n--- Exponential Model (H1) ---\n")
cat(sprintf("Rate parameter (lambda): %.6f\n", exp_rate))
cat(sprintf("Expected mean: %.1f steps\n", 1 / exp_rate))
cat(sprintf("Decay coefficient: %.6f\n", coef(exp_model)[2]))
# KS tests
geom_ks <- tryCatch({
ks.test(survdat$steps, function(x) pgeom(x, prob = geom_p))
}, error = function(e) list(p.value = 0, statistic = 1))
exp_ks <- tryCatch({
ks.test(survdat$steps, function(x) pexp(x, rate = exp_rate))
}, error = function(e) list(p.value = 0, statistic = 1))
cat("\n--- Kolmogorov-Smirnov Tests ---\n")
cat(sprintf("Geometric KS statistic: %.4f, p-value: %.4f\n",
geom_ks$statistic, geom_ks$p.value))
cat(sprintf("Exponential KS statistic: %.4f, p-value: %.4f\n",
exp_ks$statistic, exp_ks$p.value))
# Significance threshold
alpha <- 0.05
cat("\n=== CONCLUSION ===\n")
# Check if BOTH models are rejected (p < alpha)
geom_rejected <- geom_ks$p.value < alpha
exp_rejected <- exp_ks$p.value < alpha
if (geom_rejected && exp_rejected) {
cat("NEITHER model fits well (both p-values < 0.05)\n")
cat("The step distribution may be more complex than geometric or exponential.\n")
cat("Consider: mixture models, heavy-tailed distributions, or spatial effects.\n")
} else if (!geom_rejected && !exp_rejected) {
# Both acceptable - compare which is better
if (geom_ks$p.value > exp_ks$p.value) {
cat("GEOMETRIC model fits better (p=", round(geom_ks$p.value, 4), ")\n")
cat("Step counts suggest discrete process with constant termination probability.\n")
} else {
cat("EXPONENTIAL model fits better (p=", round(exp_ks$p.value, 4), ")\n")
cat("Step counts suggest continuous decay process.\n")
}
} else if (!geom_rejected) {
cat("GEOMETRIC model fits (p=", round(geom_ks$p.value, 4), ")\n")
cat("Exponential model rejected (p=", round(exp_ks$p.value, 4), ")\n")
} else {
cat("EXPONENTIAL model fits (p=", round(exp_ks$p.value, 4), ")\n")
cat("Geometric model rejected (p=", round(geom_ks$p.value, 4), ")\n")
}
})
# PAGE 6: Debug - Versions
output$debug_versions <- renderText({
paste0(
"=== PACKAGE VERSIONS ===\n",
"randomwalk: ", packageVersion("randomwalk"), "\n",
"shiny: ", packageVersion("shiny"), "\n",
"ggplot2: ", packageVersion("ggplot2"), "\n",
if (requireNamespace("mirai", quietly = TRUE)) {
paste0("mirai: ", packageVersion("mirai"), "\n")
} else "",
if (requireNamespace("nanonext", quietly = TRUE)) {
paste0("nanonext: ", packageVersion("nanonext"), "\n")
} else "",
if (requireNamespace("crew", quietly = TRUE)) {
paste0("crew: ", packageVersion("crew"), "\n")
} else ""
)
})
# PAGE 5: Debug - Inputs
output$debug_inputs <- renderText({
paste0(
"=== CURRENT INPUTS ===\n",
"Grid Size: ", input$grid_size, "\n",
"Number of Walkers: ", input$n_walkers, "\n",
"Mode: sync sequential (workers=0)\n",
"Neighborhood: ", input$neighborhood, "\n",
"Boundary: ", input$boundary, "\n",
"Max Steps: ", input$max_steps
)
})
# PAGE 5: Debug - Backend
output$debug_backend <- renderText({
env_is_webr <- is_webr()
paste0(
"=== BACKEND INFORMATION ===\n",
"Environment: WebR/Browser\n",
"R Version: ", R.version.string, "\n",
"Backend: synchronous (no parallelization)\n",
"Sync Mode: static (grid snapshots)\n",
"\n=== ENVIRONMENT DETECTION ===\n",
".webr_env exists: ", exists(".webr_env", envir = .GlobalEnv), "\n",
"WEBR env var: ", Sys.getenv("WEBR"), "\n",
"is_webr() result: ", env_is_webr
)
})
# PAGE 5: Debug - Periodic Updates
output$debug_periodic <- renderText({
# Depend on timer to update every 100ms
autoInvalidate()
# Enhanced debug info during execution
state_info <- switch(sim_state(),
"idle" = "💤 Idle - waiting for user to click 'Run Simulation'",
"running" = {
elapsed <- if (!is.null(sim_start_time())) {
as.numeric(difftime(Sys.time(), sim_start_time(), units = "secs"))
} else {
0
}
paste0("Running - ", format_duration(elapsed), " elapsed")
},
"complete" = {
elapsed <- if (!is.null(sim_end_time()) && !is.null(sim_start_time())) {
as.numeric(difftime(sim_end_time(), sim_start_time(), units = "secs"))
} else {
0
}
paste0("Complete - took ", format_duration(elapsed))
},
"error" = "❌ Error - see status tab for details",
"Unknown state"
)
paste0(
"=== PERIODIC UPDATES ===\n",
"Current Time: ", format(Sys.time(), "%Y-%m-%d %H:%M:%S"), "\n",
"Simulation State: ", state_info, "\n",
"Simulation Count: ", sim_count(), "\n",
"Result Available: ", if (is.null(sim_result())) "No" else "Yes", "\n",
"Active Tab: ", input$tabs, "\n",
if (sim_state() == "running" && !is.null(sim_start_time())) {
paste0(
"\n=== EXECUTION INFO ===\n",
"Grid Size: ", input$grid_size, "×", input$grid_size, "\n",
"Walkers: ", input$n_walkers, "\n",
"Workers: 0 (sync sequential)\n",
"Backend: synchronous"
)
} else {
""
}
)
})
}
# ============================================================================
# Run Shiny App
# ============================================================================
shinyApp(ui = ui, server = server)Tips
- Larger grid (200-400): More exploration space
- Fewer walkers (100-500): Faster simulation
- 8-hood: Diagonal movement for different patterns
- Wrap boundary: Toroidal grid (no edge termination)
Build Info
randomwalk 2.1.1 | Git a38b203 | R 4.5.2 | Built 2026-02-18 17:34:52