# =====================================================================
# inventory_voi.jl
#
# The Value of a Forecast: Inventory Management under Demand Uncertainty
# with an Imperfect (Machine-Learning) Demand Signal.
#
# Companion script for the "Value of Information" technical showcase.
# Dynamic Optimization Series.
#
# ---------------------------------------------------------------------
# Model
# ---------------------------------------------------------------------
# A monopolist lives forever, discounting the future at rate delta.
# Each period t:
#   (1) starts with inventory  S_t  (the continuous STATE),
#   (2) observes an ML signal  s_hat in {H, L}  about demand,
#   (3) chooses production  q_t >= 0  at constant unit cost c,
#   (4) the true multiplier  theta_t in {theta_H, theta_L}  is realised
#       (i.i.d. across periods),
#   (5) chooses sales  y_t  with  0 <= y_t <= S_t + q_t,
#   (6) carries leftover  S_{t+1} = S_t + q_t - y_t,  paying a CONVEX
#       holding cost (h/2) S^2.
#
#   inverse demand : p(y, theta) = a*theta - b*y
#   production cost: c * q                 (constant MC: CRS Cobb-Douglas)
#   holding cost   : (h/2) * S^2           (convex -> inventory matters)
#
# Information regimes:
#   :noinfo   q chosen under the prior  pi_H
#   :signal   q chosen under the Bayesian posterior implied by a
#             classifier with TPR = alpha, TNR = beta
#   :perfect  q chosen knowing theta exactly  (upper bound)
#
#   VOI_signal  = V_signal  - V_noinfo
#   VOI_perfect = V_perfect - V_noinfo
#   efficiency  = VOI_signal / VOI_perfect  in [0, 1]
#
# ---------------------------------------------------------------------
# Numerical method: continuous-state dynamic programming
# ---------------------------------------------------------------------
# Standard approach for a continuous state:
#   * the value function V is stored on a grid over inventory S;
#   * V(S') at off-grid points is recovered by LINEAR INTERPOLATION
#     (Interpolations.jl), with flat extrapolation at the boundaries;
#   * production q and sales y are CONTINUOUS choices, solved by bounded
#     scalar optimization (Optim.jl, Brent).
#
# Within one Bellman sweep, the sales subproblem
#   W(A, theta) = max_{0<=y<=A} (a*theta - b*y) y + delta * V(A - y)
# is solved DIRECTLY at A = S + q whenever the production search needs it.
# There is no grid over available stock A and no W cache: A is only an
# intermediate quantity (S + q), never a state variable. The single state
# is inventory S. (Caching W on an A-grid would trade a little accuracy for
# speed; we keep the direct call for simplicity and exactness.)
#
# Requires: Interpolations, Optim, Printf, Plots
# Run:      include("inventory_voi.jl")
# =====================================================================

# =====================================================================
# PACKAGES
# =====================================================================
using Interpolations   # builds the linear interpolant of the value function V(S)
using Optim            # 1-D bounded optimizer (Brent) for the continuous q and y choices
using Printf           # @printf for clean, aligned console output
using Plots            # savefig / plot / heatmap for the two output figures

# =====================================================================
# 1. PARAMETERS
#
# A single immutable struct holds every primitive. `Base.@kwdef` lets us
# construct it with keyword defaults, e.g. `Params(delta = 0.97)` to
# override just one field while keeping the rest.
# =====================================================================
Base.@kwdef struct Params
    # --- demand side ---
    a::Float64       = 10.0     # demand intercept scale: p = a*theta - b*y
    b::Float64       = 1.0      # demand slope (price sensitivity)
    theta_H::Float64 = 1.3      # demand multiplier in the GOOD state (+30%)
    theta_L::Float64 = 0.8      # demand multiplier in the BAD  state (-20%)
    pi_H::Float64    = 0.5      # prior probability of the good state (50/50) -- Well, just my simple choice 

    # --- cost side ---
    c::Float64       = 2.0      # constant unit production cost (marginal cost)
    h::Float64       = 0.5      # curvature of the convex holding cost (h/2)*S^2

    # --- dynamics ---
    delta::Float64   = 0.95     # per-period discount factor (must be < 1 for convergence)

    # --- numerical grid ---
    Smax::Float64    = 12.0     # largest inventory level the STATE grid covers
    nS::Int          = 61       # number of points on the inventory (state) grid

    # --- value function iteration controls ---
    tol::Float64     = 1e-8     # convergence threshold on the sup-norm of V updates
    maxiter::Int     = 5_000    # safety cap on the number of VFI sweeps
