---
title: "bGMYC4 Interactive Workflow"
author: "Dmitry Karabanov & Qwen"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{bGMYC4 Interactive Workflow}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
# EN: Disable execution for CRAN checks to prevent timeouts on interactive MCMC.
# RU: Отключено выполнение для проверки CRAN, чтобы избежать таймаутов.
# Для локального запуска: измените eval = FALSE на eval = TRUE или запускайте блоки вручную.
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7, fig.height = 5,
  warning = FALSE, message = FALSE,
  eval = FALSE
)
```

## 🌍 Introduction / Введение    
**EN:** This vignette documents the complete interactive workflow for Bayesian species delimitation using bGMYC4. It guides you through data selection, tree preprocessing (BEAST2 annotations), interactive MCMC diagnostics with parameter auto-tuning, multi-tree uncertainty pooling via parallel computing, and custom interactive visualization.
**RU:** Эта виньетка описывает полный интерактивный рабочий процесс байесовской делимитации видов с помощью bGMYC4. Она проведёт вас через выбор данных, предобработку деревьев (аннотации BEAST2), интерактивную диагностику MCMC с автонастройкой параметров, усреднение филогенетической неопределённости через параллельные вычисления и кастомную интерактивную визуализацию.    
---

## 1. Setup & Helper Functions / Настройка окружения и вспомогательные функции    
**EN:** Load required libraries and define safe input handlers. The workflow automatically cleans BEAST2 annotations, checks ultrametricity, and validates parameter bounds.    
**RU:** Загрузка библиотек и определение безопасных обработчиков ввода. Рабочий процесс автоматически очищает аннотации BEAST2, проверяет ультраметричность и валидирует границы параметров.    

```{r setup}
compiler::enableJIT(3) # JIT acceleration / Ускорение JIT
library(ape)
library(bGMYC4)
library(treeio)
library(ggtree)
library(dplyr)
library(plotly)
library(htmlwidgets)
library(future)
library(future.apply)
library(mcmcse)

# EN: Safe numeric scalar input / Безопасный ввод скалярных значений
input_scalar <- function(label, default_val) {
  if (!interactive()) return(default_val)
  val <- readline(sprintf("  %s (default / по умолчанию: %s): ", label, default_val))
  if (nchar(trimws(val)) == 0) return(default_val)
  num <- suppressWarnings(as.numeric(val))
  if (is.na(num)) return(default_val)
  return(num)
}

# EN: Safe numeric vector input / Безопасный ввод векторов
input_vector <- function(label, default_vec) {
  if (!interactive()) return(default_vec)
  val <- readline(sprintf("  %s (comma-separated, default / через запятую, по умолчанию: %s): ", 
                          label, paste(default_vec, collapse = ", ")))
  if (nchar(trimws(val)) == 0) return(default_vec)
  nums <- suppressWarnings(as.numeric(unlist(strsplit(val, "[,;\\s]+"))))
  nums <- nums[!is.na(nums)]
  if (length(nums) != length(default_vec)) return(default_vec)
  return(nums)
}

# EN: Enforce ultrametricity / Обеспечение ультраметричности
fix_ultrametric <- function(tr) {
  if (!is.ultrametric(tr)) {
    tr$edge.length <- round(tr$edge.length, 8)
    if (!is.ultrametric(tr)) stop("❌ Tree remains non-ultrametric / Дерево остаётся неультраметричным.")
  }
  return(tr)
}
```

---

## 2. Data Selection & Tree Preprocessing / Выбор данных и предобработка деревьев    
**EN:** Load consensus tree with BEAST2 annotations (for posterior visualization).  
**RU:** Загрузка консенсусного дерева с аннотациями BEAST2 (для визуализации posterior).    

```{r data}
consensus_path <- if (interactive()) readline("📥 Path to consensus tree / Путь к консенсусному дереву (.tree): ") else "dummy.tree"
posterior_path <- if (interactive()) readline("📥 Path to posterior trees / Путь к набору деревьев (.trees): ") else "dummy.trees"

