--- Numerical backend for the PGFPlots `autonode` library.
--
-- Copyright (C) 2026 Christophe Jorssen.
-- Distributed under LPPL 1.3c or later.
--
-- This module is the numerical backend of pgfplots-autonode.  TeX samples the
-- current pgfplots path and sends candidate positions to Lua in physical canvas
-- coordinates, measured in TeX points.  Lua then chooses one candidate per label
-- under hard geometric constraints and emits TeX callbacks for the final nodes.
--
-- The code deliberately avoids TeX-specific logic.  Node contents and TikZ
-- options remain on the TeX side; Lua only sees geometry, priorities, costs and
-- solver configuration.
--
-- @module pgfplots-autonode
-- @author Christophe Jorssen
-- @license LPPL-1.3c-or-later

local M = {
  _NAME = "pgfplots-autonode",
  _VERSION = "1.1.0",
}

local labels = {}
local config = {}

local DEFAULT_CONFIG = {
  debug_level = 0,
  max_iterations = 40,
  algorithm = "repair",
  bbox_mode = "axis-aligned",
  failure_mode = "warn",
  border_margin = 2.0,
  overlap_tolerance = 0.2,
  allow_outside = false,
  show_bounding_boxes = false,
  show_candidates = false,
  show_rejected_candidates = false,
  exact_max_labels = 10,
  exact_max_states = 50000,
}

local axis_rect = nil

local function reset_config()
  config = {}
  for key, value in pairs(DEFAULT_CONFIG) do
    config[key] = value
  end
end

reset_config()

local function log(level, message)
  if (config.debug_level or 0) >= level then
    texio.write_nl("term and log", "pgfplots-autonode: " .. message)
  end
end

local function warn(message)
  texio.write_nl("term and log", "pgfplots-autonode warning: " .. message)
end

local function package_error(message)
  if tex and tex.error then
    texio.write_nl("term and log", "pgfplots-autonode error: " .. message)
    tex.error("pgfplots-autonode error", { message })
  else
    error("pgfplots-autonode error: " .. message, 2)
  end
end

local function clamp(value, lo, hi)
  if value < lo then return lo end
  if value > hi then return hi end
  return value
end

local function validated_number(value, fallback, name, minimum, maximum,
    integer)
  local parsed = tonumber(value)
  local valid = parsed ~= nil and parsed == parsed
    and parsed ~= math.huge and parsed ~= -math.huge
  if valid and minimum ~= nil and parsed < minimum then valid = false end
  if valid and maximum ~= nil and parsed > maximum then valid = false end
  if valid and integer and parsed ~= math.floor(parsed) then valid = false end
  if valid then return parsed end
  warn("invalid " .. tostring(name) .. " `" .. tostring(value)
    .. "'; using " .. tostring(fallback))
  return fallback
end

local function tex_number(value)
  return string.format("%.6f", value or 0)
end

local function fmt(value)
  return string.format("%.4f", value or 0)
end

local function bool(value)
  return value and true or false
end

local function radians_to_degrees(angle)
  return angle * 180.0 / math.pi
end

-- Return an equivalent rotation angle that keeps text upright.  PGF/TikZ's
-- `sloped` nodes avoid upside-down text unless explicitly requested; the
-- autonode backend draws rotated nodes itself, therefore it must reproduce this
-- convention in Lua.
local function upright_angle(angle)
  angle = tonumber(angle) or 0.0
  while angle > 180.0 do angle = angle - 360.0 end
  while angle <= -180.0 do angle = angle + 360.0 end
  if angle > 90.0 then
    angle = angle - 180.0
  elseif angle < -90.0 then
    angle = angle + 180.0
  end
  return angle
end

local function atan2(y, x)
  -- LuaTeX provides the two-argument form math.atan(y, x).  The wrapper makes
  -- the intended operation explicit and keeps potential compatibility changes in
  -- one place.
  return math.atan(y, x)
end