end

# Inventory (state) grid: nS evenly spaced points from 0 to Smax.
state_grid(P::Params) = range(0.0, P.Smax, length = P.nS)

# Upper bound on production q. We never need available stock S + q to exceed
# what could possibly be sold in one period (demand intercept over slope in
# the good state), so this caps the q-search range.
qmax(P::Params) = P.Smax + P.a * P.theta_H / P.b

# =====================================================================
# 2. CONFUSION MATRIX  ->  BAYESIAN POSTERIOR
#
# The ML classifier is summarised by two rates:
#   alpha = TPR = Pr(signal = H | true state = H)   "true positive rate"
#   beta  = TNR = Pr(signal = L | true state = L)   "true negative rate"
#
# The firm acts on the POSTERIOR: given the signal it just saw, how likely
# is the good state? Bayes' rule converts (alpha, beta) + prior into that.
#
# Returns a 3-tuple:
#   piH_sH = Pr(H | signal = H)   belief after an optimistic signal
#   piH_sL = Pr(H | signal = L)   belief after a pessimistic signal
#   pr_sH  = Pr(signal = H)       how often the optimistic signal fires
# =====================================================================
function posterior(P::Params, alpha::Float64, beta::Float64)
    piH, piL = P.pi_H, 1.0 - P.pi_H          # prior probabilities of H and L

    # Joint probabilities Pr(true state, observed signal):
    pH_sH = alpha       * piH                # state H AND signal H  (correct hit)
    pL_sH = (1 - beta)  * piL                # state L AND signal H  (false positive)
    pH_sL = (1 - alpha) * piH                # state H AND signal L  (false negative)
    pL_sL = beta        * piL                # state L AND signal L  (correct rejection)

    # Marginal probabilities of each signal (denominators in Bayes' rule):
    pr_sH = pH_sH + pL_sH                     # Pr(signal = H)
    pr_sL = pH_sL + pL_sL                     # Pr(signal = L)

    # Posteriors = joint / marginal. The guards avoid 0/0 at degenerate
    # classifiers (e.g. one that never emits a given signal); there we
    # simply fall back to the prior, which is the correct limiting belief.
    piH_sH = pr_sH > 0 ? pH_sH / pr_sH : piH
    piH_sL = pr_sL > 0 ? pH_sL / pr_sL : piH

    return piH_sH, piH_sL, pr_sH
end

# =====================================================================
# 3. LINEAR INTERPOLANT OF THE VALUE FUNCTION
#
# V is known only at the grid points Sg. To evaluate V(S') at an arbitrary
# (off-grid) S' inside the optimizer, we wrap it in a linear interpolant.
# `Flat()` extrapolation clamps to the boundary value outside [0, Smax],
# which keeps the optimizer well-behaved if it probes slightly past an edge.
# =====================================================================
function make_V_interp(P::Params, Sg, V::Vector{Float64})
    itp = interpolate((collect(Sg),), V, Gridded(Linear()))  # piecewise-linear V
    return extrapolate(itp, Interpolations.Flat())           # clamp outside the grid
end