tree_beast <- treeio::read.beast(consensus_path)
tree_consensus <- tree_beast@phylo
all_trees <- read.nexus(posterior_path)
class(all_trees) <- "multiPhylo"

# EN: Analysis Mode Selection / Выбор режима анализа
mode_str <- if (interactive()) readline("Mode 1 (Single) or 2 (Multi)? / Режим 1 или 2? [2]: ") else "2"
analysis_mode <- ifelse(nchar(trimws(mode_str)) == 0, 2, as.integer(mode_str))

# EN: Tree sampling with burnin / Выборка деревьев с учетом burnin
if (analysis_mode == 2) {
  n_total <- length(all_trees)
  burnin_idx <- floor(n_total * 0.10) # 10% burnin
  n_sample <- input_scalar("Number of trees to sample / Кол-во деревьев", 10)
  set.seed(42)
  trees_sample <- all_trees[sample((burnin_idx + 1):n_total, n_sample)]
  class(trees_sample) <- "multiPhylo"
} else {
  trees_sample <- NULL
}

# EN: Outgroup removal via regex patterns / Удаление аутгрупп по паттернам
if (interactive() && tolower(readline("Drop outgroups? / Удалять аутгруппы? (y/n): ")) == "y") {
  og_str <- readline("Enter patterns (comma-separated) / Введите паттерны: ")
  og_patterns <- trimws(unlist(strsplit(og_str, "[,;]+")))
  matching_tips <- unique(unlist(lapply(og_patterns, function(p) {
    pattern <- paste0("(^|[^A-Za-z0-9])", p, "($|[^A-Za-z0-9])")
    tree_consensus$tip.label[grepl(pattern, tree_consensus$tip.label, ignore.case = TRUE, perl = TRUE)]
  })))
  
  if (length(matching_tips) > 0) {
    tree_consensus <- drop.tip(tree_consensus, matching_tips)
    if (analysis_mode == 2) {
      trees_sample <- lapply(trees_sample, drop.tip, tip = matching_tips)
      class(trees_sample) <- "multiPhylo"
    }
    tree_beast <- treeio::drop.tip(tree_beast, matching_tips)
  }
}

tree_consensus <- fix_ultrametric(tree_consensus)
ntips <- length(tree_consensus$tip.label)
```

---

## 3. Interactive Diagnostics & Parameter Tuning / Интерактивная диагностика и настройка
**EN:** The core diagnostic loop. It runs bgmyc.singlephy, evaluates acceptance rates and logposterior stationarity, plots trace graphs, and prompts you to adjust MCMC parameters before scaling to multiple trees.    
**RU:** Основной диагностический цикл. Запускает bgmyc.singlephy, оценивает acceptance rates и стационарность логарифма правдоподобия, строит графики и предлагает настроить параметры MCMC перед многодеревным анализом.

```{r diagnostics}
params <- list(
  mcmc = 10000, burnin = 1000, thinning = 10,
  py1 = 0, py2 = 1.5, pc1 = 0, pc2 = 2,
  t1 = 2, t2 = min(35, ntips - 1),
  scale = c(20, 10, 5), start = c(1, 1, floor(ntips/3))
)

