API Reference

ContinuousDP{Tf,Tg,TR,TA}

Type representing a continuous-state dynamic program. A ContinuousDP holds the primitives of the problem; the interpolation scheme used to solve it is supplied separately through a solver object (see CollocationSolver and LQASolver).

Fields

• f::Tf: Reward function f(s, x).
• g::Tg: State transition function g(s, x, e).
• discount::Float64: Discount factor.
• shocks::TR<:AbstractVecOrMat: Discretized shock nodes.
• weights::TW: Probability weights for the shock nodes: a Vector{Float64} with one weight per shock node (length(weights) == size(shocks, 1)) for a fixed distribution, or a callable weights(s) / weights(s, x) for a state- or state-action-dependent distribution (see the constructor).
• actions::TA<:ActionSpace: Action space (see ContinuousActions and DiscreteActions).

ContinuousDP(; f, g, discount, x_lb, x_ub, shocks, weights)
ContinuousDP(; f, g, discount, actions, shocks, weights)
ContinuousDP(f, g, discount, shocks, weights, actions)
ContinuousDP(f, g, discount, shocks, weights, x_lb, x_ub)

Constructor for ContinuousDP. The problem is specified by its primitives only; the interpolation basis is supplied separately through the solver passed to solve (see CollocationSolver and LQASolver). Give either actions, or both x_lb and x_ub (equivalent to actions = ContinuousActions(x_lb, x_ub)).

Arguments

• f: Reward function f(s, x).
• g: State transition function g(s, x, e).
• discount::Real: Discount factor.
• actions::ActionSpace: Action space; alternatively give x_lb, x_ub (lower and upper bound of the action as functions of the state) for a one-dimensional continuous action space.
• shocks::AbstractVecOrMat: Discretized shock nodes.
• weights: Probability weights for the shock nodes. Either a fixed probability vector with one weight per shock node (length(weights) == size(shocks, 1)), or a callable returning such a collection: weights(s) for a state-dependent distribution, or weights(s, x) for a state-action-dependent one (if both arities apply, the state-action form is used). A callable returning a Tuple or a statically-sized vector (e.g. a StaticArrays.SVector) keeps the solver sweeps allocation-free; returning a freshly allocated Vector is supported but allocation-lean. The callable is evaluated on every objective evaluation of the inner maximization (several times per visited state), so it should be cheap; precompute anything expensive outside it. At solve initialization, callable weights are probed once for an indexable, real-valued return with one entry per shock node; because the probe state-action pair need not be feasible, an exception from that probe is tolerated. Regardless of the probe result, every callable-weight evaluation checks that its length matches the shock support — this per-fetch check is the correctness guarantee. The solve operators do not check that weights sum to one (sub-stochastic weights are permitted and act as additional discounting). Simulation, by contrast, requires a proper probability vector at every visited (s, x) — simulate throws an ArgumentError for weights that are negative, non-finite, or do not sum to one, since the missing mass has no path-wise interpretation. With action-dependent weights the first-order-condition inner solver does not apply and solve automatically falls back to Brent (see solve).

ContinuousDP(cdp::ContinuousDP; f=cdp.f, g=cdp.g, discount=cdp.discount,
             shocks=cdp.shocks, weights=cdp.weights, actions=cdp.actions,
             x_lb=nothing, x_ub=nothing)

Construct a copy of cdp, optionally replacing selected model components. The x_lb/x_ub keywords replace the corresponding bound of a continuous action space.

ActionSpace

Abstract type for action spaces of a ContinuousDP. Concrete subtypes: ContinuousActions, DiscreteActions.

ContinuousActions{M,Tlb,Tub}

Continuous action space: an M-dimensional box [x_lb(s), x_ub(s)] that may depend on the state. Construct with ContinuousActions(x_lb, x_ub) for M == 1 (with x_lb(s), x_ub(s) returning scalars) or ContinuousActions{M}(x_lb, x_ub) for M > 1 (with the bound functions returning length-M tuples or vectors).

For M > 1, actions are passed to the reward and transition functions as length-M collections indexable by x[1], ..., x[M] (tuples or views), and policy functions are stored as n x M matrices (one row per state node).