# =====================================================================
# 4. SALES SUBPROBLEM (continuous y):  W(A, theta)
# Why subproblem?: Because the firm needs to decide how much to produce before the true stochastic state is realised.
# At saleing time, the firm knows the true state, but not the signal. However when production is decided, the firm knows the signal, but not the true state.
# So we need to solve the sales subproblem for both cases.
# AFTER the true state theta is realised and the firm holds available
# stock A = S + q, it picks sales y in [0, A] to solve
#
#     W(A, theta) = max_{0 <= y <= A}  (a*theta - b*y)*y  +  delta * V(A - y)
#                    \__________________/    \____________/
#                       this period's revenue   discounted value of
#                                               the leftover stock (A - y)
#
# Revenue is concave in y and the interpolated continuation value is concave
# here, so a 1-D bounded search (Brent) reliably finds the global optimum.
# The production step (Section 5) calls this DIRECTLY at A = S + q. There is
# no available-stock grid and no cache: A is just an intermediate quantity,
# never a state variable. (One could cache W on a grid of A to trade a little
# accuracy for speed, but keeping the call direct is simpler and exact.)
# =====================================================================
function sales_value(P::Params, A::Float64, theta::Float64, Vitp)
    # With no stock there is nothing to sell; value is just the discounted
    # continuation at S' = 0. (Also avoids calling optimize on an empty range.)
    A <= 0.0 && return P.delta * Vitp(0.0) #Make sure the value function is defined at 0.

    # Optim MINIMISES, so we minimise the NEGATIVE of the objective.
    neg(y) = -((P.a * theta - P.b * y) * y + P.delta * Vitp(A - y))

    res = optimize(neg, 0.0, A, Brent())     # search y over [0, A]
    return -Optim.minimum(res)               # flip the sign back to a value found as the value of the sales subproblem at the given A and theta
end

# =====================================================================
# 5. BELLMAN OPERATOR  (one full sweep T[V] -> Vnew)
#
# For each starting inventory S, the firm chooses production q >= 0 BEFORE
# the state is realised, paying c*q, to maximise
#
#     -c*q  +  E_{theta | signal}  W(S + q, theta)
#
# where the expectation uses the belief implied by the information regime.
# W (the sales subproblem, Section 4) is solved DIRECTLY at A = S + q each
# time it is needed -- there is no available-stock grid and no cache. The
# current-stock holding cost (h/2)*S^2 depends on S alone (not on q or y),
# so it is added once, outside.
## THis bellman is solved by backward induction and recursion until the value function converges.
# =====================================================================
function bellman(P::Params, Sg, V::Vector{Float64},
                 regime::Symbol, alpha::Float64, beta::Float64)

    # Build this sweep's value-function interpolant from the current V.
    Vitp = make_V_interp(P, Sg, V) #At each Sg, we have a value function V(Sg) from the previous iteration!

    # Beliefs implied by the classifier (used only in the :signal regime).
    piH_sH, piH_sL, pr_sH = posterior(P, alpha, beta)

    qcap = qmax(P)   # upper bound on the production search

    # --- production value given a belief bH on the high state ---
    # Used for :noinfo (bH = prior) and :signal (bH = posterior).
    # The sales value W is computed on the spot via sales_value(S+q, theta).
    function best_belief(S::Float64, bH::Float64)
        qhi = max(qcap - S, 1e-8)             # cap the q-search range
        neg(q) = -(-P.c * q
                   + bH       * sales_value(P, S + q, P.theta_H, Vitp)
                   + (1 - bH) * sales_value(P, S + q, P.theta_L, Vitp))
        return -Optim.minimum(optimize(neg, 0.0, qhi, Brent()))
    end

    # --- production value under PERFECT information ---
    # Here q is chosen AFTER theta is known, so we optimise q separately in
    # each state and then average with the prior (the upper bound on VOI).
    function best_perfect(S::Float64)
        qhi = max(qcap - S, 1e-8)
        negH(q) = -(-P.c * q + sales_value(P, S + q, P.theta_H, Vitp)) #why negatvie?: because we are minimizing the objective function by Optim
        negL(q) = -(-P.c * q + sales_value(P, S + q, P.theta_L, Vitp))
        vH = -Optim.minimum(optimize(negH, 0.0, qhi, Brent()))
        vL = -Optim.minimum(optimize(negL, 0.0, qhi, Brent()))
        return P.pi_H * vH + (1 - P.pi_H) * vL
    end

    # Sweep every state and apply the operator for the chosen regime.
    Vnew = Vector{Float64}(undef, P.nS)
    for (i, S) in enumerate(Sg)
        holding = 0.5 * P.h * S^2             # convex holding cost on current stock
        if regime == :noinfo
            # No signal: choose q under the prior belief.
            Vnew[i] = best_belief(S, P.pi_H) - holding
        elseif regime == :signal
            # See a signal first, THEN choose q under the matching posterior;
            # average over which signal fires, weighting by Pr(signal).
            Vnew[i] = pr_sH * best_belief(S, piH_sH) +
                      (1 - pr_sH) * best_belief(S, piH_sL) - holding
        elseif regime == :perfect
            # Know theta exactly before choosing q (the ceiling).
            Vnew[i] = best_perfect(S) - holding
        else
            error("unknown regime $regime")
        end
    end
    return Vnew
