Tier 3: Production-Scale Algorithms

UCB(a, t) = f(context, a; theta) + beta * sqrt(g^T * Sigma_inv * g)

sample ~ Normal(f(context, a; theta), nu^2 * g^T * Sigma_inv * g)

class LastLayerBayesianBandit
  def initialize(num_arms, hidden_dim:, learning_rate: 0.001, retrain_every: 100)
    @network = NeuralNetwork.new(hidden_dim: hidden_dim)
    @retrain_every = retrain_every
    @replay_buffer = []
    @step = 0

    # Bayesian linear regression on last-layer features
    # One set of parameters per arm
    @precision = Array.new(num_arms) do
      IdentityMatrix.new(hidden_dim)  # Sigma_inv = I (prior)
    end
    @theta = Array.new(num_arms) do
      ZeroVector.new(hidden_dim)
    end
    @b = Array.new(num_arms) do
      ZeroVector.new(hidden_dim)
    end
  end

  def select_arm(context)
    features = @network.extract_last_layer(context)

    samples = @theta.each_index.map do |a|
      # Posterior mean and variance from Bayesian linear regression
      mu = features.dot(@theta[a])
      sigma = Math.sqrt(features.dot(@precision[a].solve(features)))

      # Thompson sample
      NormalDistribution.sample(mu, sigma)
    end

    max_sample = samples.max
    best_arms = samples.each_index.select do |i|
      samples[i] == max_sample
    end
    best_arms.sample
  end

  def update(arm, reward, context)
    @replay_buffer.push({ context: context, arm: arm, reward: reward })
    @step += 1

    # Update Bayesian linear regression (cheap, per-step)
    features = @network.extract_last_layer(context)
    # Sherman-Morrison rank-1 update for Sigma_inv
    @precision[arm] = sherman_morrison_update(@precision[arm], features)
    @b[arm] = @b[arm] + features * reward
    @theta[arm] = @precision[arm].solve(@b[arm])

    # Periodic neural network retraining (expensive, batched)
    if @step % @retrain_every == 0
      retrain_network
      reset_bayesian_heads  # re-extract features, refit linear heads
    end
  end

  private

  def sherman_morrison_update(sigma_inv, x)
    # O(d^2) update instead of O(d^3) inversion
    # (A + xx^T)^{-1} = A^{-1} - A^{-1}xx^T A^{-1} / (1 + x^T A^{-1} x)
    a_inv_x = sigma_inv.dot(x)
    sigma_inv - a_inv_x.outer(a_inv_x) / (1.0 + x.dot(a_inv_x))
  end

  def retrain_network
    @network.train_on(@replay_buffer)
  end

  def reset_bayesian_heads
    # Re-extract features with updated network, refit linear regression
    @precision.each_index do |a|
      @precision[a] = IdentityMatrix.new(@precision[a].size)
      @b[a] = ZeroVector.new(@b[a].size)
    end
    @replay_buffer.each do |sample|
      features = @network.extract_last_layer(sample[:context])
      if sample[:arm] == a
        @precision[a] = sherman_morrison_update(@precision[a], features)
        @b[a] = @b[a] + features * sample[:reward]
      end
    end
    @theta.each_index do |a|
      @theta[a] = @precision[a].solve(@b[a])
    end
  end
end

# 1. Online per-step (research only — can be unstable)
network.train_step(context, reward)

# 2. Periodic batch retraining every L steps (avoids catastrophic forgetting)
if step % retrain_interval == 0
  network.train_on(replay_buffer)
end

# 3. Replay buffer + mini-batch SGD (most production-oriented, used by AgileRL)
replay_buffer.push(transition)
if replay_buffer.size >= batch_size
  batch = replay_buffer.sample(batch_size)
  network.train_on_batch(batch)
end

# Semi-bandit feedback: observe reward for each selected arm (preferred)
selected_arms = [3, 7, 12]
rewards = { 3 => 0.8, 7 => 0.2, 12 => 0.6 }  # one reward per arm

# Full-bandit feedback: observe only aggregate reward (harder)
selected_arms = [3, 7, 12]
aggregate_reward = 1.6  # no per-arm decomposition

# Top-K (m-set): select K arms with highest estimates
# O(n log n) — sort and take top K
def top_k_oracle(estimates, k)
  estimates.each_with_index
    .sort_by { |value, _index| -value }
    .first(k)
    .map { |_value, index| index }
end

# Matching / assignment: bipartite matching
# O(n^3) — Hungarian algorithm
def matching_oracle(estimates_matrix)
  HungarianAlgorithm.solve(estimates_matrix)
end

# Submodular maximization: greedy with (1-1/e) approximation guarantee
def submodular_oracle(ground_set, k, submodular_fn)
  selected = []
  k.times do
    best = ground_set.max_by do |item|
      submodular_fn.marginal_gain(selected, item)
    end
    selected.push(best)
    ground_set.delete(best)
  end
  selected
end