local function ordered_labels()
  local list = {}
  for _, label in pairs(labels) do
    list[#list + 1] = label
  end
  table.sort(list, function(a, b)
    if a.priority ~= b.priority then
      return a.priority > b.priority
    end
    return a.id < b.id
  end)
  return list
end

local function ordered_labels_by_id()
  local list = {}
  for _, label in pairs(labels) do
    list[#list + 1] = label
  end
  table.sort(list, function(a, b) return a.id < b.id end)
  return list
end

local function aabb_overlap_area(a, b)
  local dx = math.min(a.aabb.xmax, b.aabb.xmax) - math.max(a.aabb.xmin, b.aabb.xmin)
  local dy = math.min(a.aabb.ymax, b.aabb.ymax) - math.max(a.aabb.ymin, b.aabb.ymin)
  if dx > 0 and dy > 0 then
    return dx * dy
  end
  return 0.0
end

local function aabb_overlap_dims(a, b)
  local dx = math.min(a.aabb.xmax, b.aabb.xmax) - math.max(a.aabb.xmin, b.aabb.xmin)
  local dy = math.min(a.aabb.ymax, b.aabb.ymax) - math.max(a.aabb.ymin, b.aabb.ymin)
  if dx > 0 and dy > 0 then
    return dx, dy
  end
  return 0.0, 0.0
end

local function dot(ax, ay, bx, by)
  return ax * bx + ay * by
end

local function projected_interval(corners, ax, ay)
  local first = corners[1]
  local min_value = dot(first.x, first.y, ax, ay)
  local max_value = min_value
  for i = 2, #corners do
    local value = dot(corners[i].x, corners[i].y, ax, ay)
    if value < min_value then min_value = value end
    if value > max_value then max_value = value end
  end
  return min_value, max_value
end

local function oriented_overlap(a, b, tolerance)
  -- Separating axis theorem for two oriented rectangles.  The axes to test are
  -- the local edge normals of both rectangles.  This is still O(1) per pair but
  -- has a larger constant than an axis-aligned test.
  local axes = {
    { x = a.ux, y = a.uy },
    { x = a.vx, y = a.vy },
    { x = b.ux, y = b.uy },
    { x = b.vx, y = b.vy },
  }
  for _, axis in ipairs(axes) do
    local amin, amax = projected_interval(a.corners, axis.x, axis.y)
    local bmin, bmax = projected_interval(b.corners, axis.x, axis.y)
    local overlap = math.min(amax, bmax) - math.max(amin, bmin)
    if overlap <= tolerance then
      return false
    end
  end
  return true
end

local function boxes_overlap(a, b)
  local tolerance = config.overlap_tolerance or 0.0
  if config.bbox_mode == "oriented" then
    return oriented_overlap(a, b, tolerance)
  end
  local dx, dy = aabb_overlap_dims(a, b)
  return dx > tolerance and dy > tolerance
end

local function overlap_penalty(a, b)
  if not boxes_overlap(a, b) then
    return 0.0
  end
  local area = aabb_overlap_area(a, b)
  if area <= 0 then
    -- Oriented rectangles may overlap while their AABB intersection is tiny due
    -- to numerical tolerance.  Return a small positive penalty so that the
    -- optimizer still treats the pair as undesirable.
    return 1.0
  end
  return area
end

local function make_box(label, candidate)
  local angle = label.sloped and candidate.angle or 0.0
  local rad = math.rad(angle)
  local ux, uy = math.cos(rad), math.sin(rad)
  local vx, vy = -math.sin(rad), math.cos(rad)

  local cx, cy = candidate.x, candidate.y
  local xmin = label.local_xmin - label.clearance
  local xmax = label.local_xmax + label.clearance
  local ymin = label.local_ymin - label.clearance
  local ymax = label.local_ymax + label.clearance
  local function corner(x, y)
    return {x = cx + x * ux + y * vx, y = cy + x * uy + y * vy}
  end

  local corners = {
    corner(xmin, ymin),
    corner(xmax, ymin),
    corner(xmax, ymax),
    corner(xmin, ymax),
  }

  local xmin, xmax = corners[1].x, corners[1].x
  local ymin, ymax = corners[1].y, corners[1].y
  for i = 2, #corners do
    local c = corners[i]
    if c.x < xmin then xmin = c.x end
    if c.x > xmax then xmax = c.x end
    if c.y < ymin then ymin = c.y end
    if c.y > ymax then ymax = c.y end
  end

  return {
    aabb = { xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax },
    corners = corners,
    ux = ux, uy = uy,
    vx = vx, vy = vy,
  }
end

local function candidate_inside_visible_box(candidate)
  if config.allow_outside then
    return true, nil
  end
  if not axis_rect then
    return false, "visible axis rectangle is unavailable"
  end

  local margin = config.border_margin or 0.0
  local left = axis_rect.left + margin
  local right = axis_rect.right - margin
  local bottom = axis_rect.bottom + margin
  local top = axis_rect.top - margin

  if config.bbox_mode == "oriented" then
    for _, corner in ipairs(candidate.box.corners) do
      if corner.x < left then
        return false, string.format("outside visible box: left violation %.4fpt", left - corner.x)
      end
      if corner.x > right then
        return false, string.format("outside visible box: right violation %.4fpt", corner.x - right)
      end
      if corner.y < bottom then
        return false, string.format("outside visible box: bottom violation %.4fpt", bottom - corner.y)
      end
      if corner.y > top then
        return false, string.format("outside visible box: top violation %.4fpt", corner.y - top)
      end
    end
    return true, nil
  end

  local box = candidate.box.aabb
  if box.xmin < left then
    return false, string.format("outside visible box: left violation %.4fpt", left - box.xmin)
  end
  if box.xmax > right then
    return false, string.format("outside visible box: right violation %.4fpt", box.xmax - right)
  end
  if box.ymin < bottom then
    return false, string.format("outside visible box: bottom violation %.4fpt", bottom - box.ymin)
  end
  if box.ymax > top then
    return false, string.format("outside visible box: top violation %.4fpt", box.ymax - top)
  end
  return true, nil
end

local function local_cost(label, candidate)
  return label.preferred_weight * (candidate.pos - label.preferred) ^ 2
end

local function count_potential_conflicts(label)
  local count = 0
  for _, other in pairs(labels) do
    if other.id ~= label.id then
      for _, c1 in ipairs(label.valid_candidates) do
        for _, c2 in ipairs(other.valid_candidates or {}) do
          if boxes_overlap(c1.box, c2.box) then
            count = count + 1
            break
          end
        end
      end
    end
  end
  return count
end

local function sort_candidates_for_label(label)
  table.sort(label.candidates, function(a, b)
    if a.valid ~= b.valid then
      return a.valid and not b.valid
    end
    if math.abs(a.local_cost - b.local_cost) > 1e-12 then
      return a.local_cost < b.local_cost
    end
    return a.pos < b.pos
  end)

  label.valid_candidates = {}
  label.invalid_candidates = {}
  for _, candidate in ipairs(label.candidates) do
    if candidate.valid then
      label.valid_candidates[#label.valid_candidates + 1] = candidate
    else
      label.invalid_candidates[#label.invalid_candidates + 1] = candidate
    end
  end
end

local function prepare_candidates(label)
  local candidates = label.candidates
  local n = #candidates
  if n == 0 then
    log(1, string.format("label #%d has no candidate", label.id))
    return
  end

  table.sort(candidates, function(a, b) return a.pos < b.pos end)

  for i, candidate in ipairs(candidates) do
    -- Tangents are estimated in physical canvas coordinates.  This is essential
    -- for anisotropic axes: a data slope is not the same thing as a visual slope.
    local previous_candidate = candidates[math.max(1, i - 1)]
    local next_candidate = candidates[math.min(n, i + 1)]
    if previous_candidate.component ~= candidate.component then
      previous_candidate = candidate
    end
    if next_candidate.component ~= candidate.component then
      next_candidate = candidate
    end
    local dx = next_candidate.base_x - previous_candidate.base_x
    local dy = next_candidate.base_y - previous_candidate.base_y
    local angle = 0.0

    if math.abs(dx) + math.abs(dy) > 1e-9 then
      angle = radians_to_degrees(atan2(dy, dx))
    end
    if label.sloped and not label.allow_upside_down then
      angle = upright_angle(angle)
    end

    candidate.angle = angle

    -- Signed normal shift.  Positive values use the left normal of the local
    -- tangent, matching the user-level aliases `above` and `below`.
    local rad = math.rad(angle)
    candidate.x = candidate.base_x - math.sin(rad) * label.normal_shift
    candidate.y = candidate.base_y + math.cos(rad) * label.normal_shift
    candidate.box = make_box(label, candidate)
    candidate.local_cost = local_cost(label, candidate)

    local inside, reason = candidate_inside_visible_box(candidate)
    candidate.valid = inside
    candidate.rejection_reason = reason

    if candidate.valid then
      log(3, string.format(
        "candidate label #%d pos=%s center=(%s,%s) angle=%s valid",
        label.id, fmt(candidate.pos), fmt(candidate.x), fmt(candidate.y), fmt(angle)
      ))
    else
      log(2, string.format(
        "reject label #%d candidate pos=%s: %s",
        label.id, fmt(candidate.pos), reason or "invalid"
      ))
    end
  end

  sort_candidates_for_label(label)
end

local function compute_difficulty(label)
  local area = (label.width + 2 * label.clearance) * (label.height + 2 * label.clearance)
  local candidates = #(label.valid_candidates or {})
  local scarcity = candidates > 0 and 1.0 / candidates or 1e6
  return 1000.0 * label.priority + 0.01 * area + 100.0 * scarcity + count_potential_conflicts(label)
end

local function ordered_by_difficulty(list)
  local sorted = {}
  for i, label in ipairs(list) do sorted[i] = label end
  table.sort(sorted, function(a, b)
    local da, db = compute_difficulty(a), compute_difficulty(b)
    if math.abs(da - db) > 1e-9 then
      return da > db
    end
    return a.id < b.id
  end)
  return sorted
end

local function candidate_conflicts_with_placed(candidate, placed)
  local penalty = 0.0
  local details = {}
  for _, other in ipairs(placed) do
    if other.choice then
      local p = overlap_penalty(candidate.box, other.choice.box)
      if p > 0 then
        penalty = penalty + p
        details[#details + 1] = string.format("overlap with #%d, approximate area=%s", other.id, fmt(p))
      end
    end
  end
  return penalty > 0, penalty, details
end

local function total_overlap_penalty_for_label(label, candidate, list)
  local penalty = 0.0
  local details = {}
  for _, other in ipairs(list) do
    if other.id ~= label.id and other.choice and not other.hidden then
      local area = overlap_penalty(candidate.box, other.choice.box)
      if area > 0 then
        -- A pair belongs to both labels.  Averaging their declared weights
        -- gives one symmetric contribution to the assignment objective and
        -- makes both labels' declared preferences participate.
        local pair_weight = 0.5 * ((label.overlap_weight or 1000.0)
          + (other.overlap_weight or 1000.0))
        local contribution = pair_weight * area
        penalty = penalty + contribution
        details[#details + 1] = string.format(
          "overlap with #%d, approximate area=%s weighted cost=%s",
          other.id, fmt(area), fmt(contribution))
      end
    end
  end
  return penalty, details
end

local function final_status(list)
  local overlaps = 0
  local hidden = 0
  local outside = 0
  for _, label in ipairs(list) do
    if label.hidden then hidden = hidden + 1 end
    if label.choice and not candidate_inside_visible_box(label.choice) then outside = outside + 1 end
  end
  for i = 1, #list do
    local a = list[i]
    if a.choice and not a.hidden then
      for j = i + 1, #list do
        local b = list[j]
        if b.choice and not b.hidden and boxes_overlap(a.choice.box, b.choice.box) then
          overlaps = overlaps + 1
        end
      end
    end
  end
  return { overlaps = overlaps, hidden = hidden, outside = outside }
end

local function greedy_place(list)
  local order = ordered_by_difficulty(list)
  local placed = {}

  for _, label in ipairs(order) do
    label.choice = nil
    label.hidden = false

    for _, candidate in ipairs(label.valid_candidates or {}) do
      local conflicts, penalty, details = candidate_conflicts_with_placed(candidate, placed)
      log(2, string.format(
        "greedy try label #%d pos=%s cost=%s conflicts=%s",
        label.id, fmt(candidate.pos), fmt(candidate.local_cost), tostring(conflicts)
      ))
      if conflicts then
        for _, detail in ipairs(details) do
          log(2, "  reject: " .. detail)
        end
      else
        label.choice = candidate
        placed[#placed + 1] = label
        log(1, string.format("greedy place label #%d at pos=%s", label.id, fmt(candidate.pos)))
        break
      end
    end

    if not label.choice then
      log(1, string.format("greedy could not place label #%d without conflicts", label.id))
      if config.failure_mode == "hide-low-priority" then
        label.hidden = true
        log(1, string.format("hide label #%d according to failure mode", label.id))
      else
        -- Keep a fallback candidate so that `allow-minimal-overlap` and `warn`
        -- can still draw something.  The final failure handler decides whether
        -- this is acceptable.
        label.choice = label.valid_candidates[1]
        placed[#placed + 1] = label
      end
    end
  end
end

local function local_search(list, only_conflicting)
  -- Discrete coordinate descent under hard visibility constraints.  Candidates
  -- that overlap existing choices are allowed during the search only when no
  -- non-overlapping alternative exists; the final status still reports conflicts.
  for iteration = 1, config.max_iterations do
    local changed = false
    log(2, string.format("local-search iteration %d", iteration))

    for _, label in ipairs(list) do
      if not label.hidden and #(label.valid_candidates or {}) > 0 then
        if not only_conflicting or not label.choice then
          -- proceed
        elseif label.choice then
          local current_penalty = total_overlap_penalty_for_label(label, label.choice, list)
          if current_penalty <= 0 then
            goto continue_label
          end
        end

        local current = label.choice
        local current_total = math.huge
        if current then
          local overlap = total_overlap_penalty_for_label(label, current, list)
          current_total = current.local_cost + overlap
        end

        local best_candidate = current
        local best_total = current_total
        local best_overlap = current and total_overlap_penalty_for_label(label, current, list) or math.huge
        local best_details = {}

        for _, candidate in ipairs(label.valid_candidates) do
          local overlap, details = total_overlap_penalty_for_label(label, candidate, list)
          local total = candidate.local_cost + overlap
          log(2, string.format(
            "local-search try label #%d pos=%s total=%s local=%s overlap=%s",
            label.id, fmt(candidate.pos), fmt(total), fmt(candidate.local_cost), fmt(overlap)
          ))
          if #details > 0 then
            for _, detail in ipairs(details) do log(2, "  conflict: " .. detail) end
          end
          if total < best_total - 1e-9 then
            best_candidate = candidate
            best_total = total
            best_overlap = overlap
            best_details = details
          end
        end

        if best_candidate and best_candidate ~= current then
          log(1, string.format(
            "move label #%d: pos %s -> %s, score %s -> %s",
            label.id,
            fmt(current and current.pos or -1),
            fmt(best_candidate.pos),
            fmt(current_total),
            fmt(best_total)
          ))
          if best_overlap > 0 then
            for _, detail in ipairs(best_details) do log(1, "  remaining " .. detail) end
          end
          label.choice = best_candidate
          changed = true
        end
      end
      ::continue_label::
    end

    if not changed then
      log(1, string.format("local-search converged after %d iteration(s)", iteration))
      break
    end
  end
end

local function solve_greedy(list)
  greedy_place(list)
end

local function solve_repair(list)
  greedy_place(list)
  local_search(list, true)
end

local function solve_local_search(list)
  -- Initialise all labels independently at the closest valid candidate to the
  -- preferred position, then run coordinate descent on the global objective.
  for _, label in ipairs(list) do
    label.hidden = false
    label.choice = label.valid_candidates[1]
    if label.choice then
      log(1, string.format("initial label #%d at pos=%s", label.id, fmt(label.choice.pos)))
    end
  end
  local_search(list, false)
end

local function exact_search(list)
  local order = ordered_by_difficulty(list)
  local best_cost = math.huge
  local best_assignment = {}
  local states = 0
  local max_states = config.exact_max_states
  local feasible_found = false

  local function recurse(index, cost)
    if states >= max_states then return end
    states = states + 1
    if cost >= best_cost then return end

    if index > #order then
      feasible_found = true
      best_cost = cost
      best_assignment = {}
      for _, label in ipairs(order) do
        best_assignment[label.id] = label.choice
      end
      log(1, string.format("exact-small new best cost=%s after %d state(s)", fmt(best_cost), states))
      return
    end

    local label = order[index]
    label.hidden = false
    local previous = label.choice

    for _, candidate in ipairs(label.valid_candidates or {}) do
      local conflicts = false
      for j = 1, index - 1 do
        local other = order[j]
        if other.choice and boxes_overlap(candidate.box, other.choice.box) then
          conflicts = true
          log(2, string.format(
            "exact-small reject label #%d pos=%s: overlaps label #%d",
            label.id, fmt(candidate.pos), other.id
          ))
          break
        end
      end
      if not conflicts then
        label.choice = candidate
        recurse(index + 1, cost + candidate.local_cost)
        if states >= max_states then break end
      end
    end

    label.choice = previous
  end

  if #list > config.exact_max_labels then
    log(1, string.format(
      "exact-small skipped: %d labels exceed exact max labels=%d; using repair",
      #list, config.exact_max_labels
    ))
    solve_repair(list)
    return
  end

  recurse(1, 0.0)

  if not feasible_found then
    log(1, string.format(
      "exact-small found no feasible assignment after %d state(s); using repair",
      states
    ))
    solve_repair(list)
    return
  end

  if states >= max_states then
    log(1, string.format("exact-small reached max states=%d", max_states))
  end

  for _, label in ipairs(list) do
    label.choice = best_assignment[label.id]
  end
  log(1, string.format("exact-small finished: states=%d best cost=%s", states, fmt(best_cost)))
end

local function handle_final_failures(list)
  local status = final_status(list)
  if status.overlaps == 0 and status.outside == 0 then
    return status
  end

  local message = string.format(
    "placement has %d overlapping pair(s), %d outside label(s), %d hidden label(s)",
    status.overlaps, status.outside, status.hidden
  )

  if config.failure_mode == "error" then
    package_error(message)
  elseif config.failure_mode == "hide-low-priority" then
    -- Greedy may already have hidden labels.  Hide the lowest-priority label of
    -- each remaining conflicting pair until the status becomes feasible or no
    -- more action is possible.
    local changed = true
    while changed do
      changed = false
      for i = 1, #list do
        local a = list[i]
        if a.choice and not a.hidden then
          for j = i + 1, #list do
            local b = list[j]
            if b.choice and not b.hidden and boxes_overlap(a.choice.box, b.choice.box) then
              local victim = a
              if b.priority < a.priority or (b.priority == a.priority and b.id > a.id) then
                victim = b
              end
              victim.hidden = true
              changed = true
              log(1, string.format("hide label #%d to remove a remaining conflict", victim.id))
              break
            end
          end
        end
        if changed then break end
      end
    end
  elseif config.failure_mode == "allow-minimal-overlap" then
    warn(message .. "; drawing least-cost assignment because failure mode is allow-minimal-overlap")
  else
    warn(message)
  end

  return final_status(list)
end

local function prepare_all_labels(list)
  local active = {}
  for _, label in ipairs(list) do
    prepare_candidates(label)
    log(1, string.format(
      "label #%d: %d candidate(s), %d valid inside visible box",
      label.id, #label.candidates, #(label.valid_candidates or {})
    ))
    if #(label.valid_candidates or {}) == 0 then
      label.choice = nil
      label.hidden = true
      local message = string.format(
        "label #%d has no valid candidate inside the visible box", label.id)
      if config.failure_mode == "error" then
        package_error(message)
      else
        warn(message .. "; hiding label")
      end
    else
      active[#active + 1] = label
    end
  end
  return active
end

local function enforce_visibility(list)
  for _, label in ipairs(list) do
    if label.choice and not label.hidden then
      local inside, reason = candidate_inside_visible_box(label.choice)
      if not inside then
        label.choice = nil
        label.hidden = true
        local message = string.format(
          "label #%d selected an outside candidate (%s)",
          label.id, reason or "unknown reason")
        if config.failure_mode == "error" then
          package_error(message)
        else
          warn(message .. "; hiding label")
        end
      end
    end
  end
end

--- Reset all axis-local labels, candidates, and geometry.
-- Configuration values are preserved for the next axis.
function M.reset()
  labels = {}
  axis_rect = nil
  log(2, "reset axis state")
end

--- Set the diagnostic verbosity.
-- @tparam number level Non-negative debug level.
function M.set_debug(level)
  config.debug_level = validated_number(level, DEFAULT_CONFIG.debug_level,
    "debug level", 0, nil, true)
  if config.debug_level > 0 then
    texio.write_nl("term and log", "pgfplots-autonode: debug level=" .. tostring(config.debug_level))
  end
end

--- Configure the placement solver for the current axis.
-- Unknown table fields are ignored. Invalid enumerated or numerical values
-- produce a warning and fall back to their documented defaults.
-- @tparam[opt={}] table options Solver configuration.
function M.configure(options)
  options = options or {}
  for key, value in pairs(options) do
    if key == "debug_level" then
      config.debug_level = validated_number(value,
        DEFAULT_CONFIG.debug_level, "debug level", 0, nil, true)
    elseif key == "max_iterations" then
      config.max_iterations = validated_number(value,
        DEFAULT_CONFIG.max_iterations, "max iterations", 1, nil, true)
    elseif key == "algorithm" then
      if value == "greedy" or value == "repair" or value == "local-search" or value == "exact-small" then
        config.algorithm = value
      else
        warn("unknown algorithm `" .. tostring(value) .. "'; using repair")
        config.algorithm = "repair"
      end
    elseif key == "bbox_mode" then
      if value == "axis-aligned" or value == "oriented" then
        config.bbox_mode = value
      else
        warn("unknown bbox mode `" .. tostring(value) .. "'; using axis-aligned")
        config.bbox_mode = "axis-aligned"
      end
    elseif key == "failure_mode" then
      if value == "warn" or value == "error" or value == "hide-low-priority" or value == "allow-minimal-overlap" then
        config.failure_mode = value
      else
        warn("unknown failure mode `" .. tostring(value) .. "'; using warn")
        config.failure_mode = "warn"
      end
    elseif key == "border_margin" then
      config.border_margin = validated_number(value,
        DEFAULT_CONFIG.border_margin, "border margin", 0)
    elseif key == "overlap_tolerance" then
      config.overlap_tolerance = validated_number(value,
        DEFAULT_CONFIG.overlap_tolerance, "overlap tolerance", 0)
    elseif key == "allow_outside" then
      config.allow_outside = bool(value)
    elseif key == "show_bounding_boxes" then
      config.show_bounding_boxes = bool(value)
    elseif key == "show_candidates" then
      config.show_candidates = bool(value)
    elseif key == "show_rejected_candidates" then
      config.show_rejected_candidates = bool(value)
    elseif key == "exact_max_labels" then
      config.exact_max_labels = validated_number(value,
        DEFAULT_CONFIG.exact_max_labels, "exact max labels", 1, nil, true)
    elseif key == "exact_max_states" then
      config.exact_max_states = validated_number(value,
        DEFAULT_CONFIG.exact_max_states, "exact max states", 1, nil, true)
    end
  end

  log(1, string.format(
    "configuration: algorithm=%s bbox=%s failure=%s margin=%spt tolerance=%spt",
    config.algorithm, config.bbox_mode, config.failure_mode,
    fmt(config.border_margin), fmt(config.overlap_tolerance)
  ))
end

--- Set the visible axis rectangle in TeX points.
-- @tparam number left Left canvas coordinate.
-- @tparam number right Right canvas coordinate.
-- @tparam number bottom Bottom canvas coordinate.
-- @tparam number top Top canvas coordinate.
function M.set_axis_rect(left, right, bottom, top)
  left = validated_number(left, 0.0, "axis left")
  right = validated_number(right, 0.0, "axis right")
  bottom = validated_number(bottom, 0.0, "axis bottom")
  top = validated_number(top, 0.0, "axis top")
  axis_rect = {
    left = math.min(left, right),
    right = math.max(left, right),
    bottom = math.min(bottom, top),
    top = math.max(bottom, top),
  }
  log(1, string.format(
    "visible rectangle: left=%s right=%s bottom=%s top=%s",
    fmt(axis_rect.left), fmt(axis_rect.right), fmt(axis_rect.bottom), fmt(axis_rect.top)
  ))
end

--- Register one measured label.
-- The positional and size values use TeX points, except `preferred`, which is
-- a normalized path position in the interval [0,1].
-- @tparam number id Positive label identifier.
-- @tparam number width Measured label width.
-- @tparam number height Measured label height.
-- @tparam number preferred Preferred normalized path position.
-- @tparam number normal_shift Shift perpendicular to the curve.
-- @tparam number clearance Extra collision clearance.
-- @tparam boolean sloped Whether to rotate the label with the curve.
-- @tparam boolean allow_upside_down Whether raw path rotation may be inverted.
-- @tparam number preferred_weight Cost of moving from the preferred position.
-- @tparam number overlap_weight Nonnegative overlap-area weight. A pair uses
-- the arithmetic mean of the two labels' weights during local optimization.
-- @tparam number priority Placement priority; higher values are handled first.
function M.add_label(id, width, height, preferred, normal_shift, clearance, sloped, allow_upside_down, preferred_weight, overlap_weight, priority, local_xmin, local_xmax, local_ymin, local_ymax)
  id = validated_number(id, 0, "label id", 1, nil, true)
  if id <= 0 then
    warn("ignored label with invalid id")
    return
  end

  width = validated_number(width, 0.0, "label width", 0)
  height = validated_number(height, 0.0, "label height", 0)
  local_xmin = local_xmin == nil and -width / 2
    or validated_number(local_xmin, -width / 2, "node left offset")
  local_xmax = local_xmax == nil and width / 2
    or validated_number(local_xmax, width / 2, "node right offset")
  local_ymin = local_ymin == nil and -height / 2
    or validated_number(local_ymin, -height / 2, "node bottom offset")
  local_ymax = local_ymax == nil and height / 2
    or validated_number(local_ymax, height / 2, "node top offset")
  if local_xmax < local_xmin or local_ymax < local_ymin then
    package_error("invalid measured node bounds for label #" .. id)
  end

  labels[id] = {
    id = id,
    width = width,
    height = height,
    local_xmin = local_xmin,
    local_xmax = local_xmax,
    local_ymin = local_ymin,
    local_ymax = local_ymax,
    preferred = validated_number(preferred, 0.5,
      "preferred position", 0, 1),
    normal_shift = validated_number(normal_shift, 0.0, "normal shift"),
    clearance = validated_number(clearance, 0.0, "label clearance", 0),
    sloped = bool(sloped),
    allow_upside_down = bool(allow_upside_down),
    preferred_weight = validated_number(preferred_weight, 8.0,
      "preferred weight", 0),
    overlap_weight = validated_number(overlap_weight, 1000.0,
      "overlap weight", 0),
    priority = validated_number(priority, 0.0, "priority"),
    candidates = {},
    path_points = {},
    valid_candidates = {},
    invalid_candidates = {},
    choice = nil,
    hidden = false,
  }

  log(1, string.format(
    "registered label #%d: width=%spt height=%spt preferred=%s normal_shift=%spt sloped=%s allow_upside_down=%s priority=%s",
    id,
    fmt(labels[id].width),
    fmt(labels[id].height),
    fmt(labels[id].preferred),
    fmt(labels[id].normal_shift),
    tostring(labels[id].sloped),
    tostring(labels[id].allow_upside_down),
    fmt(labels[id].priority)
  ))
end

--- Replace provisional dimensions with the actual unrotated TikZ node bounds.
-- Bounds are relative to the final drawing anchor, not to a presumed center.
function M.set_label_geometry(id, xmin, xmax, ymin, ymax)
  local label = labels[validated_number(id, 0, "label id", 1, nil, true)]
  if not label then
    package_error("cannot measure an unknown label #" .. tostring(id))
    return
  end
  xmin = validated_number(xmin, 0, "node left offset")
  xmax = validated_number(xmax, 0, "node right offset")
  ymin = validated_number(ymin, 0, "node bottom offset")
  ymax = validated_number(ymax, 0, "node top offset")
  if xmax <= xmin or ymax <= ymin then
    package_error("invalid measured node bounds for label #" .. id)
    return
  end
  label.local_xmin, label.local_xmax = xmin, xmax
  label.local_ymin, label.local_ymax = ymin, ymax
  label.width, label.height = xmax - xmin, ymax - ymin
  log(2, string.format(
    "measured label #%d: offsets=(%s,%s,%s,%s)pt",
    id, fmt(xmin), fmt(xmax), fmt(ymin), fmt(ymax)))
end

--- Add one sampled position for a registered label.
-- @tparam number id Label identifier.
-- @tparam number pos Normalized position on the PGFPlots path.
-- @tparam number x Canvas x coordinate in TeX points.
-- @tparam number y Canvas y coordinate in TeX points.
function M.add_candidate(id, pos, x, y, component)
  id = validated_number(id, 0, "candidate label id", 1, nil, true)
  local label = labels[id]
  if not label then
    warn(string.format("ignored candidate for unknown label #%s", tostring(id)))
    return
  end
  local candidate_x = validated_number(x, 0.0, "candidate x coordinate")
  local candidate_y = validated_number(y, 0.0, "candidate y coordinate")

  label.candidates[#label.candidates + 1] = {
    pos = validated_number(pos, 0.5, "candidate position", 0, 1),
    base_x = candidate_x,
    base_y = candidate_y,
    x = candidate_x,
    y = candidate_y,
    component = component == nil and 0
      or validated_number(component, 0, "candidate component", 0, nil, true),
    valid = true,
    rejection_reason = nil,
  }
end

local candidate_positions_uniform
local candidate_positions_around_preferred
local candidate_positions_adaptive

--- Begin a new disconnected component of the surveyed PGFPlots path.
-- No candidate interval may cross this break.
function M.add_path_break(id)
  local label = labels[validated_number(id, 0, "path label id", 1, nil, true)]
  if label then
    label.path_points[#label.path_points + 1] = {break_before = true}
  end
end

--- Record one surveyed path vertex and PGFPlots' preceding segment length.
-- The latter is the metric used by \pgfplotspointplotattime, which preserves
-- the documented normalized path-position convention.
function M.add_path_point(id, x, y, segment_length)
  local label = labels[validated_number(id, 0, "path label id", 1, nil, true)]
  if not label then return end
  local points = label.path_points
  local previous = points[#points]
  local length = validated_number(segment_length, 0, "path segment length", 0)
  if not previous or previous.break_before then length = 0 end
  label.path_length = (label.path_length or 0) + length
  points[#points + 1] = {
    x = validated_number(x, 0, "path x coordinate"),
    y = validated_number(y, 0, "path y coordinate"),
    distance = label.path_length,
  }
end

local function clip_segment_to_axis(a, b)
  if config.allow_outside then return 0, 1 end
  if not axis_rect then return nil end
  local dx, dy = b.x - a.x, b.y - a.y
  local start, finish = 0, 1
  local function cut(p, q)
    if math.abs(p) < 1e-12 then return q >= 0 end
    local fraction = q / p
    if p < 0 then
      if fraction > finish then return false end
      start = math.max(start, fraction)
    else
      if fraction < start then return false end
      finish = math.min(finish, fraction)
    end
    return true
  end
  if not cut(-dx, a.x - axis_rect.left)
      or not cut(dx, axis_rect.right - a.x)
      or not cut(-dy, a.y - axis_rect.bottom)
      or not cut(dy, axis_rect.top - a.y)
      or finish <= start then
    return nil
  end
  return start, finish
end

local function visible_path_intervals(label, min_pos, max_pos)
  local intervals = {}
  local total = label.path_length or 0
  if total <= 0 then return intervals end
  local previous
  for _, point in ipairs(label.path_points) do
    if point.break_before then
      previous = nil
    else
      if not previous then point.continues_visible = false end
      if previous and point.distance > previous.distance then
        local first, last = clip_segment_to_axis(previous, point)
        if first then
          local delta = point.distance - previous.distance
          local lo = math.max(min_pos, (previous.distance + first * delta) / total)
          local hi = math.min(max_pos, (previous.distance + last * delta) / total)
          if hi > lo + 1e-12 then
            local prior = intervals[#intervals]
            if prior and previous.continues_visible
                and lo <= prior[2] + 1e-10 then
              prior[2] = math.max(prior[2], hi)
            else
              intervals[#intervals + 1] = {lo, hi}
            end
            point.continues_visible = true
          else
            point.continues_visible = false
          end
        else
          point.continues_visible = false
        end
      end
      previous = point
    end
  end
  return intervals
end

local function positions_on_visible_path(intervals, preferred, n, strategy)
  local total = 0
  for _, interval in ipairs(intervals) do
    total = total + interval[2] - interval[1]
  end
  if total <= 0 then return {} end

  local preferred_distance, elapsed = 0, 0
  local best_distance = math.huge
  for _, interval in ipairs(intervals) do
    local projected = clamp(preferred, interval[1], interval[2])
    local distance = math.abs(projected - preferred)
    if distance < best_distance then
      best_distance = distance
      preferred_distance = elapsed + projected - interval[1]
    end
    elapsed = elapsed + interval[2] - interval[1]
  end
  local visible_preferred = preferred_distance / total
  local fractions
  if strategy == "uniform" then
    fractions = candidate_positions_uniform(visible_preferred, 0, 1, n)
  elseif strategy == "adaptive" then
    fractions = candidate_positions_adaptive(visible_preferred, 0, 1, n)
  else
    fractions = candidate_positions_around_preferred(visible_preferred, 0, 1, n)
  end

  local positions = {}
  for _, fraction in ipairs(fractions) do
    local target = fraction * total
    local offset = 0
    for index, interval in ipairs(intervals) do
      local span = interval[2] - interval[1]
      if target <= offset + span or index == #intervals then
        positions[#positions + 1] = {
          pos = interval[1] + clamp(target - offset, 0, span),
          component = index,
        }
        break
      end
      offset = offset + span
    end
  end
  return positions
end

candidate_positions_uniform = function(preferred, min_pos, max_pos, n)
  local positions = {}
  if n <= 1 then
    return { clamp(preferred, min_pos, max_pos) }
  end
  for i = 0, n - 1 do
    positions[#positions + 1] = min_pos + (max_pos - min_pos) * i / (n - 1)
  end
  return positions
end

candidate_positions_around_preferred = function(preferred, min_pos, max_pos, n)
  local positions = {}
  local seen = {}
  local function add(p)
    p = clamp(p, min_pos, max_pos)
    local key = string.format("%.8f", p)
    if not seen[key] and #positions < n then
      positions[#positions + 1] = p
      seen[key] = true
    end
  end

  preferred = clamp(preferred, min_pos, max_pos)
  add(preferred)
  if n <= 1 then return positions end

  local step = (max_pos - min_pos) / math.max(1, n - 1)
  for k = 1, 4 * n do
    add(preferred - k * step)
    add(preferred + k * step)
    if #positions >= n then break end
  end
  return positions
end

candidate_positions_adaptive = function(preferred, min_pos, max_pos, n)
  -- The adaptive strategy keeps a dense cluster around the preferred position
  -- and a coarse coverage of the whole admissible interval.  It is still a
  -- deterministic finite set, which is important for reproducible TeX builds.
  if n <= 3 then
    return candidate_positions_around_preferred(preferred, min_pos, max_pos, n)
  end
  local positions = {}
  local seen = {}
  local function add(p)
    p = clamp(p, min_pos, max_pos)
    local key = string.format("%.8f", p)
    if not seen[key] and #positions < n then
      positions[#positions + 1] = p
      seen[key] = true
    end
  end

  preferred = clamp(preferred, min_pos, max_pos)
  add(preferred)
  local local_count = math.floor(0.6 * n)
  local global_count = n - local_count
  local local_span = 0.25 * (max_pos - min_pos)
  local local_min = clamp(preferred - local_span, min_pos, max_pos)
  local local_max = clamp(preferred + local_span, min_pos, max_pos)

  for _, p in ipairs(candidate_positions_around_preferred(preferred, local_min, local_max, local_count)) do add(p) end
  for _, p in ipairs(candidate_positions_uniform(preferred, min_pos, max_pos, global_count + 2)) do add(p) end
  for _, p in ipairs(candidate_positions_around_preferred(preferred, min_pos, max_pos, n)) do add(p) end
  return positions
end

--- Emit TeX callbacks that sample candidate positions on the current path.
-- @tparam number id Label identifier.
-- @tparam number preferred Preferred normalized path position.
-- @tparam number min_pos Smallest admissible normalized position.
-- @tparam number max_pos Largest admissible normalized position.
-- @tparam number n Requested number of candidates.
-- @tparam string strategy `uniform`, `around-preferred`, or `adaptive`.
function M.emit_candidate_positions(id, preferred, min_pos, max_pos, n, strategy)
  id = validated_number(id, 0, "candidate label id", 1, nil, true)
  preferred = validated_number(preferred, 0.5,
    "preferred position", 0, 1)
  min_pos = validated_number(min_pos, 0.05, "minimum position", 0, 1)
  max_pos = validated_number(max_pos, 0.95, "maximum position", 0, 1)
  if max_pos < min_pos then
    min_pos, max_pos = max_pos, min_pos
  end
  n = validated_number(n, 31, "candidate count", 2, nil, true)
  strategy = strategy or "around-preferred"

  local positions
  if strategy == "uniform" then
    positions = candidate_positions_uniform(preferred, min_pos, max_pos, n)
  elseif strategy == "adaptive" then
    positions = candidate_positions_adaptive(preferred, min_pos, max_pos, n)
  else
    if strategy ~= "around-preferred" and strategy ~= "around preferred" then
      warn("unknown candidate strategy `" .. tostring(strategy)
        .. "'; using around-preferred")
    end
    positions = candidate_positions_around_preferred(preferred, min_pos, max_pos, n)
  end

  log(2, string.format("emit %d candidate sample(s) for label #%d using strategy=%s", #positions, id, strategy))
  for _, pos in ipairs(positions) do
    tex.sprint(string.format(
      "\\csname pgfplotsautonode@samplecandidate\\endcsname{%s}{0}",
      tex_number(pos)))
  end
end

--- Sample only visible portions of the surveyed plot polyline.
-- The path's normalized positions retain PGFPlots' full-path meaning.  If the
-- surveyed stream is unavailable, fall back to the original point-at-time
-- sampler; the hard axis-border check still applies to every candidate.
function M.emit_visible_candidate_positions(id, preferred, min_pos, max_pos,
    n, strategy)
  id = validated_number(id, 0, "candidate label id", 1, nil, true)
  preferred = validated_number(preferred, 0.5, "preferred position", 0, 1)
  min_pos = validated_number(min_pos, 0, "minimum position", 0, 1)
  max_pos = validated_number(max_pos, 1, "maximum position", 0, 1)
  if max_pos < min_pos then min_pos, max_pos = max_pos, min_pos end
  n = validated_number(n, 51, "candidate count", 2, nil, true)
  strategy = strategy or "around-preferred"
  if strategy == "around preferred" then strategy = "around-preferred" end
  if strategy ~= "uniform" and strategy ~= "adaptive"
      and strategy ~= "around-preferred" then
    warn("unknown candidate strategy `" .. tostring(strategy)
      .. "'; using around-preferred")
    strategy = "around-preferred"
  end

  local label = labels[id]
  if not label then return end
  if not label.path_length or label.path_length <= 0 then
    warn(string.format(
      "label #%d has no surveyed path metric; using full-path candidate sampling",
      id))
    M.emit_candidate_positions(id, preferred, min_pos, max_pos, n, strategy)
    return
  end
  local intervals, positions
  if max_pos == min_pos then
    -- A fixed label has one admissible path position, not a zero-length
    -- interval to be discarded by the segment clipper.
    intervals = visible_path_intervals(label, 0, 1)
    positions = {}
    for index, interval in ipairs(intervals) do
      if min_pos >= interval[1] - 1e-10
          and min_pos <= interval[2] + 1e-10 then
        positions[1] = {pos = min_pos, component = index}
        break
      end
    end
  else
    intervals = visible_path_intervals(label, min_pos, max_pos)
    positions = positions_on_visible_path(intervals, preferred, n, strategy)
  end
  log(1, string.format(
    "label #%d: %d visible path component(s), %d candidate sample(s)",
    id, #intervals, #positions))
  for _, item in ipairs(positions) do
    -- PGFPlots' point-at-time interpolator can fail at exactly 0 or 1 on a
    -- disconnected path. The surveyed endpoint is already in canvas
    -- coordinates, so use it directly without changing the requested position.
    local endpoint
    if item.pos <= 1e-12 then
      for _, point in ipairs(label.path_points) do
        if not point.break_before then
          endpoint = point
          break
        end
      end
    elseif item.pos >= 1 - 1e-12 then
      for index = #label.path_points, 1, -1 do
        local point = label.path_points[index]
        if not point.break_before then
          endpoint = point
          break
        end
      end
    end
    if endpoint then
      M.add_candidate(id, item.pos, endpoint.x, endpoint.y, item.component)
    else
      tex.sprint(string.format(
        "\\csname pgfplotsautonode@samplecandidate\\endcsname{%s}{%d}",
        tex_number(item.pos), item.component))
    end
  end
end

--- Solve the current label assignment without emitting drawing commands.
-- @treturn table Labels ordered by identifier, including their chosen
-- candidates and visibility state.
function M.solve()
  local list = ordered_labels()
  log(1, string.format("solve: %d label(s), algorithm=%s", #list, config.algorithm))

  if #list > 0 and not axis_rect and not config.allow_outside then
    package_error("visible axis rectangle is unavailable; automatic labels cannot be placed")
  end

  -- A label with no inside candidate is not part of the assignment problem.
  -- Overlap policies may relax label-to-label collisions, never the axis border.
  local active = prepare_all_labels(list)

  if config.algorithm == "greedy" then
    solve_greedy(active)
  elseif config.algorithm == "local-search" then
    solve_local_search(active)
  elseif config.algorithm == "exact-small" then
    exact_search(active)
  else
    solve_repair(active)
  end

  enforce_visibility(list)
  local status = handle_final_failures(list)
  log(1, string.format(
    "final status: overlaps=%d outside=%d hidden=%d",
    status.overlaps, status.outside, status.hidden
  ))

  for _, label in ipairs(ordered_labels_by_id()) do
    if label.hidden then
      log(1, string.format("final label #%d: hidden", label.id))
    elseif label.choice then
      log(1, string.format(
        "final label #%d: pos=%s center=(%s,%s) angle=%s",
        label.id,
        fmt(label.choice.pos),
        fmt(label.choice.x),
        fmt(label.choice.y),
        fmt(label.choice.angle)
      ))
    end
  end

  return ordered_labels_by_id()
end

local function emit_candidate_debug(list)
  if not config.show_candidates and not config.show_rejected_candidates then
    return
  end
  for _, label in ipairs(list) do
    for _, candidate in ipairs(label.candidates or {}) do
      if candidate.valid or config.show_rejected_candidates then
        local state = candidate.valid and "valid" or "rejected"
        tex.sprint(string.format(
          "\\pgfplotsautonodedrawcandidate{%s}{%s}{%s}",
          state,
          tex_number(candidate.x or candidate.base_x),
          tex_number(candidate.y or candidate.base_y)
        ))
      end
    end
  end
end

local function emit_box_debug(list)
  if not config.show_bounding_boxes then return end
  for _, label in ipairs(list) do
    if label.choice and not label.hidden then
      local box = label.choice.box.aabb
      tex.sprint(string.format(
        "\\pgfplotsautonodedrawbbox{%d}{%s}{%s}{%s}{%s}",
        label.id,
        tex_number(box.xmin),
        tex_number(box.ymin),
        tex_number(box.xmax),
        tex_number(box.ymax)
      ))
    end
  end
end

--- Solve the current label assignment and emit final TeX drawing callbacks.
-- Debug candidates and bounding boxes are emitted when enabled. Axis-local
-- state is reset after emission.
function M.solve_and_emit()
  local list = M.solve()

  emit_candidate_debug(list)

  for _, label in ipairs(list) do
    local choice = label.choice
    if choice and not label.hidden then
      local angle = label.sloped and choice.angle or 0.0
      tex.sprint(string.format(
        "\\pgfplotsautonodedraw{%d}{%s}{%s}{%s}",
        label.id,
        tex_number(choice.x),
        tex_number(choice.y),
        tex_number(angle)
      ))
    end
  end

  emit_box_debug(list)
  M.reset()
end

return M