end

# =====================================================================
# 6. VALUE FUNCTION ITERATION
#
# Because delta < 1 the Bellman operator is a contraction mapping, so
# iterating T[V] from any starting guess converges to the unique fixed
# point V*. We stop when consecutive iterates differ by less than `tol`
# in the sup-norm (the largest absolute change across all grid points).
# =====================================================================
function solve_vfi(P::Params, regime::Symbol;
                   alpha::Float64 = 1.0, beta::Float64 = 1.0,
                   verbose::Bool = false)
    Sg = state_grid(P)                        # build the state grid once
    V = zeros(P.nS)                           # initial guess: V = 0 everywhere
    diff = Inf
    for it in 1:P.maxiter
        Vnew = bellman(P, Sg, V, regime, alpha, beta)       # one operator sweep
        diff = maximum(abs.(Vnew .- V))                     # sup-norm change
        V = Vnew                                            # accept the update
        if verbose && it % 25 == 0
            @printf("  [%-7s] iter %4d  sup-norm = %.3e\n", regime, it, diff)
        end
        diff < P.tol && break                  # converged: stop early
    end
    verbose && @printf("  [%-7s] done  sup-norm = %.2e\n", regime, diff)
    return Sg, V                               # return the grid and the solved V
end

# =====================================================================
# 7. VALUE OF INFORMATION FOR A GIVEN CLASSIFIER (alpha, beta)
#
# Solve all three regimes and read each value at S = 0 (an empty warehouse,
# a clean common reference point). VOI is the value the signal adds over the
# no-information benchmark; perfect information is the ceiling; efficiency is
# the share of that ceiling the classifier captures.
# =====================================================================
"""Solve V(S) on the inventory grid for all three information regimes."""
function solve_values_all_regimes(P::Params, alpha::Float64, beta::Float64)
    Sg, Vno  = solve_vfi(P, :noinfo)
    _, Vsig = solve_vfi(P, :signal; alpha = alpha, beta = beta)
    _, Vper = solve_vfi(P, :perfect)
    Sg = collect(Sg)
    return (; Sg, Vno, Vsig, Vper,
            voi_sig = Vsig .- Vno,
            voi_per = Vper .- Vno)
end

function compute_voi(P::Params, alpha::Float64, beta::Float64)
    r = solve_values_all_regimes(P, alpha, beta)
    v_no, v_sig, v_per = r.Vno[1], r.Vsig[1], r.Vper[1]
    voi_sig = v_sig - v_no
    voi_per = v_per - v_no
    eff = voi_per > 1e-9 ? voi_sig / voi_per : 0.0
    return (; v_no, v_sig, v_per, voi_sig, voi_per, eff)
end

# Parameters that widen V(S) gaps across regimes (larger demand spread, lower h).
high_voi_params() = Params(
    theta_H = 1.6, theta_L = 0.6,   # wider demand states → more to learn
    h       = 0.15,                 # cheaper to hold inventory → policies respond more
    a       = 12.0,                 # larger demand scale → higher stakes
)