class CombinatorialUCB
  def initialize(num_base_arms, oracle:, alpha: 2.0)
    @oracle = oracle
    @alpha = alpha
    @counts = Array.new(num_base_arms, 0)
    @values = Array.new(num_base_arms, 0.0)
    @total_rounds = 0
  end

  def select_arms
    @total_rounds += 1

    # Compute UCB for each base arm
    ucb_values = @counts.each_index.map do |i|
      if @counts[i] == 0
        Float::INFINITY  # ensure unplayed arms are selected
      else
        exploration_bonus = Math.sqrt(@alpha * Math.log(@total_rounds) / @counts[i])
        @values[i] + exploration_bonus
      end
    end

    # Feed UCB estimates to the combinatorial oracle
    @oracle.solve(ucb_values)
  end

  def update(arms, rewards)
    # Semi-bandit feedback: update each selected arm individually
    arms.each_with_index do |arm, i|
      @counts[arm] += 1
      @values[arm] += (rewards[i] - @values[arm]) / @counts[arm]
    end
  end
end

class SlidingWindowUCB
  def initialize(num_arms, window_size:, alpha: 2.0)
    @window_size = window_size
    @alpha = alpha
    @num_arms = num_arms
    @history = []  # list of { arm:, reward:, time: }
    @total_rounds = 0
  end

  def select_arm
    @total_rounds += 1

    # Get recent history within the window
    recent = @history.last(@window_size)

    ucb_values = (0...@num_arms).map do |arm|
      arm_history = recent.select { |h| h[:arm] == arm }
      if arm_history.empty?
        Float::INFINITY  # unplayed in window — explore
      else
        mean = arm_history.sum { |h| h[:reward] } / arm_history.size.to_f
        exploration = Math.sqrt(@alpha * Math.log(@total_rounds) / arm_history.size)
        mean + exploration
      end
    end

    max_ucb = ucb_values.max
    best_arms = ucb_values.each_index.select { |i| ucb_values[i] == max_ucb }
    best_arms.sample
  end

  def update(arm, reward)
    @history.push({ arm: arm, reward: reward, time: @total_rounds })
  end
end

# Window size heuristic (Garivier & Moulines):
# tau = 2 * sqrt(T * log(T) / (1 + num_change_points))

class DiscountedUCB
  def initialize(num_arms, gamma: 0.99, alpha: 2.0)
    @gamma = gamma
    @alpha = alpha
    @num_arms = num_arms
    @discounted_counts = Array.new(num_arms, 0.0)
    @discounted_sums   = Array.new(num_arms, 0.0)
    @total_discounted   = 0.0
  end

  def select_arm
    ucb_values = (0...@num_arms).map do |arm|
      if @discounted_counts[arm] < 1e-10
        Float::INFINITY
      else
        mean = @discounted_sums[arm] / @discounted_counts[arm]
        exploration = Math.sqrt(@alpha * Math.log(@total_discounted) / @discounted_counts[arm])
        mean + exploration
      end
    end

    max_ucb = ucb_values.max
    best_arms = ucb_values.each_index.select { |i| ucb_values[i] == max_ucb }
    best_arms.sample
  end

  def update(arm, reward)
    # Discount all existing counts and sums
    @discounted_counts.each_index do |i|
      @discounted_counts[i] *= @gamma
      @discounted_sums[i]   *= @gamma
    end
    @total_discounted = @total_discounted * @gamma + 1.0

    # Add new observation
    @discounted_counts[arm] += 1.0
    @discounted_sums[arm]   += reward
  end
end

# gamma typically 0.95 to 0.99
# Lower gamma = faster forgetting = better for rapidly changing environments

class CusumUCB
  def initialize(num_arms, base_algorithm:, threshold: 10.0, window: 50)
    @base = base_algorithm
    @threshold = threshold
    @window = window
    @reward_history = Array.new(num_arms) { [] }
  end

  def select_arm
    @base.select_arm
  end

  def update(arm, reward)
    @base.update(arm, reward)
    @reward_history[arm].push(reward)

    # CUSUM change-point detection
    if detect_change(arm)
      reset_arm(arm)
    end
  end

  private

  def detect_change(arm)
    history = @reward_history[arm]
    return false if history.size < @window

    recent = history.last(@window)
    reference_mean = history[0...-@window].sum / (history.size - @window).to_f

    # Cumulative sum of deviations from reference mean
    cusum_pos = 0.0
    cusum_neg = 0.0
    recent.each do |r|
      cusum_pos = [0.0, cusum_pos + (r - reference_mean)].max
      cusum_neg = [0.0, cusum_neg - (r - reference_mean)].max
      return true if cusum_pos > @threshold or cusum_neg > @threshold
    end

    false
  end

  def reset_arm(arm)
    @reward_history[arm] = []
    @base.reset(arm)  # base algorithm clears statistics for this arm
  end
end