Fields

• x_lb::Tlb: Lower bound of the action as a function of the state.
• x_ub::Tub: Upper bound of the action as a function of the state.

DiscreteActions{TA}

Discrete (finite) action space, given by a vector of action values of any homogeneous type TA — numbers, tuples, labels, etc.; the values are passed opaquely to the reward and transition functions. Internally the solvers work with the indices into vals, and the solution exposes both the values (res.X) and the indices (res.X_ind), following the MarkovChain/state_values convention of QuantEcon.

Infeasible state-action pairs are expressed by the reward function returning -Inf; the enumeration over actions then skips the transition evaluation for that candidate. A well-posed model should have at least one feasible action at every state the solver evaluates. If every action is infeasible at some state, the first action is retained as a fallback, and subsequent policy evaluation may call g for that pair — g should tolerate such calls.

Fields

• vals::Vector{TA}: Action values.

CollocationSolver(basis; algorithm=PFI, inner_solver=:foc,
                  tol=sqrt(eps()), max_iter=500)
CollocationSolver(; basis, algorithm=PFI, inner_solver=:foc,
                  tol=sqrt(eps()), max_iter=500)

Bellman equation collocation solver configuration for solve: the interpolation basis together with the algorithm parameters. The problem itself (a ContinuousDP) holds only the model primitives.

Arguments

• basis::Basis: Object that contains the interpolation basis information; its domain is the approximation domain of the value function.
• algorithm::Type{<:DPAlgorithm}: PFI for policy function iteration or VFI for value function iteration (for linear-quadratic approximation see LQASolver).
• inner_solver::Symbol: How to solve the inner maximization over actions. :foc (default) solves the first-order condition by safeguarded root-finding, warm-started across iterations, using the exact gradient of the fitted value function and finite differences of f and g. It is intended for smooth, effectively concave inner problems where the first-order condition identifies the maximizing action: it falls back to Brent maximization state-by-state when derivative evaluation is unavailable or non-finite (and is used only for continuously differentiable bases: any Chebyshev, or splines of degree >= 2), but it does not attempt to detect nonconcavity or multiple stationary points. Use inner_solver=:brent for the derivative-free path. The choice has no effect for discrete action spaces (solved by exact enumeration), but the value is still validated. With action-dependent shock weights weights(s, x) (see ContinuousDP) the first-order condition acquires a term the FOC solver does not compute, so :foc automatically behaves as :brent.
• tol::Real: Convergence tolerance.
• max_iter::Integer: Maximum number of iterations.

LQASolver(basis; point)
LQASolver(; basis, point)

Linear-quadratic approximation solver configuration for solve: the model is approximated around the reference point point = (s, x, e) and the resulting LQ problem is solved exactly; the value function of the LQ solution is then represented in the interpolation basis (so the result has the same interface as a collocation solution).

Arguments

• basis::Basis: Object that contains the interpolation basis information.
• point::Tuple: The reference point (s, x, e) (typically a steady state) around which the LQ approximation is constructed.

solve(cdp, solver)

Solve the continuous-state dynamic program cdp with the given solver configuration (CollocationSolver or LQASolver).

Arguments

• cdp::ContinuousDP: The dynamic program to solve.
• solver: Solver configuration.
• v_init::Vector{Float64}: Optional keyword; initial value function values at the interpolation nodes of solver.basis.
• verbose::Integer: Optional keyword; level of feedback (0 for no output, 1 for warnings only, 2 for warnings and convergence messages during iteration).
• print_skip::Integer: Optional keyword; if verbose == 2, how many iterations between print messages.

Returns

• res::CDPSolveResult: Solution object of the dynamic program.

VFI

Value function iteration algorithm for solve.

PFI

Policy function iteration algorithm for solve.

LQA

Linear-quadratic approximation algorithm for solve.

Used via solve(cdp, LQASolver(basis; point=(s, x, e))) to approximate the model around a reference point and solve the resulting LQ problem.

set_eval_nodes!(res, s_nodes_coord)

Set the evaluation nodes and recompute the value/policy functions.

Arguments

• res::CDPSolveResult: Solution object to update in place.
• s_nodes_coord::NTuple{N,AbstractVector}: Coordinate vectors of the new evaluation nodes.