# EN: Diagnostic tuning loop / Цикл диагностики и настройки
repeat {
  res_single <- bgmyc.singlephy(
    phylo = tree_consensus, mcmc = params$mcmc, burnin = params$burnin, 
    thinning = params$thinning, py1 = params$py1, py2 = params$py2, 
    pc1 = params$pc1, pc2 = params$pc2, t1 = params$t1, t2 = params$t2, 
    scale = params$scale, start = params$start
  )
  
  # EN: Convergence checks / Проверка сходимости
  ar <- res_single$accept
  cat(sprintf("Acceptance rates: py=%.3f | pc=%.3f | th=%.3f\n", ar[1], ar[2], ar[3]))
  if (requireNamespace("mcmcse", quietly = TRUE)) {
    ess_vals <- sapply(1:4, function(col) round(mcmcse::ess(res_single$par[, col])))
    cat(sprintf("ESS: py=%d | pc=%d | th=%d | logL=%d [Target > 200]\n", 
                ess_vals[1], ess_vals[2], ess_vals[3], ess_vals[4]))
  }
  
  plot(res_single)
  
  if (!interactive() || tolower(readline("Accept parameters? / Принять параметры? (y/n): ")) != "n") break
  
  # EN: Update parameters / Обновление параметров
  params$mcmc <- input_scalar("mcmc", params$mcmc)
  params$burnin <- input_scalar("burnin", params$burnin)
  params$thinning <- input_scalar("thinning", params$thinning)
  params$scale <- input_vector("scale", params$scale)
  params$start <- input_vector("start", params$start)
}
```

---

## 4. Multi-Tree Analysis / Многодеревный анализ    
**EN:** Once diagnostics are stable, bgmyc.multiphylo runs on the sampled trees.     
**RU:** После стабилизации диагностики bgmyc.multiphylo запускается на выбранных деревьях.     

```{r multi}
if (analysis_mode == 2) {
  # EN: Parallel execution / Параллельное выполнение
  n_workers <- min(parallel::detectCores(logical = FALSE) - 1, length(trees_sample))
  plan(multisession, workers = max(1, n_workers))
  
  final_res <- future_lapply(seq_along(trees_sample), function(i) {
    bgmyc.singlephy(
      phylo = trees_sample[[i]], mcmc = params$mcmc, burnin = params$burnin, 
      thinning = params$thinning, py1 = params$py1, py2 = params$py2, 
      pc1 = params$pc1, pc2 = params$pc2, t1 = params$t1, t2 = params$t2, 
      scale = params$scale, start = params$start
    )
  }, future.seed = TRUE)
  class(final_res) <- "multibgmyc"
  
  # EN: Gelman-Rubin R-hat / Статистика Гелмана-Рубина
  if (requireNamespace("mcmcse", quietly = TRUE) && length(final_res) > 1) {
    chains_list <- lapply(final_res, function(res) res$par[, 3])
    names(chains_list) <- paste0("Tree_", seq_along(final_res))
    rhat <- round(mcmcse::gelman(chains_list)$Rhat, 3)
    cat(sprintf("Gelman-Rubin R-hat: %.3f [Target < 1.05]\n", rhat))
  }
} else {
  final_res <- list(res_single)
  class(final_res) <- "multibgmyc"
}
```
---


## 5. Custom Interactive Visualization / Кастомная визуализация (Plotly)    
**EN:** Custom interactive visualization.     
**RU:** Визуализация результатов.     

```{r plotly_viz}
probmat <- spec.probmat(final_res)

# EN: Synchronize tip order between tree and matrix / Синхронизация порядка таксонов
p_tree <- suppressWarnings(ggtree(tree_beast, layout = "rectangular"))
tips_data <- p_tree$data %>% filter(isTip) %>% arrange(y)
tip_order <- tips_data$label
probmat <- probmat[tip_order, tip_order]

# EN: Extract posterior probabilities for branches / Извлечение posterior для ветвей
post_col <- intersect(c("posterior", "prob", "Posterior"), colnames(p_tree$data))[1]
node_posterior <- p_tree$data[[post_col]]
names(node_posterior) <- p_tree$data$node

