phase3_display_labels <- c(
  mae = "MAE",
  rmse = "RMSE",
  pbias = "PBIAS (%)",
  nse = "NSE",
  kge_2009 = "KGE (2009)"
)

test_that("phase 3 approved metrics remain live and keep finalized registry linkage", {
  public_metrics <- .hydroeval_registry_view(
    include_non_public = FALSE,
    status = "planned"
  )
  phase3_metrics <- c("mae", "rmse", "pbias", "nse", "kge_2009")

  expect_true(all(phase3_metrics %in% public_metrics$canonical_id))

  for (metric_name in phase3_metrics) {
    index <- match(metric_name, public_metrics$canonical_id)
    provenance <- .hydroeval_get_metric_provenance(
      public_metrics$provenance_key[[index]]
    )

    expect_true(length(provenance$primary_refs) >= 1L)
    expect_true(
      exists(
        public_metrics$implementation_target[[index]],
        mode = "function",
        inherits = TRUE
      )
    )
  }
})

test_that("public metric provenance keys resolve to bibliography entries", {
  public_metrics <- .hydroeval_registry_view(
    include_non_public = FALSE,
    status = "planned"
  )
  bibliography <- readLines(
    .hydroeval_test_reference_path(),
    warn = FALSE
  )
  bib_keys <- sub(
    "^@(?:article|book|manual|misc|incollection|inproceedings)\\{([^,]+),.*$",
    "\\1",
    grep(
      "^@(?:article|book|manual|misc|incollection|inproceedings)\\{",
      bibliography,
      value = TRUE
    )
  )

  for (index in seq_len(nrow(public_metrics))) {
    provenance <- .hydroeval_get_metric_provenance(
      public_metrics$provenance_key[[index]]
    )

    expect_true(all(provenance$primary_refs %in% bib_keys))
    expect_true(all(provenance$interpretive_refs %in% bib_keys))
  }

  pbias_provenance <- .hydroeval_get_metric_provenance("gupta_1999_pbias")

  expect_identical(pbias_provenance$primary_refs, "gupta_1999")
  expect_identical(pbias_provenance$interpretive_refs, "moriasi_2007")
})

test_that("phase 3 approved metrics keep locked scientific display labels", {
  public_metrics <- .hydroeval_registry_view(
    include_non_public = FALSE,
    status = "planned"
  )
  index <- match(names(phase3_display_labels), public_metrics$canonical_id)

  expect_false(anyNA(index))
  expect_false(anyNA(public_metrics$display_label[index]))
  expect_true(all(nzchar(public_metrics$display_label[index])))
  expect_identical(
    stats::setNames(public_metrics$display_label[index], public_metrics$canonical_id[index]),
    phase3_display_labels
  )
})

test_that("public wrappers return unnamed scalar doubles on success", {
  observed <- c(1, 2, 3, 4)
  simulated <- c(2, 3, 4, 5)
  wrappers <- list(
    mae = mae,
    rmse = rmse,
    pbias = pbias,
    nse = nse,
    kge_2009 = kge_2009
  )

  for (metric_name in names(wrappers)) {
    value <- wrappers[[metric_name]](observed, simulated)

    expect_type(value, "double")
    expect_length(value, 1L)
    expect_null(names(value))
    expect_false(is.na(value))
  }
})

test_that("approved wrappers compute expected values", {
  observed <- c(1, 2, 3, 4)
  simulated <- c(1, 3, 2, 5)
  correlation <- stats::cor(observed, simulated)
  alpha <- stats::sd(simulated) / stats::sd(observed)
  beta <- mean(simulated) / mean(observed)

  expect_equal(mae(observed, simulated), 0.75, tolerance = 1e-12)
  expect_equal(rmse(observed, simulated), sqrt(0.75), tolerance = 1e-12)
  expect_equal(pbias(observed, simulated), 10, tolerance = 1e-12)
  expect_equal(nse(observed, simulated), 0.4, tolerance = 1e-12)
  expect_equal(
    kge_2009(observed, simulated),
    1 - sqrt((correlation - 1)^2 + (alpha - 1)^2 + (beta - 1)^2),
    tolerance = 1e-12
  )
})