class WhittleIndexPolicy
  def initialize(arms, budget:)
    @arms = arms       # each arm has transition_active, transition_passive matrices
    @budget = budget
  end

  def allocate_budget(states)
    # Compute Whittle index for each arm given its current state
    indices = @arms.each_with_index.map do |arm, i|
      compute_whittle_index(arm, states[i])
    end

    # Activate the B arms with highest indices
    ranked = indices.each_with_index
      .sort_by { |index_value, _arm_id| -index_value }

    active_arms = ranked.first(@budget).map { |_value, arm_id| arm_id }
    active_arms
  end

  def update(new_states)
    # States evolve according to transition matrices
    # The agent observes new states after activation decisions
    # No explicit update needed if transitions are known
  end

  private

  def compute_whittle_index(arm, state)
    # Method depends on problem structure:

    # 1. Closed-form (only for simple 1D-state models)
    #    e.g., for two-state arms: analytic formula from transition probabilities

    # 2. Adaptive-greedy (Nino-Mora): O(K^3) for K states
    #    Iteratively computes indices by solving a sequence of LP relaxations

    # 3. Value iteration: solve the single-arm problem for a range of subsidies
    #    Binary search over subsidy lambda until indifference condition is met
    binary_search_whittle_index(arm, state)
  end

  def binary_search_whittle_index(arm, state)
    lo = -1.0
    hi = 1.0

    100.times do
      mid = (lo + hi) / 2.0
      # Solve single-arm problem: is it better to activate or passivate
      # when receiving subsidy mid for passivating?
      active_value = arm.reward(state) + arm.transition_active[state].expected_future_value
      passive_value = mid + arm.transition_passive[state].expected_future_value

      if active_value > passive_value
        lo = mid   # need higher subsidy for indifference
      else
        hi = mid
      end
    end

    (lo + hi) / 2.0
  end
end

class UCWhittle
  def initialize(num_arms, num_states:, budget:)
    @budget = budget
    # Track transition counts per arm, per state, per action (active/passive)
    @transition_counts = Array.new(num_arms) do
      Array.new(num_states) do
        { active: Array.new(num_states, 0), passive: Array.new(num_states, 0) }
      end
    end
  end

  def allocate_budget(states)
    # Estimate transition probabilities with optimistic confidence bounds
    optimistic_arms = @transition_counts.each_with_index.map do |arm_counts, i|
      build_optimistic_transitions(arm_counts, states[i])
    end

    # Compute Whittle indices using optimistic transition estimates
    indices = optimistic_arms.each_with_index.map do |arm, i|
      compute_whittle_index(arm, states[i])
    end

    ranked = indices.each_with_index
      .sort_by { |index_value, _arm_id| -index_value }
    ranked.first(@budget).map { |_value, arm_id| arm_id }
  end

  def update(states, actions, new_states)
    states.each_with_index do |state, i|
      action = actions.include?(i) ? :active : :passive
      @transition_counts[i][state][action][new_states[i]] += 1
    end
  end
end

class CascadeUCB1
  def initialize(num_items, slate_size:, alpha: 2.0)
    @num_items = num_items
    @slate_size = slate_size
    @alpha = alpha
    @counts = Array.new(num_items, 0)
    @values = Array.new(num_items, 0.0)  # estimated attraction probability
    @total_rounds = 0
  end

  def select_ranking
    @total_rounds += 1

    # Compute UCB for each item's attraction probability
    ucb_values = (0...@num_items).map do |item|
      if @counts[item] == 0
        Float::INFINITY
      else
        exploration = Math.sqrt(@alpha * Math.log(@total_rounds) / @counts[item])
        @values[item] + exploration
      end
    end

    # Rank by UCB score, return top slate_size items
    ucb_values.each_with_index
      .sort_by { |value, _item| -value }
      .first(@slate_size)
      .map { |_value, item| item }
  end

  def update(ranking, click_position)
    # click_position = index within ranking where user clicked
    # nil if user examined all items without clicking

    if click_position.nil?
      # User examined all items, clicked none — all are negative feedback
      ranking.each do |item|
        @counts[item] += 1
        @values[item] += (0.0 - @values[item]) / @counts[item]
      end
    else
      # Items BEFORE click: examined but not attractive (reward = 0)
      ranking[0...click_position].each do |item|
        @counts[item] += 1
        @values[item] += (0.0 - @values[item]) / @counts[item]
      end

      # Clicked item: attractive (reward = 1)
      clicked_item = ranking[click_position]
      @counts[clicked_item] += 1
      @values[clicked_item] += (1.0 - @values[clicked_item]) / @counts[clicked_item]

      # Items AFTER click: NOT updated — user never saw them
      # This is the key property of the cascade model
    end
  end
end

# Replace the UCB computation with KL-UCB:
def kl_ucb(mean, count, total_rounds)
  # Find the largest q such that KL(mean, q) <= log(total_rounds) / count
  # KL(p, q) = p * log(p/q) + (1-p) * log((1-p)/(1-q))  (Bernoulli KL)
  threshold = Math.log(total_rounds) / count
  binary_search_kl_upper_bound(mean, threshold)
end