# EN: Custom color gradient function / Функция градиента цвета
get_pp_color <- function(pp) {
  if(is.na(pp)) return("#CCCCCC")
  pp <- max(0, min(1, pp))
  colors <- list(c(0.0, 1.0, 0.0, 0.0), c(0.25, 1.0, 0.5, 0.0), 
                 c(0.5, 1.0, 1.0, 0.0), c(0.75, 0.5, 1.0, 0.0), c(1.0, 0.0, 0.7, 0.0))
  for(i in 1:(length(colors)-1)) {
    if(pp >= colors[[i]][1] && pp <= colors[[i+1]][1]) {
      t <- (pp - colors[[i]][1]) / (colors[[i+1]][1] - colors[[i]][1])
      r <- colors[[i]][2] + t * (colors[[i+1]][2] - colors[[i]][2])
      g <- colors[[i]][3] + t * (colors[[i+1]][3] - colors[[i]][3])
      b <- colors[[i]][4] + t * (colors[[i+1]][4] - colors[[i]][4])
      return(sprintf("#%02X%02X%02X", round(r*255), round(g*255), round(b*255)))
    }
  }
  return("#CCCCCC")
}

# EN: Build Plotly Tree / Построение дерева в Plotly
edges <- p_tree$data %>% filter(!is.na(parent))
fig_tree <- plot_ly()
for(i in 1:nrow(edges)) {
  child <- edges[i, ]
  parent <- p_tree$data %>% filter(node == child$parent)
  branch_color <- get_pp_color(node_posterior[as.character(child$node)])
  fig_tree <- fig_tree %>% add_segments(
    x = parent$x, xend = child$x, y = child$y, yend = child$y,
    line = list(color = branch_color, width = 5), showlegend = FALSE)
  fig_tree <- fig_tree %>% add_segments(
    x = parent$x, xend = parent$x, y = parent$y, yend = child$y,
    line = list(color = branch_color, width = 5), showlegend = FALSE)
}

# EN: Build Plotly Heatmap / Построение тепловой карты
fig_heat <- plot_ly(
  z = probmat, x = 1:nrow(probmat), y = 1:ncol(probmat), type = "heatmap",
  colorscale = list(list(0.0, "#F7FCF5"), list(0.5, "#41AB5D"), list(1.0, "#00441B")),
  zmin = 0, zmax = 1, showscale = TRUE
)

# EN: Combine 1:1 / Объединение 1:1
fig_combined <- subplot(fig_tree, fig_heat, nrows = 1, widths = c(0.5, 0.5), shareY = TRUE) %>%
  layout(yaxis = list(autorange = "reversed", showticklabels = FALSE),
         xaxis = list(showticklabels = FALSE),
         xaxis2 = list(showticklabels = FALSE),
         yaxis2 = list(autorange = "reversed", showticklabels = FALSE))

htmlwidgets::saveWidget(fig_combined, "bGMYC_interactive_heatmap.html", selfcontained = TRUE)
```
---

## 6. Delimitation & Export / Делимитация и экспорт    
**EN:** Delimitation.     
**RU:** Делимитация.     

```{r export}
# EN: Export clusters at different thresholds / Экспорт кластеров на разных порогах
for (p in c(0.05, 0.01)) {
  out <- bgmyc.point(probmat, ppcutoff = p)
  df <- data.frame(
    Sequence = unlist(out),
    MOTU_bGMYC = rep(seq_along(out), lengths(out)),
    stringsAsFactors = FALSE
  )
  write.table(df, file = sprintf("Delimitation_bGMYC_%.2f.csv", p),
              row.names = FALSE, sep = ";", dec = ".", quote = FALSE, fileEncoding = "UTF-8")
}

