From 05cc1192f581e289e0bb6ecce5b314bdad2a5840 Mon Sep 17 00:00:00 2001 From: Daisuke Oyama Date: Sun, 26 Jul 2026 15:11:06 +0900 Subject: [PATCH 1/2] ENH: Make DPSolveResult immutable with solve metadata Implements issue #390 with decisions D1(a), D2(c), D3(b), D4 skip, D5, plus the ddp back-reference pulled forward from the #398/#120 plan. DPSolveResult becomes an immutable struct DPSolveResult{Algo,Tval,TD<:DiscreteDP,TMC<:MarkovChain}, constructed complete at the end of solve, with new fields ddp (held by reference), epsilon, max_iter, and k (the options passed to or defaulted by solve). The solver internals go array-level: _solve! now has the signature _solve!(ddp, ::Type{Algo}, v, Tv, sigma, max_iter, epsilon, k) -> num_iter, and the five (ddp, ddpr)-taking methods (bellman_operator!, compute_greedy!, evaluate_policy, RQ_sigma, MarkovChain) are removed, replaced by their existing array-level twins plus the new MarkovChain(ddp, sigma). compute_greedy! is removed entirely (its only method was the (ddp, ddpr) one). New exported accessor controlled_mc(res) = res.mc returns the controlled Markov chain type-stably. Behavior-preserving for solution values, iteration counts, and convergence criteria. Co-Authored-By: Claude Fable 5 --- src/QuantEcon.jl | 4 +- src/markov/ddp.jl | 295 +++++++++++++++++++++------------------------- test/test_ddp.jl | 39 +++++- 3 files changed, 170 insertions(+), 168 deletions(-) diff --git a/src/QuantEcon.jl b/src/QuantEcon.jl index 32e73f39..e6ad3ecf 100644 --- a/src/QuantEcon.jl +++ b/src/QuantEcon.jl @@ -130,8 +130,8 @@ export # ddp DiscreteDP, VFI, PFI, MPFI, solve, RQ_sigma, evaluate_policy, bellman_operator, compute_greedy, - bellman_operator!, compute_greedy!, num_states, backward_induction, - num_sa_pairs, to_sa_pair_form, to_product_form, + bellman_operator!, num_states, backward_induction, + num_sa_pairs, to_sa_pair_form, to_product_form, controlled_mc, # zeros / optimization bisect, brenth, brent, ridder, expand_bracket, divide_bracket, diff --git a/src/markov/ddp.jl b/src/markov/ddp.jl index 80b1abf8..092948a1 100644 --- a/src/markov/ddp.jl +++ b/src/markov/ddp.jl @@ -448,50 +448,52 @@ struct MPFI <: DDPAlgorithm end DPSolveResult Object for retaining results and associated metadata after solving the model. +Returned by `solve`; immutable and constructed complete (the contents of the +array fields may still be mutated). # Fields +- `ddp::DiscreteDP`: The model solved. Held by reference, not copied: if the + model's `R`, `Q`, or `beta` attributes are mutated after `solve`, this field + reflects the mutated model while the solution fields below remain those of + the model as solved. - `v::Vector{Tval}`: Value function vector. -- `Tv::Vector{Tval}`: Temporary value function vector. +- `Tv::Vector{Tval}`: Value function vector from the last iteration step + (internal work buffer). - `num_iter::Int`: Number of iterations. - `sigma::Vector{Int}`: Policy function vector. - `mc::MarkovChain`: Controlled Markov chain. +- `epsilon::Float64`: Value for epsilon-optimality passed to (or defaulted by) + `solve`, whether or not the solution method uses it. +- `max_iter::Int`: Maximum number of iterations passed to (or defaulted by) + `solve`. +- `k::Int`: Number of iterations for partial policy evaluation passed to (or + defaulted by) `solve`, whether or not the solution method uses it. """ -mutable struct DPSolveResult{Algo<:DDPAlgorithm,Tval<:Real} +struct DPSolveResult{Algo<:DDPAlgorithm,Tval<:Real,TD<:DiscreteDP, + TMC<:MarkovChain} + ddp::TD v::Vector{Tval} Tv::Vector{Tval} num_iter::Int sigma::Vector{Int} - mc::MarkovChain - - function DPSolveResult{Algo,Tval}( - ddp::DiscreteDP - ) where {Algo,Tval} - v = s_wise_max(ddp, ddp.R) # Initialise v with v_init - ddpr = new{Algo,Tval}(v, similar(v), 0, similar(v, Int)) - - # fill in sigma with proper values - compute_greedy!(ddp, ddpr) - ddpr - end - - # method to pass initial value function (skip the s-wise max) - function DPSolveResult{Algo,Tval}( - ddp::DiscreteDP, v::AbstractVector - ) where {Algo,Tval} - # copy the input so that the caller's array is not overwritten by - # the solution methods - v_own = Vector{Tval}(undef, length(v)) - copyto!(v_own, v) - ddpr = new{Algo,Tval}(v_own, similar(v_own), 0, similar(v_own, Int)) - - # fill in sigma with proper values - compute_greedy!(ddp, ddpr) - ddpr - end + mc::TMC + epsilon::Float64 + max_iter::Int + k::Int end +# metadata as required keywords, so that nothing can be silently omitted; +# the only constructor `solve` calls +DPSolveResult{Algo}( + ddp::TD, v::Vector{Tval}, Tv::Vector{Tval}, num_iter::Integer, + sigma::Vector{Int}, mc::TMC; + epsilon::Real, max_iter::Integer, k::Integer + ) where {Algo<:DDPAlgorithm,Tval,TD<:DiscreteDP,TMC<:MarkovChain} = + DPSolveResult{Algo,Tval,TD,TMC}(ddp, v, Tv, num_iter, sigma, mc, + epsilon, max_iter, k) + # ------------------------ # # Bellman operator methods # # ------------------------ # @@ -539,29 +541,6 @@ function bellman_operator!( Tv, sigma end -@doc doc""" - bellman_operator!(ddp, ddpr) - -Apply the Bellman operator using `v=ddpr.v`, `Tv=ddpr.Tv`, and `sigma=ddpr.sigma`. - -# Arguments - -- `ddp::DiscreteDP`: Object that contains the model parameters. -- `ddpr::DPSolveResult`: Object that contains result variables. - -# Returns - -- `Tv::typeof(ddpr.Tv)`: Updated value function vector. -- `sigma::typeof(ddpr.sigma)`: Updated policy function vector. - -# Notes - -Updates `ddpr.Tv` and `ddpr.sigma` inplace. - -""" -bellman_operator!(ddp::DiscreteDP, ddpr::DPSolveResult) = - bellman_operator!(ddp, ddpr.v, ddpr.Tv, ddpr.sigma) - """ bellman_operator!(ddp, v, sigma) @@ -619,28 +598,6 @@ bellman_operator(ddp::DiscreteDP, v::AbstractVector) = # Compute greedy methods # # ---------------------- # -@doc doc""" - compute_greedy!(ddp, ddpr) - -Compute the ``v``-greedy policy. - -# Arguments - -- `ddp::DiscreteDP`: Object that contains the model parameters. -- `ddpr::DPSolveResult`: Object that contains result variables. - -# Returns - -- `sigma::Vector{Int}`: Array containing `v`-greedy policy rule. - -# Notes - -Modifies ddpr.sigma and ddpr.Tv in place. - -""" -compute_greedy!(ddp::DiscreteDP, ddpr::DPSolveResult) = - (bellman_operator!(ddp, ddpr); ddpr.sigma) - @doc doc""" compute_greedy(ddp, v) @@ -667,26 +624,6 @@ end # Evaluate policy methods # # ----------------------- # -""" - evaluate_policy(ddp, ddpr) - -Method of `evaluate_policy` that extracts sigma from a `DPSolveResult`. - -See other docstring for details. - -# Arguments - -- `ddp::DiscreteDP`: Object that contains the model parameters. -- `ddpr::DPSolveResult`: Object that contains result variables. - -# Returns - -- `v_sigma::Array{Float64}`: Value vector of `sigma`, of length `n`. - -""" -evaluate_policy(ddp::DiscreteDP, ddpr::DPSolveResult) = - evaluate_policy(ddp, ddpr.sigma) - """ evaluate_policy(ddp, sigma) @@ -751,19 +688,32 @@ convergence. function solve(ddp::DiscreteDP{T}, method::Type{Algo}=VFI; max_iter::Integer=250, epsilon::Real=1e-3, k::Integer=20) where {Algo<:DDPAlgorithm,T} - ddpr = DPSolveResult{Algo,T}(ddp, _default_v_init(ddp, Algo)) - _solve!(ddp, ddpr, max_iter, epsilon, k) - ddpr.mc = MarkovChain(ddp, ddpr) - ddpr + _solve(ddp, _default_v_init(ddp, Algo), Algo, max_iter, epsilon, k) end function solve(ddp::DiscreteDP{T}, v_init::AbstractVector{T}, method::Type{Algo}=VFI; max_iter::Integer=250, epsilon::Real=1e-3, k::Integer=20) where {Algo<:DDPAlgorithm,T} - ddpr = DPSolveResult{Algo,T}(ddp, v_init) - _solve!(ddp, ddpr, max_iter, epsilon, k) - ddpr.mc = MarkovChain(ddp, ddpr) - ddpr + _solve(ddp, v_init, Algo, max_iter, epsilon, k) +end + +# allocate the solution arrays, seed the v-greedy policy, run the +# iteration loop, and construct the result once, complete +function _solve(ddp::DiscreteDP{T}, v_init::AbstractVector, + ::Type{Algo}, max_iter::Integer, epsilon::Real, + k::Integer) where {Algo<:DDPAlgorithm,T} + # copy the input so that the caller's array is not overwritten by + # the solution methods + v = Vector{T}(undef, length(v_init)) + copyto!(v, v_init) + Tv = similar(v) + sigma = similar(v, Int) + bellman_operator!(ddp, v, Tv, sigma) # initial v-greedy sigma + + num_iter = _solve!(ddp, Algo, v, Tv, sigma, max_iter, epsilon, k) + + DPSolveResult{Algo}(ddp, v, Tv, num_iter, sigma, MarkovChain(ddp, sigma); + epsilon=epsilon, max_iter=max_iter, k=k) end # default initial value function: v(s) = max_a r(s, a) in general, and @@ -839,44 +789,39 @@ end # --------- # """ - MarkovChain(ddp, ddpr) + MarkovChain(ddp, sigma) Returns the controlled Markov chain for a given policy `sigma`. # Arguments - `ddp::DiscreteDP`: Object that contains the model parameters. -- `ddpr::DPSolveResult`: Object that contains result variables. +- `sigma::AbstractVector{<:Integer}`: Policy rule vector. # Returns - `mc::MarkovChain`: Controlled Markov chain. """ -QuantEcon.MarkovChain(ddp::DiscreteDP, ddpr::DPSolveResult) = - MarkovChain(RQ_sigma(ddp, ddpr)[2]) +QuantEcon.MarkovChain(ddp::DiscreteDP, sigma::AbstractVector{<:Integer}) = + MarkovChain(RQ_sigma(ddp, sigma)[2]) """ - RQ_sigma(ddp, ddpr) + controlled_mc(res) -Method of `RQ_sigma` that extracts sigma from a `DPSolveResult`. - -See other docstring for details. +Return the Markov chain controlled by the optimal policy of a solved +model, `res.mc`. The chain is returned by reference, not copied. # Arguments -- `ddp::DiscreteDP`: Object that contains the model parameters. -- `ddpr::DPSolveResult`: Object that contains result variables. +- `res::DPSolveResult`: Object that contains result variables. # Returns -- `R_sigma::Vector`: Reward vector for `sigma`, of length `n`. -- `Q_sigma::AbstractMatrix`: Transition probability matrix for `sigma`, of - shape `(n, n)`; dense for the product formulation, and sparse if `ddp.Q` - is sparse. +- `mc::MarkovChain`: Controlled Markov chain. """ -RQ_sigma(ddp::DiscreteDP, ddpr::DPSolveResult) = RQ_sigma(ddp, ddpr.sigma) +controlled_mc(res::DPSolveResult) = res.mc """ RQ_sigma(ddp, sigma) @@ -1273,7 +1218,7 @@ function _max_abs_diff(x::AbstractVector, y::AbstractVector) end """ - _solve!(ddp, ddpr, max_iter, epsilon, k) + _solve!(ddp, ::Type{VFI}, v, Tv, sigma, max_iter, epsilon, k) Implements Value Iteration. @@ -1282,19 +1227,26 @@ NOTE: See `solve` for further details. # Arguments - `ddp::DiscreteDP`: Object that contains the model parameters. -- `ddpr::DPSolveResult{VFI}`: Object that contains result variables. +- `::Type{VFI}`: Algorithm type used for dispatch. +- `v::AbstractVector`: Value function vector, seeded with the initial + value function; updated in place with the solution. +- `Tv::AbstractVector`: Work buffer for the updated value function; + updated in place. +- `sigma::AbstractVector{Int}`: Policy function vector, seeded with the + initial `v`-greedy policy; updated in place with the solution. - `max_iter::Integer`: Maximum number of iterations. - `epsilon::Real`: Value for epsilon-optimality. - `k::Integer`: Number of iterations (not used for VFI). # Returns -- `ddpr::DPSolveResult{VFI}`: Updated result object. +- `num_iter::Int`: Number of iterations performed. """ function _solve!( - ddp::DiscreteDP, ddpr::DPSolveResult{VFI}, max_iter::Integer, - epsilon::Real, k::Integer + ddp::DiscreteDP, ::Type{VFI}, v::AbstractVector, Tv::AbstractVector, + sigma::AbstractVector{Int}, max_iter::Integer, epsilon::Real, + k::Integer ) if ddp.beta == 0.0 tol = Inf @@ -1304,27 +1256,28 @@ function _solve!( tol = epsilon * (1-ddp.beta) / (2*ddp.beta) end - vals = similar(ddp.R, eltype(ddpr.v)) # buffer for state-action values + vals = similar(ddp.R, eltype(v)) # buffer for state-action values + num_iter = 0 for i in 1:max_iter # updates Tv in place - bellman_operator!(ddp, ddpr.v, ddpr.Tv, ddpr.sigma, vals) + bellman_operator!(ddp, v, Tv, sigma, vals) - # compute error and update the v inside ddpr - err = _max_abs_diff(ddpr.Tv, ddpr.v) - copyto!(ddpr.v, ddpr.Tv) - ddpr.num_iter += 1 + # compute error and update v + err = _max_abs_diff(Tv, v) + copyto!(v, Tv) + num_iter += 1 if err < tol break end end - ddpr + num_iter end """ - _solve!(ddp, ddpr, max_iter, epsilon, k) + _solve!(ddp, ::Type{PFI}, v, Tv, sigma, max_iter, epsilon, k) Policy Function Iteration. @@ -1335,107 +1288,125 @@ NOTE: The epsilon is ignored in this method. It is only here so dispatch can # Arguments - `ddp::DiscreteDP`: Object that contains the model parameters. -- `ddpr::DPSolveResult{PFI}`: Object that contains result variables. +- `::Type{PFI}`: Algorithm type used for dispatch. +- `v::AbstractVector`: Value function vector; updated in place with the + solution. +- `Tv::AbstractVector`: Work buffer for the updated value function; + updated in place. +- `sigma::AbstractVector{Int}`: Policy function vector, seeded with the + initial `v`-greedy policy; updated in place with the solution. - `max_iter::Integer`: Maximum number of iterations. - `epsilon::Real`: Value for epsilon-optimality (not used for PFI). - `k::Integer`: Number of iterations (not used for PFI). # Returns -- `ddpr::DPSolveResult{PFI}`: Updated result object. +- `num_iter::Int`: Number of iterations performed. """ function _solve!( - ddp::DiscreteDP, ddpr::DPSolveResult{PFI}, max_iter::Integer, - epsilon::Real, k::Integer + ddp::DiscreteDP, ::Type{PFI}, v::AbstractVector, Tv::AbstractVector, + sigma::AbstractVector{Int}, max_iter::Integer, epsilon::Real, + k::Integer ) - old_sigma = copy(ddpr.sigma) + old_sigma = copy(sigma) if ddp.beta == 1.0 throw(ArgumentError("method invalid for beta = 1")) end + vals = similar(ddp.R, eltype(v)) # buffer for state-action values + + num_iter = 0 for i in 1:max_iter - ddpr.v = evaluate_policy(ddp, ddpr) - compute_greedy!(ddp, ddpr) + copyto!(v, evaluate_policy(ddp, sigma)) + bellman_operator!(ddp, v, Tv, sigma, vals) # v-greedy sigma - ddpr.num_iter += 1 - if old_sigma == ddpr.sigma + num_iter += 1 + if old_sigma == sigma break end - copyto!(old_sigma, ddpr.sigma) + copyto!(old_sigma, sigma) end - ddpr + num_iter end span(x::AbstractVector) = maximum(x) - minimum(x) midrange(x::AbstractVector) = mean(extrema(x)) """ - _solve!(ddp, ddpr, max_iter, epsilon, k) + _solve!(ddp, ::Type{MPFI}, v, Tv, sigma, max_iter, epsilon, k) Modified Policy Function Iteration. # Arguments - `ddp::DiscreteDP`: Object that contains the model parameters. -- `ddpr::DPSolveResult{MPFI}`: Object that contains result variables. +- `::Type{MPFI}`: Algorithm type used for dispatch. +- `v::AbstractVector`: Value function vector, seeded with the initial + value function; updated in place with the solution. +- `Tv::AbstractVector`: Work buffer for the updated value function; + updated in place. +- `sigma::AbstractVector{Int}`: Policy function vector, seeded with the + initial `v`-greedy policy; updated in place with the solution. - `max_iter::Integer`: Maximum number of iterations. - `epsilon::Real`: Value for epsilon-optimality. - `k::Integer`: Number of iterations for partial policy evaluation. # Returns -- `ddpr::DPSolveResult{MPFI}`: Updated result object. +- `num_iter::Int`: Number of iterations performed. """ function _solve!( - ddp::DiscreteDP, ddpr::DPSolveResult{MPFI}, max_iter::Integer, - epsilon::Real, k::Integer + ddp::DiscreteDP, ::Type{MPFI}, v::AbstractVector, Tv::AbstractVector, + sigma::AbstractVector{Int}, max_iter::Integer, epsilon::Real, + k::Integer ) if ddp.beta == 1.0 throw(ArgumentError("method invalid for beta = 1")) end - # NOTE: ddpr.v is used as is; when called through `solve` without - # v_init, it has been initialized with min r / (1 - beta), which - # guarantees convergence + # NOTE: v is used as is; when called through `solve` without v_init, + # it has been initialized with min r / (1 - beta), which guarantees + # convergence beta = ddp.beta - old_sigma = copy(ddpr.sigma) + old_sigma = copy(sigma) tol = beta > 0 ? epsilon * (1-beta) / beta : Inf - vals = similar(ddp.R, eltype(ddpr.v)) # buffer for state-action values - dif = similar(ddpr.v) + vals = similar(ddp.R, eltype(v)) # buffer for state-action values + dif = similar(v) + num_iter = 0 for i in 1:max_iter # updates Tv, sigma inplace - bellman_operator!(ddp, ddpr.v, ddpr.Tv, ddpr.sigma, vals) - dif .= ddpr.Tv .- ddpr.v + bellman_operator!(ddp, v, Tv, sigma, vals) + dif .= Tv .- v - ddpr.num_iter += 1 + num_iter += 1 # check convergence if span(dif) < tol - ddpr.v .= ddpr.Tv .+ midrange(dif) * beta / (1-beta) + v .= Tv .+ midrange(dif) * beta / (1-beta) break end # now update v to use the output of the bellman step when entering # policy loop - copyto!(ddpr.v, ddpr.Tv) + copyto!(v, Tv) # now do k iterations of policy iteration - R_sigma, Q_sigma = RQ_sigma(ddp, ddpr) + R_sigma, Q_sigma = RQ_sigma(ddp, sigma) for i in 1:k - mul!(ddpr.Tv, Q_sigma, ddpr.v) - ddpr.Tv .= R_sigma .+ beta .* ddpr.Tv - copyto!(ddpr.v, ddpr.Tv) + mul!(Tv, Q_sigma, v) + Tv .= R_sigma .+ beta .* Tv + copyto!(v, Tv) end end - ddpr + num_iter end diff --git a/test/test_ddp.jl b/test/test_ddp.jl index 36568174..354ad35d 100644 --- a/test/test_ddp.jl +++ b/test/test_ddp.jl @@ -151,11 +151,42 @@ Tests for markov/ddp.jl end end - @testset "compute_greedy! changes ddpr.v" begin + @testset "bellman_operator! overwrites the Tv buffer" begin res = solve(ddp0, VFI) - res.Tv[:] .= 500.0 - compute_greedy!(ddp0, res) - @test maximum(abs, res.Tv .- 500.0) > 0 + _Tv = fill(500.0, length(res.v)) + _sigma = similar(res.sigma) + bellman_operator!(ddp0, res.v, _Tv, _sigma) + @test maximum(abs, _Tv .- 500.0) > 0 + @test _sigma == res.sigma + end + + @testset "DPSolveResult metadata and construction" begin + for _ddp in ddp0_collection + for Algo in (VFI, PFI, MPFI) + # explicitly passed options are recorded + _res = solve(_ddp, Algo; max_iter=137, epsilon=1e-4, k=17) + @test _res.epsilon == 1e-4 + @test _res.max_iter == 137 + @test _res.k == 17 + + # defaulted options are recorded + _res_d = solve(_ddp, Algo) + @test _res_d.epsilon == 1e-3 + @test _res_d.max_iter == 250 + @test _res_d.k == 20 + + # the model is held by reference, not copied + @test _res_d.ddp === _ddp + + # the result is immutable and concretely typed + @test !ismutable(_res_d) + @test isconcretetype(typeof(_res_d)) + + # controlled_mc returns the stored chain, type-stably + @test controlled_mc(_res_d) === _res_d.mc + @test (@inferred controlled_mc(_res_d)) === _res_d.mc + end + end end @testset "value_iteration" begin From fe5aff579553942917fe799a758232273d094e24 Mon Sep 17 00:00:00 2001 From: Daisuke Oyama Date: Sun, 26 Jul 2026 16:49:07 +0900 Subject: [PATCH 2/2] FIX: Preserve promoted PFI workspace eltype; remove controlled_mc Addresses review of the DPSolveResult rebuild. The hoisted PFI state-action buffer is now allocated at promote_type(eltype(R), eltype(v), typeof(beta)), matching what the 4-arg bellman_operator! the greedy step previously went through allocates each iteration. A buffer in eltype(v) rounds near-tied action values to equality before the argmax for mixed-precision models (e.g. Float32 rewards with a Float64 beta) and can flip the selected policy via the first-action tie-break. Regression test added; verified failing before this fix and matching master after it. VFI and MPFI keep their eltype(v) buffers, which is pre-existing behavior not changed by this PR. controlled_mc is removed before ever being released: it duplicated the documented and now concretely-typed res.mc field, and the name misreads as computing the chain from the current res.sigma, which is what MarkovChain(ddp, sigma) does. No replacement accessor is added here; a uniform accessor over results and the planned POMDPs.jl policy wrapper is deferred to the extension PR for #398. The `@inferred` pin on the concrete mc field remains via a test-local helper. Co-Authored-By: Claude Fable 5 --- src/QuantEcon.jl | 2 +- src/markov/ddp.jl | 25 +++++++------------------ test/test_ddp.jl | 26 +++++++++++++++++++++++--- 3 files changed, 31 insertions(+), 22 deletions(-) diff --git a/src/QuantEcon.jl b/src/QuantEcon.jl index e6ad3ecf..169594fa 100644 --- a/src/QuantEcon.jl +++ b/src/QuantEcon.jl @@ -131,7 +131,7 @@ export DiscreteDP, VFI, PFI, MPFI, solve, RQ_sigma, evaluate_policy, bellman_operator, compute_greedy, bellman_operator!, num_states, backward_induction, - num_sa_pairs, to_sa_pair_form, to_product_form, controlled_mc, + num_sa_pairs, to_sa_pair_form, to_product_form, # zeros / optimization bisect, brenth, brent, ridder, expand_bracket, divide_bracket, diff --git a/src/markov/ddp.jl b/src/markov/ddp.jl index 092948a1..b8d5e567 100644 --- a/src/markov/ddp.jl +++ b/src/markov/ddp.jl @@ -806,23 +806,6 @@ Returns the controlled Markov chain for a given policy `sigma`. QuantEcon.MarkovChain(ddp::DiscreteDP, sigma::AbstractVector{<:Integer}) = MarkovChain(RQ_sigma(ddp, sigma)[2]) -""" - controlled_mc(res) - -Return the Markov chain controlled by the optimal policy of a solved -model, `res.mc`. The chain is returned by reference, not copied. - -# Arguments - -- `res::DPSolveResult`: Object that contains result variables. - -# Returns - -- `mc::MarkovChain`: Controlled Markov chain. - -""" -controlled_mc(res::DPSolveResult) = res.mc - """ RQ_sigma(ddp, sigma) @@ -1315,7 +1298,13 @@ function _solve!( throw(ArgumentError("method invalid for beta = 1")) end - vals = similar(ddp.R, eltype(v)) # buffer for state-action values + # buffer for state-action values, in the same promoted element type + # as allocated by the 4-arg bellman_operator! method that the greedy + # step used to go through: with mixed-precision models (e.g. Float32 + # rewards with a Float64 beta), a buffer in eltype(v) would round + # near-tied action values to equality before the argmax + vals = similar(ddp.R, promote_type(eltype(ddp.R), eltype(v), + typeof(ddp.beta))) num_iter = 0 for i in 1:max_iter diff --git a/test/test_ddp.jl b/test/test_ddp.jl index 354ad35d..b7e0d0ab 100644 --- a/test/test_ddp.jl +++ b/test/test_ddp.jl @@ -161,6 +161,10 @@ Tests for markov/ddp.jl end @testset "DPSolveResult metadata and construction" begin + # local helper: @inferred needs a function call to check that + # field access on the result infers concretely + _get_mc(r) = r.mc + for _ddp in ddp0_collection for Algo in (VFI, PFI, MPFI) # explicitly passed options are recorded @@ -182,9 +186,8 @@ Tests for markov/ddp.jl @test !ismutable(_res_d) @test isconcretetype(typeof(_res_d)) - # controlled_mc returns the stored chain, type-stably - @test controlled_mc(_res_d) === _res_d.mc - @test (@inferred controlled_mc(_res_d)) === _res_d.mc + # the mc field is concretely typed: access infers + @test (@inferred _get_mc(_res_d)) === _res_d.mc end end end @@ -656,6 +659,23 @@ Tests for markov/ddp.jl @test QuantEcon._max_abs_diff([1.0, 3.0], [0.0, 0.0]) == 3.0 end + @testset "PFI argmax with mixed precision" begin + # the state-action value buffer of the PFI greedy step must + # use the promoted element type, not eltype(v): with Float32 + # rewards and a Float64 beta, the two actions at state 1 are + # near-tied (~1e-11 apart at ~10), distinguishable in Float64 + # but rounded to equality in Float32, where the first-action + # tie-break would flip the policy to [1, 1] + _R_mp = Float32[1.0 -8.0; 2.0 2.0] + _Q_mp = zeros(Float32, 2, 2, 2) + _Q_mp[1, 1, 1] = 1.0 + _Q_mp[1, 2, 2] = 1.0 + _Q_mp[2, 1, 2] = 1.0 + _Q_mp[2, 2, 2] = 1.0 + _ddp_mp = DiscreteDP(_R_mp, _Q_mp, 0.900000000001) + @test solve(_ddp_mp, PFI).sigma == [2, 1] + end + @testset "s_wise_max! argmax with mixed precision" begin # the argmax must be decided at the precision of vals, not of # out: with a Float32 out, 1.0 + 4e-8 rounds back to 1.0f0,