simulate([rng=GLOBAL_RNG], res, s_init, ts_length)

Generate a sample path of state variable(s) from a solved model.

Arguments

• rng::AbstractRNG: Random number generator.
• res::CDPSolveResult: Solution object of the dynamic program.
• s_init: Initial value of state variable(s).
• ts_length::Integer: Length of simulation.

Returns

• s_path::VecOrMat: Generated sample path of state variable(s).

simulate!([rng=GLOBAL_RNG], s_path, res, s_init)

Generate a sample path of state variable(s) from a solved model.

Arguments

• rng::AbstractRNG: Random number generator.
• s_path::VecOrMat: Array to store the generated sample path.
• res::CDPSolveResult: Solution object of the dynamic program.
• s_init: Initial value of state variable(s).

Returns

• s_path::VecOrMat: Generated sample path of state variable(s).

approx_lq(s_star, x_star, f_star, Df_star, DDf_star, g_star, Dg_star,
          discount)

Construct a QuantEcon.LQ instance that approximates the dynamic program around a steady state.

Arguments

• s_star::ScalarOrArray{T}: Steady-state value of the state variable(s).
• x_star::ScalarOrArray{T}: Steady-state value of the action variable(s).
• f_star::Real: Reward function evaluated at the steady state.
• Df_star::AbstractVector{T}: Gradient of the reward function f at the steady state, Df_star = [f_s', f_x'].
• DDf_star::AbstractMatrix{T}: Hessian of the reward function f at the steady state, DDf_star = [f_ss f_sx; f_xs f_xx].
• g_star::ScalarOrArray{T}: State transition function evaluated at the steady state.
• Dg_star::AbstractMatrix{T}: Jacobian of the transition function g at the steady state, Dg_star = [g_s, g_x].
• discount::Real: Discount factor.

Returns

• lq::QuantEcon.LQ: The LQ approximation.

CDPSolveResult{Algo,N,TCDP,TI,TE,TX}

Type storing the solution of a continuous-state dynamic program obtained by algorithm Algo.

Fields

• cdp::TCDP<:ContinuousDP: The dynamic program that was solved.
• interp::TI<:Interp{N}: The interpolation scheme used by the solver.
• tol::Float64: Convergence tolerance used by the solver.
• max_iter::Int: Maximum number of iterations allowed.
• C::Vector{Float64}: Basis coefficient vector for the fitted value function.
• converged::Bool: Whether the algorithm converged.
• num_iter::Int: Number of iterations performed.
• inner_solver::Symbol: Inner solver used by solve (:foc or :brent); also used when re-evaluating the policy (e.g. by set_eval_nodes!) for multi-dimensional continuous actions.
• eval_nodes::TE<:VecOrMat: Nodes at which the solution is evaluated. Defaults to interp.S.
• eval_nodes_coord::NTuple{N,Vector{Float64}}: Coordinate vectors of the evaluation nodes along each dimension. Defaults to interp.Scoord.
• V::Vector{Float64}: Value function evaluated at eval_nodes.
• X::TX<:AbstractVecOrMat: Policy function (action values) evaluated at eval_nodes (an n x M matrix for M-dimensional continuous actions).
• X_ind::Vector{Int}: For a discrete action space, the indices into cdp.actions.vals corresponding to X; empty otherwise.
• resid::Vector{Float64}: Approximation residuals at eval_nodes.

(res::CDPSolveResult)(s_nodes)

Evaluate the solved model at user-supplied state nodes.

Returns (V, X, resid), where V is the value function, X is the greedy policy, and resid is the Bellman residual at s_nodes.

CDPWorkspace{TF,TD}

Preallocated buffers used by the solution algorithms for a dynamic program bound to an interpolation scheme. Construct with CDPWorkspace(cp; inner_solver=:foc).

Not thread-safe: use one workspace per thread.

Fields