# =====================================================================
# 8. PRODUCTION POLICY AT S = 0
#
# Recover the optimal production q* at an empty warehouse under three
# beliefs: after an optimistic signal, under the prior, and after a
# pessimistic signal. Shows the signal actually MOVES the decision.
# =====================================================================
function q_policy_at_zero(P::Params, alpha::Float64, beta::Float64)
    Sg = state_grid(P)
    _, V = solve_vfi(P, :signal; alpha = alpha, beta = beta)  # solve the signal regime
    Vitp = make_V_interp(P, Sg, V)                           # interpolant of the solved V
    piH_sH, piH_sL, _ = posterior(P, alpha, beta)
    qcap = qmax(P)

    # argmax of expected payoff over q at S = 0, for a given belief bH.
    # (We return the MINIMIZER of the negated objective, i.e. the optimal q.)
    # At S = 0 available stock equals q, so we evaluate sales_value at q.
    qstar(bH) = begin
        neg(q) = -(-P.c * q
                   + bH       * sales_value(P, q, P.theta_H, Vitp)
                   + (1 - bH) * sales_value(P, q, P.theta_L, Vitp))
        Optim.minimizer(optimize(neg, 0.0, qcap, Brent()))
    end

    return qstar(piH_sH), qstar(P.pi_H), qstar(piH_sL)   # (signal H, prior, signal L)
end

# =====================================================================
# 9. PLOTTING STYLE & PRODUCTION POLICIES  q*(S)
# =====================================================================
const PLOT_COLORS = (
    noinfo     = "#2F4B7C",
    signal     = "#7A5195",
    signal_h   = "#2CA02C",
    signal_l   = "#D62728",
    perfect    = "#555555",
    perfect_h  = "#FF7F0E",
    perfect_l  = "#17BECF",
    voi        = "#BC5090",
)

function apply_plot_style!(p; title = "", xlab = "", ylab = "")
    plot!(p;
          background_color      = :white,
          foreground_color      = "#2b2b2b",
          foreground_color_axis = "#2b2b2b",
          grid                  = true,
          gridcolor             = "#e6e6e6",
          gridwidth             = 1,
          framestyle            = :box,
          linewidth             = 2.4,
          linealpha             = 0.95,
          legend                = :best,
          legendfontsize        = 9,
          legendframealpha      = 0.92,
          titlefontsize         = 13,
          guidefontsize         = 11,
          tickfontsize          = 9,
          title                 = title,
          xlabel                = xlab,
          ylabel                = ylab,
          size                  = (820, 520),
          margin                = 6Plots.mm,
          dpi                   = 220)
    return p
end

function save_plot(p, path::String)
    d = dirname(path)
    !isempty(d) && mkpath(d)
    savefig(p, path)
end

"""Optimal production q*(S) for a fixed belief on the high demand state."""
function q_opt_belief(P::Params, S::Float64, bH::Float64, Vitp)
    qhi = max(qmax(P) - S, 1e-8)
    neg(q) = -(-P.c * q
               + bH       * sales_value(P, S + q, P.theta_H, Vitp)
               + (1 - bH) * sales_value(P, S + q, P.theta_L, Vitp))
    return Optim.minimizer(optimize(neg, 0.0, qhi, Brent()))
end

"""Optimal production when theta is known before choosing q (perfect info)."""
function q_opt_theta(P::Params, S::Float64, theta::Float64, Vitp)
    qhi = max(qmax(P) - S, 1e-8)
    neg(q) = -(-P.c * q + sales_value(P, S + q, theta, Vitp))
    return Optim.minimizer(optimize(neg, 0.0, qhi, Brent()))
end

"""
Extract production policies on the inventory grid for each information regime.
Returns state grid Sg and q*(S) curves (signal regime split by s=H / s=L).
"""
function extract_production_policies(P::Params, alpha::Float64, beta::Float64)
    Sg = collect(state_grid(P))
    _, Vno  = solve_vfi(P, :noinfo)
    _, Vsig = solve_vfi(P, :signal; alpha = alpha, beta = beta)
    _, Vper = solve_vfi(P, :perfect)

    Vitp_no  = make_V_interp(P, Sg, Vno)
    Vitp_sig = make_V_interp(P, Sg, Vsig)
    Vitp_per = make_V_interp(P, Sg, Vper)

    piH_sH, piH_sL, pr_sH = posterior(P, alpha, beta)

    q_noinfo = [q_opt_belief(P, S, P.pi_H, Vitp_no) for S in Sg]
    q_sig_h  = [q_opt_belief(P, S, piH_sH, Vitp_sig) for S in Sg]
    q_sig_l  = [q_opt_belief(P, S, piH_sL, Vitp_sig) for S in Sg]
    q_sig    = pr_sH .* q_sig_h .+ (1 - pr_sH) .* q_sig_l

    q_per_h  = [q_opt_theta(P, S, P.theta_H, Vitp_per) for S in Sg]
    q_per_l  = [q_opt_theta(P, S, P.theta_L, Vitp_per) for S in Sg]
    q_per    = P.pi_H .* q_per_h .+ (1 - P.pi_H) .* q_per_l

    return (; Sg, q_noinfo, q_sig, q_sig_h, q_sig_l, q_per, q_per_h, q_per_l)