test_that("wrappers respect shared validation contract and NA policy", {
  expect_equal(
    mae(
      observed = c(1, NA, 3),
      simulated = c(1, 9, 5),
      na_policy = "omit"
    ),
    1,
    tolerance = 1e-12
  )

  expect_error(
    mae(c(1, 2, 3), c(1, 2)),
    class = "hydroeval_validation_error",
    regexp = "equal length"
  )

  expect_error(
    rmse(c(1, NA, 3), c(1, 9, 5), na_policy = "fail"),
    class = "hydroeval_validation_error",
    regexp = "Missing `NA` values are not allowed"
  )
})

test_that("implemented undefined degeneracies fail explicitly", {
  expect_error(
    pbias(c(-1, 1), c(0, 0)),
    class = "hydroeval_metric_degeneracy",
    regexp = "zero_observed_sum"
  )

  expect_error(
    nse(c(2, 2, 2), c(1, 2, 3)),
    class = "hydroeval_metric_degeneracy",
    regexp = "constant_observed_series"
  )

  expect_error(
    kge_2009(c(1, 2, 3), c(5, 5, 5)),
    class = "hydroeval_metric_degeneracy",
    regexp = "constant_simulated_series"
  )
})

test_that("kge_2009 allows zero simulated mean when components stay defined", {
  observed <- c(1, 2, 3)
  simulated <- c(-1, 0, 1)

  expect_equal(kge_2009(observed, simulated), 0, tolerance = 1e-12)
})

test_that("approved metrics return ideal values under perfect agreement", {
  observed <- c(1, 2, 3, 4)
  simulated <- c(1, 2, 3, 4)

  expect_equal(mae(observed, simulated), 0, tolerance = 1e-12)
  expect_equal(rmse(observed, simulated), 0, tolerance = 1e-12)
  expect_equal(pbias(observed, simulated), 0, tolerance = 1e-12)
  expect_equal(nse(observed, simulated), 1, tolerance = 1e-12)
  expect_equal(kge_2009(observed, simulated), 1, tolerance = 1e-12)
})

test_that("approved metrics return explicit shifted-case values", {
  observed <- c(1, 2, 3, 4)
  simulated <- c(2, 3, 4, 5)

  expect_equal(mae(observed, simulated), 1, tolerance = 1e-12)
  expect_equal(rmse(observed, simulated), 1, tolerance = 1e-12)
  expect_equal(pbias(observed, simulated), 40, tolerance = 1e-12)
  expect_equal(nse(observed, simulated), 0.2, tolerance = 1e-12)
  expect_equal(kge_2009(observed, simulated), 0.6, tolerance = 1e-12)
})

test_that("nse returns zero for mean-predictor baseline", {
  observed <- c(1, 2, 3, 4)
  simulated <- c(2.5, 2.5, 2.5, 2.5)

  expect_equal(nse(observed, simulated), 0, tolerance = 1e-12)
})

test_that("constant observed degeneracy is explicit for efficiency metrics", {
  observed <- c(3, 3, 3, 3)
  simulated <- c(3, 4, 5, 6)

  expect_error(
    nse(observed, simulated),
    class = "hydroeval_metric_degeneracy",
    regexp = "constant_observed_series"
  )

  expect_error(
    kge_2009(observed, simulated),
    class = "hydroeval_metric_degeneracy",
    regexp = "constant_observed_series"
  )
})

test_that("pbias fails when observed sum is zero", {
  observed <- c(-1, 1, -2, 2)
  simulated <- c(0, 0, 0, 0)

  expect_error(
    pbias(observed, simulated),
    class = "hydroeval_metric_degeneracy",
    regexp = "zero_observed_sum"
  )
})

test_that("pairwise NA handling preserves complete-pair shifted behavior", {
  observed <- c(1, NA, 2, 3, 4)
  simulated <- c(2, 99, 3, 4, 5)

  expect_equal(mae(observed, simulated, na_policy = "omit"), 1, tolerance = 1e-12)
  expect_equal(rmse(observed, simulated, na_policy = "omit"), 1, tolerance = 1e-12)
  expect_equal(pbias(observed, simulated, na_policy = "omit"), 40, tolerance = 1e-12)
  expect_equal(nse(observed, simulated, na_policy = "omit"), 0.2, tolerance = 1e-12)
  expect_equal(kge_2009(observed, simulated, na_policy = "omit"), 0.6, tolerance = 1e-12)
})