• fec::TF<:FunEvalCache: Workspace for point evaluation of the value function.
• dfecs::TD: Tuple of DerivFunEvalCaches for the gradient of the value function (one per state dimension), or nothing if the first-order-condition solver does not apply (basis not continuously differentiable, or action-dependent shock weights; see below).
• Tv::Vector{Float64}: Buffer for updated values at the interpolation nodes.
• X::TX<:VecOrMat{Float64}: Buffer for updated actions at the interpolation nodes for a continuous action space, an n x M matrix for M-dimensional actions (initialized to NaN; also serves as the warm start for the inner maximization in the next sweep).
• X_ind::Vector{Int}: Buffer for updated action indices at the interpolation nodes for a discrete action space; empty otherwise.
• inner_solver::Symbol: :foc to solve the inner maximization via its first-order condition (with Brent as automatic fallback), or :brent to always use derivative-free Brent maximization. Has no effect for a discrete action space (solved by enumeration), but is still validated.

CDPWorkspace(cp::_CollocationProblem; inner_solver=:foc)

Construct a CDPWorkspace for the bound problem cp. See solve for the meaning of inner_solver.

Interp{N,TB,TS,TM,TL}

Type representing an interpolation scheme on an N-dimensional domain.

Fields

• basis::TB<:Basis{N}: Object that contains the interpolation basis information.
• S::TS<:VecOrMat: Vector or Matrix that contains interpolation nodes (collocation points).
• Scoord::NTuple{N,Vector{Float64}}: Coordinate vectors of the interpolation nodes along each dimension.
• length::Int: Total number of interpolation nodes on the tensor grid.
• size::NTuple{N,Int}: Number of interpolation nodes along each dimension.
• lb::NTuple{N,Float64}: Lower bounds of the domain.
• ub::NTuple{N,Float64}: Upper bounds of the domain.
• Phi::TM<:AbstractMatrix: Basis matrix evaluated at the interpolation nodes.
• Phi_lu::TL<:Factorization: LU factorization of Phi.

Interp(basis)

Construct an Interp from a Basis.

Arguments

• basis::Basis: Object that contains the interpolation basis information.

PolicyFunction(res::CDPSolveResult)

Callable object evaluating the policy function at a single state point: pf = PolicyFunction(res); pf(s). For a discrete action space, the greedy action is recomputed exactly at s (a discrete policy is never interpolated); for a continuous action space, the policy values res.X are interpolated piecewise linearly over the evaluation nodes and clamped into [x_lb(s), x_ub(s)]. Returns an action value for scalar actions and a length-M tuple for M-dimensional continuous actions.

The policy data are shared with res at construction time: for continuous actions, construct after the final set_eval_nodes! call. The evaluator machinery is allocation-free; total per-call allocations depend on the user-supplied functions invoked during evaluation (the action-bound functions for continuous actions; the reward and transition functions for discrete ones). Not thread-safe: use one instance per thread.

ValueFunction(res::CDPSolveResult)

Callable object evaluating the fitted value function at a single state point: vf = ValueFunction(res); vf(s). Evaluation is non-allocating; the basis coefficients are shared with res (not copied).

Not thread-safe: use one instance per thread (as with CDPWorkspace).

_solve!(cp, res, ws, verbose, print_skip; point)

Implement linear quadratic approximation. See solve for further details.

_solve!(cp, res, ws, verbose, print_skip)

Implement policy iteration. See solve for further details.

_solve!(cp, res, ws, verbose, print_skip)

Implement value iteration. See solve for further details.

compute_greedy!(cp, C, X)
compute_greedy!(cp, ss, C, X[, fec])

Compute the greedy policy for the given basis coefficients.

Arguments

• cp::_CollocationProblem: The dynamic program bound to its interpolation scheme.
• ss::AbstractArray{Float64}: Interpolation nodes.
• C::Vector{Float64}: Basis coefficient vector for the value function.
• X::Vector{Float64}: A buffer array to hold the updated policy function.
• fec::FunEvalCache: Workspace for point evaluation of the value function. Constructed internally if not given.

Returns

• X::Vector{Float64}: Updated policy function vector.

evaluate!(res[, fec])

Evaluate the value function and the policy function at the evaluation nodes.

The result arrays res.V, res.X, and res.resid are updated in place when their lengths already match the number of evaluation nodes, and reallocated otherwise.

Arguments