end

"""Ex-ante production policies q*(S) across information regimes."""
function plot_production_policies(P::Params, alpha::Float64, beta::Float64;
                                  outfile = "policy_production_by_regime.png")
    pol = extract_production_policies(P, alpha, beta)
    Sg = pol.Sg

    fig = plot(Sg, pol.q_noinfo;
               label = "No information (prior)",
               color = PLOT_COLORS.noinfo,
               linestyle = :solid)
    plot!(fig, Sg, pol.q_sig;
          label = "Signal (ex ante)",
          color = PLOT_COLORS.signal,
          linestyle = :solid)
    plot!(fig, Sg, pol.q_per;
          label = "Perfect info (ex ante)",
          color = PLOT_COLORS.perfect,
          linestyle = :dash)
    apply_plot_style!(fig;
        title = "Production policy q*(S) by information regime",
        xlab  = "Beginning-of-period inventory  S",
        ylab  = "Optimal production  q")
    save_plot(fig, outfile)
    return fig
end

"""Value functions V(S) by regime (top) and state-wise VOI (bottom)."""
function plot_value_and_voi_by_state(P::Params, alpha::Float64, beta::Float64;
                                     outfile = "value_and_voi_by_state.png")
    r = solve_values_all_regimes(P, alpha, beta)
    Sg = r.Sg

    p_v = plot(Sg, r.Vno;
               label = "No information",
               color = PLOT_COLORS.noinfo,
               linestyle = :solid)
    plot!(p_v, Sg, r.Vsig;
          label = "Signal",
          color = PLOT_COLORS.signal,
          linestyle = :solid)
    plot!(p_v, Sg, r.Vper;
          label = "Perfect information",
          color = PLOT_COLORS.perfect,
          linestyle = :dash)
    apply_plot_style!(p_v;
        title = "Value function V(S) by information regime",
        xlab  = "Inventory  S",
        ylab  = "Value  V(S)")

    p_d = plot(Sg, r.voi_sig;
               label = "VOI from signal  (V_signal − V_no info)",
               color = PLOT_COLORS.signal,
               linestyle = :solid)
    plot!(p_d, Sg, r.voi_per;
          label = "VOI from perfect info  (V_perfect − V_no info)",
          color = PLOT_COLORS.perfect,
          linestyle = :dash)
    apply_plot_style!(p_d;
        title = "Value of information at each inventory level",
        xlab  = "Inventory  S",
        ylab  = "ΔV(S)")

    fig = plot(p_v, p_d; layout = (2, 1), size = (860, 760), link = :x)
    save_plot(fig, outfile)
    return fig, r
end

function print_voi_summary(P::Params, alpha::Float64, beta::Float64; label = "")
    r = solve_values_all_regimes(P, alpha, beta)
    !isempty(label) && println("\n--- $label ---")
    @printf("  Classifier TPR=%.2f  TNR=%.2f\n", alpha, beta)
    @printf("  V(no info) @ S=0     = %9.4f\n", r.Vno[1])
    @printf("  V(signal)  @ S=0     = %9.4f\n", r.Vsig[1])
    @printf("  V(perfect) @ S=0     = %9.4f\n", r.Vper[1])
    @printf("  VOI signal @ S=0     = %9.4f  (max over S: %.4f)\n",
            r.voi_sig[1], maximum(r.voi_sig))
    @printf("  VOI perfect @ S=0    = %9.4f  (max over S: %.4f)\n",
            r.voi_per[1], maximum(r.voi_per))
    eff = r.voi_per[1] > 1e-9 ? r.voi_sig[1] / r.voi_per[1] : 0.0
    @printf("  Efficiency @ S=0     = %8.1f%%\n", 100eff)
    return r