# EN: Export full probability matrix / Экспорт полной матрицы вероятностей
spec_out <- bgmyc.spec(final_res)
write.csv(spec_out$specprobs, "bGMYC_full_probs.csv", row.names = FALSE)
```
---

## ⚡ Performance & Export / Производительность и экспорт    
**EN:** bgmyc.multiphylo() automatically parallelizes across physical CPU cores. To limit cores (e.g., to 6), run before analysis: options(mc.cores = 6).    
**RU:** bgmyc.multiphylo() автоматически использует все физические ядра. Чтобы ограничить число ядер (например, до 6), выполните перед запуском: options(mc.cores = 6).    
---

## 📊 Parameter Reference Guide / Справочник по параметрам
| Parameter / Параметр | Role / Роль | Recommended Range | Biological/Statistical Notes / Примечания |
|:---|:---|:---|:---|
| `mcmc` | Chain length / Длина цепи | `10k` (test), `50k+` (final) | Longer chains improve posterior resolution. / Длинные цепи улучшают апостериорную оценку. |
| `burnin` | Warm-up / Разогрев | `20–30%` of `mcmc` | Discards non-stationary start. Increase if logposterior drifts. / Отбрасывает неустановившуюся фазу. |
| `thinning` | Sampling interval / Интервал выборки | `10–50` | Reduces autocorrelation & RAM. Higher for long chains. / Снижает автокорреляцию и нагрузку на память. |
| `py1`, `py2` | Yule rate prior / Априор видообразования | `0`, `0.5–1.5` | Model: `λ ∝ n^py`. `py > 1.5` blurs Yule/Coalescent boundary. / >1.5 размывает границу модели. |
| `pc1`, `pc2` | Coalescent prior / Априор коалесценции | `0`, `1.0–2.0` | Models `Ne` change. `pc < 1` → decline, `pc > 1` → growth. / Моделирует динамику эффективного размера популяции. |
| `t1`, `t2` | Threshold prior (species count) / Априор числа видов | `2`, `min(35, ntips-5)` | Must be `< ntips`. Auto-capped to prevent crashes. / Строго `< числа таксонов`. |
| `scale` | Proposal step widths / Ширина предложений MCMC | `c(20–30, 10–15, 3–7)` | Higher = more conservative. Tune via acceptance rates. / Выше = консервативнее. Настраивается по acceptance rates. |
| `ppcutoff` | Lumping threshold / Порог объединения | `0.05` (variable), `0.95` (strict) | Low = captures high intraspecific variation/ILS. / Низкий = учитывает высокую внутривидовую изменчивость. |

---

## 🔍 Convergence Diagnostics Guide / Руководство по диагностике сходимости    
**EN:** After the diagnostic run, evaluate two metrics:    
1. **Acceptance Rates**: `0.20–0.40` is optimal. `<0.15` → decrease `scale`. `>0.50` → increase `scale`.    
2. **Trace Plots ("Fuzzy Caterpillars")**: After burnin, parameters should oscillate horizontally around a stable mean. Upward/downward trends → increase burnin or mcmc. Sticky steps → scale too high. Smooth/lazy curves → scale too low. Hitting bounds → widen priors.    

**RU:** После диагностического запуска проверьте два показателя:    
1. **Acceptance Rates**: Оптимум `0.20–0.40`. `<0.15` → уменьшите `scale`. `>0.50` → увеличьте `scale`.    
2. **Графики ("пушистые гусеницы")**: После burnin параметры должны колебаться горизонтально вокруг стабильного среднего. Тренды вверх/вниз → увеличьте `burnin` или `mcmc`. Ступеньки → `scale` слишком высок. Гладкие кривые → `scale` слишком низок. Прилипание к границам → расширьте априоры.    

---

## 📎 Notes for CRAN & Local Use / Примечания для CRAN и локального запуска    
**EN:**
- `eval = FALSE` in code chunks prevents CRAN from timing out during automated checks.    
- To run interactively: change `eval = FALSE` to `eval = TRUE` in the first chunk, or simply copy-paste chunks into RStudio console and execute sequentially.    
- Always run diagnostics on 1 tree before scaling to `multiPhylo`.    

**RU:**
- `eval = FALSE` предотвращает таймауты при автоматической проверке CRAN.    
- Для локального запуска: измените `eval = FALSE` на `eval = TRUE` в первом блоке, или копируйте блоки в консоль RStudio и запускайте последовательно.    
- Всегда запускайте диагностику на 1 дереве перед переходом к `multiPhylo`.    

---
**References / Ссылки:** Pons et al. 2006 *Syst. Biol.* 55:595; Reid & Carstens 2012 *Mol. Ecol. Res.* 12:446.

`vignette("bGMYC4-interactive", package = "bGMYC4")`