• res::CDPSolveResult: Solution object to update in place.
• fec::FunEvalCache: Workspace for point evaluation of the value function. Constructed internally if not given.

evaluate_policy!(cp, X, C[, fec])

Compute the value function for a given policy and update the basis coefficients: solve (Phi - beta * E[Phi(g(S, X, e))]) C = f(S, X), where the expected-basis matrix is assembled row by row with the non-allocating point-evaluation kernels. When the collocation matrix Phi is sparse (spline and piecewise linear bases), the system matrix is assembled and factorized in sparse form.

Arguments

• cp::_CollocationProblem: The dynamic program bound to its interpolation scheme.
• X::AbstractVecOrMat: Policy function (action values); an n x M matrix for M-dimensional continuous actions.
• C::Vector{Float64}: A buffer array to hold the basis coefficients.
• fec::FunEvalCache: Workspace whose per-dimension caches are used for basis evaluation at the next states. Constructed internally if not given.

Returns

• C::Vector{Float64}: Updated basis coefficient vector.

operator_iteration!(T, C, tol, max_iter; verbose=2, print_skip=50)

Iterate an operator on the basis coefficients until convergence.

Arguments

• T::Function: Operator that updates basis coefficients (one step of VFI or PFI).
• C::Vector{Float64}: Initial basis coefficient vector.
• tol::Float64: Convergence tolerance.
• max_iter::Integer: Maximum number of iterations.
• verbose::Integer: Level of feedback (0 for no output, 1 for warnings only, 2 for warning and convergence messages during iteration).
• print_skip::Integer: If verbose == 2, how many iterations between print messages.

Returns

• converged::Bool: Whether the iteration converged.
• i::Int: Number of iterations performed.

policy_iteration_operator!(cp, C, X)
policy_iteration_operator!(cp, C, ws)

Perform one step of policy function iteration and update the basis coefficients.

Arguments

• cp::_CollocationProblem: The dynamic program bound to its interpolation scheme.
• C::Vector{Float64}: Basis coefficient vector for the value function.
• X::Vector{Float64}: A buffer array to hold the updated policy function.
• ws::CDPWorkspace: Workspace for the solution algorithms.

Returns

• C::Vector{Float64}: Updated basis coefficient vector.

s_wise_max!(cp, ss, C, Tv[, fec])

Find optimal value for each grid point. These helpers apply to one-dimensional continuous action spaces; discrete and multi-dimensional action spaces are handled internally by the operators.

Arguments

• cp::_CollocationProblem: The dynamic program bound to its interpolation scheme.
• ss::AbstractArray{Float64}: Interpolation nodes.
• C::Vector{Float64}: Basis coefficient vector for the value function.
• Tv::Vector{Float64}: A buffer array to hold the updated value function.
• fec::FunEvalCache: Workspace for point evaluation of the value function. Constructed internally if not given.

Returns

• Tv::Vector{Float64}: Updated value function vector.

s_wise_max!(cp, ss, C, Tv, X[, fec])

Find optimal value and action for each grid point.

Arguments

• cp::_CollocationProblem: The dynamic program bound to its interpolation scheme.
• ss::AbstractArray{Float64}: Interpolation nodes.
• C::Vector{Float64}: Basis coefficient vector for the value function.
• Tv::Vector{Float64}: A buffer array to hold the updated value function.
• X::Vector{Float64}: A buffer array to hold the updated policy function.
• fec::FunEvalCache: Workspace for point evaluation of the value function. Constructed internally if not given.

Returns

• Tv::Vector{Float64}: Updated value function vector.
• X::Vector{Float64}: Updated policy function vector.

s_wise_max(cp, ss, C[, fec])

Find optimal value and action for each grid point.

Arguments

• cp::_CollocationProblem: The dynamic program bound to its interpolation scheme.
• ss::AbstractArray{Float64}: Interpolation nodes.
• C::Vector{Float64}: Basis coefficient vector for the value function.
• fec::FunEvalCache: Workspace for point evaluation of the value function. Constructed internally if not given.

Returns

• Tv::Vector{Float64}: Value function vector.
• X::Vector{Float64}: Policy function vector.

bellman_operator!(cp, C, Tv)
bellman_operator!(cp, C, ws)