end

# =====================================================================
# 10. MAIN  -  print headline numbers and save figures
# =====================================================================
function main(; P::Params = high_voi_params(),
              alpha_plot::Float64 = 0.85,
              heatmap_n::Int = 15,
              make_figures::Bool = true)
    println("="^66)
    println("Inventory management & the value of an ML demand forecast")
    println("(continuous-state DP: linear interpolation + Optim)")
    println("="^66)

    # --- (A) baseline vs high-VOI parameters (80% classifier) ---
    print_voi_summary(Params(), 0.80, 0.80; label = "Default parameters")
    print_voi_summary(P, 0.80, 0.80; label = "High-VOI parameters (current run)")

    # --- (B) production policy at S = 0 for an 85% classifier ---
    qH, qN, qL = q_policy_at_zero(P, 0.85, 0.85)
    @printf("\nProduction at S=0 (85%% classifier): qL=%.2f  prior=%.2f  qH=%.2f\n",
            qL, qN, qH)

    if !make_figures
        println("\n(Skipping figures: make_figures = false)")
        return nothing
    end

    # --- (C) value functions & state-wise VOI ---
    @printf("\nPlotting V(S) and VOI(S) (classifier TPR=TNR=%.2f)...\n", alpha_plot)
    plot_value_and_voi_by_state(P, alpha_plot, alpha_plot)

    # --- (D) production policies by information regime ---
    @printf("Plotting production policies...\n")
    plot_production_policies(P, alpha_plot, alpha_plot;
                             outfile = "policy_production_by_regime.png")

    # --- (E) Figure: VOI vs symmetric accuracy (TPR = TNR) ---
    accs = collect(range(0.5, 1.0, length = 21))
    _, Vno  = solve_vfi(P, :noinfo)
    _, Vper = solve_vfi(P, :perfect)
    vno = Vno[1]; voi_p = Vper[1] - vno
    voi_sym = Float64[]
    for q in accs
        _, Vs = solve_vfi(P, :signal; alpha = q, beta = q)
        push!(voi_sym, Vs[1] - vno)
    end
    p1 = plot(accs, voi_sym;
              label = false,
              color = PLOT_COLORS.voi,
              linestyle = :solid)
    apply_plot_style!(p1;
        title = "Value of the demand forecast",
        xlab  = "Classifier accuracy  (TPR = TNR)",
        ylab  = "VOI  (V_signal − V_no info)")
    hline!(p1, [voi_p]; label = "Perfect-info ceiling",
           color = PLOT_COLORS.perfect, linestyle = :dash, linewidth = 1.8)
    save_plot(p1, "voi_vs_accuracy.png")

    # --- (F) efficiency heatmap over (TPR, TNR) ---
    gg = collect(range(0.5, 1.0, length = heatmap_n))
    Z = zeros(length(gg), length(gg))
    for (i, tpr) in enumerate(gg), (j, tnr) in enumerate(gg)
        _, Vs = solve_vfi(P, :signal; alpha = tpr, beta = tnr)
        Z[j, i] = voi_p > 1e-9 ? (Vs[1] - vno) / voi_p : 0.0
    end
    p2 = heatmap(gg, gg, Z;
                 clims = (0, 1),
                 color = :viridis,
                 colorbar_title = "Efficiency",
                 xlabel = "TPR  Pr(ŝ = H | H)",
                 ylabel = "TNR  Pr(ŝ = L | L)",
                 title  = "Share of perfect-information value captured",
                 size = (640, 560),
                 dpi = 220,
                 framestyle = :box,
                 tickfontsize = 9,
                 titlefontsize = 13,
                 guidefontsize = 11)
    save_plot(p2, "voi_efficiency_heatmap.png")

    println("\nSaved:")
    println("  value_and_voi_by_state.png")
    println("  policy_production_by_regime.png")
    println("  voi_vs_accuracy.png")
    println("  voi_efficiency_heatmap.png")
    return nothing
end

# Run when executed as a script (not when `include`d from the REPL).
if abspath(PROGRAM_FILE) == @__FILE__
    main()
end