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

## ----setup--------------------------------------------------------------------
library(taxdiv)

## ----data---------------------------------------------------------------------
# Species abundances (individuals per plot)
community <- c(
  Quercus_coccifera    = 25,
  Quercus_infectoria   = 18,
  Pinus_brutia         = 30,
  Pinus_nigra          = 12,
  Juniperus_excelsa    = 8,
  Juniperus_oxycedrus  = 6,
  Arbutus_andrachne    = 15,
  Styrax_officinalis   = 4,
  Cercis_siliquastrum  = 3,
  Olea_europaea        = 10
)

# Taxonomic hierarchy
tax_tree <- build_tax_tree(
  species = names(community),
  Genus   = c("Quercus", "Quercus", "Pinus", "Pinus",
              "Juniperus", "Juniperus", "Arbutus", "Styrax",
              "Cercis", "Olea"),
  Family  = c("Fagaceae", "Fagaceae", "Pinaceae", "Pinaceae",
              "Cupressaceae", "Cupressaceae", "Ericaceae", "Styracaceae",
              "Fabaceae", "Oleaceae"),
  Order   = c("Fagales", "Fagales", "Pinales", "Pinales",
              "Pinales", "Pinales", "Ericales", "Ericales",
              "Fabales", "Lamiales")
)

tax_tree

## ----classical----------------------------------------------------------------
# Shannon diversity (natural log)
H <- shannon(community)
cat("Shannon H':", round(H, 4), "\n")

# Simpson indices
D <- simpson(community, type = "dominance")
GS <- simpson(community, type = "gini_simpson")
inv_D <- simpson(community, type = "inverse")

cat("Simpson dominance (D):", round(D, 4), "\n")
cat("Gini-Simpson (1-D):", round(GS, 4), "\n")
cat("Inverse Simpson (1/D):", round(inv_D, 4), "\n")

## ----distance-----------------------------------------------------------------
dist_mat <- tax_distance_matrix(tax_tree)
round(dist_mat, 2)

## ----cw-----------------------------------------------------------------------
# Delta: average taxonomic diversity (abundance-weighted)
d <- delta(community, tax_tree)
cat("Delta (taxonomic diversity):", round(d, 4), "\n")

# Delta*: taxonomic distinctness (abundance-weighted, excludes same-species)
ds <- delta_star(community, tax_tree)
cat("Delta* (taxonomic distinctness):", round(ds, 4), "\n")

# AvTD (Delta+): average taxonomic distinctness (presence/absence)
spp <- names(community)
avg_td <- avtd(spp, tax_tree)
cat("AvTD (Delta+):", round(avg_td, 4), "\n")

# VarTD (Lambda+): variation in taxonomic distinctness
var_td <- vartd(spp, tax_tree)
cat("VarTD (Lambda+):", round(var_td, 4), "\n")

## ----pto----------------------------------------------------------------------
result <- ozkan_pto(community, tax_tree)

cat("uTO  (unweighted diversity):", round(result$uTO, 4), "\n")
cat("TO   (weighted diversity):", round(result$TO, 4), "\n")
cat("uTO+ (unweighted distance):", round(result$uTO_plus, 4), "\n")
cat("TO+  (weighted distance):", round(result$TO_plus, 4), "\n")

## ----ed_levels----------------------------------------------------------------
cat("Deng entropy at each level:\n")
for (i in seq_along(result$Ed_levels)) {
  cat("  ", names(result$Ed_levels)[i], ":",
      round(result$Ed_levels[i], 4), "\n")
}

## ----components---------------------------------------------------------------
pto_components(community, tax_tree)

## ----compare------------------------------------------------------------------
# Degraded community (fewer species, less taxonomic breadth)
degraded <- c(
  Quercus_coccifera   = 40,
  Pinus_brutia        = 35,
  Juniperus_oxycedrus = 10
)

tax_degraded <- tax_tree[tax_tree$Species %in% names(degraded), ]

cat("=== Original community (10 species) ===\n")
cat("Shannon:", round(shannon(community), 4), "\n")
r1 <- ozkan_pto(community, tax_tree)
cat("uTO+:", round(r1$uTO_plus, 4), "\n")
cat("TO+:", round(r1$TO_plus, 4), "\n\n")

cat("=== Degraded community (3 species) ===\n")
cat("Shannon:", round(shannon(degraded), 4), "\n")
r2 <- ozkan_pto(degraded, tax_degraded)
cat("uTO+:", round(r2$uTO_plus, 4), "\n")
cat("TO+:", round(r2$TO_plus, 4), "\n")

## ----run2---------------------------------------------------------------------
set.seed(42)
run2 <- ozkan_pto_resample(community, tax_tree, n_iter = 101, seed = 42)

cat("=== Run 1 (deterministic) ===\n")
cat("uTO+:", round(run2$uTO_plus_det, 4), "\n")
cat("TO+: ", round(run2$TO_plus_det, 4), "\n")
cat("uTO: ", round(run2$uTO_det, 4), "\n")
cat("TO:  ", round(run2$TO_det, 4), "\n\n")