Apply the Bellman operator and update the basis coefficients. Values are stored in Tv (or ws.Tv).

Arguments

• cp::_CollocationProblem: The dynamic program bound to its interpolation scheme.
• C::Vector{Float64}: Basis coefficient vector for the value function.
• Tv::Vector{Float64}: Vector to store values.
• ws::CDPWorkspace: Workspace for the solution algorithms.

Returns

• C::Vector{Float64}: Updated basis coefficient vector.

_QuadratureKernel{TG,TR,TW} <: _TransitionKernel

The structured transition kernel of a ContinuousDP: the transition function g paired with fixed quadrature nodes and a weights carrier. Branch j has next state g(s, x, shocks[j]) and weight w[j], where w is the weight container at (s, x): the fixed vector itself, or the value of a user-supplied callable weights(s) / weights(s, x) (wrapped in _StateWeights / _StateActionWeights by _build_kernel).

Besides the general _TransitionKernel contract, this kernel provides indexed access — _branch_weights(ker, s, x) (the weight container, fetched once per (s, x)) and _branch_state(ker, s, x, j) — on which the FOC solvers rely: derivative-based consumers re-evaluate _branch_state at perturbed actions holding the branch index fixed. Hence _forces_brent is false for fixed or state-only weights, true for action-dependent weights (whose derivative term the FOC solvers do not compute). The sampling helper _draw_branch_index(rng, ker, s, x) (used by simulate! for callable weights) also lives at this tier: it returns an index into the fixed shock nodes.

Allocation contract: a callable weights function returning a Tuple or a statically-sized vector (e.g. a StaticArrays.SVector) keeps the sweeps allocation-free; returning a freshly allocated Vector is supported but makes the path allocation-lean instead. Cost contract: the callable is invoked on every objective evaluation of the inner maximization — several times per state per sweep, not once per state — so it should be cheap; anything expensive belongs in a table or interpolant precomputed outside the callable.

_TransitionKernel