cat("=== Run 2 (max across", run2$n_iter, "iterations) ===\n")
cat("uTO+:", round(run2$uTO_plus_max, 4), "\n")
cat("TO+: ", round(run2$TO_plus_max, 4), "\n")
cat("uTO: ", round(run2$uTO_max, 4), "\n")
cat("TO:  ", round(run2$TO_max, 4), "\n")

## ----run2_summary-------------------------------------------------------------
iter_df <- run2$iteration_results
cat("uTO+ range:", round(min(iter_df$uTO_plus), 4), "to",
    round(max(iter_df$uTO_plus), 4), "\n")
cat("TO+  range:", round(min(iter_df$TO_plus), 4), "to",
    round(max(iter_df$TO_plus), 4), "\n")

## ----run3---------------------------------------------------------------------
run3 <- ozkan_pto_sensitivity(community, tax_tree, run2, seed = 123)

cat("=== Run 3 (sensitivity analysis) ===\n")
cat("Run 3 max uTO+:", round(run3$run3_uTO_plus_max, 4), "\n")
cat("Run 3 max TO+: ", round(run3$run3_TO_plus_max, 4), "\n\n")

cat("=== Overall max (across Run 1 + 2 + 3) ===\n")
cat("uTO+:", round(run3$uTO_plus_max, 4), "\n")
cat("TO+: ", round(run3$TO_plus_max, 4), "\n")
cat("uTO: ", round(run3$uTO_max, 4), "\n")
cat("TO:  ", round(run3$TO_max, 4), "\n")

## ----species_probs------------------------------------------------------------
probs <- run3$species_probs
prob_df <- data.frame(
  Species = names(probs),
  Probability = round(probs, 4)
)
print(prob_df, row.names = FALSE)

## ----pipeline_summary---------------------------------------------------------
cat("Pipeline: Run 1 -> Run 2 -> Run 3\n\n")

cat("         uTO+      TO+       uTO       TO\n")
cat("Run 1:", sprintf("%9.4f %9.4f %9.4f %9.4f",
    run2$uTO_plus_det, run2$TO_plus_det,
    run2$uTO_det, run2$TO_det), "\n")
cat("Run 2:", sprintf("%9.4f %9.4f %9.4f %9.4f",
    run2$uTO_plus_max, run2$TO_plus_max,
    run2$uTO_max, run2$TO_max), "\n")
cat("Run 3:", sprintf("%9.4f %9.4f %9.4f %9.4f",
    run3$uTO_plus_max, run3$TO_plus_max,
    run3$uTO_max, run3$TO_max), "\n")

## ----rarefaction--------------------------------------------------------------
rare <- rarefaction_taxonomic(community, tax_tree,
                               index = "shannon",
                               steps = 10, n_boot = 100, seed = 42)
cat("Rarefaction results (Shannon):\n")
print(round(rare, 4))

## ----dendrogram, fig.width=9, fig.height=5.5, fig.alt="Dendrogram showing species grouped by family with abundance labels"----
plot_taxonomic_tree(tax_tree, community = community,
                    color_by = "Family", label_size = 3.5,
                    title = "Mediterranean Forest - Taxonomic Tree")

## ----heatmap, fig.width=9, fig.height=8, fig.alt="Heatmap showing pairwise taxonomic distances between species"----
plot_heatmap(tax_tree, label_size = 2.8,
             title = "Taxonomic Distance Heatmap")

## ----comparison_data----------------------------------------------------------
# Dominant community: same species, uneven abundances
dominant_community <- c(
  Quercus_coccifera   = 80, Quercus_infectoria  = 5,
  Pinus_brutia        = 3,  Pinus_nigra         = 2,
  Juniperus_excelsa   = 2,  Juniperus_oxycedrus = 1,
  Arbutus_andrachne   = 3,  Styrax_officinalis  = 1,
  Cercis_siliquastrum = 2,  Olea_europaea       = 1
)

communities <- list(
  Diverse  = community,
  Dominant = dominant_community
)

## ----comparison_plot, fig.width=10, fig.height=6, fig.alt="Bar chart comparing 14 diversity indices between diverse and dominant communities"----
comparison <- compare_indices(communities, tax_tree, plot = TRUE)
comparison$plot

## ----comparison_table---------------------------------------------------------
comparison$table

## ----iteration, fig.width=9, fig.height=5, fig.alt="Scatter plot showing TO values across 101 stochastic resampling iterations"----
plot_iteration(run2, component = "TO",
               title = "Run 2 - TO Values Across Iterations")

## ----bubble, fig.width=10, fig.height=7, fig.alt="Bubble chart showing species contributions to diversity colored by family"----
plot_bubble(community, tax_tree, color_by = "Family",
            title = "Species Contributions to Diversity")

## ----radar, fig.width=8, fig.height=8, fig.alt="Radar chart comparing normalized diversity index values between diverse and dominant communities"----
plot_radar(communities, tax_tree,
           title = "Diverse vs Dominant - Radar Comparison")

## ----rarefaction_plot, fig.width=8, fig.height=5, fig.alt="Rarefaction curve for Shannon index with bootstrap confidence interval"----
plot_rarefaction(rare)