Abstract supertype of internal transition-kernel representations: the distribution of the next state given (s, x), as finitely many weighted branches (s'_k, w_k).

The general contract is:

• _branch_sum(f, ker, s, x, args...): return the sum of f(s', w, args...) over the branches at (s, x);
• _foreach_branch(f, ker, s, x, args...): call f(s', w, args...) on each branch at (s, x) (for side-effecting consumers such as policy-system row assembly);
• _draw_next_state(rng, ker, s, x): draw a next state from the branch distribution — required only when simulate! support is desired;
• _forces_brent(ker): whether the first-order-condition inner solver must fall back to Brent. Defaults to true: the FOC paths additionally require the indexed access described under _QuadratureKernel, which a general kernel need not provide.

These operations cover the generic Bellman expectation (_expected_value), Brent maximization, discrete-action enumeration, policy-evaluation system assembly (_append_state_rows! via _foreach_branch), and — when _draw_next_state is implemented — simulation. A general kernel need not expose stable branch indices: _draw_branch_index is a helper of the structured _QuadratureKernel, not part of this contract.

In the traversals, f is a top-level function and args its explicit payload. This avoids a capturing closure (which would materialize on the heap here), but the generic traversal indirection is nevertheless not allocation-free under the current implementation: measurements show roughly 50-80 bytes per call. Allocation-sensitive structured consumers therefore retain direct _QuadratureKernel specializations.

Opting into the FOC solvers by overriding _forces_brent(ker) = false declares more than the traversal contract: such a kernel must additionally provide the structured tier's indexed access (_branch_weights, _branch_state) with a stable branch count and stable branch identity under action perturbations, and branch probabilities that do not depend on the action (the FOC solvers do not compute probability derivatives).

_build_kernel(cdp, s_probe, x_probe)
_build_kernel(cp::_CollocationProblem)

Construct the _QuadratureKernel of cdp. For a fixed weight vector this is free (a non-escaping bundle of references). For callable weights the arity is detected by hasmethod (weights(s, x) takes precedence if both apply); the user callable is never evaluated here — construction sits on the sweep path, so return-shape validation happens once per solve instead, in _validate_weights at workspace creation. Sweep-level entry points construct the kernel once per call and pass it to the per-state solvers.

_draw_next_state(rng, ker, s, x)

Draw a next state from the branch distribution at (s, x). Part of the general _TransitionKernel contract for kernels used with simulate!; the structured kernel draws a branch index by inverse CDF and evaluates its next state. There is no generic fallback: a general kernel that should support simulation defines its own method.

_validate_weights(cp::_CollocationProblem)

One-time validation of callable weights, called at workspace creation (once per solve) so that nothing here rides the sweep path: classify the arity (throws for a callable of the wrong arity) and evaluate one probe call at the first interpolation node and its probe action to validate the returned container. The probe point is not guaranteed feasible, and a model may legitimately error there (e.g. weights undefined for an action with f(s, x) == -Inf, whose transition the solvers never evaluate): if the probe call throws, this validation is skipped — correctness then rests on the per-fetch length check built into the weight carriers, so a malformed return still errors at its first actual use instead of silently truncating the branch loop.

_objective_and_deriv(cdp, ker, s, C, fec, dfecs, x, x_lb, x_ub)

Evaluate the inner objective H(x) = f(s, x) + beta * E[V^(g(s, x, e))] and its derivative H'(x) = f_x + beta * E[grad V^(g(s, x, e)) . g_x] at the action x. The gradient of the fitted value function is evaluated exactly via dfecs (whose coefficients must be set with set_coefs! beforehand), while f_x and g_x are computed by central finite differences with the step shrunk so that x ± h stay within [x_lb, x_ub] (so that f and g are never called at infeasible actions); if no adequate step exists, (NaN, NaN) is returned to trigger the Brent fallback.

Returns

• H::Float64, Hp::Float64: Objective value and derivative (may be non-finite, in which case the caller should fall back to Brent).

_s_wise_max!(cdp, ker, s, C, fec)

Find the optimal value and action at a given state s.

Arguments

• cdp::ContinuousDP: The dynamic program.
• ker::_TransitionKernel: Its transition kernel (see _build_kernel).
• s: State point at which to maximize.
• C: Basis coefficient vector for the value function.
• fec::FunEvalCache: Workspace for evaluating the value function at the next states.

Returns

• v::Float64: Optimal value at s.
• x::Float64: Optimal action at s.

_s_wise_max_discrete!(cdp, ker, s, C, fec)

Find the optimal value and action at state s by enumeration over the discrete action set.

Returns

• v::Float64: Optimal value at s (-Inf if every action is infeasible).
• k::Int: Index of the optimal action in cdp.actions.vals (the first index if every action is infeasible; ties go to the lowest index).

_s_wise_max_foc!(cdp, ker, s, C, fec, dfecs, x_prev)

Find the optimal value and action at state s by solving the first-order condition H'(x) = 0 with safeguarded bracketing root-finding (regula falsi with the Illinois modification), warm-started at x_prev (NaN for a cold start). Falls back to the Brent-based _s_wise_max! whenever the objective or its derivative is non-finite at a required point. Corner solutions are detected from the sign of H' during the bracketing expansion.

The coefficients of dfecs must have been set with set_coefs!(., C).

_s_wise_max_foc_sweep!(cp, C, Tv, X, fec, dfecs)

Run the FOC-based inner maximization over all interpolation nodes, storing values in Tv and maximizers in X. The previous contents of X serve as warm starts (NaN entries mean cold start). Sets the coefficients of dfecs from C. Falls back to Brent state-by-state on exceptions from the model functions (e.g. a DomainError at a finite-difference point).

_s_wise_max_multi!(cdp, ker, s, C, fec, dfecs, xout, use_foc)

Find the optimal value and action at state s for an M-dimensional continuous action space by box-constrained maximization, warm-started at xout (NaN entries mean cold start from the box center, with a coarse feasible-start probe if that is infeasible). With use_foc = true (and dfecs available), Fminbox(LBFGS) with the analytic objective gradient is tried first, falling back to derivative-free cyclic coordinate-wise Brent maximization on any failure or non-finite outcome; with use_foc = false, coordinate-wise Brent is used directly.

Writes the maximizer into xout and returns the maximized value.