From c6c2be840dc03a68b49e4ee65242595b7175851c Mon Sep 17 00:00:00 2001 From: Jason Jensen Date: Fri, 17 Feb 2023 09:48:00 -0500 Subject: [PATCH 01/32] Accound for equations in ordered dicts --- src/stackedtime/solverdata.jl | 4 ++-- src/steadystate/global.jl | 4 ++-- src/steadystate/presolve.jl | 45 +++++++++++++++++++---------------- src/steadystate/solverdata.jl | 2 +- 4 files changed, 29 insertions(+), 26 deletions(-) diff --git a/src/stackedtime/solverdata.jl b/src/stackedtime/solverdata.jl index eb2edc7..d30a3ce 100644 --- a/src/stackedtime/solverdata.jl +++ b/src/stackedtime/solverdata.jl @@ -366,8 +366,8 @@ function StackedTimeSolverData(model::Model, plan::Plan, fctype::AbstractVector{ # Prep the Jacobian matrix neq = 0 # running counter of equations added to matrix # Model equations are the same for each sim period, just shifted according to t - Jblock = [ti + NT * (ModelBaseEcon._index_of_var(var, unknowns) - 1) for eqn in equations for (var, ti) in keys(eqn.tsrefs)] - Iblock = [i for (i, eqn) in enumerate(equations) for _ in eqn.tsrefs] + Jblock = [ti + NT * (ModelBaseEcon._index_of_var(var, unknowns) - 1) for eqn_pair in equations for (var, ti) in keys(eqn_pair[2].tsrefs)] + Iblock = [i for (i, eqn_pair) in enumerate(equations) for _ in eqn_pair[2].tsrefs] Tblock = -model.maxlag:model.maxlead for t in sim push!(TT, t .+ Tblock) diff --git a/src/steadystate/global.jl b/src/steadystate/global.jl index 2883ba7..c732979 100644 --- a/src/steadystate/global.jl +++ b/src/steadystate/global.jl @@ -59,7 +59,7 @@ end When a model is given, we compute the residual of the entire steady state system. """ -global_SS_R!(resid::AbstractVector{Float64}, point::AbstractVector{Float64}, model::Model) = global_SS_R!(resid, point, ModelBaseEcon.alleqns(model.sstate)) +global_SS_R!(resid::AbstractVector{Float64}, point::AbstractVector{Float64}, model::Model) = global_SS_R!(resid, point, values(ModelBaseEcon.alleqns(model.sstate))) @assert precompile(global_SS_R!, (Vector{Float64}, Vector{Float64}, Model)) """ @@ -110,5 +110,5 @@ end When a model is given, we compute the residual of the entire steady state system. """ -global_SS_RJ(point::Vector{Float64}, model::Model) = global_SS_RJ(point, ModelBaseEcon.alleqns(model.sstate)) +global_SS_RJ(point::Vector{Float64}, model::Model) = global_SS_RJ(point, values(ModelBaseEcon.alleqns(model.sstate))) @assert precompile(global_SS_RJ, (Vector{Float64}, Model)) diff --git a/src/steadystate/presolve.jl b/src/steadystate/presolve.jl index 427575c..99cb74a 100644 --- a/src/steadystate/presolve.jl +++ b/src/steadystate/presolve.jl @@ -5,6 +5,8 @@ # All rights reserved. ################################################################################## +using ModelBaseEcon.OrderedCollections + """ presolve_sstate!(model; ) presolve_sstate!(model, mask, values; ) @@ -18,12 +20,12 @@ loop. ### Arguments * `model` - the Model instance * `mask` - a vector of Bool. Defaults to `model.sstate.mask` - * `values` - a vector of numbers. Defaults to `model.sstate.values` Caller - must specify either both `mask` and `values` or neither of them. `mask[i]` + * `vals` - a vector of numbers. Defaults to `model.sstate.values` Caller + must specify either both `mask` and `vals` or neither of them. `mask[i]` equals `true` if and only if the i-th steady state value has alredy been solved for. -`mask` and `values` are both input and output data to the algorithm. Any values +`mask` and `vals` are both input and output data to the algorithm. Any vals that are successfully presolved are updated in place, and their `mask` enties are set to `true`. @@ -50,44 +52,45 @@ function solve1d(sseqn, vals, ind, ::Val{:newton}, tol = 1e-12, maxiter = 5) return newton1!(sseqn.eval_RJ, vals, ind; tol, maxiter) end -function _presolve_equations!(eqns, mask, values, method, verbose, tol) - eqns_solved = falses(size(eqns)) - eqns_resid = zeros(size(eqns)) +function _presolve_equations!(eqns, mask, vals, method, verbose, tol) + eqns_solved = OrderedDict{Symbol, Bool}(key => false for key in keys(eqns)) + eqns_resid = OrderedDict{Symbol, Float64}(key => 0.0 for key in keys(eqns)) + # eqns_resid = zeros(length(eqns)) retval = false # return value: true if any mask changed, false otherwise updated = true # keep track if any mask changed within the loop - while updated && !all(eqns_solved) + while updated && !all(values(eqns_solved)) updated = false # keep track if anything changes this outer iteration - for (eqn_idx, sseqn) in enumerate(eqns) - eqns_solved[eqn_idx] && continue + for (eqn_key, sseqn) in eqns + eqns_solved[eqn_key] && continue # mask is true for variables that are already solved unsolved = .!mask[sseqn.vinds] nunsolved = sum(unsolved) if nunsolved == 0 # all variables are solved, yet equation is not marked solved. # check if equation is satisfied - eqns_resid[eqn_idx] = R = sseqn.eval_resid(values[sseqn.vinds]) + eqns_resid[eqn_key] = R = sseqn.eval_resid(vals[sseqn.vinds]) # mark it solved either way, but issue a warning if residual is not zero - eqns_solved[eqn_idx] = true + eqns_solved[eqn_key] = true if verbose && abs(R) > 100tol - @warn "Equation $eqn_idx has residual $R:\n $sseqn" + @warn "Equation $eqn_key has residual $R:\n $sseqn" end elseif nunsolved == 1 # only one variable left unsolved. call the 1d solver ind = findall(unsolved)[1] - vals = values[sseqn.vinds] - success = solve1d(sseqn, vals, ind, Val(method), 0.1tol) + _vals = vals[sseqn.vinds] + success = solve1d(sseqn, _vals, ind, Val(method), 0.1tol) if success - if abs(vals[ind]) < 1e-10 - vals[ind] = 0.0 + if abs(_vals[ind]) < 1e-10 + _vals[ind] = 0.0 end - values[sseqn.vinds[ind]] = vals[ind] + vals[sseqn.vinds[ind]] = _vals[ind] retval = updated = mask[sseqn.vinds[ind]] = true if verbose - @info "Presolved equation $eqn_idx for $(sseqn.vsyms[ind]) = $(vals[ind])\n $sseqn" + @info "Presolved equation $eqn_key for $(sseqn.vsyms[ind]) = $(_vals[ind])\n $sseqn" end else if verbose - @info "Failed to presolve $eqn_idx for $(sseqn.vsyms[ind])\n $sseqn" + @info "Failed to presolve $eqn_key for $(sseqn.vsyms[ind])\n $sseqn" end end else @@ -106,10 +109,10 @@ presolve_sstate!(model::Model; kwargs...) = presolve_sstate!(model::Model, mask::AbstractVector{Bool}, values::AbstractVector{Float64}; verbose = model.verbose, tol = model.tol, _1dsolver = :bisect, method = _1dsolver) = _presolve_equations!(ModelBaseEcon.alleqns(model.sstate), mask, values, method, verbose, tol) -presolve_sstate!(eqns::Vector{SteadyStateEquation}, mask::AbstractVector{Bool}, values::AbstractVector{Float64}; +presolve_sstate!(eqns::LittleDict{Symbol,SteadyStateEquation}, mask::AbstractVector{Bool}, values::AbstractVector{Float64}; verbose = false, tol = 1e-12, _1dsolver = :bisect, method = _1dsolver) = _presolve_equations!(eqns, mask, values, method, verbose, tol) @assert precompile(presolve_sstate!, (Model, Vector{Bool}, Vector{Float64})) -@assert precompile(presolve_sstate!, (Vector{SteadyStateEquation}, Vector{Bool}, Vector{Float64})) +# @assert precompile(presolve_sstate!, (LittleDict{Symbol,SteadyStateEquation}, Vector{Bool}, Vector{Float64})) diff --git a/src/steadystate/solverdata.jl b/src/steadystate/solverdata.jl index e716af8..5f8f01d 100644 --- a/src/steadystate/solverdata.jl +++ b/src/steadystate/solverdata.jl @@ -102,7 +102,7 @@ function SolverData(model::Model; presolve::Bool=true, ) # constructor where we're solving all equations for all variables local sstate = model.sstate - local alleqns = ModelBaseEcon.alleqns(sstate) + local alleqns = collect(values(ModelBaseEcon.alleqns(sstate))) local neqns = length(alleqns) local nvars = length(sstate.values) sd = SolverData(copy(model.sstate.values), Vector{Float64}(undef, neqns), # point and resid From 2e1594d66e6407e59d2befaeada1d2c41919f778 Mon Sep 17 00:00:00 2001 From: KristofferC Date: Thu, 9 Mar 2023 17:02:07 +0100 Subject: [PATCH 02/32] proof of concept: allow using Pardiso as a linear solver MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Pardiso is a sparse linear solver that comes with MKL and can have a significant speed improvement over the default solver in julia (UMFPACK). This PR is a proof of concept of using Pardiso in the stacked time solver. Running the FRB US model with the noise it performs around 2x better on my machine compared to the default: ``` Section ncalls time %tot avg alloc %tot avg ───────────────────────────────────────────────────────────────────────────── factorize 8 3.66s 57.9% 458ms 34.7MiB 1.8% 4.34MiB StateSpaceEcon.StackedTimeSolver.USE_PARDISO_FOR_LU = false: Section ncalls time %tot avg alloc %tot avg ───────────────────────────────────────────────────────────────────────────── factorize 8 8.43s 82.3% 1.05s 5.01GiB 73.4% 642MiB ``` --- Project.toml | 2 + src/stackedtime/shockdecomp.jl | 9 ++++- src/stackedtime/simulate.jl | 11 +++++- src/stackedtime/solverdata.jl | 67 ++++++++++++++++++++++++++++------ test/simtests.jl | 3 +- 5 files changed, 76 insertions(+), 16 deletions(-) diff --git a/Project.toml b/Project.toml index 7472c53..f0a8674 100644 --- a/Project.toml +++ b/Project.toml @@ -10,6 +10,7 @@ ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" ModelBaseEcon = "1d16192e-b65e-11ea-11ed-0789cee22d2f" +Pardiso = "46dd5b70-b6fb-5a00-ae2d-e8fea33afaf2" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" @@ -21,6 +22,7 @@ DiffResults = "1.0" ForwardDiff = "0.10" JLD2 = "0.4" ModelBaseEcon = "0.5" +Pardiso = "0.5" Suppressor = "0.2" TimeSeriesEcon = "0.5" julia = "1.7" diff --git a/src/stackedtime/shockdecomp.jl b/src/stackedtime/shockdecomp.jl index 749e077..5163d7c 100644 --- a/src/stackedtime/shockdecomp.jl +++ b/src/stackedtime/shockdecomp.jl @@ -37,7 +37,7 @@ function shockdecomp(m::Model, p::Plan, exog_data::SimData; verbose::Bool=getoption(m, :verbose, false), maxiter::Int=getoption(m, :maxiter, 20), tol::Float64=getoption(m, :tol, 1e-12), - sparse_solver::Function=(A, b) -> A \ b, + sparse_solver::Function=(A, b) -> sparse_solve(A, b), linesearch::Bool=getoption(m, :linesearch, false), fctype::FinalCondition=getoption(m, :fctype, fcgiven), _debug=false @@ -196,7 +196,12 @@ function shockdecomp(m::Model, p::Plan, exog_data::SimData; end SDMAT = Matrix(gdata.J * sparse(SDI, SDJ, -SDV, size(gdata.J, 2), length(exog_inds))) - ldiv!(gdata.J_factorized[], SDMAT) + factor = gdata.J_factorized[] + if factor isa PardisoFactorization + pardiso_solve!(factor, SDMAT) + else + ldiv!(factor, SDMAT) + end end diff --git a/src/stackedtime/simulate.jl b/src/stackedtime/simulate.jl index f4edfc4..4e597f0 100644 --- a/src/stackedtime/simulate.jl +++ b/src/stackedtime/simulate.jl @@ -9,6 +9,13 @@ return model end +function sparse_solve(A, b) + if A isa PardisoFactorization + pardiso_solve!(A, copy(b)) + else + A \ b + end +end """ sim_nr!(x, sd, maxiter, tol, verbose [, sparse_solver [, linesearch]]) @@ -29,7 +36,7 @@ Solve the simulation problem. """ function sim_nr!(x::AbstractArray{Float64}, sd::StackedTimeSolverData, maxiter::Int64, tol::Float64, verbose::Bool, - sparse_solver::Function=(A, b) -> A \ b, linesearch::Bool=false) + sparse_solver::Function=(A, b) -> sparse_solve(A, b), linesearch::Bool=false) for it = 1:maxiter Fx, Jx = stackedtime_RJ(x, x, sd) nFx = norm(Fx, Inf) @@ -106,7 +113,7 @@ function simulate(m::Model, fctype=getoption(m, :fctype, fcgiven), expectation_horizon::Union{Nothing,Int64}=nothing, #= Newton-Raphson options =# - sparse_solver::Function=(A, b) -> A \ b, + sparse_solver::Function=(A, b) -> sparse_solve(A, b), linesearch::Bool=getoption(m, :linesearch, false), warn_maxiter::Bool=getoption(getoption(m, :warn, Options()), :maxiter, false) ) diff --git a/src/stackedtime/solverdata.jl b/src/stackedtime/solverdata.jl index eb2edc7..50b2de6 100644 --- a/src/stackedtime/solverdata.jl +++ b/src/stackedtime/solverdata.jl @@ -5,6 +5,8 @@ # All rights reserved. ################################################################################## +USE_PARDISO_FOR_LU = true + """ StackedTimeSolverData @@ -68,7 +70,7 @@ function var_CiSc end """ assign_fc!(x::Vector, exog::Vector, vind::Int, sd::StackedTimeSolverData, fc::FinalCondition) - + Applying the final condition `fc` for variable with index `vind`. Exogenous data is provided in `exog` and stacked time solver data in `sd`. This function updates the solution vector `x` in place and returns `x`. @@ -84,7 +86,7 @@ function assign_fc! end ####### @inline function var_CiSc(::StackedTimeSolverData, ::ModelVariable, ::FCNone) - # No Jacobian correction + # No Jacobian correction return Dict{Int,Vector{Float64}}() end @@ -96,7 +98,7 @@ end ####### @inline function var_CiSc(::StackedTimeSolverData, ::ModelVariable, ::FCGiven) - # No Jacobian correction + # No Jacobian correction return Dict{Int,Vector{Float64}}() end @@ -111,7 +113,7 @@ end ####### @inline function var_CiSc(::StackedTimeSolverData, ::ModelVariable, ::FCMatchSSLevel) - # No Jacobian correction + # No Jacobian correction return Dict{Int,Vector{Float64}}() end @@ -392,7 +394,7 @@ function StackedTimeSolverData(model::Model, plan::Plan, fctype::AbstractVector{ # We also need the times of the final conditions push!(TT, term) - # + # exog_mask = falses(nunknowns * NT) # Initial conditions are set as exogenous values exog_mask[vec(LI[init, :])] .= true @@ -418,7 +420,7 @@ end """ assign_exog_data!(x::Matrix, exog::Matrix, sd::StackedTimeSolverData) -Assign the exogenous points into `x` according to the plan with which `sd` was created using +Assign the exogenous points into `x` according to the plan with which `sd` was created using exogenous data from `exog`. Also call [`assign_final_condition!`](@ref). !!! warning @@ -434,7 +436,7 @@ end """ assign_final_condition!(x::Matrix, exog::Matrix, sd::StackedTimeSolver) -Assign the final conditions into `x`. The final condition types for the different variables of the model +Assign the final conditions into `x`. The final condition types for the different variables of the model are stored in the the solver data `sd`. `exog` is used for [`fcgiven`](@ref). !!! warning @@ -464,7 +466,7 @@ function stackedtime_R!(res::AbstractArray{Float64,1}, point::AbstractArray{Floa return res end -#= disable +#= disable # this is not necessary with log-transformed variables # we keep it because it might be necessary for other transformations in the future. @@ -571,7 +573,7 @@ function update_CiSc!(x, sd, ::Val{fcnatural}) end end return nothing -end +end =# @@ -615,9 +617,9 @@ function stackedtime_RJ(point::AbstractArray{Float64}, exog_data::AbstractArray{ end try # compute lu decomposition of the active part of J and cache it. - sd.J_factorized[] = sd.factorization == :qr ? qr(JJ) : lu(JJ) + sd.J_factorized[] = sd.factorization == :qr ? qr(JJ) : _factorize(JJ) catch e - if e isa SingularException + if e isa SingularException || e isa Pardiso.PardisoException || e isa Pardiso.PardisoPosDefException @error("The system is underdetermined with the given set of equations and final conditions.") end debugging && return RES, deepcopy(JJ) @@ -627,6 +629,49 @@ function stackedtime_RJ(point::AbstractArray{Float64}, exog_data::AbstractArray{ return RES, sd.J_factorized[] end +using Pardiso + +mutable struct PardisoFactorization + ps::MKLPardisoSolver + J::SparseMatrixCSC +end + +# See https://github.com/JuliaSparse/Pardiso.jl/blob/master/examples/exampleunsym.jl +function pardiso_factorize(JJ) + ps = MKLPardisoSolver() + set_matrixtype!(ps, Pardiso.REAL_NONSYM) + pardisoinit(ps) + fix_iparm!(ps, :N) + # See https://www.intel.com/content/www/us/en/develop/documentation/onemkl-developer-reference-c/top/sparse-solver-routines/onemkl-pardiso-parallel-direct-sparse-solver-iface/pardiso-iparm-parameter.html + set_iparm!(ps, 2, 2) # The parallel (OpenMP) version of the nested dissection algorithm. + JJ_pardiso = get_matrix(ps, JJ, :N) + set_phase!(ps, Pardiso.ANALYSIS) + pardiso(ps, JJ_pardiso, Float64[]) + set_phase!(ps, Pardiso.NUM_FACT) + pardiso(ps, JJ, Float64[]) + pf = PardisoFactorization(ps, JJ_pardiso) + finalizer(pf) do x + set_phase!(x.ps, Pardiso.RELEASE_ALL) + pardiso(x.ps) + end +end + +function pardiso_solve!(ps::PardisoFactorization, x::AbstractArray) + ps, J = ps.ps, ps.J + set_phase!(ps, Pardiso.SOLVE_ITERATIVE_REFINE) + X = similar(x) # Solution is stored in X + pardiso(ps, X, J, x) + copy!(x, X) +end + + +function _factorize(JJ) + if USE_PARDISO_FOR_LU + return pardiso_factorize(JJ) + else + return lu(JJ) + end +end """ assign_update_step!(x::Array, lambda, dx, sd::StackedTimeSolverData) diff --git a/test/simtests.jl b/test/simtests.jl index c97fd27..4103433 100644 --- a/test/simtests.jl +++ b/test/simtests.jl @@ -352,7 +352,8 @@ end p = Plan(m, 3:17) ed = zerodata(m, p) ed .= rand(Float64, size(ed)) - @test_logs (:error, r"The system is underdetermined.*") (@test_throws SingularException simulate(m, p, ed, fctype=fcnatural)) + # @test_throws "The system is underdetermined.*" simulate(m, p, ed, fctype=fcnatural) + sd = StateSpaceEcon.StackedTimeSolver.StackedTimeSolverData(m, p, fcgiven) @test_throws ErrorException StateSpaceEcon.StackedTimeSolver.update_plan!(sd, m, Plan(m, 3:8)) x = rand(Float64, size(ed)) From 68581dcd925be5aa8461ea7e6fb8490d4de9af70 Mon Sep 17 00:00:00 2001 From: KristofferC Date: Fri, 10 Mar 2023 11:37:14 +0100 Subject: [PATCH 03/32] reuse the symbolic factorization if the sparsity pattern has not changed --- src/stackedtime/solverdata.jl | 50 ++++++++++++++++++++++++++--------- 1 file changed, 38 insertions(+), 12 deletions(-) diff --git a/src/stackedtime/solverdata.jl b/src/stackedtime/solverdata.jl index 50b2de6..33540d5 100644 --- a/src/stackedtime/solverdata.jl +++ b/src/stackedtime/solverdata.jl @@ -616,8 +616,16 @@ function stackedtime_RJ(point::AbstractArray{Float64}, exog_data::AbstractArray{ JJ = JAC[:, sd.solve_mask] end try + # Check if sparsity pattern has changed, if not, we can reuse the symbolic factorization + # if Pardiso is used + psf = nothing + if sd.J_factorized[] isa PardisoFactorization + J = sd.J_factorized[].J + same_sparsity_pattern = (J.colptr == JJ.colptr && J.rowval == JJ.rowval) + same_sparsity_pattern && (psf = sd.J_factorized[]) + end # compute lu decomposition of the active part of J and cache it. - sd.J_factorized[] = sd.factorization == :qr ? qr(JJ) : _factorize(JJ) + sd.J_factorized[] = sd.factorization == :qr ? qr(JJ) : _factorize(JJ; psf) catch e if e isa SingularException || e isa Pardiso.PardisoException || e isa Pardiso.PardisoPosDefException @error("The system is underdetermined with the given set of equations and final conditions.") @@ -634,26 +642,44 @@ using Pardiso mutable struct PardisoFactorization ps::MKLPardisoSolver J::SparseMatrixCSC + PardisoFactorization(ps::MKLPardisoSolver) = new(ps) end -# See https://github.com/JuliaSparse/Pardiso.jl/blob/master/examples/exampleunsym.jl -function pardiso_factorize(JJ) +function pardiso_init() ps = MKLPardisoSolver() set_matrixtype!(ps, Pardiso.REAL_NONSYM) pardisoinit(ps) fix_iparm!(ps, :N) # See https://www.intel.com/content/www/us/en/develop/documentation/onemkl-developer-reference-c/top/sparse-solver-routines/onemkl-pardiso-parallel-direct-sparse-solver-iface/pardiso-iparm-parameter.html set_iparm!(ps, 2, 2) # The parallel (OpenMP) version of the nested dissection algorithm. - JJ_pardiso = get_matrix(ps, JJ, :N) - set_phase!(ps, Pardiso.ANALYSIS) - pardiso(ps, JJ_pardiso, Float64[]) - set_phase!(ps, Pardiso.NUM_FACT) - pardiso(ps, JJ, Float64[]) - pf = PardisoFactorization(ps, JJ_pardiso) - finalizer(pf) do x + psf = PardisoFactorization(ps) + finalizer(psf) do x set_phase!(x.ps, Pardiso.RELEASE_ALL) pardiso(x.ps) end + return psf +end + + +# See https://github.com/JuliaSparse/Pardiso.jl/blob/master/examples/exampleunsym.jl +function pardiso_factorize(JJ; psf::Union{Nothing, PardisoFactorization}) + reuse_ps = psf !== nothing + if psf === nothing + psf = pardiso_init() + end + + ps = psf.ps + psf.J = get_matrix(ps, JJ, :N) + + # No need to run the analysis phase if the sparsity pattern is unchanged + # and we reuse the same solver object + if !reuse_ps + set_phase!(ps, Pardiso.ANALYSIS) + pardiso(ps, psf.J, Float64[]) + end + set_phase!(ps, Pardiso.NUM_FACT) + pardiso(ps, psf.J, Float64[]) + return psf end function pardiso_solve!(ps::PardisoFactorization, x::AbstractArray) @@ -665,9 +691,9 @@ function pardiso_solve!(ps::PardisoFactorization, x::AbstractArray) end -function _factorize(JJ) +function _factorize(JJ; psf) if USE_PARDISO_FOR_LU - return pardiso_factorize(JJ) + return pardiso_factorize(JJ; psf) else return lu(JJ) end From 36ecc5f3ab32d1567d15ef28391303bf95b17328 Mon Sep 17 00:00:00 2001 From: Boyan Bejanov Date: Wed, 5 Apr 2023 09:37:42 -0400 Subject: [PATCH 04/32] version 0.4.2 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 9174697..8e9bb84 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "StateSpaceEcon" uuid = "e4c825b0-b65c-11ea-0b5a-6176b64e7b7f" authors = ["Atai Akunov ", "Boyan Bejanov ", "Nicholas Labelle St-Pierre "] -version = "0.4.1" +version = "0.4.2" [deps] DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab" From 1cb740e5d36f45e531ae324d4dfdc99654f5fe0a Mon Sep 17 00:00:00 2001 From: Jason Jensen Date: Wed, 5 Apr 2023 15:09:31 -0400 Subject: [PATCH 05/32] handle equation dictionaries --- src/steadystate/diagnose.jl | 2 +- src/steadystate/global.jl | 8 ++++---- src/steadystate/initial.jl | 2 +- src/steadystate/sssolve.jl | 2 +- 4 files changed, 7 insertions(+), 7 deletions(-) diff --git a/src/steadystate/diagnose.jl b/src/steadystate/diagnose.jl index ac1d07b..ba96b5f 100644 --- a/src/steadystate/diagnose.jl +++ b/src/steadystate/diagnose.jl @@ -27,7 +27,7 @@ Standard options (default values from `model.options`) """ function check_sstate(model::Model; verbose::Bool = model.options.verbose, tol::Float64 = model.options.tol) # if parameters changed, see that constraints use the latest - for eqn in model.sstate.constraints + for eqn in values(model.sstate.constraints) ModelBaseEcon._update_eqn_params!(eqn.eval_resid, model.parameters) end R, J = global_SS_RJ(model.sstate.values, model) diff --git a/src/steadystate/global.jl b/src/steadystate/global.jl index c732979..33e8923 100644 --- a/src/steadystate/global.jl +++ b/src/steadystate/global.jl @@ -14,11 +14,11 @@ function inadmissible_error(eqind, eqn, point, val) varsargs = (; (v => a for (v, a) in zip(vars, args))...) params = (; _getparams(eqn)...) if typeof(val) <: AbstractArray - @error "Singular gradient in steady state equation $eqind: $eqn" params varsargs val - error("Singular gradient in steady state equation $eqind: $eqn") + @error "Singular gradient in steady state equation $(eqn.name) => $eqn" params varsargs val + error("Singular gradient in steady state equation $(eqn.name) => $eqn") else - @error "Inadmissible point in steady state equation $eqind: $eqn" params varsargs val - error("Inadmissible point in steady state equation $eqind: $eqn") + @error "Inadmissible point in steady state equation $(eqn.name)) => $eqn" params varsargs val + error("Inadmissible point in steady state equation $(eqn.name) => $eqn") end end diff --git a/src/steadystate/initial.jl b/src/steadystate/initial.jl index f38bea0..4dc104b 100644 --- a/src/steadystate/initial.jl +++ b/src/steadystate/initial.jl @@ -75,7 +75,7 @@ function _do_update_auxvars_presolve!(model::Model, verbose::Bool, method::Symbo end # presolve only the steady state constraints (to apply the values of exog variables) if !isempty(ss.constraints) - for eqn in ss.constraints + for eqn in values(ss.constraints) ModelBaseEcon._update_eqn_params!(eqn.eval_resid, model.parameters) end presolve_sstate!(ss.constraints, ss.mask, ss.values; model.tol, verbose, method) diff --git a/src/steadystate/sssolve.jl b/src/steadystate/sssolve.jl index 38c25cf..aebae15 100644 --- a/src/steadystate/sssolve.jl +++ b/src/steadystate/sssolve.jl @@ -43,7 +43,7 @@ function sssolve!(model::Model; # !! make sure constraints use the latest parameter values # !! before calling SolverData constructor, because it computes the initial residual - for eqn in ss.constraints + for eqn in values(ss.constraints) ModelBaseEcon._update_eqn_params!(eqn.eval_resid, model.parameters) end From 50f27d167fab11058e1c11385af960fd48a8fbad Mon Sep 17 00:00:00 2001 From: Jason Jensen Date: Thu, 6 Apr 2023 11:38:17 -0400 Subject: [PATCH 06/32] Use OrderedDict for steadystate --- src/firstorder/solve.jl | 2 +- src/steadystate/presolve.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/firstorder/solve.jl b/src/firstorder/solve.jl index cdfd407..e4f45d9 100644 --- a/src/firstorder/solve.jl +++ b/src/firstorder/solve.jl @@ -50,7 +50,7 @@ function VarMaps(model::Model) ex_vars = Tuple{Symbol,Int}[] for mvar in model.allvars var = mvar.name - lags, leads = extrema(lag for eqn in model.equations + lags, leads = extrema(lag for eqn in values(model.equations) for (name, lag) in keys(eqn.tsrefs) if name == var) if isexog(mvar) || isshock(mvar) diff --git a/src/steadystate/presolve.jl b/src/steadystate/presolve.jl index 99cb74a..746b9e2 100644 --- a/src/steadystate/presolve.jl +++ b/src/steadystate/presolve.jl @@ -109,7 +109,7 @@ presolve_sstate!(model::Model; kwargs...) = presolve_sstate!(model::Model, mask::AbstractVector{Bool}, values::AbstractVector{Float64}; verbose = model.verbose, tol = model.tol, _1dsolver = :bisect, method = _1dsolver) = _presolve_equations!(ModelBaseEcon.alleqns(model.sstate), mask, values, method, verbose, tol) -presolve_sstate!(eqns::LittleDict{Symbol,SteadyStateEquation}, mask::AbstractVector{Bool}, values::AbstractVector{Float64}; +presolve_sstate!(eqns::OrderedDict{Symbol,SteadyStateEquation}, mask::AbstractVector{Bool}, values::AbstractVector{Float64}; verbose = false, tol = 1e-12, _1dsolver = :bisect, method = _1dsolver) = _presolve_equations!(eqns, mask, values, method, verbose, tol) From 4a10ffb3e10241685cb13290e89fa856b38d2a20 Mon Sep 17 00:00:00 2001 From: Jason Jensen Date: Thu, 6 Apr 2023 11:38:34 -0400 Subject: [PATCH 07/32] run tests on model copies --- test/dynss.jl | 2 +- test/sim_fo.jl | 28 ++++++++++++++-------------- test/simtests.jl | 2 +- test/sstests.jl | 15 ++++++++------- 4 files changed, 24 insertions(+), 23 deletions(-) diff --git a/test/dynss.jl b/test/dynss.jl index bd9745a..356f636 100644 --- a/test/dynss.jl +++ b/test/dynss.jl @@ -95,7 +95,7 @@ end @testset "dynss+" begin # tests to make sure that changes in parameter values are applied to steady state constraints correctly - local m = PC.model + local m = deepcopy(PC.model) m.sstate.values .= 0.0 @test begin sssolve!(m, presolve=false) diff --git a/test/sim_fo.jl b/test/sim_fo.jl index 9961a00..87a845d 100644 --- a/test/sim_fo.jl +++ b/test/sim_fo.jl @@ -58,7 +58,7 @@ end end @testset "M.fo.unant" begin - run_fo_unant_tests(M.model) + run_fo_unant_tests(deepcopy(M.model)) end module R @@ -79,23 +79,23 @@ end @initialize model end @testset "R.fo.unant" begin - run_fo_unant_tests(R.model) + run_fo_unant_tests(deepcopy(R.model)) end @using_example E2 @testset "E2.fo.unant" begin - run_fo_unant_tests(E2.model) + run_fo_unant_tests(deepcopy(E2.model)) end @using_example E3 @testset "E3.fo.unant" begin - run_fo_unant_tests(E3.model) + run_fo_unant_tests(deepcopy(E3.model)) end @using_example E6 @testset "E6.fo.unant" begin - let m = E6.model + let m = deepcopy(E6.model) # set slopes to 0, otherwise we're not allowed to linearize m.p_dly = 0 m.p_dlp = 0 @@ -130,7 +130,7 @@ function test_shockdecomp_firstorder(m, rng=1U:20U, fctype=fcslope) end @testset "shkdcmp.fo" begin - for m in (M.model, R.model) + for m in (deepcopy(M.model), deepcopy(R.model)) m.rho = 0.6 empty!(m.sstate.constraints) test_shockdecomp_firstorder(m) @@ -138,11 +138,11 @@ end @steadystate m x = 6 test_shockdecomp_firstorder(m, 1U:20U, fclevel) end - for m in (E2.model, E3.model) + for m in (deepcopy(E2.model), deepcopy(E3.model)) empty!(m.sstate.constraints) test_shockdecomp_firstorder(m, 1U:500U) end - let m = E6.model + let m = deepcopy(E6.model) # set slopes to 0, otherwise we're not allowed to linearize m.p_dly = 0 m.p_dlp = 0 @@ -188,7 +188,7 @@ end @testset "inidcmp.fo" begin - for m in (M.model, R.model) + for m in (deepcopy(M.model), deepcopy(R.model)) m.rho = 0.6 empty!(m.sstate.constraints) test_initdecomp_firstorder(m) @@ -196,11 +196,11 @@ end @steadystate m x = 6 test_initdecomp_firstorder(m) end - for m in (E2.model, E3.model) + for m in (deepcopy(E2.model), deepcopy(E3.model)) empty!(m.sstate.constraints) test_initdecomp_firstorder(m) end - let m = E6.model + let m = deepcopy(E6.model) # set slopes to 0, otherwise we're not allowed to linearize m.p_dly = 0 m.p_dlp = 0 @@ -313,7 +313,7 @@ function test_initdecomp_stackedtime(m; nonlin=!m.linear, rng=2001Q1:2024Q4, fct end @testset "inidcmp.st" begin - for m in (M.model, R.model) + for m in (deepcopy(M.model), deepcopy(R.model)) m.rho = 0.6 empty!(m.sstate.constraints) test_initdecomp_stackedtime(m) @@ -321,11 +321,11 @@ end @steadystate m x = 6 test_initdecomp_stackedtime(m, fctype=fclevel) end - for m in (E2.model, E3.model) + for m in (deepcopy(E2.model), deepcopy(E3.model)) empty!(m.sstate.constraints) test_initdecomp_stackedtime(m) end - let m = E6.model + let m = deepcopy(E6.model) # set slopes to 0, otherwise we're not allowed to linearize m.p_dly = 0 m.p_dlp = 0 diff --git a/test/simtests.jl b/test/simtests.jl index c97fd27..ac00e36 100644 --- a/test/simtests.jl +++ b/test/simtests.jl @@ -409,7 +409,7 @@ end end @testset "new.unant" begin - m = E1.model + m = deepcopy(E1.model) m.α = m.β = 0.5 test_rng = 20Q1:22Q1 diff --git a/test/sstests.jl b/test/sstests.jl index 4a598c4..43553ba 100644 --- a/test/sstests.jl +++ b/test/sstests.jl @@ -5,9 +5,10 @@ # All rights reserved. ################################################################################## -empty!(E1.model.sstate.constraints) + @testset "E1.sstate" begin let m = deepcopy(E1.model) + empty!(m.sstate.constraints) m.α = 0.5 m.β = 1.0 - m.α clear_sstate!(m) @@ -51,9 +52,9 @@ empty!(E1.model.sstate.constraints) end end -empty!(E1.model.sstate.constraints) @testset "E1.sstate, auto" begin let m = deepcopy(E1.model) + empty!(m.sstate.constraints) m.α = 0.5 m.β = 1.0 - m.α clear_sstate!(m) @@ -91,9 +92,9 @@ empty!(E1.model.sstate.constraints) end end -empty!(E1.model.sstate.constraints) @testset "E1.sstate, lm" begin let m = deepcopy(E1.model) + empty!(m.sstate.constraints) m.α = 0.5 m.β = 1.0 - m.α clear_sstate!(m) @@ -131,9 +132,9 @@ empty!(E1.model.sstate.constraints) end end -empty!(E2.model.sstate.constraints) @testset "E2.sstate" begin let m = deepcopy(E2.model) + empty!(m.sstate.constraints) clear_sstate!(m) sssolve!(m) @test check_sstate(m) == 0 @@ -241,7 +242,7 @@ end end @testset "SSTEST" begin - let m = SSTEST.model + let m = deepcopy(SSTEST.model) @test sum(issteady.(m.allvars)) == 2 clear_sstate!(m) @test issssolved(m) @@ -318,8 +319,8 @@ end end @testset "presolve, ssZeroSlope" begin - empty!(E1.model.sstate.constraints) let m = deepcopy(E1.model) + empty!(m.sstate.constraints) m.α = 0.5 m.β = 1.0 - m.α clear_sstate!(m) @@ -328,8 +329,8 @@ end @test m.sstate.mask == [false, false, true, true] end - empty!(E1.model.sstate.constraints) let m = deepcopy(E1.model) + empty!(m.sstate.constraints) m.α = 0.5 m.β = 1.0 - m.α m.flags.ssZeroSlope = true From 984eacf7c906491a979deaf2dd862a82cc8521a1 Mon Sep 17 00:00:00 2001 From: Boyan Bejanov Date: Thu, 6 Apr 2023 15:56:13 -0400 Subject: [PATCH 08/32] using TimerOutputs for profiling --- Project.toml | 3 ++- src/StackedTimeSolver.jl | 5 +++++ src/stackedtime/shockdecomp.jl | 2 +- src/stackedtime/simulate.jl | 10 +++++----- src/stackedtime/solverdata.jl | 24 ++++++++++++------------ src/stackedtime/sparse.jl | 15 +++++++++++++++ 6 files changed, 40 insertions(+), 19 deletions(-) create mode 100644 src/stackedtime/sparse.jl diff --git a/Project.toml b/Project.toml index fd908d0..bde8608 100644 --- a/Project.toml +++ b/Project.toml @@ -16,13 +16,14 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb" TimeSeriesEcon = "8b6756d2-c55c-11ea-2998-5f67ea17da60" +TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" [compat] DiffResults = "1.0" ForwardDiff = "0.10" JLD2 = "0.4" -Pardiso = "0.5" ModelBaseEcon = "0.5.2" +Pardiso = "0.5" Suppressor = "0.2" TimeSeriesEcon = "0.5.1" julia = "1.7" diff --git a/src/StackedTimeSolver.jl b/src/StackedTimeSolver.jl index b0433c4..0a9f5d2 100644 --- a/src/StackedTimeSolver.jl +++ b/src/StackedTimeSolver.jl @@ -23,6 +23,9 @@ using Printf using ModelBaseEcon using TimeSeriesEcon +using TimerOutputs +const timer = TimerOutput() + using ..Plans import ..steadystatearray @@ -40,6 +43,8 @@ import ModelBaseEcon.hassolverdata import ModelBaseEcon.getsolverdata import ModelBaseEcon.setsolverdata! +include("stackedtime/sparse.jl") + include("stackedtime/fctypes.jl") include("stackedtime/misc.jl") include("stackedtime/solverdata.jl") diff --git a/src/stackedtime/shockdecomp.jl b/src/stackedtime/shockdecomp.jl index 5163d7c..dd8be3c 100644 --- a/src/stackedtime/shockdecomp.jl +++ b/src/stackedtime/shockdecomp.jl @@ -67,7 +67,7 @@ function shockdecomp(m::Model, p::Plan, exog_data::SimData; end # prepare the stacked-time - gdata = StackedTimeSolverData(m, p, fctype, variant) + gdata = @timeit_debug timer "StackedTimeSolverData" StackedTimeSolverData(m, p, fctype, variant) # check the residual. why are we doing this? we know it's 0! shocked = copy(exog_data) diff --git a/src/stackedtime/simulate.jl b/src/stackedtime/simulate.jl index 9c39f2e..7fd3a66 100644 --- a/src/stackedtime/simulate.jl +++ b/src/stackedtime/simulate.jl @@ -9,7 +9,7 @@ return model end -function sparse_solve(A, b) +@timeit_debug timer function sparse_solve(A, b) if A isa PardisoFactorization pardiso_solve!(A, copy(b)) else @@ -34,7 +34,7 @@ Solve the simulation problem. do this at each iteration. Default is `false`. """ -function sim_nr!(x::AbstractArray{Float64}, sd::StackedTimeSolverData, +@timeit_debug timer function sim_nr!(x::AbstractArray{Float64}, sd::StackedTimeSolverData, maxiter::Int64, tol::Float64, verbose::Bool, sparse_solver::Function=(A, b) -> sparse_solve(A, b), linesearch::Bool=false) for it = 1:maxiter @@ -201,14 +201,14 @@ function simulate(m::Model, if p_unant.range != p_ant.range error("Anticipated and unanticipated ranges don't match.") end - if size(exog_unant) != size(exog_ant) - error("Anticipated and unanticipated data don't match.") - end if deviation_unant exog_unant[:, logvars] .*= baseline[:, logvars] exog_unant[:, .!logvars] .+= baseline[:, .!logvars] end exog_unant = ModelBaseEcon.update_auxvars(transform(exog_unant, m), m) + if size(exog_unant) != size(exog_ant) + error("Anticipated and unanticipated data don't match.") + end else #=== prepare unanticipated data and plan (backward compatibility) ===# p_unant = p_ant diff --git a/src/stackedtime/solverdata.jl b/src/stackedtime/solverdata.jl index 071bc4d..4cd8ee8 100644 --- a/src/stackedtime/solverdata.jl +++ b/src/stackedtime/solverdata.jl @@ -5,7 +5,7 @@ # All rights reserved. ################################################################################## -USE_PARDISO_FOR_LU = true +USE_PARDISO_FOR_LU = false """ StackedTimeSolverData @@ -321,7 +321,7 @@ function make_BI(JMAT::SparseMatrixCSC{Float64,Int}, II::AbstractVector{<:Abstra end StackedTimeSolverData(model::Model, plan::Plan, fctype::FinalCondition, variant::Symbol=model.options.variant) = StackedTimeSolverData(model, plan, setfc(model, fctype), variant) -function StackedTimeSolverData(model::Model, plan::Plan, fctype::AbstractVector{FinalCondition}, variant::Symbol=model.options.variant) +@timeit_debug timer function StackedTimeSolverData(model::Model, plan::Plan, fctype::AbstractVector{FinalCondition}, variant::Symbol=model.options.variant) evaldata = getevaldata(model, variant) var_to_idx = ModelBaseEcon.get_var_to_idx(model) @@ -457,13 +457,13 @@ end Compute the residual of the stacked time system at the given `point`. R is updated in place and returned. """ -function stackedtime_R!(res::AbstractArray{Float64,1}, point::AbstractArray{Float64}, exog_data::AbstractArray{Float64}, sd::StackedTimeSolverData) +@timeit_debug timer function stackedtime_R!(res::AbstractArray{Float64,1}, point::AbstractArray{Float64}, exog_data::AbstractArray{Float64}, sd::StackedTimeSolverData) @assert size(point) == size(exog_data) == (sd.NT, sd.NU) # point = reshape(point, sd.NT, sd.NU) # exog_data = reshape(exog_data, sd.NT, sd.NU) @assert(length(res) == size(sd.J, 1), "Length of residual vector doesn't match.") for (ii, tt) in zip(sd.II, sd.TT) - eval_R!(view(res, ii), point[tt, :], sd.evaldata) + @timeit_debug timer "eval_R!" eval_R!(view(res, ii), point[tt, :], sd.evaldata) end return res end @@ -585,7 +585,7 @@ end Compute the residual and Jacobian of the stacked time system at the given `point`. """ -function stackedtime_RJ(point::AbstractArray{Float64}, exog_data::AbstractArray{Float64}, sd::StackedTimeSolverData; +@timeit_debug timer function stackedtime_RJ(point::AbstractArray{Float64}, exog_data::AbstractArray{Float64}, sd::StackedTimeSolverData; debugging=false) nunknowns = sd.NU @assert size(point) == size(exog_data) == (sd.NT, nunknowns) @@ -599,12 +599,12 @@ function stackedtime_RJ(point::AbstractArray{Float64}, exog_data::AbstractArray{ if haveJ # update only RES for i = 1:length(sd.BI) - eval_R!(view(RES, sd.II[i]), point[sd.TT[i], :], sd.evaldata) + @timeit_debug timer "eval_R!" eval_R!(view(RES, sd.II[i]), point[sd.TT[i], :], sd.evaldata) end else # update both RES and JAC for i = 1:length(sd.BI) - R, J = eval_RJ(point[sd.TT[i], :], sd.evaldata) + @timeit_debug timer "eval_RJ" R, J = eval_RJ(point[sd.TT[i], :], sd.evaldata) RES[sd.II[i]] .= R JAC.nzval[sd.BI[i]] .= J.nzval end @@ -647,7 +647,7 @@ mutable struct PardisoFactorization PardisoFactorization(ps::MKLPardisoSolver) = new(ps) end -function pardiso_init() +@timeit_debug timer function pardiso_init() ps = MKLPardisoSolver() set_matrixtype!(ps, Pardiso.REAL_NONSYM) pardisoinit(ps) @@ -664,7 +664,7 @@ end # See https://github.com/JuliaSparse/Pardiso.jl/blob/master/examples/exampleunsym.jl -function pardiso_factorize(JJ; psf::Union{Nothing, PardisoFactorization}) +@timeit_debug timer function pardiso_factorize(JJ; psf::Union{Nothing,PardisoFactorization}) reuse_ps = psf !== nothing if psf === nothing psf = pardiso_init() @@ -684,7 +684,7 @@ function pardiso_factorize(JJ; psf::Union{Nothing, PardisoFactorization}) return psf end -function pardiso_solve!(ps::PardisoFactorization, x::AbstractArray) +@timeit_debug timer function pardiso_solve!(ps::PardisoFactorization, x::AbstractArray) ps, J = ps.ps, ps.J set_phase!(ps, Pardiso.SOLVE_ITERATIVE_REFINE) X = similar(x) # Solution is stored in X @@ -693,8 +693,8 @@ function pardiso_solve!(ps::PardisoFactorization, x::AbstractArray) end -function _factorize(JJ; psf) - if USE_PARDISO_FOR_LU +@timeit_debug timer function _factorize(JJ; psf) + if USE_PARDISO_FOR_LU return pardiso_factorize(JJ; psf) else return lu(JJ) diff --git a/src/stackedtime/sparse.jl b/src/stackedtime/sparse.jl new file mode 100644 index 0000000..f94ffc4 --- /dev/null +++ b/src/stackedtime/sparse.jl @@ -0,0 +1,15 @@ +################################################################################## +# This file is part of StateSpaceEcon.jl +# BSD 3-Clause License +# Copyright (c) 2020-2023, Bank of Canada +# All rights reserved. +################################################################################## + + +### API + +sf_factorize(A::SparseMatrixCSC)::Factorization = lu(A)::Factorization +sf_update!(F::Factorization, A::SparseMatrixCSC) = lu!(F, A) +sf_solve!(F::Factorization, b) = ldiv!(F, b) +sf_solve(F::Factorization, b) = (cb = copy(b); ldiv!(F, cb); cb) + From 0b9f3d8fbb272e872e1c2ac707b3f3f3f283b95a Mon Sep 17 00:00:00 2001 From: Boyan Bejanov Date: Thu, 6 Apr 2023 21:12:47 -0400 Subject: [PATCH 09/32] sparse factorization [wip] simulate works shockdecomp broken stoch_simulate broken --- Project.toml | 1 + src/StackedTimeSolver.jl | 6 ++ src/stackedtime/shockdecomp.jl | 3 +- src/stackedtime/simulate.jl | 29 ++++---- src/stackedtime/solverdata.jl | 36 +++++----- src/stackedtime/sparse.jl | 110 ++++++++++++++++++++++++++++-- src/stackedtime/stoch_simulate.jl | 3 +- 7 files changed, 150 insertions(+), 38 deletions(-) diff --git a/Project.toml b/Project.toml index bde8608..c2dce00 100644 --- a/Project.toml +++ b/Project.toml @@ -14,6 +14,7 @@ Pardiso = "46dd5b70-b6fb-5a00-ae2d-e8fea33afaf2" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" +SuiteSparse = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9" Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb" TimeSeriesEcon = "8b6756d2-c55c-11ea-2998-5f67ea17da60" TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" diff --git a/src/StackedTimeSolver.jl b/src/StackedTimeSolver.jl index 0a9f5d2..6e61edf 100644 --- a/src/StackedTimeSolver.jl +++ b/src/StackedTimeSolver.jl @@ -43,6 +43,12 @@ import ModelBaseEcon.hassolverdata import ModelBaseEcon.getsolverdata import ModelBaseEcon.setsolverdata! +using Pardiso + + +using SuiteSparse +using SuiteSparse.UMFPACK + include("stackedtime/sparse.jl") include("stackedtime/fctypes.jl") diff --git a/src/stackedtime/shockdecomp.jl b/src/stackedtime/shockdecomp.jl index dd8be3c..43450bf 100644 --- a/src/stackedtime/shockdecomp.jl +++ b/src/stackedtime/shockdecomp.jl @@ -37,7 +37,6 @@ function shockdecomp(m::Model, p::Plan, exog_data::SimData; verbose::Bool=getoption(m, :verbose, false), maxiter::Int=getoption(m, :maxiter, 20), tol::Float64=getoption(m, :tol, 1e-12), - sparse_solver::Function=(A, b) -> sparse_solve(A, b), linesearch::Bool=getoption(m, :linesearch, false), fctype::FinalCondition=getoption(m, :fctype, fcgiven), _debug=false @@ -75,7 +74,7 @@ function shockdecomp(m::Model, p::Plan, exog_data::SimData; stackedtime_R!(res_shocked, shocked, shocked, gdata) # Run the "shocked" simulation with the given exogenous data. if norm(res_shocked, Inf) > tol assign_exog_data!(shocked, exog_data, gdata) - sim_nr!(shocked, gdata, maxiter, tol, verbose, sparse_solver, linesearch) + sim_nr!(shocked, gdata, maxiter, tol, verbose, linesearch) end # We need the Jacobian matrix evaluated at the control solution. diff --git a/src/stackedtime/simulate.jl b/src/stackedtime/simulate.jl index 7fd3a66..6013aac 100644 --- a/src/stackedtime/simulate.jl +++ b/src/stackedtime/simulate.jl @@ -9,6 +9,7 @@ return model end +#= @timeit_debug timer function sparse_solve(A, b) if A isa PardisoFactorization pardiso_solve!(A, copy(b)) @@ -16,8 +17,10 @@ end A \ b end end +=# + """ - sim_nr!(x, sd, maxiter, tol, verbose [, sparse_solver [, linesearch]]) + sim_nr!(x, sd, maxiter, tol, verbose [, linesearch]) Solve the simulation problem. * `x` - the array of data for the simulation. All initial, final and exogenous @@ -27,18 +30,15 @@ Solve the simulation problem. * `maxiter` - maximum number of iterations. * `tol` - desired accuracy. * `verbose` - whether or not to print progress information. - * `sparse_solver` (optional) - a function called to solve the linear system A - x = b for x. Defaults to A\\b * `linesearch::Bool` - When `true` the Newton-Raphson is modified to include a search along the descent direction for a sufficient decrease in f. It will do this at each iteration. Default is `false`. """ @timeit_debug timer function sim_nr!(x::AbstractArray{Float64}, sd::StackedTimeSolverData, - maxiter::Int64, tol::Float64, verbose::Bool, - sparse_solver::Function=(A, b) -> sparse_solve(A, b), linesearch::Bool=false) + maxiter::Int64, tol::Float64, verbose::Bool, linesearch::Bool=false) for it = 1:maxiter - Fx, Jx = stackedtime_RJ(x, x, sd) + Fx = stackedtime_R!(Vector{Float64}(undef, size(sd.J, 1)), x, x, sd) nFx = norm(Fx, Inf) if nFx < tol if verbose @@ -46,7 +46,9 @@ Solve the simulation problem. end return true end - Δx = sparse_solver(Jx, Fx) + Δx, Jx = stackedtime_RJ(x, x, sd) + Δx = sf_solve!(Jx, Δx) + λ = 1.0 if linesearch nf = norm(Fx) @@ -113,7 +115,6 @@ function simulate(m::Model, fctype=getoption(m, :fctype, fcgiven), expectation_horizon::Union{Nothing,Int64}=nothing, #= Newton-Raphson options =# - sparse_solver::Function=(A, b) -> sparse_solve(A, b), linesearch::Bool=getoption(m, :linesearch, false), warn_maxiter::Bool=getoption(getoption(m, :warn, Options()), :maxiter, false) ) @@ -168,7 +169,7 @@ function simulate(m::Model, if verbose @info "Simulating $(p_ant.range[1 + m.maxlag:NT - m.maxlead])" # anticipate gdata.FC end - converged = sim_nr!(x, gdata, maxiter, tol, verbose, sparse_solver, linesearch) + converged = sim_nr!(x, gdata, maxiter, tol, verbose, linesearch) if warn_maxiter && !converged @warn("Newton-Raphson reached maximum number of iterations (`maxiter`).") end @@ -247,7 +248,7 @@ function simulate(m::Model, if verbose @info "Simulating $(p_ant.range[t:T]) with $((itol, imaxiter))" # anticipate expectation_horizon gdata.FC end - converged = sim_nr!(xx, gdata, imaxiter, itol, verbose, sparse_solver, linesearch) + converged = sim_nr!(xx, gdata, imaxiter, itol, verbose, linesearch) if warn_maxiter && !converged @warn("Newton-Raphson reached maximum number of iterations (`imaxiter`).") end @@ -260,7 +261,7 @@ function simulate(m::Model, if verbose @info "Simulating $(p_ant.range[t:T]) with $((tol, maxiter))" # anticipate expectation_horizon gdata.FC end - converged = sim_nr!(xx, gdata, maxiter, tol, verbose, sparse_solver, linesearch) + converged = sim_nr!(xx, gdata, maxiter, tol, verbose, linesearch) if warn_maxiter && !converged @warn("Newton-Raphson reached maximum number of iterations (`maxiter`).") end @@ -294,7 +295,7 @@ function simulate(m::Model, if verbose @info "Simulating $(p_ant.range[t:T])" # anticipate expectation_horizon sdata.FC end - converged = sim_nr!(xx, sdata, maxiter, tol, verbose, sparse_solver, linesearch) + converged = sim_nr!(xx, sdata, maxiter, tol, verbose, linesearch) if warn_maxiter && !converged @warn("Newton-Raphson reached maximum number of iterations (`maxiter`).") end @@ -342,7 +343,7 @@ function simulate(m::Model, if verbose @info("Simulating $(p_ant.range[t] .+ (0:expectation_horizon - 1))") # anticipate expectation_horizon sdata.FC end - converged = sim_nr!(xx, sdata, 5, sqrt(tol), verbose, sparse_solver, linesearch) + converged = sim_nr!(xx, sdata, 5, sqrt(tol), verbose, linesearch) if warn_maxiter && !converged @warn("Newton-Raphson reached maximum number of iterations (`maxiter`).") end @@ -366,7 +367,7 @@ function simulate(m::Model, if verbose @info "Simulating $(p_ant.range[last_t + 1:T])" # anticipate expectation_horizon sdata.FC end - converged = sim_nr!(xx, sdata, maxiter, tol, verbose, sparse_solver, linesearch) + converged = sim_nr!(xx, sdata, maxiter, tol, verbose, linesearch) if warn_maxiter && !converged @warn("Newton-Raphson reached maximum number of iterations (`maxiter`).") end diff --git a/src/stackedtime/solverdata.jl b/src/stackedtime/solverdata.jl index 4cd8ee8..62dd6d3 100644 --- a/src/stackedtime/solverdata.jl +++ b/src/stackedtime/solverdata.jl @@ -5,8 +5,6 @@ # All rights reserved. ################################################################################## -USE_PARDISO_FOR_LU = false - """ StackedTimeSolverData @@ -321,7 +319,7 @@ function make_BI(JMAT::SparseMatrixCSC{Float64,Int}, II::AbstractVector{<:Abstra end StackedTimeSolverData(model::Model, plan::Plan, fctype::FinalCondition, variant::Symbol=model.options.variant) = StackedTimeSolverData(model, plan, setfc(model, fctype), variant) -@timeit_debug timer function StackedTimeSolverData(model::Model, plan::Plan, fctype::AbstractVector{FinalCondition}, variant::Symbol=model.options.variant) +@timeit_debug timer function StackedTimeSolverData(model::Model, plan::Plan, fctype::AbstractVector{FinalCondition}, variant::Symbol=model.options.variant) evaldata = getevaldata(model, variant) var_to_idx = ModelBaseEcon.get_var_to_idx(model) @@ -414,7 +412,7 @@ StackedTimeSolverData(model::Model, plan::Plan, fctype::FinalCondition, variant: need_SS ? transform(steadystatearray(model, plan), model) : zeros(0, 0), islog.(model.allvars) .| isneglog.(model.allvars), islin.(model.allvars), evaldata, exog_mask, fc_mask, solve_mask, - Ref{Any}(nothing), getoption(model, :factorization, :lu)) + Ref{Any}(nothing), getoption(model, :factorization, :default)) return update_plan!(sd, model, plan; changed=true) end @@ -457,7 +455,7 @@ end Compute the residual of the stacked time system at the given `point`. R is updated in place and returned. """ -@timeit_debug timer function stackedtime_R!(res::AbstractArray{Float64,1}, point::AbstractArray{Float64}, exog_data::AbstractArray{Float64}, sd::StackedTimeSolverData) +@timeit_debug timer function stackedtime_R!(res::AbstractArray{Float64,1}, point::AbstractArray{Float64}, exog_data::AbstractArray{Float64}, sd::StackedTimeSolverData) @assert size(point) == size(exog_data) == (sd.NT, sd.NU) # point = reshape(point, sd.NT, sd.NU) # exog_data = reshape(exog_data, sd.NT, sd.NU) @@ -585,7 +583,7 @@ end Compute the residual and Jacobian of the stacked time system at the given `point`. """ -@timeit_debug timer function stackedtime_RJ(point::AbstractArray{Float64}, exog_data::AbstractArray{Float64}, sd::StackedTimeSolverData; +@timeit_debug timer function stackedtime_RJ(point::AbstractArray{Float64}, exog_data::AbstractArray{Float64}, sd::StackedTimeSolverData; debugging=false) nunknowns = sd.NU @assert size(point) == size(exog_data) == (sd.NT, nunknowns) @@ -618,16 +616,8 @@ Compute the residual and Jacobian of the stacked time system at the given JJ = JAC[:, sd.solve_mask] end try - # Check if sparsity pattern has changed, if not, we can reuse the symbolic factorization - # if Pardiso is used - psf = nothing - if sd.J_factorized[] isa PardisoFactorization - J = sd.J_factorized[].J - same_sparsity_pattern = (J.colptr == JJ.colptr && J.rowval == JJ.rowval) - same_sparsity_pattern && (psf = sd.J_factorized[]) - end - # compute lu decomposition of the active part of J and cache it. - sd.J_factorized[] = sd.factorization == :qr ? qr(JJ) : _factorize(JJ; psf) + # compute factorization of the active part of J and cache it. + sf_factorize!(Val(sd.factorization), sd.J_factorized, JJ) catch e if e isa SingularException || e isa Pardiso.PardisoException || e isa Pardiso.PardisoPosDefException @error("The system is underdetermined with the given set of equations and final conditions.") @@ -639,6 +629,17 @@ Compute the residual and Jacobian of the stacked time system at the given return RES, sd.J_factorized[] end +function sf_factorize!(v::Val, Rf::Ref{Any}, A::SparseMatrixCSC) + if isnothing(Rf[]) + Rf[] = sf_prepare(v, A) + else + Rf[] = sf_factor!(Rf[], A) + end +end + + +#= KristofferC' implemetation (delete eventually, but keep it for now, just in case) + using Pardiso mutable struct PardisoFactorization @@ -700,6 +701,9 @@ end return lu(JJ) end end + +=# + """ assign_update_step!(x::Array, lambda, dx, sd::StackedTimeSolverData) diff --git a/src/stackedtime/sparse.jl b/src/stackedtime/sparse.jl index f94ffc4..be4dfc7 100644 --- a/src/stackedtime/sparse.jl +++ b/src/stackedtime/sparse.jl @@ -6,10 +6,112 @@ ################################################################################## +# using LinearAlgebra +# using SparseArrays +# using SuiteSparse +# using SuiteSparse.UMFPACK +# using Pardiso + + ### API -sf_factorize(A::SparseMatrixCSC)::Factorization = lu(A)::Factorization -sf_update!(F::Factorization, A::SparseMatrixCSC) = lu!(F, A) -sf_solve!(F::Factorization, b) = ldiv!(F, b) -sf_solve(F::Factorization, b) = (cb = copy(b); ldiv!(F, cb); cb) +# selected sparse linear algebra library is a Symbol +const sf_libs = ( + :default, # use Julia's standard library (UMFPACK) + :lu, # same as :default + # :qr, # use qr-decomposition, rather than lu (still UMFPACK) + :pardiso, # use Pardiso - the one included with MKL +) + +global sf_default = :lu + +# a function to initialize a Factorization instance +# this is also a good place to do the symbolic analysis +sf_prepare(A::SparseMatrixCSC, sparse_lib::Symbol=:default) = sf_prepare(Val(sparse_lib), A) +sf_prepare(::Val{S}, args...) where {S} = throw(ArgumentError("Unknown sparse library $S. Try one of $(sf_libs).")) + +# a function to calculate the numerical factors +sf_factor!(f::Factorization, A::SparseMatrixCSC) = throw(ArgumentError("Unknown factorization type $(typeof(f)).")) + +# a function to solve the linear system +sf_solve!(f::Factorization, x::AbstractArray) = throw(ArgumentError("Unknown factorization type $(typeof(f)).")) + + +### Default (UMFPACK) +sf_prepare(::Val{:default}, A::SparseMatrixCSC) = sf_prepare(Val(sf_default), A) + +@timeit_debug timer "sf_prepare_lu" sf_prepare(::Val{:lu}, A::SparseMatrixCSC) = lu(A) +@timeit_debug timer "sf_factor!_lu" sf_factor!(f::SuiteSparse.UMFPACK.UmfpackLU, A::SparseMatrixCSC) = (lu!(f, A); f) +@timeit_debug timer "sf_solve!_lu" sf_solve!(f::SuiteSparse.UMFPACK.UmfpackLU, x::AbstractArray) = (ldiv!(f, x); x) + +# @timeit_debug timer "sf_prepare_qr" sf_prepare(::Val{:qr}, A::SparseMatrixCSC) = qr(A) +# @timeit_debug timer "sf_factor!_qr" sf_factor!(f::SuiteSparse.SPQR.QRSparse, A::SparseMatrixCSC) = (qr!(f, A); f) +# @timeit_debug timer "sf_solve!_qr" sf_solve!(f::SuiteSparse.SPQR.QRSparse, x::AbstractArray) = (ldiv!(f, x); x) + +### Pardiso (thanks to @KristofferC) + +# See https://github.com/JuliaSparse/Pardiso.jl/blob/master/examples/exampleunsym.jl + +mutable struct PardisoFactorization{Tv<:Real} <: Factorization{Tv} + ps::MKLPardisoSolver + A::SparseMatrixCSC{Tv,Int} +end + +@timeit_debug timer "sf_prepare_par" function sf_prepare(::Val{:pardiso}, A::SparseMatrixCSC) + Tv = eltype(A) + ps = MKLPardisoSolver() + set_matrixtype!(ps, Pardiso.REAL_NONSYM) + pardisoinit(ps) + # See https://www.intel.com/content/www/us/en/develop/documentation/onemkl-developer-reference-c/top/sparse-solver-routines/onemkl-pardiso-parallel-direct-sparse-solver-iface/pardiso-iparm-parameter.html + fix_iparm!(ps, :N) + set_iparm!(ps, 2, 2) # The parallel (OpenMP) version of the nested dissection algorithm. + pf = PardisoFactorization{Tv}(ps, get_matrix(ps, A, :N)) + finalizer(pf) do x + set_phase!(x.ps, Pardiso.RELEASE_ALL) + pardiso(x.ps) + end + _pardiso_full!(pf) + return pf +end + +@timeit_debug timer "_pardso_full" function _pardiso_full!(pf::PardisoFactorization) + # run the analysis phase + ps = pf.ps + set_phase!(ps, Pardiso.ANALYSIS_NUM_FACT) + pardiso(ps, pf.A, Float64[]) + return pf +end + +@timeit_debug timer "_pardso_num" function _pardiso_numeric!(pf::PardisoFactorization) + # run the analysis phase + ps = pf.ps + set_phase!(ps, Pardiso.NUM_FACT) + pardiso(ps, pf.A, Float64[]) + return pf +end + +@timeit_debug timer "sf_factor!_par" function sf_factor!(pf::PardisoFactorization, A::SparseMatrixCSC) + A = get_matrix(pf.ps, A, :N)::typeof(A) + # can we reuse the + _A = pf.A + if A.m == _A.m && A.n == _A.n && A.colptr == _A.colptr && A.rowval == _A.rowval + if A.nzval ≈ _A.nzval + nothing + else + pf.A = A + _pardiso_numeric!(pf) + end + else + pf.A = A + _pardiso_full!(pf) + end + return pf +end + +@timeit_debug timer "sf_solve!_par" function sf_solve!(pf::PardisoFactorization, x::AbstractArray) + ps = pf.ps + set_phase!(ps, Pardiso.SOLVE_ITERATIVE_REFINE) + pardiso(ps, x, pf.A, copy(x)) + return x +end diff --git a/src/stackedtime/stoch_simulate.jl b/src/stackedtime/stoch_simulate.jl index bea08b8..b9a794a 100644 --- a/src/stackedtime/stoch_simulate.jl +++ b/src/stackedtime/stoch_simulate.jl @@ -19,7 +19,6 @@ function stoch_simulate(m::Model, p::Plan, baseline::SimData, shocks; maxiter::Int=m.options.maxiter, fctype=getoption(m, :fctype, fcgiven), #= Newton-Raphson options =# - sparse_solver::Function=\, linesearch::Bool=getoption(m, :linesearch, false), warn_maxiter::Bool=getoption(getoption(m, :warn, Options()), :maxiter, false) ) @@ -122,7 +121,7 @@ function stoch_simulate(m::Model, p::Plan, baseline::SimData, shocks; # solve try - converged = sim_nr!(eₜ, dₜ, maxiter, tol, verbose, sparse_solver, linesearch) + converged = sim_nr!(eₜ, dₜ, maxiter, tol, verbose, linesearch) if warn_maxiter && !converged @warn("Newton-Raphson reached maximum number of iterations for $skey at $(t:sim_end)") end From a47e239b47b0deb0571558207ef917468002d7f0 Mon Sep 17 00:00:00 2001 From: Boyan Bejanov Date: Mon, 10 Apr 2023 16:48:12 -0400 Subject: [PATCH 10/32] renamed :lu case to :umfpack and refactored it to have the same structure as pardiso --- src/stackedtime/solverdata.jl | 6 ++-- src/stackedtime/sparse.jl | 67 ++++++++++++++++++++++++----------- 2 files changed, 49 insertions(+), 24 deletions(-) diff --git a/src/stackedtime/solverdata.jl b/src/stackedtime/solverdata.jl index 62dd6d3..9e40208 100644 --- a/src/stackedtime/solverdata.jl +++ b/src/stackedtime/solverdata.jl @@ -461,7 +461,7 @@ updated in place and returned. # exog_data = reshape(exog_data, sd.NT, sd.NU) @assert(length(res) == size(sd.J, 1), "Length of residual vector doesn't match.") for (ii, tt) in zip(sd.II, sd.TT) - @timeit_debug timer "eval_R!" eval_R!(view(res, ii), point[tt, :], sd.evaldata) + eval_R!(view(res, ii), point[tt, :], sd.evaldata) end return res end @@ -597,12 +597,12 @@ Compute the residual and Jacobian of the stacked time system at the given if haveJ # update only RES for i = 1:length(sd.BI) - @timeit_debug timer "eval_R!" eval_R!(view(RES, sd.II[i]), point[sd.TT[i], :], sd.evaldata) + eval_R!(view(RES, sd.II[i]), point[sd.TT[i], :], sd.evaldata) end else # update both RES and JAC for i = 1:length(sd.BI) - @timeit_debug timer "eval_RJ" R, J = eval_RJ(point[sd.TT[i], :], sd.evaldata) + R, J = eval_RJ(point[sd.TT[i], :], sd.evaldata) RES[sd.II[i]] .= R JAC.nzval[sd.BI[i]] .= J.nzval end diff --git a/src/stackedtime/sparse.jl b/src/stackedtime/sparse.jl index be4dfc7..21b1ee1 100644 --- a/src/stackedtime/sparse.jl +++ b/src/stackedtime/sparse.jl @@ -6,27 +6,19 @@ ################################################################################## -# using LinearAlgebra -# using SparseArrays -# using SuiteSparse -# using SuiteSparse.UMFPACK -# using Pardiso - - ### API # selected sparse linear algebra library is a Symbol const sf_libs = ( :default, # use Julia's standard library (UMFPACK) - :lu, # same as :default - # :qr, # use qr-decomposition, rather than lu (still UMFPACK) + :umfpack, # same as :default :pardiso, # use Pardiso - the one included with MKL ) -global sf_default = :lu +global sf_default = :umfpack # a function to initialize a Factorization instance -# this is also a good place to do the symbolic analysis +# this is also a good place to do the symbolic analysis sf_prepare(A::SparseMatrixCSC, sparse_lib::Symbol=:default) = sf_prepare(Val(sparse_lib), A) sf_prepare(::Val{S}, args...) where {S} = throw(ArgumentError("Unknown sparse library $S. Try one of $(sf_libs).")) @@ -37,20 +29,51 @@ sf_factor!(f::Factorization, A::SparseMatrixCSC) = throw(ArgumentError("Unknown sf_solve!(f::Factorization, x::AbstractArray) = throw(ArgumentError("Unknown factorization type $(typeof(f)).")) +########################################################################### ### Default (UMFPACK) sf_prepare(::Val{:default}, A::SparseMatrixCSC) = sf_prepare(Val(sf_default), A) -@timeit_debug timer "sf_prepare_lu" sf_prepare(::Val{:lu}, A::SparseMatrixCSC) = lu(A) -@timeit_debug timer "sf_factor!_lu" sf_factor!(f::SuiteSparse.UMFPACK.UmfpackLU, A::SparseMatrixCSC) = (lu!(f, A); f) -@timeit_debug timer "sf_solve!_lu" sf_solve!(f::SuiteSparse.UMFPACK.UmfpackLU, x::AbstractArray) = (ldiv!(f, x); x) +function _sf_same_sparse_pattern(A::SparseMatrixCSC, B::SparseMatrixCSC) + return (A.m == B.m) && (A.n == B.n) && (A.colptr == B.colptr) && (A.rowval == B.rowval) +end + +mutable struct LUFactorization{Tv<:Real} <: Factorization{Tv} + F::SuiteSparse.UMFPACK.UmfpackLU{Tv,Int} + A::SparseMatrixCSC{Tv,Int} +end + +@timeit_debug timer "sf_prepare_lu" function sf_prepare(::Val{:umfpack}, A::SparseMatrixCSC) + Tv = eltype(A) + @timeit_debug timer "_lu_full" F = lu(A) + return LUFactorization{Tv}(F, A) +end + +@timeit_debug timer "sf_factor!_lu" function sf_factor!(f::LUFactorization, A::SparseMatrixCSC) + _A = f.A + if _sf_same_sparse_pattern(A, _A) + if A.nzval ≈ _A.nzval + # matrix hasn't changed significantly + nothing + else + # sparse pattern is the same, different numbers + f.A = A + @timeit_debug timer "_lu_num" lu!(f.F, A) + end + else + # totally new matrix, start over + f.A = A + @timeit_debug timer "_lu_full" f.F = lu(A) + end + return f +end -# @timeit_debug timer "sf_prepare_qr" sf_prepare(::Val{:qr}, A::SparseMatrixCSC) = qr(A) -# @timeit_debug timer "sf_factor!_qr" sf_factor!(f::SuiteSparse.SPQR.QRSparse, A::SparseMatrixCSC) = (qr!(f, A); f) -# @timeit_debug timer "sf_solve!_qr" sf_solve!(f::SuiteSparse.SPQR.QRSparse, x::AbstractArray) = (ldiv!(f, x); x) +@timeit_debug timer "sf_solve!_lu" sf_solve!(f::LUFactorization, x::AbstractArray) = (ldiv!(f.F, x); x) +########################################################################### ### Pardiso (thanks to @KristofferC) # See https://github.com/JuliaSparse/Pardiso.jl/blob/master/examples/exampleunsym.jl +# See https://www.intel.com/content/www/us/en/develop/documentation/onemkl-developer-reference-c/top/sparse-solver-routines/onemkl-pardiso-parallel-direct-sparse-solver-iface/pardiso-iparm-parameter.html mutable struct PardisoFactorization{Tv<:Real} <: Factorization{Tv} ps::MKLPardisoSolver @@ -62,9 +85,9 @@ end ps = MKLPardisoSolver() set_matrixtype!(ps, Pardiso.REAL_NONSYM) pardisoinit(ps) - # See https://www.intel.com/content/www/us/en/develop/documentation/onemkl-developer-reference-c/top/sparse-solver-routines/onemkl-pardiso-parallel-direct-sparse-solver-iface/pardiso-iparm-parameter.html fix_iparm!(ps, :N) - set_iparm!(ps, 2, 2) # The parallel (OpenMP) version of the nested dissection algorithm. + # set_iparm!(ps, 1, 1) # Override defaults + # set_iparm!(ps, 2, 2) # Select algorithm pf = PardisoFactorization{Tv}(ps, get_matrix(ps, A, :N)) finalizer(pf) do x set_phase!(x.ps, Pardiso.RELEASE_ALL) @@ -92,16 +115,18 @@ end @timeit_debug timer "sf_factor!_par" function sf_factor!(pf::PardisoFactorization, A::SparseMatrixCSC) A = get_matrix(pf.ps, A, :N)::typeof(A) - # can we reuse the _A = pf.A - if A.m == _A.m && A.n == _A.n && A.colptr == _A.colptr && A.rowval == _A.rowval + if _sf_same_sparse_pattern(A, _A) if A.nzval ≈ _A.nzval + # same matrix, factorization hasn't changed nothing else + # same sparsity pattern, but different numbers pf.A = A _pardiso_numeric!(pf) end else + # totally new matrix, start over pf.A = A _pardiso_full!(pf) end From abdb2ed073b543438a45410d2d9b3bf050575ff6 Mon Sep 17 00:00:00 2001 From: Boyan Bejanov Date: Mon, 10 Apr 2023 19:58:24 -0400 Subject: [PATCH 11/32] unit tests run with all sparse factorization libraries --- src/stackedtime/shockdecomp.jl | 20 +-- src/stackedtime/solverdata.jl | 12 +- src/stackedtime/sparse.jl | 62 ++++++-- test/dynss.jl | 4 +- test/logsimtests.jl | 24 +-- test/misc.jl | 4 +- test/runtests.jl | 24 ++- test/shockdecomp.jl | 17 ++- test/sim_fo.jl | 30 ++-- test/simdatatests.jl | 8 +- test/simtests.jl | 260 +++++++++++++++++---------------- test/sstests.jl | 34 ++--- test/stochsims.jl | 10 +- 13 files changed, 284 insertions(+), 225 deletions(-) diff --git a/src/stackedtime/shockdecomp.jl b/src/stackedtime/shockdecomp.jl index 43450bf..860b3b9 100644 --- a/src/stackedtime/shockdecomp.jl +++ b/src/stackedtime/shockdecomp.jl @@ -19,7 +19,7 @@ dynamic model for the given plan range and final condition. We verify the residual and issue a warning, but do not enforce this. See [`steadystatedata`](@ref). -As part of the algorithm we run a simulation with the given `plan`, +As part of the algorithm we run a simulation with the given `plan`, `exog_data` and `fctype`. See [`simulate`](@ref) for other options. !!! note @@ -66,7 +66,7 @@ function shockdecomp(m::Model, p::Plan, exog_data::SimData; end # prepare the stacked-time - gdata = @timeit_debug timer "StackedTimeSolverData" StackedTimeSolverData(m, p, fctype, variant) + gdata = StackedTimeSolverData(m, p, fctype, variant) # check the residual. why are we doing this? we know it's 0! shocked = copy(exog_data) @@ -78,7 +78,7 @@ function shockdecomp(m::Model, p::Plan, exog_data::SimData; end # We need the Jacobian matrix evaluated at the control solution. - res_control, _ = stackedtime_RJ(control, control, gdata) + res_control, Jac_control = stackedtime_RJ(control, control, gdata) if norm(res_control, Inf) > tol # What to do if it's not a solution? We just give a warning for now, but maybe we should error()?! @warn "Control is not a solution:" norm(res_control, Inf) @@ -128,7 +128,7 @@ function shockdecomp(m::Model, p::Plan, exog_data::SimData; # now do the decomposition begin - # we build a sparse matrix + # we build a sparse matrix SDI = Int[] # row indexes SDJ = Int[] # column indexes SDV = Float64[] # values @@ -160,7 +160,7 @@ function shockdecomp(m::Model, p::Plan, exog_data::SimData; if haskey(id, v.name) # have initial decomposition idv = id[v.name] - if norm(delta[pinit, v.name]-sum(idv[pinit, :],dims=2),Inf)>tol + if norm(delta[pinit, v.name] - sum(idv[pinit, :], dims=2), Inf) > tol error("The given `initdecomp` does not add up for $(v.name).") end rsdv[pinit] .= idv @@ -195,13 +195,7 @@ function shockdecomp(m::Model, p::Plan, exog_data::SimData; end SDMAT = Matrix(gdata.J * sparse(SDI, SDJ, -SDV, size(gdata.J, 2), length(exog_inds))) - factor = gdata.J_factorized[] - if factor isa PardisoFactorization - pardiso_solve!(factor, SDMAT) - else - ldiv!(factor, SDMAT) - end - + sf_solve!(Jac_control, SDMAT) end # now split and decorate the shock-decomposition matrix @@ -216,7 +210,7 @@ function shockdecomp(m::Model, p::Plan, exog_data::SimData; continue end v_inv_endo_inds = inv_endo_inds[v_inds[v_endo_mask]] - v_data = result.sd[v] + v_data = result.sd[v] # assign contributions to endogenous points v_data[v_endo_mask, :] = SDMAT[v_inv_endo_inds, :] # contributions of initial conditions already assigned diff --git a/src/stackedtime/solverdata.jl b/src/stackedtime/solverdata.jl index 9e40208..3f03a8c 100644 --- a/src/stackedtime/solverdata.jl +++ b/src/stackedtime/solverdata.jl @@ -615,16 +615,8 @@ Compute the residual and Jacobian of the stacked time system at the given else JJ = JAC[:, sd.solve_mask] end - try - # compute factorization of the active part of J and cache it. - sf_factorize!(Val(sd.factorization), sd.J_factorized, JJ) - catch e - if e isa SingularException || e isa Pardiso.PardisoException || e isa Pardiso.PardisoPosDefException - @error("The system is underdetermined with the given set of equations and final conditions.") - end - debugging && return RES, deepcopy(JJ) - rethrow() - end + # compute factorization of the active part of J and cache it. + sf_factorize!(Val(sd.factorization), sd.J_factorized, JJ) end return RES, sd.J_factorized[] end diff --git a/src/stackedtime/sparse.jl b/src/stackedtime/sparse.jl index 21b1ee1..65b5cd6 100644 --- a/src/stackedtime/sparse.jl +++ b/src/stackedtime/sparse.jl @@ -44,7 +44,14 @@ end @timeit_debug timer "sf_prepare_lu" function sf_prepare(::Val{:umfpack}, A::SparseMatrixCSC) Tv = eltype(A) - @timeit_debug timer "_lu_full" F = lu(A) + F = try + @timeit_debug timer "_lu_full" lu(A) + catch error + if error isa SingularException + @error("The system is underdetermined with the given set of equations and final conditions.") + end + rethrow() + end return LUFactorization{Tv}(F, A) end @@ -57,12 +64,26 @@ end else # sparse pattern is the same, different numbers f.A = A - @timeit_debug timer "_lu_num" lu!(f.F, A) + try + @timeit_debug timer "_lu_num" lu!(f.F, A) + catch error + if error isa SingularException + @error("The system is underdetermined with the given set of equations and final conditions.") + end + rethrow() + end end else # totally new matrix, start over f.A = A - @timeit_debug timer "_lu_full" f.F = lu(A) + f.F = try + @timeit_debug timer "_lu_full" lu(A) + catch error + if error isa SingularException + @error("The system is underdetermined with the given set of equations and final conditions.") + end + rethrow() + end end return f end @@ -87,6 +108,8 @@ end pardisoinit(ps) fix_iparm!(ps, :N) # set_iparm!(ps, 1, 1) # Override defaults + # set_iparm!(ps, 11, 1) # disable automatic scaling for non-symmetric matrices + # set_iparm!(ps, 13, 1) # disable automatic scaling for non-symmetric matrices # set_iparm!(ps, 2, 2) # Select algorithm pf = PardisoFactorization{Tv}(ps, get_matrix(ps, A, :N)) finalizer(pf) do x @@ -100,16 +123,30 @@ end @timeit_debug timer "_pardso_full" function _pardiso_full!(pf::PardisoFactorization) # run the analysis phase ps = pf.ps - set_phase!(ps, Pardiso.ANALYSIS_NUM_FACT) - pardiso(ps, pf.A, Float64[]) + try + set_phase!(ps, Pardiso.ANALYSIS_NUM_FACT) + pardiso(ps, pf.A, Float64[]) + catch error + if error isa Pardiso.PardisoException || error isa Pardiso.PardisoPosDefException + @error("The system is underdetermined with the given set of equations and final conditions.") + end + rethrow() + end return pf end @timeit_debug timer "_pardso_num" function _pardiso_numeric!(pf::PardisoFactorization) # run the analysis phase ps = pf.ps - set_phase!(ps, Pardiso.NUM_FACT) - pardiso(ps, pf.A, Float64[]) + try + set_phase!(ps, Pardiso.NUM_FACT) + pardiso(ps, pf.A, Float64[]) + catch error + if error isa Pardiso.PardisoException || error isa Pardiso.PardisoPosDefException + @error("The system is underdetermined with the given set of equations and final conditions.") + end + rethrow() + end return pf end @@ -135,8 +172,15 @@ end @timeit_debug timer "sf_solve!_par" function sf_solve!(pf::PardisoFactorization, x::AbstractArray) ps = pf.ps - set_phase!(ps, Pardiso.SOLVE_ITERATIVE_REFINE) - pardiso(ps, x, pf.A, copy(x)) + try + set_phase!(ps, Pardiso.SOLVE_ITERATIVE_REFINE) + pardiso(ps, x, pf.A, copy(x)) + catch error + if error isa Pardiso.PardisoException || error isa Pardiso.PardisoPosDefException + @error("The system is underdetermined with the given set of equations and final conditions.") + end + rethrow() + end return x end diff --git a/test/dynss.jl b/test/dynss.jl index bd9745a..203d33e 100644 --- a/test/dynss.jl +++ b/test/dynss.jl @@ -1,7 +1,7 @@ ################################################################################## # This file is part of StateSpaceEcon.jl # BSD 3-Clause License -# Copyright (c) 2020-2022, Bank of Canada +# Copyright (c) 2020-2023, Bank of Canada # All rights reserved. ################################################################################## @@ -95,7 +95,7 @@ end @testset "dynss+" begin # tests to make sure that changes in parameter values are applied to steady state constraints correctly - local m = PC.model + local m = deepcopy(PC.model) m.sstate.values .= 0.0 @test begin sssolve!(m, presolve=false) diff --git a/test/logsimtests.jl b/test/logsimtests.jl index 06986a0..478ff26 100644 --- a/test/logsimtests.jl +++ b/test/logsimtests.jl @@ -1,7 +1,7 @@ ################################################################################## # This file is part of StateSpaceEcon.jl # BSD 3-Clause License -# Copyright (c) 2020-2022, Bank of Canada +# Copyright (c) 2020-2023, Bank of Canada # All rights reserved. ################################################################################## @@ -12,8 +12,8 @@ @shocks m Ex Ey @parameters m rate = 1.015 @equations m begin - @log X[t] = rate * X[t - 1] + Ex[t] - Y[t] = Y[t - 1] + rate + Ey[t] + @log X[t] = rate * X[t-1] + Ex[t] + Y[t] = Y[t-1] + rate + Ey[t] end @initialize m @@ -22,10 +22,10 @@ clear_sstate!(m) sssolve!(m) - @test m.sstate.X.level ≈ 1.0 + @test m.sstate.X.level ≈ 1.0 @test m.sstate.X.slope ≈ m.rate - @test m.sstate.Y.level + 1 ≈ 0.0 + 1 - @test m.sstate.Y.slope + 1 ≈ m.rate + 1 + @test m.sstate.Y.level + 1 ≈ 0.0 + 1 + @test m.sstate.Y.slope + 1 ≈ m.rate + 1 p = Plan(m, 1U:10U) @@ -39,7 +39,7 @@ @test k[s] == d[s] end @test k.X[1U] ≈ 1.0 - @test k.Y[1U] + 1 ≈ 1.0 + @test k.Y[1U] + 1 ≈ 1.0 @test pct(k.X, -1) ≈ TSeries(p.range, (m.rate - 1) * 100) @test diff(k.Y) ≈ TSeries(p.range, m.rate) @@ -65,8 +65,8 @@ @shocks m Ex Ey @parameters m rate = 1.015 @equations m begin - @log X[t + 1] = rate * X[t] + Ex[t] - Y[t + 1] = Y[t] + rate + Ey[t] + @log X[t+1] = rate * X[t] + Ex[t] + Y[t+1] = Y[t] + rate + Ey[t] end @initialize m @@ -77,8 +77,8 @@ sssolve!(m) @test m.sstate.X.level + 1 ≈ 1.0 + 1 @test m.sstate.X.slope + 1 ≈ m.rate + 1 - @test m.sstate.Y.level + 1 ≈ 0.0 + 1 - @test m.sstate.Y.slope + 1 ≈ m.rate + 1 + @test m.sstate.Y.level + 1 ≈ 0.0 + 1 + @test m.sstate.Y.slope + 1 ≈ m.rate + 1 p = Plan(m, 1U:10U) @@ -119,7 +119,7 @@ # update exogenous data accordingly ed[firstdate(p)] = k[firstdate(p)] - ed[lastdate(p),1:2] .= rand(2) + ed[lastdate(p), 1:2] .= rand(2) @test simulate(m, p, ed, fctype=fcslope) ≈ k @test simulate(m, p, ed, fctype=fcnatural) ≈ k diff --git a/test/misc.jl b/test/misc.jl index 9e4e1fe..0483e03 100644 --- a/test/misc.jl +++ b/test/misc.jl @@ -1,7 +1,7 @@ ################################################################################## # This file is part of StateSpaceEcon.jl # BSD 3-Clause License -# Copyright (c) 2020-2022, Bank of Canada +# Copyright (c) 2020-2023, Bank of Canada # All rights reserved. ################################################################################## @@ -14,7 +14,7 @@ end @testset "fcset" begin m = deepcopy(E7A.model) - + empty!(m.sstate.constraints) @steadystate m lc = 1 @steadystate m linv = 1 diff --git a/test/runtests.jl b/test/runtests.jl index b2d2210..36d914c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -23,6 +23,7 @@ using ModelBaseEcon using Test using Random using Suppressor +import Pardiso @testset "1dsolvers" begin # f(x) = (x-2)*(x-3) = a x^2 + b x + c with vals = [a, x, b, c] @@ -215,11 +216,22 @@ end @test t1 == TSeries(1U, [1, 1, 3, 5, 5, 5, 5, 5]) end -include("simtests.jl") -include("logsimtests.jl") -include("sim_fo.jl") -include("shockdecomp.jl") -include("dynss.jl") include("misc.jl") -include("stochsims.jl") +for sfdef = QuoteNode.(StateSpaceEcon.StackedTimeSolver.sf_libs) + + sfdef.value == :default && continue + + @info "Using $(sfdef)" + + Core.eval(StateSpaceEcon.StackedTimeSolver, :(sf_default = $(sfdef))) + + include("simtests.jl") + include("logsimtests.jl") + include("sim_fo.jl") + include("shockdecomp.jl") + include("dynss.jl") + + include("stochsims.jl") + +end diff --git a/test/shockdecomp.jl b/test/shockdecomp.jl index ebbcf0f..82a6d68 100644 --- a/test/shockdecomp.jl +++ b/test/shockdecomp.jl @@ -1,10 +1,15 @@ +################################################################################## +# This file is part of StateSpaceEcon.jl +# BSD 3-Clause License +# Copyright (c) 2020-2023, Bank of Canada +# All rights reserved. +################################################################################## using LinearAlgebra using JLD2 @testset "shkdcmp" begin m = deepcopy(E7A.model) expected = Workspace(load("data/shkdcmp_E7A.jld2")) - TimeSeriesEcon.clean_old_frequencies!(expected) empty!(m.sstate.constraints) @steadystate m lc = 1 @@ -23,10 +28,10 @@ using JLD2 # shock on dly is half inventory half consumption dlinv_shk = result.sd.dlinv.dlinv_shk[result.sd.dlinv.dlinv_shk.>1e-9] # positive shock values dly_dlinv_shk = result.sd.dly.dlinv_shk[result.sd.dlinv.dlinv_shk.>1e-9] # same rows for dly - @test all(dly_dlinv_shk ./ dlinv_shk .≈ 0.5) + @test norm(dly_dlinv_shk ./ dlinv_shk .- 0.5, Inf) < 1e-8 dlc_shk = result.sd.dlc.dlc_shk[result.sd.dlc.dlc_shk.>1e-9] # positive shock values dly_dlc_shk = result.sd.dly.dlc_shk[result.sd.dlc.dlc_shk.>1e-9] # same rows for dly - @test all(dly_dlc_shk ./ dlc_shk .≈ 0.5) + @test norm(dly_dlc_shk ./ dlc_shk .- 0.5, Inf) < 1e-8 sd_zeros = fill!(similar(result.sd.dlc), 0) @@ -56,10 +61,10 @@ using JLD2 # and the same relation between shocks dlinv_shk = result_dev.sd.dlinv.dlinv_shk[result_dev.sd.dlinv.dlinv_shk.>1e-9] # positive shock values dly_dlinv_shk = result_dev.sd.dly.dlinv_shk[result_dev.sd.dlinv.dlinv_shk.>1e-9] # same rows for dly - @test all(dly_dlinv_shk ./ dlinv_shk .≈ 0.5) + @test norm(dly_dlinv_shk ./ dlinv_shk .- 0.5, Inf) < 1e-8 dlc_shk = result_dev.sd.dlc.dlc_shk[result_dev.sd.dlc.dlc_shk.>1e-9] # positive shock values dly_dlc_shk = result_dev.sd.dly.dlc_shk[result_dev.sd.dlc.dlc_shk.>1e-9] # same rows for dly - @test all(dly_dlc_shk ./ dlc_shk .≈ 0.5) + @test norm(dly_dlc_shk ./ dlc_shk .- 0.5, Inf) < 1e-8 # for consumption, term, init and nonlinear are zero @test norm(result_dev.sd.dlc.init, Inf) < 1e-12 @@ -106,7 +111,7 @@ using JLD2 result3fo = shockdecomp(m, p, exog_data3; solver=:firstorder) @test compare(result3, result3fo; ignoremissing=true, atol=2^10 * eps(), quiet=true) - # additive for linearized model + # additive for linearized model exog_data4 = copy(exog_data3) exog_data4[first(rng), m.shocks] .+= 0.1 result4 = shockdecomp(m, p, exog_data4; control=result3.s, variant=:linearize, fctype=fcnatural) diff --git a/test/sim_fo.jl b/test/sim_fo.jl index 9961a00..7593785 100644 --- a/test/sim_fo.jl +++ b/test/sim_fo.jl @@ -58,7 +58,7 @@ end end @testset "M.fo.unant" begin - run_fo_unant_tests(M.model) + run_fo_unant_tests(deepcopy(M.model)) end module R @@ -79,23 +79,23 @@ end @initialize model end @testset "R.fo.unant" begin - run_fo_unant_tests(R.model) + run_fo_unant_tests(deepcopy(R.model)) end @using_example E2 @testset "E2.fo.unant" begin - run_fo_unant_tests(E2.model) + run_fo_unant_tests(deepcopy(E2.model)) end @using_example E3 @testset "E3.fo.unant" begin - run_fo_unant_tests(E3.model) + run_fo_unant_tests(deepcopy(E3.model)) end @using_example E6 @testset "E6.fo.unant" begin - let m = E6.model + let m = deepcopy(E6.model) # set slopes to 0, otherwise we're not allowed to linearize m.p_dly = 0 m.p_dlp = 0 @@ -131,6 +131,7 @@ end @testset "shkdcmp.fo" begin for m in (M.model, R.model) + m = deepcopy(m) m.rho = 0.6 empty!(m.sstate.constraints) test_shockdecomp_firstorder(m) @@ -139,10 +140,11 @@ end test_shockdecomp_firstorder(m, 1U:20U, fclevel) end for m in (E2.model, E3.model) + m = deepcopy(m) empty!(m.sstate.constraints) test_shockdecomp_firstorder(m, 1U:500U) end - let m = E6.model + let m = deepcopy(E6.model) # set slopes to 0, otherwise we're not allowed to linearize m.p_dly = 0 m.p_dlp = 0 @@ -189,6 +191,7 @@ end @testset "inidcmp.fo" begin for m in (M.model, R.model) + m = deepcopy(m) m.rho = 0.6 empty!(m.sstate.constraints) test_initdecomp_firstorder(m) @@ -197,10 +200,11 @@ end test_initdecomp_firstorder(m) end for m in (E2.model, E3.model) + m = deepcopy(m) empty!(m.sstate.constraints) test_initdecomp_firstorder(m) end - let m = E6.model + let m = deepcopy(E6.model) # set slopes to 0, otherwise we're not allowed to linearize m.p_dly = 0 m.p_dlp = 0 @@ -229,7 +233,7 @@ function test_initdecomp_stackedtime(m; nonlin=!m.linear, rng=2001Q1:2024Q4, fct p = Plan(m, rng) exog = steadystatedata(m, p) exog[begin:first(rng)-1, vars] = rand(m.maxlag, m.nvars) - # zero shocks + # zero shocks ref = shockdecomp(m, p, exog; solver, fctype) for v in vars, s in shks @test norm(ref.sd[v][:, s], Inf) < atol @@ -241,7 +245,7 @@ function test_initdecomp_stackedtime(m; nonlin=!m.linear, rng=2001Q1:2024Q4, fct exog1[begin:first(rng1)-1, vars] = ref.s exog1[rng1, shks] = ref.s[rng1, shks] if fctype === fclevel - exog1[last(rng1)+1:end,vars] = ref.s + exog1[last(rng1)+1:end, vars] = ref.s end res1 = shockdecomp(m, p1, exog1; solver, fctype, initdecomp=ref) res1a = shockdecomp(m, p1, exog1; solver, fctype) @@ -273,7 +277,7 @@ function test_initdecomp_stackedtime(m; nonlin=!m.linear, rng=2001Q1:2024Q4, fct exog1[begin:first(rng1)-1, shks] .= NaN exog1[rng1, shks] = ref.s[rng1, shks] if fctype === fclevel - exog1[last(rng1)+1:end,vars] = ref.s + exog1[last(rng1)+1:end, vars] = ref.s end res1 = shockdecomp(m, p1, exog1; solver, fctype, initdecomp=ref) res1a = shockdecomp(m, p1, exog1; solver, fctype) @@ -304,7 +308,7 @@ function test_initdecomp_stackedtime(m; nonlin=!m.linear, rng=2001Q1:2024Q4, fct exog1[begin:first(rng1)-1, shks] .= NaN # mask shocks in initial conditions exog1[rng1, shks] = ref.s[rng1, shks] if fctype === fclevel - exog1[last(rng1)+1:end,vars] = ref.s + exog1[last(rng1)+1:end, vars] = ref.s end res1 = shockdecomp(m, p1, exog1; solver, fctype, initdecomp=ref) @test rangeof(ref) == rangeof(res1) @@ -314,6 +318,7 @@ end @testset "inidcmp.st" begin for m in (M.model, R.model) + m = deepcopy(m) m.rho = 0.6 empty!(m.sstate.constraints) test_initdecomp_stackedtime(m) @@ -322,10 +327,11 @@ end test_initdecomp_stackedtime(m, fctype=fclevel) end for m in (E2.model, E3.model) + m = deepcopy(m) empty!(m.sstate.constraints) test_initdecomp_stackedtime(m) end - let m = E6.model + let m = deepcopy(E6.model) # set slopes to 0, otherwise we're not allowed to linearize m.p_dly = 0 m.p_dlp = 0 diff --git a/test/simdatatests.jl b/test/simdatatests.jl index 34982b2..33e0136 100644 --- a/test/simdatatests.jl +++ b/test/simdatatests.jl @@ -47,11 +47,11 @@ @test sd[:] == dta[:] @test sd[:a] isa TSeries && sd[:a].values == dta[:, 1] @test sd["b"] isa TSeries && sd["b"].values == dta[:, 2] - # + # sd.a = dta[:, 1] sd.b = dta[:, 2] @test sd[:] == dta[:] - # + # sd[:] = dta[:] sd.a = sd.b @test sd[:, 1] == sd[:, 2] @@ -70,7 +70,7 @@ @test_throws DimensionMismatch sd[2000Q3] = zeros(size(dta, 2) + 1) @test_throws ArgumentError sd[2000Y] @test_throws BoundsError sd[2010Q4] = rand(2) - # + # @test sd[2001Q1:2002Q4] isa SimData @test sd[2001Q1:2002Q4] ≈ sd[5:12, :] sd[2001Q1:2002Q4] = 1:16 @@ -95,7 +95,7 @@ @test similar(sd) isa typeof(sd) @test_throws BoundsError sd[1999Q1:2000Q4] = zeros(8, 2) @test_throws BoundsError sd[2004Q1:2013Q4, [:a, :b]] - # setindex with two + # setindex with two sd[2001Q1, (:b,)] = [3.5] @test sd[5, 2] == 3.5 sd[2000Q1:2001Q4, (:b,)] .= 3.7 diff --git a/test/simtests.jl b/test/simtests.jl index 4103433..4e48cec 100644 --- a/test/simtests.jl +++ b/test/simtests.jl @@ -1,7 +1,7 @@ ################################################################################## # This file is part of StateSpaceEcon.jl # BSD 3-Clause License -# Copyright (c) 2020-2022, Bank of Canada +# Copyright (c) 2020-2023, Bank of Canada # All rights reserved. ################################################################################## @@ -13,21 +13,21 @@ using DelimitedFiles m.α = 0.5 m.β = 1.0 - m.α for T = 6:10 - p = Plan(m, 2:T - 1) + p = Plan(m, 2:T-1) data = zerodata(m, p) data[1, :] = [1 0] # initial condition data[end, :] = [5 0] # final condition sim00 = simulate(m, p, data) - exp00 = hcat(1.0:(5.0 - 1.0) / (T - 1):5.0, zeros(T)) + exp00 = hcat(1.0:(5.0-1.0)/(T-1):5.0, zeros(T)) # @info "T=$T" sim00 exp00 data @test sim00 ≈ exp00 # solution piecewise linear y - one line before and another line after the shock data1 = copy(data) - shk = .1 + shk = 0.1 y2val = ((T - 2) * (1 + 2 * shk) + 5.0) / (T - 1) - data1[2U,:y_shk] = shk # shock at time 2 + data1[2U, :y_shk] = shk # shock at time 2 sim01 = simulate(m, p, data1) - exp01 = vcat([1.0 y2val:(5.0 - y2val) / (T - 2):5.0...], [0 shk zeros(T - 2)...])' + exp01 = vcat([1.0 y2val:(5.0-y2val)/(T-2):5.0...], [0 shk zeros(T - 2)...])' @test sim01 ≈ exp01 # exogenous-endogenous swap # - replicate the solution above backing out the shock that produces it @@ -43,36 +43,36 @@ using DelimitedFiles # test_simulation(E1.model, "data/M1_TestMatrix.csv") end -function test_simulation(m_in, path; atol = 1.0e-9) +function test_simulation(m_in, path; atol=1.0e-9) m = deepcopy(m_in) nvars = ModelBaseEcon.nvariables(m) nshks = ModelBaseEcon.nshocks(m) - # + # data_chk = readdlm(path, ',', Float64) IDX = size(data_chk, 1) - plan = Plan(m, 1 + m.maxlag:IDX - m.maxlead) - # + plan = Plan(m, 1+m.maxlag:IDX-m.maxlead) + # p01 = deepcopy(plan) - autoexogenize!(p01, m, m.maxlag + 1:IDX - m.maxlead) + autoexogenize!(p01, m, m.maxlag+1:IDX-m.maxlead) data01 = zeroarray(m, p01) - data01[1:m.maxlag,:] = data_chk[1:m.maxlag,:] - data01[end - m.maxlead + 1:end,:] = data_chk[end - m.maxlead + 1:end,:] - data01[m.maxlag + 1:end - m.maxlead,1:nvars] = data_chk[m.maxlag + 1:end - m.maxlead,1:nvars] + data01[1:m.maxlag, :] = data_chk[1:m.maxlag, :] + data01[end-m.maxlead+1:end, :] = data_chk[end-m.maxlead+1:end, :] + data01[m.maxlag+1:end-m.maxlead, 1:nvars] = data_chk[m.maxlag+1:end-m.maxlead, 1:nvars] res01 = simulate(m, p01, data01) - @test isapprox(res01, data_chk; atol = atol) - + @test isapprox(res01, data_chk; atol=atol) + p02 = deepcopy(plan) data02 = zeroarray(m, p02) - data02[1:m.maxlag,:] = data_chk[1:m.maxlag,:] - data02[end - m.maxlead + 1:end,:] = data_chk[end - m.maxlead + 1:end,:] - data02[m.maxlag + 1:end - m.maxlead,nvars .+ (1:nshks)] = data_chk[m.maxlag + 1:end - m.maxlead,nvars .+ (1:nshks)] + data02[1:m.maxlag, :] = data_chk[1:m.maxlag, :] + data02[end-m.maxlead+1:end, :] = data_chk[end-m.maxlead+1:end, :] + data02[m.maxlag+1:end-m.maxlead, nvars.+(1:nshks)] = data_chk[m.maxlag+1:end-m.maxlead, nvars.+(1:nshks)] res02 = simulate(m, p02, data02) - @test isapprox(res02, data_chk; atol = atol) + @test isapprox(res02, data_chk; atol=atol) # linesearch m.options.linesearch = true res01_line = simulate(m, p01, data01) - @test isapprox(res01_line, data_chk; atol = atol) + @test isapprox(res01_line, data_chk; atol=atol) m.options.linesearch = false # deviation @@ -80,19 +80,19 @@ function test_simulation(m_in, path; atol = 1.0e-9) data_chk_dev = deepcopy(data_chk) sssolve!(m) for (i, var) in enumerate(m.allvars) - data01_dev[:,i] .-= m.sstate[var].level - data_chk_dev[:,i] .-= m.sstate[var].level + data01_dev[:, i] .-= m.sstate[var].level + data_chk_dev[:, i] .-= m.sstate[var].level if m.sstate[var].slope != 0 period = 0 for j in p01.range - data01_dev[j,i] -= m.sstate[var].slope*period - data_chk_dev[j,i] -= m.sstate[var].slope*period + data01_dev[j, i] -= m.sstate[var].slope * period + data_chk_dev[j, i] -= m.sstate[var].slope * period period += 1 end end end res01_dev = simulate(m, p01, data01_dev; deviation=true) - @test isapprox(res01_dev, data_chk_dev; atol = atol) + @test isapprox(res01_dev, data_chk_dev; atol=atol) end @testset "E1.sim" begin @@ -115,7 +115,7 @@ end # linearization tests @testset "linearize" begin - m3 = deepcopy(E3.model) + m3 = deepcopy(E3.model) clear_sstate!(m3) @test_throws ModelBaseEcon.LinearizationError linearize!(m3) @@ -144,10 +144,10 @@ end m7 = deepcopy(E7.model) if !issssolved(m7) empty!(m7.sstate.constraints) - @steadystate m7 linv = lc - 7; - @steadystate m7 lc = 14; + @steadystate m7 linv = lc - 7 + @steadystate m7 lc = 14 clear_sstate!(m7) - sssolve!(m7); + sssolve!(m7) end # the aux variables are present @test_throws ModelBaseEcon.LinearizationError linearize!(m7) @@ -167,13 +167,13 @@ end @testset "E1.unant" begin m = deepcopy(E1.model) m.α = m.β = 0.5 - p = Plan(m, 1:3); - data = zeroarray(m, p); - data[1,1] = 1.0; - data[end,1] = 5.0; - data[3,2] = 0.1; + p = Plan(m, 1:3) + data = zeroarray(m, p) + data[1, 1] = 1.0 + data[end, 1] = 5.0 + data[3, 2] = 0.1 res_u = simulate(m, p, data; anticipate=false) - true_u = Float64[1 0; 2 0; 47 / 15 0.1; 61 / 15 0; 5 0] + true_u = Float64[1 0; 2 0; 47/15 0.1; 61/15 0; 5 0] @test res_u ≈ true_u atol = 1e-12 end @@ -184,23 +184,24 @@ end var_inds = 1:nvars shk_inds = nvars .+ (1:nshks) # Steady state - m.sstate.values .= 0; - m.sstate.mask .= true; + m.sstate.values .= 0 + m.sstate.mask .= true let atrue = readdlm("./data/M2_Ant.csv", ',', Float64), utrue = readdlm("./data/M2_Unant.csv", ',', Float64) # Set simulation ranges and plan IDX = size(atrue, 1) init = 1:m.maxlag - sim = 1 + m.maxlag:IDX - m.maxlead - term = IDX .+ (1 - m.maxlead:0) + sim = 1+m.maxlag:IDX-m.maxlead + term = IDX .+ (1-m.maxlead:0) p = Plan(m, sim) # Make up the exogenous data for the experiment let adata = zeroarray(m, p), udata = zeroarray(m, p) - adata[sim,shk_inds] .= atrue[sim,shk_inds] - udata[sim,shk_inds] .= utrue[sim,shk_inds] - adata[term,:] .= atrue[term,:] - udata[term,:] .= utrue[term,:] + + adata[sim, shk_inds] .= atrue[sim, shk_inds] + udata[sim, shk_inds] .= utrue[sim, shk_inds] + adata[term, :] .= atrue[term, :] + udata[term, :] .= utrue[term, :] # Run the simulations and test ares = simulate(m, p, adata) @test ares ≈ atrue atol = m.options.tol @@ -209,11 +210,12 @@ end end let adata = zeroarray(m, p), udata = zeroarray(m, p) + autoexogenize!(p, m, sim) - adata[sim,var_inds] .= atrue[sim,var_inds] - udata[sim,var_inds] .= utrue[sim,var_inds] - adata[term,:] .= atrue[term,:] - udata[term,:] .= utrue[term,:] + adata[sim, var_inds] .= atrue[sim, var_inds] + udata[sim, var_inds] .= utrue[sim, var_inds] + adata[term, :] .= atrue[term, :] + udata[term, :] .= utrue[term, :] # Run the simulations and test ares = simulate(m, p, adata) @test ares ≈ atrue atol = m.options.tol @@ -227,17 +229,17 @@ end # Set simulation ranges and plan IDX = size(atrue, 1) init = 1:m.maxlag - sim = 1 + m.maxlag:IDX - m.maxlead - term = IDX .+ (1 - m.maxlead:0) + sim = 1+m.maxlag:IDX-m.maxlead + term = IDX .+ (1-m.maxlead:0) p = Plan(m, sim) # Make up the exogenous data for the experiment - adata = zeroarray(m, p); - udata = zeroarray(m, p); + adata = zeroarray(m, p) + udata = zeroarray(m, p) ygap_ind = indexin([:ygap], m.variables)[1] adata[sim[1:s], ygap_ind] .= -0.10 udata[sim[1:s], ygap_ind] .= -0.10 - adata[term,:] .= atrue[term,:] - udata[term,:] .= utrue[term,:] + adata[term, :] .= atrue[term, :] + udata[term, :] .= utrue[term, :] exogenize!(p, :ygap, sim[1:s]) endogenize!(p, :ygap_shk, sim[1:s]) # Run the simulations and test @@ -255,23 +257,24 @@ end var_inds = 1:nvars shk_inds = nvars .+ (1:nshks) # Steady state - m.sstate.values .= 0; - m.sstate.mask .= true; + m.sstate.values .= 0 + m.sstate.mask .= true let atrue = readdlm("./data/M3_Ant.csv", ',', Float64), utrue = readdlm("./data/M3_Unant.csv", ',', Float64) # Set simulation ranges and plan IDX = size(atrue, 1) init = 1:m.maxlag - sim = 1 + m.maxlag:IDX - m.maxlead - term = IDX .+ (1 - m.maxlead:0) + sim = 1+m.maxlag:IDX-m.maxlead + term = IDX .+ (1-m.maxlead:0) p = Plan(m, sim) # Make up the exogenous data for the experiment let adata = zeroarray(m, p), udata = zeroarray(m, p) - adata[sim,shk_inds] .= atrue[sim,shk_inds] - udata[sim,shk_inds] .= utrue[sim,shk_inds] - adata[term,:] .= atrue[term,:] - udata[term,:] .= utrue[term,:] + + adata[sim, shk_inds] .= atrue[sim, shk_inds] + udata[sim, shk_inds] .= utrue[sim, shk_inds] + adata[term, :] .= atrue[term, :] + udata[term, :] .= utrue[term, :] # Run the simulations and test ares = simulate(m, p, adata) @test ares ≈ atrue atol = m.options.tol @@ -280,11 +283,12 @@ end end let adata = zeroarray(m, p), udata = zeroarray(m, p) + autoexogenize!(p, m, sim) - adata[sim,var_inds] .= atrue[sim,var_inds] - udata[sim,var_inds] .= utrue[sim,var_inds] - adata[term,:] .= atrue[term,:] - udata[term,:] .= utrue[term,:] + adata[sim, var_inds] .= atrue[sim, var_inds] + udata[sim, var_inds] .= utrue[sim, var_inds] + adata[term, :] .= atrue[term, :] + udata[term, :] .= utrue[term, :] # Run the simulations and test ares = simulate(m, p, adata) @test ares ≈ atrue atol = m.options.tol @@ -292,24 +296,24 @@ end @test ures ≈ utrue atol = m.options.tol end # tests with different final conditions and expectation_horizon values. - for (fc, eh) = Iterators.product((fcgiven, fclevel, fcslope, fcnatural), ( nothing, 30, 50, 100)) + for (fc, eh) = Iterators.product((fcgiven, fclevel, fcslope, fcnatural), (nothing, 30, 50, 100)) # Make up the exogenous data for the experiment p = Plan(m, sim) udata = zeroarray(m, p) - udata[sim,shk_inds] .= utrue[sim,shk_inds] - udata[term,:] .= utrue[term,:] + udata[sim, shk_inds] .= utrue[sim, shk_inds] + udata[term, :] .= utrue[term, :] # Run the simulations and test ures = simulate(m, p, udata; anticipate=false, expectation_horizon=eh, fctype=fc) - @test ures ≈ utrue atol = m.options.tol rtol = maximum(abs, utrue[term,:]) + @test ures ≈ utrue atol = m.options.tol rtol = maximum(abs, utrue[term, :]) # Repeat with autoexogenize p = Plan(m, sim) udata = zeroarray(m, p) autoexogenize!(p, m, sim) - udata[sim,var_inds] .= utrue[sim,var_inds] - udata[term,:] .= utrue[term,:] + udata[sim, var_inds] .= utrue[sim, var_inds] + udata[term, :] .= utrue[term, :] # Run the simulations and test ures = simulate(m, p, udata; anticipate=false, expectation_horizon=eh, fctype=fc) - @test ures ≈ utrue atol = m.options.tol rtol = maximum(abs, utrue[term,:]) + @test ures ≈ utrue atol = m.options.tol rtol = maximum(abs, utrue[term, :]) end end for s in (1, 4) @@ -318,17 +322,17 @@ end # Set simulation ranges and plan IDX = size(atrue, 1) init = 1:m.maxlag - sim = 1 + m.maxlag:IDX - m.maxlead - term = IDX .+ (1 - m.maxlead:0) + sim = 1+m.maxlag:IDX-m.maxlead + term = IDX .+ (1-m.maxlead:0) p = Plan(m, sim) # Make up the exogenous data for the experiment - adata = zeroarray(m, p); - udata = zeroarray(m, p); + adata = zeroarray(m, p) + udata = zeroarray(m, p) ygap_ind = indexin([:ygap], m.variables)[1] adata[sim[1:s], ygap_ind] .= -0.10 udata[sim[1:s], ygap_ind] .= -0.10 - adata[term,:] .= atrue[term,:] - udata[term,:] .= utrue[term,:] + adata[term, :] .= atrue[term, :] + udata[term, :] .= utrue[term, :] exogenize!(p, :ygap, sim[1:s]) endogenize!(p, :ygap_shk, sim[1:s]) # Run the simulations and test @@ -342,7 +346,9 @@ end @testset "FCTypes" begin let m = Model() @variables m x - @equations m begin x[t] = x[t + 3] end + @equations m begin + x[t] = x[t+3] + end @initialize m p = Plan(m, 1:1) ed = zerodata(m, p) @@ -352,7 +358,7 @@ end p = Plan(m, 3:17) ed = zerodata(m, p) ed .= rand(Float64, size(ed)) - # @test_throws "The system is underdetermined.*" simulate(m, p, ed, fctype=fcnatural) + @test_logs (:error, r"The system is underdetermined.*") (@test_throws Union{SingularException,Pardiso.PardisoPosDefException,Pardiso.PardisoException} simulate(m, p, ed, fctype=fcnatural)) sd = StateSpaceEcon.StackedTimeSolver.StackedTimeSolverData(m, p, fcgiven) @test_throws ErrorException StateSpaceEcon.StackedTimeSolver.update_plan!(sd, m, Plan(m, 3:8)) @@ -364,7 +370,7 @@ end let m = deepcopy(E3.model) p = Plan(m, 3:177) ed = zerodata(m, p) - ed[3U, m.shocks] .= 0.2 * rand(1,3) + ed[3U, m.shocks] .= 0.2 * rand(1, 3) ed[3U:end, m.variables] .= rand(178, 3) empty!(m.sstate.constraints) clear_sstate!(m) @@ -377,9 +383,9 @@ end s1 = simulate(m, p, ed, fctype=fclevel) s2 = simulate(m, p, ed, fctype=fcslope) s3 = simulate(m, p, ed, fctype=fcnatural) - @test maximum(abs, s1 - s2; dims = :) < 1e-8 - @test maximum(abs, s1 - s3; dims = :) < 1e-8 - @test maximum(abs, s3 - s2; dims = :) < 1e-8 + @test maximum(abs, s1 - s2; dims=:) < 1e-8 + @test maximum(abs, s1 - s3; dims=:) < 1e-8 + @test maximum(abs, s3 - s2; dims=:) < 1e-8 sd = StateSpaceEcon.StackedTimeSolver.StackedTimeSolverData(m, p, fcgiven) x = rand(Float64, size(ed)) @@ -389,7 +395,7 @@ end end let m = deepcopy(E6.model) empty!(m.sstate.constraints) - @steadystate m lp = 1.0 + @steadystate m lp = 1.0 @steadystate m lp + lyn = 1.5 @steadystate m lp + ly = 1.2 clear_sstate!(m) @@ -400,7 +406,7 @@ end ed.dly_shk[3] = rand() s2 = simulate(m, p, ed, fctype=fcslope) s3 = simulate(m, p, ed, fctype=fcnatural) - @test maximum(abs, s3 - s2; dims = :) < 1e-12 + @test maximum(abs, s3 - s2; dims=:) < 1e-12 sd = StateSpaceEcon.StackedTimeSolver.StackedTimeSolverData(m, p, fcgiven) x = rand(Float64, size(ed)) @@ -410,7 +416,7 @@ end end @testset "new.unant" begin - m = E1.model + m = deepcopy(E1.model) m.α = m.β = 0.5 test_rng = 20Q1:22Q1 @@ -440,7 +446,7 @@ end # now attempt to recover shocks from solutions p_a2 = autoexogenize!(copy(p_a), m, test_rng) p_u2 = copy(p_a2) - # recover + # recover chk_a = let data = copy(res_a) data.y_shk[test_rng] .= 100rand(length(test_rng)) simulate(m, p_a2, data) @@ -474,95 +480,95 @@ end clear_sstate!(m) nvars = ModelBaseEcon.nvariables(m) nshks = ModelBaseEcon.nshocks(m) - # + # data_chk = readdlm("data/M3_TestData.csv", ',', Float64) IDX = size(data_chk, 1) - plan = Plan(m, 1 + m.maxlag:IDX - m.maxlead) - # + plan = Plan(m, 1+m.maxlag:IDX-m.maxlead) + # p01 = deepcopy(plan) - autoexogenize!(p01, m, m.maxlag + 1:IDX - m.maxlead) + autoexogenize!(p01, m, m.maxlag+1:IDX-m.maxlead) data01 = zeroarray(m, p01) - data01[1:m.maxlag,:] = data_chk[1:m.maxlag,:] - data01[end - m.maxlead + 1:end,:] = data_chk[end - m.maxlead + 1:end,:] - data01[m.maxlag + 1:end - m.maxlead,1:nvars] = data_chk[m.maxlag + 1:end - m.maxlead,1:nvars] + data01[1:m.maxlag, :] = data_chk[1:m.maxlag, :] + data01[end-m.maxlead+1:end, :] = data_chk[end-m.maxlead+1:end, :] + data01[m.maxlag+1:end-m.maxlead, 1:nvars] = data_chk[m.maxlag+1:end-m.maxlead, 1:nvars] initial_guess = similar(data01) initial_guess .= 1 - initial_guess[end-2:end,:] .= 0 #at the very end + initial_guess[end-2:end, :] .= 0 #at the very end res01 = simulate(m, p01, data01) res01_line = nothing m.options.linesearch = true out = @capture_err begin res01_line = simulate(m, p01, data01; initial_guess=initial_guess, verbose=true) end - + # line search should give the same results @test occursin("Linesearch success", out) - @test isapprox(res01, data_chk; atol = 1.0e-9) - @test isapprox(res01_line, data_chk; atol = 1.0e-9) - @test isapprox(res01, res01_line; atol = 1.0e-9) + @test isapprox(res01, data_chk; atol=1.0e-9) + @test isapprox(res01_line, data_chk; atol=1.0e-9) + @test isapprox(res01, res01_line; atol=1.0e-9) end let m = deepcopy(E1.model) clear_sstate!(m) nvars = ModelBaseEcon.nvariables(m) nshks = ModelBaseEcon.nshocks(m) - # + # data_chk = readdlm("data/M1_TestData.csv", ',', Float64) IDX = size(data_chk, 1) - plan = Plan(m, 1 + m.maxlag:IDX - m.maxlead) - # + plan = Plan(m, 1+m.maxlag:IDX-m.maxlead) + # p01 = deepcopy(plan) - autoexogenize!(p01, m, m.maxlag + 1:IDX - m.maxlead) + autoexogenize!(p01, m, m.maxlag+1:IDX-m.maxlead) data01 = zeroarray(m, p01) - data01[1:m.maxlag,:] = data_chk[1:m.maxlag,:] - data01[end - m.maxlead + 1:end,:] = data_chk[end - m.maxlead + 1:end,:] - data01[m.maxlag + 1:end - m.maxlead,1:nvars] = data_chk[m.maxlag + 1:end - m.maxlead,1:nvars] + data01[1:m.maxlag, :] = data_chk[1:m.maxlag, :] + data01[end-m.maxlead+1:end, :] = data_chk[end-m.maxlead+1:end, :] + data01[m.maxlag+1:end-m.maxlead, 1:nvars] = data_chk[m.maxlag+1:end-m.maxlead, 1:nvars] # deviation gives the same results res01 = simulate(m, p01, data01) - res01_dev = simulate(m, p01, data01; deviation=true) - @test isapprox(res01_dev, data_chk; atol = 1.0e-9) - @test isapprox(res01_dev, res01; atol = 1.0e-9) + res01_dev = simulate(m, p01, data01; deviation=true) + @test isapprox(res01_dev, data_chk; atol=1.0e-9) + @test isapprox(res01_dev, res01; atol=1.0e-9) end let m = deepcopy(E3.model) clear_sstate!(m) nvars = ModelBaseEcon.nvariables(m) nshks = ModelBaseEcon.nshocks(m) - # + # data_chk = readdlm("data/M3_TestData.csv", ',', Float64) IDX = size(data_chk, 1) - plan = Plan(m, 1 + m.maxlag:IDX - m.maxlead) - # + plan = Plan(m, 1+m.maxlag:IDX-m.maxlead) + # p01 = deepcopy(plan) - autoexogenize!(p01, m, m.maxlag + 1:IDX - m.maxlead) + autoexogenize!(p01, m, m.maxlag+1:IDX-m.maxlead) data01 = zeroarray(m, p01) - data01[1:m.maxlag,:] = data_chk[1:m.maxlag,:] - data01[end - m.maxlead + 1:end,:] = data_chk[end - m.maxlead + 1:end,:] - data01[m.maxlag + 1:end - m.maxlead,1:nvars] = data_chk[m.maxlag + 1:end - m.maxlead,1:nvars] + data01[1:m.maxlag, :] = data_chk[1:m.maxlag, :] + data01[end-m.maxlead+1:end, :] = data_chk[end-m.maxlead+1:end, :] + data01[m.maxlag+1:end-m.maxlead, 1:nvars] = data_chk[m.maxlag+1:end-m.maxlead, 1:nvars] # deviation do not give the same result res01 = simulate(m, p01, data01) res01_dev = simulate(m, p01, data01; deviation=true) - @test isapprox(res01, data_chk; atol = 1.0e-9) - @test !isapprox(res01_dev, data_chk; atol = 1.0e-9) - @test !isapprox(res01_dev, res01; atol = 1.0e-9) + @test isapprox(res01, data_chk; atol=1.0e-9) + @test !isapprox(res01_dev, data_chk; atol=1.0e-9) + @test !isapprox(res01_dev, res01; atol=1.0e-9) # subtract steady-state from data and results data02 = deepcopy(data01) data_chk02 = deepcopy(data_chk) sssolve!(m) for (i, var) in enumerate(m.allvars) - data02[:,i] .-= m.sstate[var].level - data_chk02[:,i] .-= m.sstate[var].level + data02[:, i] .-= m.sstate[var].level + data_chk02[:, i] .-= m.sstate[var].level end res02 = simulate(m, p01, data02) res02_dev = simulate(m, p01, data02; deviation=true) - # results equal new data_chk02 - @test isapprox(res01, data_chk02; atol = 1.0e-9) - @test isapprox(res02_dev, data_chk02; atol = 1.0e-9) - @test isapprox(res02_dev, res01; atol = 1.0e-9) + # results equal new data_chk02 + @test isapprox(res01, data_chk02; atol=1.0e-9) + @test isapprox(res02_dev, data_chk02; atol=1.0e-9) + @test isapprox(res02_dev, res01; atol=1.0e-9) end end diff --git a/test/sstests.jl b/test/sstests.jl index 4a598c4..9bf104a 100644 --- a/test/sstests.jl +++ b/test/sstests.jl @@ -1,7 +1,7 @@ ################################################################################## # This file is part of StateSpaceEcon.jl # BSD 3-Clause License -# Copyright (c) 2020-2022, Bank of Canada +# Copyright (c) 2020-2023, Bank of Canada # All rights reserved. ################################################################################## @@ -14,7 +14,7 @@ empty!(E1.model.sstate.constraints) sssolve!(m) @test check_sstate(m) == 0 @test m.sstate.mask == [false, false, true, true] - # + # @steadystate m y = 1.2 m.α = 0.5 m.β = 1.0 - m.α @@ -23,7 +23,7 @@ empty!(E1.model.sstate.constraints) @test check_sstate(m) == 0 @test m.sstate.mask == [true, false, true, true] @test m.sstate.values[1] == 1.2 - # + # m.α = 0.4 m.β = 1.0 - m.α clear_sstate!(m) @@ -32,17 +32,17 @@ empty!(E1.model.sstate.constraints) @test check_sstate(m) == 0 @test m.sstate.mask == trues(4) @test m.sstate.values == [1.2, 0.0, 0.0, 0.0] - # + # empty!(m.sstate.constraints) m.α = 0.3 m.β = 1.0 - m.α clear_sstate!(m) @suppress begin - sssolve!(m; verbose=true, nropts = Options(linesearch = true)) + sssolve!(m; verbose=true, nropts=Options(linesearch=true)) end @test check_sstate(m) == 0 @test m.sstate.values[2] == 0.0 - # + # empty!(m.sstate.constraints) m.α = 0.3 m.β = 1.0 - m.α @@ -60,7 +60,7 @@ empty!(E1.model.sstate.constraints) sssolve!(m; method=:auto) @test check_sstate(m) == 0 @test m.sstate.mask == [true, true, true, true] # different from default - # + # @steadystate m y = 1.2 m.α = 0.5 m.β = 1.0 - m.α @@ -69,7 +69,7 @@ empty!(E1.model.sstate.constraints) @test check_sstate(m) == 0 @test m.sstate.mask == [true, true, true, true] # different from default @test m.sstate.values[1] == 1.2 - # + # m.α = 0.4 m.β = 1.0 - m.α clear_sstate!(m) @@ -78,7 +78,7 @@ empty!(E1.model.sstate.constraints) @test check_sstate(m) == 0 @test m.sstate.mask == trues(4) @test m.sstate.values == [1.2, 0.0, 0.0, 0.0] - # + # empty!(m.sstate.constraints) m.α = 0.3 m.β = 1.0 - m.α @@ -100,7 +100,7 @@ empty!(E1.model.sstate.constraints) sssolve!(m; method=:lm) @test check_sstate(m) == 0 @test m.sstate.mask == [true, true, true, true] # different from default - # + # @steadystate m y = 1.2 m.α = 0.5 m.β = 1.0 - m.α @@ -109,7 +109,7 @@ empty!(E1.model.sstate.constraints) @test check_sstate(m) == 0 @test m.sstate.mask == [true, true, true, true] # different from default @test m.sstate.values[1] == 1.2 - # + # m.α = 0.4 m.β = 1.0 - m.α clear_sstate!(m) @@ -118,7 +118,7 @@ empty!(E1.model.sstate.constraints) @test check_sstate(m) == 0 @test m.sstate.mask == trues(4) @test m.sstate.values == [1.2, 0.0, 0.0, 0.0] - # + # empty!(m.sstate.constraints) m.α = 0.3 m.β = 1.0 - m.α @@ -195,7 +195,7 @@ end @steadystate m lc = 14 clear_sstate!(m) sssolve!(m) - @test check_sstate(m, tol = 10m.tol) == 0 + @test check_sstate(m, tol=10m.tol) == 0 @test all(m.sstate.mask) @test m.sstate.values ≈ [0.004, 0.0, 0.004, 0.0, 0.004, 0.0, 14.0, 0.004, 7.0, 0.004, 9.267287357063445, 0.004, 14.000911466453774, 0.004, 0, 0, 0, 0, 14.000911466453774, 0.004, 9.267287357063445, 0.004] @@ -285,7 +285,7 @@ end @test check_sstate(m; verbose=true) == 6 end @test occursin("System may be inconsistent", out) - + end end @@ -296,12 +296,12 @@ end @variables m c @equations m begin a[t] = (a[t-1] + 10)^2 - b[t] = b[t-1]/0 + b[t] = b[t-1] / 0 b[t] = (a[t-1] + 10)^2 c[t] = a[t] + b[t] end @initialize m - + # tests @test sum(issteady.(m.allvars)) == 2 clear_sstate!(m) @@ -327,7 +327,7 @@ end @test check_sstate(m) == 0 @test m.sstate.mask == [false, false, true, true] end - + empty!(E1.model.sstate.constraints) let m = deepcopy(E1.model) m.α = 0.5 diff --git a/test/stochsims.jl b/test/stochsims.jl index 11e7642..743549a 100644 --- a/test/stochsims.jl +++ b/test/stochsims.jl @@ -41,11 +41,11 @@ db_1 = copy(db_ss) shk_1 = MVTSeries(srng, m.shocks, rand) db_1 .= db_1 .+ shk_1 - - @test simulate(m, plan, db_1; anticipate=false) ≈ stoch_simulate(m, plan, db_ss, shk_1)[1] - @test simulate(m, plan, db_1; anticipate=false) ≈ stoch_simulate(m, plan, rawdata(db_ss), shk_1)[1] - @test simulate(m, plan, db_1; anticipate=false) ≈ stoch_simulate(m, plan, Workspace(; pairs(db_ss)...), shk_1)[1] - + + @test simulate(m, plan, db_1; anticipate=false) ≈ stoch_simulate(m, plan, db_ss, shk_1)[1] + @test simulate(m, plan, db_1; anticipate=false) ≈ stoch_simulate(m, plan, rawdata(db_ss), shk_1)[1] + @test simulate(m, plan, db_1; anticipate=false) ≈ stoch_simulate(m, plan, Workspace(; pairs(db_ss)...), shk_1)[1] + db_1 = simulate(m, plan, db_1; anticipate=true) res = stoch_simulate(m, plan, db_1, shk_sample; check=true) for (r, s) in zip(res, shk_sample) From c9a260f72f443223b201d97e4852982a34bd320e Mon Sep 17 00:00:00 2001 From: Boyan Bejanov Date: Mon, 10 Apr 2023 20:33:50 -0400 Subject: [PATCH 12/32] turn off pardiso tests on macOS --- test/runtests.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/test/runtests.jl b/test/runtests.jl index 36d914c..80a85fb 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -222,6 +222,9 @@ for sfdef = QuoteNode.(StateSpaceEcon.StackedTimeSolver.sf_libs) sfdef.value == :default && continue + # Pardiso in macos is giving "Not enough memory". Disable for now + sfdef.value == :pardiso && Sys.isapple() && continue + @info "Using $(sfdef)" Core.eval(StateSpaceEcon.StackedTimeSolver, :(sf_default = $(sfdef))) From 3ae961f646827855dd2bc56403f3002bd7a10a85 Mon Sep 17 00:00:00 2001 From: Boyan Bejanov Date: Mon, 10 Apr 2023 20:34:44 -0400 Subject: [PATCH 13/32] relax some test tolerances --- test/shockdecomp.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/test/shockdecomp.jl b/test/shockdecomp.jl index 82a6d68..e48eb60 100644 --- a/test/shockdecomp.jl +++ b/test/shockdecomp.jl @@ -28,10 +28,10 @@ using JLD2 # shock on dly is half inventory half consumption dlinv_shk = result.sd.dlinv.dlinv_shk[result.sd.dlinv.dlinv_shk.>1e-9] # positive shock values dly_dlinv_shk = result.sd.dly.dlinv_shk[result.sd.dlinv.dlinv_shk.>1e-9] # same rows for dly - @test norm(dly_dlinv_shk ./ dlinv_shk .- 0.5, Inf) < 1e-8 + @test norm(dly_dlinv_shk ./ dlinv_shk .- 0.5, Inf) < 1e-7 dlc_shk = result.sd.dlc.dlc_shk[result.sd.dlc.dlc_shk.>1e-9] # positive shock values dly_dlc_shk = result.sd.dly.dlc_shk[result.sd.dlc.dlc_shk.>1e-9] # same rows for dly - @test norm(dly_dlc_shk ./ dlc_shk .- 0.5, Inf) < 1e-8 + @test norm(dly_dlc_shk ./ dlc_shk .- 0.5, Inf) < 1e-7 sd_zeros = fill!(similar(result.sd.dlc), 0) @@ -61,10 +61,10 @@ using JLD2 # and the same relation between shocks dlinv_shk = result_dev.sd.dlinv.dlinv_shk[result_dev.sd.dlinv.dlinv_shk.>1e-9] # positive shock values dly_dlinv_shk = result_dev.sd.dly.dlinv_shk[result_dev.sd.dlinv.dlinv_shk.>1e-9] # same rows for dly - @test norm(dly_dlinv_shk ./ dlinv_shk .- 0.5, Inf) < 1e-8 + @test norm(dly_dlinv_shk ./ dlinv_shk .- 0.5, Inf) < 1e-7 dlc_shk = result_dev.sd.dlc.dlc_shk[result_dev.sd.dlc.dlc_shk.>1e-9] # positive shock values dly_dlc_shk = result_dev.sd.dly.dlc_shk[result_dev.sd.dlc.dlc_shk.>1e-9] # same rows for dly - @test norm(dly_dlc_shk ./ dlc_shk .- 0.5, Inf) < 1e-8 + @test norm(dly_dlc_shk ./ dlc_shk .- 0.5, Inf) < 1e-7 # for consumption, term, init and nonlinear are zero @test norm(result_dev.sd.dlc.init, Inf) < 1e-12 From d81d9f21e29c73cfdffd3e8fe72fb83d0399d6db Mon Sep 17 00:00:00 2001 From: Boyan Bejanov Date: Mon, 10 Apr 2023 20:36:47 -0400 Subject: [PATCH 14/32] 1-thread pardiso in tests --- test/runtests.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/test/runtests.jl b/test/runtests.jl index 80a85fb..4965b20 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -218,6 +218,9 @@ end include("misc.jl") +# make sure Pardiso runs deterministic (single thread) for the tests +Pardiso.set_nprocs_mkl!(1) + for sfdef = QuoteNode.(StateSpaceEcon.StackedTimeSolver.sf_libs) sfdef.value == :default && continue From da432cef8549d5489dda9f9305fce9aaadc76485 Mon Sep 17 00:00:00 2001 From: Boyan Bejanov Date: Mon, 10 Apr 2023 20:53:56 -0400 Subject: [PATCH 15/32] codecov-action v3 --- .github/workflows/main.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index 507e806..57e0ad8 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -30,7 +30,7 @@ jobs: # Steps represent a sequence of tasks that will be executed as part of the job steps: # Checks-out your repository under $GITHUB_WORKSPACE, so your job can access it - - uses: actions/checkout@v2 + - uses: actions/checkout@v3 - name: Setup Julia environment # You may pin to the exact commit or the version. @@ -64,6 +64,6 @@ jobs: # process and submit the code coverage information - uses: julia-actions/julia-processcoverage@v1 - - uses: codecov/codecov-action@v1 + - uses: codecov/codecov-action@v3 with: file: lcov.info From d35e3947345390ab1d21894e05617b502b86fd05 Mon Sep 17 00:00:00 2001 From: Boyan Bejanov Date: Mon, 10 Apr 2023 21:27:28 -0400 Subject: [PATCH 16/32] factor out try-catch blocks in sparse.jl --- Project.toml | 1 + src/stackedtime/simulate.jl | 10 ----- src/stackedtime/sparse.jl | 77 +++++++++++++------------------------ 3 files changed, 27 insertions(+), 61 deletions(-) diff --git a/Project.toml b/Project.toml index c2dce00..84bcdcc 100644 --- a/Project.toml +++ b/Project.toml @@ -26,6 +26,7 @@ JLD2 = "0.4" ModelBaseEcon = "0.5.2" Pardiso = "0.5" Suppressor = "0.2" +TimerOutputs = "0.5" TimeSeriesEcon = "0.5.1" julia = "1.7" diff --git a/src/stackedtime/simulate.jl b/src/stackedtime/simulate.jl index 6013aac..26744ed 100644 --- a/src/stackedtime/simulate.jl +++ b/src/stackedtime/simulate.jl @@ -9,16 +9,6 @@ return model end -#= -@timeit_debug timer function sparse_solve(A, b) - if A isa PardisoFactorization - pardiso_solve!(A, copy(b)) - else - A \ b - end -end -=# - """ sim_nr!(x, sd, maxiter, tol, verbose [, linesearch]) diff --git a/src/stackedtime/sparse.jl b/src/stackedtime/sparse.jl index 65b5cd6..11fc1f3 100644 --- a/src/stackedtime/sparse.jl +++ b/src/stackedtime/sparse.jl @@ -10,8 +10,8 @@ # selected sparse linear algebra library is a Symbol const sf_libs = ( - :default, # use Julia's standard library (UMFPACK) - :umfpack, # same as :default + :default, # + :umfpack, # use Julia's standard library (UMFPACK) :pardiso, # use Pardiso - the one included with MKL ) @@ -19,14 +19,14 @@ global sf_default = :umfpack # a function to initialize a Factorization instance # this is also a good place to do the symbolic analysis -sf_prepare(A::SparseMatrixCSC, sparse_lib::Symbol=:default) = sf_prepare(Val(sparse_lib), A) -sf_prepare(::Val{S}, args...) where {S} = throw(ArgumentError("Unknown sparse library $S. Try one of $(sf_libs).")) +# sf_prepare(A::SparseMatrixCSC, sparse_lib::Symbol=:default) = sf_prepare(Val(sparse_lib), A) +# sf_prepare(::Val{S}, args...) where {S} = throw(ArgumentError("Unknown sparse library $S. Try one of $(sf_libs).")) # a function to calculate the numerical factors -sf_factor!(f::Factorization, A::SparseMatrixCSC) = throw(ArgumentError("Unknown factorization type $(typeof(f)).")) +# sf_factor!(f::Factorization, A::SparseMatrixCSC) = throw(ArgumentError("Unknown factorization type $(typeof(f)).")) # a function to solve the linear system -sf_solve!(f::Factorization, x::AbstractArray) = throw(ArgumentError("Unknown factorization type $(typeof(f)).")) +# sf_solve!(f::Factorization, x::AbstractArray) = throw(ArgumentError("Unknown factorization type $(typeof(f)).")) ########################################################################### @@ -37,6 +37,20 @@ function _sf_same_sparse_pattern(A::SparseMatrixCSC, B::SparseMatrixCSC) return (A.m == B.m) && (A.n == B.n) && (A.colptr == B.colptr) && (A.rowval == B.rowval) end +macro _sf_check_factorize(exception, expression) + error = gensym("error") + return esc(quote + try + $expression + catch $error + if $error isa $exception + @error("The system is underdetermined with the given set of equations and final conditions.") + end + rethrow() + end + end) +end + mutable struct LUFactorization{Tv<:Real} <: Factorization{Tv} F::SuiteSparse.UMFPACK.UmfpackLU{Tv,Int} A::SparseMatrixCSC{Tv,Int} @@ -44,14 +58,7 @@ end @timeit_debug timer "sf_prepare_lu" function sf_prepare(::Val{:umfpack}, A::SparseMatrixCSC) Tv = eltype(A) - F = try - @timeit_debug timer "_lu_full" lu(A) - catch error - if error isa SingularException - @error("The system is underdetermined with the given set of equations and final conditions.") - end - rethrow() - end + F = @_sf_check_factorize(SingularException, @timeit_debug timer "_lu_full" lu(A)) return LUFactorization{Tv}(F, A) end @@ -64,26 +71,12 @@ end else # sparse pattern is the same, different numbers f.A = A - try - @timeit_debug timer "_lu_num" lu!(f.F, A) - catch error - if error isa SingularException - @error("The system is underdetermined with the given set of equations and final conditions.") - end - rethrow() - end + @_sf_check_factorize(SingularException, @timeit_debug timer "_lu_num" lu!(f.F, A)) end else # totally new matrix, start over f.A = A - f.F = try - @timeit_debug timer "_lu_full" lu(A) - catch error - if error isa SingularException - @error("The system is underdetermined with the given set of equations and final conditions.") - end - rethrow() - end + f.F = @_sf_check_factorize(SingularException, @timeit_debug timer "_lu_full" lu(A)) end return f end @@ -107,9 +100,6 @@ end set_matrixtype!(ps, Pardiso.REAL_NONSYM) pardisoinit(ps) fix_iparm!(ps, :N) - # set_iparm!(ps, 1, 1) # Override defaults - # set_iparm!(ps, 11, 1) # disable automatic scaling for non-symmetric matrices - # set_iparm!(ps, 13, 1) # disable automatic scaling for non-symmetric matrices # set_iparm!(ps, 2, 2) # Select algorithm pf = PardisoFactorization{Tv}(ps, get_matrix(ps, A, :N)) finalizer(pf) do x @@ -123,14 +113,9 @@ end @timeit_debug timer "_pardso_full" function _pardiso_full!(pf::PardisoFactorization) # run the analysis phase ps = pf.ps - try + @_sf_check_factorize Union{Pardiso.PardisoException,Pardiso.PardisoPosDefException} begin set_phase!(ps, Pardiso.ANALYSIS_NUM_FACT) pardiso(ps, pf.A, Float64[]) - catch error - if error isa Pardiso.PardisoException || error isa Pardiso.PardisoPosDefException - @error("The system is underdetermined with the given set of equations and final conditions.") - end - rethrow() end return pf end @@ -138,14 +123,9 @@ end @timeit_debug timer "_pardso_num" function _pardiso_numeric!(pf::PardisoFactorization) # run the analysis phase ps = pf.ps - try + @_sf_check_factorize Union{Pardiso.PardisoException,Pardiso.PardisoPosDefException} begin set_phase!(ps, Pardiso.NUM_FACT) pardiso(ps, pf.A, Float64[]) - catch error - if error isa Pardiso.PardisoException || error isa Pardiso.PardisoPosDefException - @error("The system is underdetermined with the given set of equations and final conditions.") - end - rethrow() end return pf end @@ -172,14 +152,9 @@ end @timeit_debug timer "sf_solve!_par" function sf_solve!(pf::PardisoFactorization, x::AbstractArray) ps = pf.ps - try + @_sf_check_factorize Union{Pardiso.PardisoException,Pardiso.PardisoPosDefException} begin set_phase!(ps, Pardiso.SOLVE_ITERATIVE_REFINE) pardiso(ps, x, pf.A, copy(x)) - catch error - if error isa Pardiso.PardisoException || error isa Pardiso.PardisoPosDefException - @error("The system is underdetermined with the given set of equations and final conditions.") - end - rethrow() end return x end From 25ea3d746ca8befb37c402d7a5d68ab327b868d0 Mon Sep 17 00:00:00 2001 From: Boyan Bejanov Date: Tue, 11 Apr 2023 00:06:42 -0400 Subject: [PATCH 17/32] better views --- src/stackedtime/shockdecomp.jl | 25 ++++++++++++------------- src/stackedtime/simulate.jl | 22 +++++++++++----------- src/stackedtime/stoch_simulate.jl | 3 ++- 3 files changed, 25 insertions(+), 25 deletions(-) diff --git a/src/stackedtime/shockdecomp.jl b/src/stackedtime/shockdecomp.jl index 860b3b9..064e3bb 100644 --- a/src/stackedtime/shockdecomp.jl +++ b/src/stackedtime/shockdecomp.jl @@ -108,7 +108,7 @@ function shockdecomp(m::Model, p::Plan, exog_data::SimData; end # indices for the exogenous variable-points split into groups - exog_inds = Workspace(; + exog_inds = @views @inbounds Workspace(; # contributions of initial conditions init=vec(LI[init, :]), # contributions of final conditions @@ -148,7 +148,7 @@ function shockdecomp(m::Model, p::Plan, exog_data::SimData; continue end rsdv = result.sd[v.name] = MVTSeries(p.range, keys(exog_inds), zeros) - rsdv[pinit, :init] = delta[pinit, v.name] + rsdv.init[pinit] = view(delta[v.name], pinit) end else shex = filter(v -> isshock(v) || isexog(v), m.allvars) @@ -160,13 +160,13 @@ function shockdecomp(m::Model, p::Plan, exog_data::SimData; if haskey(id, v.name) # have initial decomposition idv = id[v.name] - if norm(delta[pinit, v.name] - sum(idv[pinit, :], dims=2), Inf) > tol + if norm(view(delta[v.name], pinit) - sum(view(idv, pinit, :), dims=2), Inf) > tol error("The given `initdecomp` does not add up for $(v.name).") end rsdv[pinit] .= idv for (colind, colname) in enumerate(keys(exog_inds)) if haskey(idv, colname) || colname ∈ shex - inds = LI[init, vind] + inds = @view LI[init, vind] append!(SDI, inds) append!(SDJ, fill(colind, size(inds))) append!(SDV, idv[pinit, colname].values) @@ -174,9 +174,9 @@ function shockdecomp(m::Model, p::Plan, exog_data::SimData; end else # no initial decomposition - it all goes in init - rsdv[pinit, :init] = delta[pinit, v.name] + rsdv.init[pinit] = view(delta[v.name], pinit) colind = 1 - inds = LI[init, vind] + inds = @view LI[init, vind] append!(SDI, inds) append!(SDJ, fill(colind, size(inds))) append!(SDV, delta[inds]) @@ -202,21 +202,20 @@ function shockdecomp(m::Model, p::Plan, exog_data::SimData; begin # in order to split the rows of SDMAT by variable, we need the inverse indexing map inv_endo_inds = zeros(Int, size(gdata.solve_mask)) - inv_endo_inds[gdata.solve_mask] .= 1:sum(gdata.solve_mask) + inv_endo_inds[gdata.solve_mask] = 1:sum(gdata.solve_mask) for (index, v) in enumerate(m.allvars) - v_inds = vec(LI[:, index]) - v_endo_mask = gdata.solve_mask[v_inds] + v_endo_mask = gdata.solve_mask[@view LI[:, index]] if !any(v_endo_mask) continue end - v_inv_endo_inds = inv_endo_inds[v_inds[v_endo_mask]] + v_inv_endo_inds = inv_endo_inds[@view LI[v_endo_mask, index]] v_data = result.sd[v] # assign contributions to endogenous points v_data[v_endo_mask, :] = SDMAT[v_inv_endo_inds, :] # contributions of initial conditions already assigned # assign contributions to final conditions if fctype === fcgiven || fctype === fclevel - v_data[(begin-1).+term, :term] .= delta[v.name][term] + v_data.term[(begin-1).+term] = view(delta[v.name], term) elseif fctype === fcrate for tt in term v_data[begin-1+tt, :] = v_data[begin-2+tt, :] @@ -238,8 +237,8 @@ function shockdecomp(m::Model, p::Plan, exog_data::SimData; if deviation logvars = islog.(m.varshks) .| isneglog.(m.varshks) - result.s[:, logvars] ./= result.c[:, logvars] - result.s[:, .!logvars] .-= result.c[:, .!logvars] + @views result.s[:, logvars] ./= result.c[:, logvars] + @views result.s[:, .!logvars] .-= result.c[:, .!logvars] end if _debug diff --git a/src/stackedtime/simulate.jl b/src/stackedtime/simulate.jl index 26744ed..cc32cee 100644 --- a/src/stackedtime/simulate.jl +++ b/src/stackedtime/simulate.jl @@ -125,7 +125,7 @@ function simulate(m::Model, NT = length(p_ant.range) nauxs = length(m.auxvars) nvarshks = length(m.varshks) - logvars = [islog(var) | isneglog(var) for var in m.varshks] + logvars = islog.(m.varshks) .| isneglog.(m.varshks) if size(exog_ant) != (NT, nvarshks) error("Incorrect dimensions of exog_data. Expected $((NT, nvarshks)), got $(size(exog_ant)).") @@ -142,8 +142,8 @@ function simulate(m::Model, if size(baseline) != (NT, nvarshks) error("Incorrect dimensions of baseline. Expected $((NT, nvarshks)), got $(size(baseline)).") end - exog_ant[:, logvars] .*= baseline[:, logvars] - exog_ant[:, .!logvars] .+= baseline[:, .!logvars] + @views exog_ant[:, logvars] .*= baseline[:, logvars] + @views exog_ant[:, .!logvars] .+= baseline[:, .!logvars] end exog_ant = ModelBaseEcon.update_auxvars(transform(exog_ant, m), m) @@ -193,8 +193,8 @@ function simulate(m::Model, error("Anticipated and unanticipated ranges don't match.") end if deviation_unant - exog_unant[:, logvars] .*= baseline[:, logvars] - exog_unant[:, .!logvars] .+= baseline[:, .!logvars] + @views exog_unant[:, logvars] .*= baseline[:, logvars] + @views exog_unant[:, .!logvars] .+= baseline[:, .!logvars] end exog_unant = ModelBaseEcon.update_auxvars(transform(exog_unant, m), m) if size(exog_unant) != size(exog_ant) @@ -209,7 +209,7 @@ function simulate(m::Model, x[sim, shkinds] .= 0 end - x[init, allvarinds] .= exog_ant[init, allvarinds] + x[init, allvarinds] = exog_ant[init, allvarinds] t0 = first(sim) T = last(sim) if expectation_horizon === nothing @@ -222,7 +222,7 @@ function simulate(m::Model, if t == t0 imaxiter = maxiter itol = tol - elseif (maximum(abs, x[t, exog_inds] .- exog_unant[t, exog_inds]) < tol) #= && (psim[t0, Val(:inds)] == exog_inds) =# + elseif (maximum(abs, x[t, exog_inds] - exog_unant[t, exog_inds]) < tol) #= && (psim[t0, Val(:inds)] == exog_inds) =# continue else imaxiter = 5 @@ -300,14 +300,14 @@ function simulate(m::Model, # these intermediate simulations are always with fcnatural, # have length equal to expectation_horizon and # only the first period is imposed - if (maximum(abs, x[t, exog_inds] .- exog_unant[t, exog_inds]) < tol) #= && (exog_inds == shkinds) =# + if (maximum(abs, x[t, exog_inds] - exog_unant[t, exog_inds]) < tol) #= && (exog_inds == shkinds) =# continue end psim1 = copy(psim) # the range of psim1 might extend beyond the range of p_ant. # we copy from p_ant as far as we have and copy the last line beyond that tmp_rng = t:min(t + expectation_horizon - 1, T) - psim1.exogenous[t0.+(0:length(tmp_rng)-1), :] .= p_ant.exogenous[tmp_rng, :] + psim1.exogenous[t0.+(0:length(tmp_rng)-1), :] = p_ant.exogenous[tmp_rng, :] # ===> must leave the psim1 plan empty beyond the end of p_ant # becasue we don't have data in exog_and for any exogenized @@ -370,8 +370,8 @@ function simulate(m::Model, x = x[axes(exog_ant)...] x .= inverse_transform(x, m) if deviation - x[:, logvars] ./= baseline[:, logvars] - x[:, .!logvars] .-= baseline[:, .!logvars] + @views x[:, logvars] ./= baseline[:, logvars] + @views x[:, .!logvars] .-= baseline[:, .!logvars] end return x diff --git a/src/stackedtime/stoch_simulate.jl b/src/stackedtime/stoch_simulate.jl index b9a794a..28a0022 100644 --- a/src/stackedtime/stoch_simulate.jl +++ b/src/stackedtime/stoch_simulate.jl @@ -136,9 +136,10 @@ function stoch_simulate(m::Model, p::Plan, baseline::SimData, shocks; end # time loop # strip auxvar columns and inverse transform + have_auxs = (m.nauxs > 0) for (key, result) in pairs(results) isfailed(result) && continue - if m.nauxs > 0 + if have_auxs result = result[:, m.varshks] end results[key] = inverse_transform(result, m) From c364fb5fba3865d36ebcc3fb58f7c08a33c85ac9 Mon Sep 17 00:00:00 2001 From: Boyan Bejanov Date: Tue, 11 Apr 2023 12:00:21 -0400 Subject: [PATCH 18/32] export use_umfpack() and use_pardiso() --- src/StackedTimeSolver.jl | 3 +++ src/stackedtime/sparse.jl | 42 +++++++++++++++++++++++++++++++++++++-- test/runtests.jl | 9 +++++++++ 3 files changed, 52 insertions(+), 2 deletions(-) diff --git a/src/StackedTimeSolver.jl b/src/StackedTimeSolver.jl index 6e61edf..d5f394f 100644 --- a/src/StackedTimeSolver.jl +++ b/src/StackedTimeSolver.jl @@ -72,6 +72,9 @@ export FCMatchSSRate, fcslope, fcrate export FCConstRate, fcnatural export setfc +export use_pardiso, use_pardiso! +export use_umfpack, use_umfpack! + # the following are deprecated export seriesoverlay, dictoverlay diff --git a/src/stackedtime/sparse.jl b/src/stackedtime/sparse.jl index 11fc1f3..901f992 100644 --- a/src/stackedtime/sparse.jl +++ b/src/stackedtime/sparse.jl @@ -6,8 +6,6 @@ ################################################################################## -### API - # selected sparse linear algebra library is a Symbol const sf_libs = ( :default, # @@ -17,6 +15,46 @@ const sf_libs = ( global sf_default = :umfpack +""" + use_umfpack() + +Set the default sparse factorization library to UMFPACK (the one used in Julia's +standard library). See also [`use_pardiso`](@ref). + +""" +@inline use_umfpack() = (global sf_default = :umfpack; nothing) +# @inline use_umfpack(m::Model) = use_umfpack!(m) + +""" + use_umfpack!(model) + +Instruct the stacked-time solver to use Pardiso with this model. See also +[`use_pardiso!`](@ref). +""" +@inline use_umfpack!(m::Model) = (m.options.factorization = :umfpack; m) +# @inline use_umfpack!() = use_umfpack() +export use_umfpack, use_umfpack! + +""" + use_pardiso() + +Set the default sparse factorization library to Pardiso. See also +[`use_umfpack`](@ref). +""" +@inline use_pardiso() = (global sf_default = :pardiso; nothing) +# @inline use_pardiso(m::Model) = use_pardiso!(m) + +""" + use_pardiso!(model) + +Instruct the stacked-time solver to use Pardiso with this model. +""" +@inline use_pardiso!(m::Model) = (m.options.factorization = :pardiso; m) +# @inline use_pardiso!() = use_pardiso() +export use_pardiso, use_pardiso! + +### API + # a function to initialize a Factorization instance # this is also a good place to do the symbolic analysis # sf_prepare(A::SparseMatrixCSC, sparse_lib::Symbol=:default) = sf_prepare(Val(sparse_lib), A) diff --git a/test/runtests.jl b/test/runtests.jl index 4965b20..1642550 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -218,6 +218,15 @@ end include("misc.jl") +@testset "sparse" begin + @test (use_pardiso(); StateSpaceEcon.StackedTimeSolver.sf_default == :pardiso) + @test (use_umfpack(); StateSpaceEcon.StackedTimeSolver.sf_default == :umfpack) + let m = Model() + @test (use_pardiso!(m); m.factorization == :pardiso) + @test (use_umfpack!(m); m.factorization == :umfpack) + end +end + # make sure Pardiso runs deterministic (single thread) for the tests Pardiso.set_nprocs_mkl!(1) From 40c5991aa69009c69ca4dcaa983b3395283ac527 Mon Sep 17 00:00:00 2001 From: Boyan Bejanov Date: Wed, 12 Apr 2023 10:55:24 -0400 Subject: [PATCH 19/32] use_umfpack! and use_pardiso! return m.options rather than m --- src/stackedtime/sparse.jl | 8 ++------ test/runtests.jl | 2 +- 2 files changed, 3 insertions(+), 7 deletions(-) diff --git a/src/stackedtime/sparse.jl b/src/stackedtime/sparse.jl index 901f992..65c05f6 100644 --- a/src/stackedtime/sparse.jl +++ b/src/stackedtime/sparse.jl @@ -23,7 +23,6 @@ standard library). See also [`use_pardiso`](@ref). """ @inline use_umfpack() = (global sf_default = :umfpack; nothing) -# @inline use_umfpack(m::Model) = use_umfpack!(m) """ use_umfpack!(model) @@ -31,8 +30,7 @@ standard library). See also [`use_pardiso`](@ref). Instruct the stacked-time solver to use Pardiso with this model. See also [`use_pardiso!`](@ref). """ -@inline use_umfpack!(m::Model) = (m.options.factorization = :umfpack; m) -# @inline use_umfpack!() = use_umfpack() +@inline use_umfpack!(m::Model) = (m.options.factorization = :umfpack; m.options) export use_umfpack, use_umfpack! """ @@ -42,15 +40,13 @@ Set the default sparse factorization library to Pardiso. See also [`use_umfpack`](@ref). """ @inline use_pardiso() = (global sf_default = :pardiso; nothing) -# @inline use_pardiso(m::Model) = use_pardiso!(m) """ use_pardiso!(model) Instruct the stacked-time solver to use Pardiso with this model. """ -@inline use_pardiso!(m::Model) = (m.options.factorization = :pardiso; m) -# @inline use_pardiso!() = use_pardiso() +@inline use_pardiso!(m::Model) = (m.options.factorization = :pardiso; m.options) export use_pardiso, use_pardiso! ### API diff --git a/test/runtests.jl b/test/runtests.jl index 1642550..1c4278a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -235,7 +235,7 @@ for sfdef = QuoteNode.(StateSpaceEcon.StackedTimeSolver.sf_libs) sfdef.value == :default && continue # Pardiso in macos is giving "Not enough memory". Disable for now - sfdef.value == :pardiso && Sys.isapple() && continue + # sfdef.value == :pardiso && Sys.isapple() && continue @info "Using $(sfdef)" From 2c2849aadd1df3cdc893998b12c0c3c80b0f605c Mon Sep 17 00:00:00 2001 From: Jason Jensen Date: Wed, 12 Apr 2023 14:02:57 -0400 Subject: [PATCH 20/32] uncomment precompile line --- src/steadystate/presolve.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/steadystate/presolve.jl b/src/steadystate/presolve.jl index 746b9e2..604f4ed 100644 --- a/src/steadystate/presolve.jl +++ b/src/steadystate/presolve.jl @@ -114,5 +114,5 @@ presolve_sstate!(eqns::OrderedDict{Symbol,SteadyStateEquation}, mask::AbstractVe _presolve_equations!(eqns, mask, values, method, verbose, tol) @assert precompile(presolve_sstate!, (Model, Vector{Bool}, Vector{Float64})) -# @assert precompile(presolve_sstate!, (LittleDict{Symbol,SteadyStateEquation}, Vector{Bool}, Vector{Float64})) +@assert precompile(presolve_sstate!, (OrderedDict{Symbol,SteadyStateEquation}, Vector{Bool}, Vector{Float64})) From 200de729e45b680ab77d3de3e859b91a90c86a68 Mon Sep 17 00:00:00 2001 From: Boyan Bejanov Date: Mon, 17 Apr 2023 14:12:45 -0400 Subject: [PATCH 21/32] bugfix for stoch_simulate result container --- src/stackedtime/stoch_simulate.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/stackedtime/stoch_simulate.jl b/src/stackedtime/stoch_simulate.jl index 28a0022..06e504f 100644 --- a/src/stackedtime/stoch_simulate.jl +++ b/src/stackedtime/stoch_simulate.jl @@ -8,7 +8,8 @@ # version of simulate for stochastic simulations # NOTE: This requires TimeSeriesEcon v0.5.1 where map(func, ::Workspace) returns a Workspace. -_make_result_container(shocks, basedata) = map(_ -> copy(basedata), shocks) +_make_result_container(shocks::Workspace, basedata)::Workspace = map(_ -> copy(basedata), shocks) +_make_result_container(shocks::Vector, basedata)::Vector{MaybeSimData} = convert(Vector{MaybeSimData}, map(_ -> copy(basedata), shocks)) function stoch_simulate(m::Model, p::Plan, baseline::SimData, shocks; check::Bool=false, From f9b65fa0c83177a936dcbe71927edc8207a2ac7e Mon Sep 17 00:00:00 2001 From: Boyan Bejanov Date: Mon, 17 Apr 2023 19:07:08 -0400 Subject: [PATCH 22/32] removing empty TODO.md --- TODO.md | 2 -- 1 file changed, 2 deletions(-) delete mode 100644 TODO.md diff --git a/TODO.md b/TODO.md deleted file mode 100644 index 339d038..0000000 --- a/TODO.md +++ /dev/null @@ -1,2 +0,0 @@ - -1. From 12d7a3494ed9c276037d9a8f50921f727ed42913 Mon Sep 17 00:00:00 2001 From: Boyan Bejanov Date: Tue, 18 Apr 2023 10:23:14 -0400 Subject: [PATCH 23/32] change algorithm to 3 --- src/stackedtime/sparse.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/stackedtime/sparse.jl b/src/stackedtime/sparse.jl index 65c05f6..6cbbbde 100644 --- a/src/stackedtime/sparse.jl +++ b/src/stackedtime/sparse.jl @@ -134,7 +134,7 @@ end set_matrixtype!(ps, Pardiso.REAL_NONSYM) pardisoinit(ps) fix_iparm!(ps, :N) - # set_iparm!(ps, 2, 2) # Select algorithm + set_iparm!(ps, 2, 3) # Select algorithm pf = PardisoFactorization{Tv}(ps, get_matrix(ps, A, :N)) finalizer(pf) do x set_phase!(x.ps, Pardiso.RELEASE_ALL) From a7fe69f268bd92e2f976fb0d7813d71d34d130a7 Mon Sep 17 00:00:00 2001 From: Nic2020 Date: Fri, 21 Apr 2023 12:37:39 -0400 Subject: [PATCH 24/32] Remove itol/imaxiter (#48) * Remove itol/imaxiter --- src/stackedtime/simulate.jl | 16 +++++----------- 1 file changed, 5 insertions(+), 11 deletions(-) diff --git a/src/stackedtime/simulate.jl b/src/stackedtime/simulate.jl index cc32cee..6dc3fa9 100644 --- a/src/stackedtime/simulate.jl +++ b/src/stackedtime/simulate.jl @@ -219,14 +219,8 @@ function simulate(m::Model, exog_inds = p_unant[t, Val(:inds)] psim = Plan(m, t:T) psim.exogenous .= p_ant.exogenous[begin+Int(t - t0):end, :] - if t == t0 - imaxiter = maxiter - itol = tol - elseif (maximum(abs, x[t, exog_inds] - exog_unant[t, exog_inds]) < tol) #= && (psim[t0, Val(:inds)] == exog_inds) =# + if t !== t0 && (maximum(abs, x[t, exog_inds] - exog_unant[t, exog_inds]) < tol) #= && (psim[t0, Val(:inds)] == exog_inds) =# continue - else - imaxiter = 5 - itol = sqrt(tol) end setexog!(psim, t0, exog_inds) gdata = StackedTimeSolverData(m, psim, fctype, variant) @@ -236,11 +230,11 @@ function simulate(m::Model, xx = view(x, sim_range, :) assign_final_condition!(xx, exog_unant[sim_range, :], gdata) if verbose - @info "Simulating $(p_ant.range[t:T]) with $((itol, imaxiter))" # anticipate expectation_horizon gdata.FC + @info "Simulating $(p_ant.range[t:T]) with $((tol, maxiter))" # anticipate expectation_horizon gdata.FC end - converged = sim_nr!(xx, gdata, imaxiter, itol, verbose, linesearch) + converged = sim_nr!(xx, gdata, maxiter, tol, verbose, linesearch) if warn_maxiter && !converged - @warn("Newton-Raphson reached maximum number of iterations (`imaxiter`).") + @warn("Newton-Raphson reached maximum number of iterations (`maxiter`).") end last_run = Workspace(; t, xx, gdata) end @@ -333,7 +327,7 @@ function simulate(m::Model, if verbose @info("Simulating $(p_ant.range[t] .+ (0:expectation_horizon - 1))") # anticipate expectation_horizon sdata.FC end - converged = sim_nr!(xx, sdata, 5, sqrt(tol), verbose, linesearch) + converged = sim_nr!(xx, sdata, maxiter, tol, verbose, linesearch) if warn_maxiter && !converged @warn("Newton-Raphson reached maximum number of iterations (`maxiter`).") end From a45d4ae0d20f345b7ea0c7d7cd7d96d77f5bcb6a Mon Sep 17 00:00:00 2001 From: Boyan Bejanov Date: Thu, 27 Apr 2023 11:46:35 -0400 Subject: [PATCH 25/32] index plan with ModelVariable --- src/Plans.jl | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/src/Plans.jl b/src/Plans.jl index 650988e..ffacc80 100644 --- a/src/Plans.jl +++ b/src/Plans.jl @@ -78,7 +78,7 @@ function Plan(model::Model, range::AbstractUnitRange) range = (first(range)-model.maxlag):(last(range)+model.maxlead) local varshks = model.varshks local N = length(varshks) - local names = tuple(Symbol.(varshks)...) # force conversion of Vector{ModelSymbol} to NTuple{N,Symbol} + local names = tuple(Symbol.(varshks)...) # force conversion of Vector{ModelVariable} to NTuple{N,Symbol} p = Plan(range, NamedTuple{names}(1:N), falses(length(range), length(varshks))) for (ind, var) = enumerate(varshks) if isshock(var) || isexog(var) @@ -146,12 +146,16 @@ Base.setindex!(p::Plan, x, i...) = error("Cannot assign directly. Use `exogenize ####################################### # query the exo-end status of a variable -@inline Base.getindex(p::Plan{T}, vars::Symbol...) where {T} = begin +@inline Base.getindex(p::Plan, vars::Union{AbstractString,Symbol,ModelVariable}...) = Base.getindex(p, Symbol[vars...]) +@inline Base.getindex(p::Plan, vars::AbstractVector{<:Union{AbstractString,ModelVariable}}) = Base.getindex(p, Symbol[vars...]) +@inline Base.getindex(p::Plan{T}, vars::AbstractVector{Symbol}) where {T} = begin var_inds = Int[p.varshks[vars]...] Plan{T}(p.range, NamedTuple{(vars...,)}(eachindex(vars)), p.exogenous[:, var_inds]) end -@inline Base.getindex(p::Plan{T}, rng::AbstractUnitRange{T}, vars::Symbol...) where {T} = begin +@inline Base.getindex(p::Plan, rng::AbstractUnitRange, vars::Union{AbstractString,Symbol,ModelVariable}...) = Base.getindex(p, rng, Symbol[vars...]) +@inline Base.getindex(p::Plan, rng::AbstractUnitRange, vars::AbstractVector{<:Union{AbstractString,ModelVariable}}) = Base.getindex(p, rng, Symbol[vars...]) +@inline Base.getindex(p::Plan{T}, rng::AbstractUnitRange{T}, vars::AbstractVector{Symbol}) where {T} = begin rng.start < p.range.start && throw(BoundsError(p, rng.start)) rng.stop > p.range.stop && throw(BoundsError(p, rng.stop)) var_inds = Int[p.varshks[vars]...] From 344c527a7ae7acfe88fbe370f2c237627897b3f8 Mon Sep 17 00:00:00 2001 From: Boyan Bejanov Date: Mon, 1 May 2023 15:03:40 -0400 Subject: [PATCH 26/32] bugfix: stoch_simulate() with empty shocks - now it returns empty results instead of crashing --- src/stackedtime/stoch_simulate.jl | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/stackedtime/stoch_simulate.jl b/src/stackedtime/stoch_simulate.jl index 06e504f..fc8dfc8 100644 --- a/src/stackedtime/stoch_simulate.jl +++ b/src/stackedtime/stoch_simulate.jl @@ -24,6 +24,10 @@ function stoch_simulate(m::Model, p::Plan, baseline::SimData, shocks; warn_maxiter::Bool=getoption(getoption(m, :warn, Options()), :maxiter, false) ) + if isempty(shocks) + return _make_result_container(shocks, baseline); + end + # get the range of all shocks realizations (this will be the full simulation range) shkrng = mapreduce(rangeof, union, values(shocks)) From 1633ff48472f1ec652026677adde2543b7671a92 Mon Sep 17 00:00:00 2001 From: Boyan Bejanov Date: Mon, 1 May 2023 16:50:28 -0400 Subject: [PATCH 27/32] warn if stochastic shock is endogenous (not error) --- src/stackedtime/stoch_simulate.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/stackedtime/stoch_simulate.jl b/src/stackedtime/stoch_simulate.jl index fc8dfc8..6da7a1e 100644 --- a/src/stackedtime/stoch_simulate.jl +++ b/src/stackedtime/stoch_simulate.jl @@ -51,7 +51,7 @@ function stoch_simulate(m::Model, p::Plan, baseline::SimData, shocks; tinds = Plans._offset(p, rangeof(val)) vind = p.varshks[shk] if !all(p.exogenous[tinds, vind]) - throw(ArgumentError("$shk in shocks realization $key is endogenous in the given plan.")) + @warn "$shk in shocks[$key] is endogenous in the given plan." end end end From fcce6632458e6963f6c5d4c6ac485c587c2e722f Mon Sep 17 00:00:00 2001 From: Jason Jensen Date: Mon, 29 May 2023 10:29:20 -0400 Subject: [PATCH 28/32] add equation change tests --- test/simtests.jl | 91 ++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 91 insertions(+) diff --git a/test/simtests.jl b/test/simtests.jl index 4e48cec..3c240ad 100644 --- a/test/simtests.jl +++ b/test/simtests.jl @@ -572,3 +572,94 @@ end end end + +############################################################# +# Tests of equation changes +@testset "E1 to E2.sim" begin + m = E1.newmodel() + @parameters m begin + cp = [0.5, 0.02] + cr = [0.75, 1.5, 0.5] + cy = [0.5, -0.02] + end + @variables m begin + @delete y + pinf + rate + ygap + end + @shocks m begin + @delete y_shk + pinf_shk + rate_shk + ygap_shk + end + @autoexogenize m begin + @delete y => y_shk + pinf = pinf_shk + rate = rate_shk + ygap = ygap_shk + end + @equations m begin + @delete _EQ1; + pinf[t]=cp[1]*pinf[t-1]+(.98-cp[1])*pinf[t+1]+cp[2]*ygap[t]+pinf_shk[t] + rate[t]=cr[1]*rate[t-1]+(1-cr[1])*(cr[2]*pinf[t]+cr[3]*ygap[t])+rate_shk[t] + ygap[t]=cy[1]*ygap[t-1]+(.98-cy[1])*ygap[t+1]+cy[2]*(rate[t]-pinf[t+1])+ygap_shk[t] + end + + @reinitialize m + + test_simulation(m, "data/M2_TestData.csv") +end + +@testset "E2 to E3.sim" begin + m = E2.newmodel() + @equations m begin + :_EQ1 => pinf[t]=cp[1]*pinf[t-1]+0.3*pinf[t+1]+0.05*pinf[t+2]+0.05*pinf[t+3]+cp[2]*ygap[t]+pinf_shk[t] + :_EQ2 => rate[t]=cr[1]*rate[t-1]+(1-cr[1])*(cr[2]*pinf[t]+cr[3]*ygap[t])+rate_shk[t] + :_EQ3 => ygap[t]=cy[1]/2*ygap[t-2]+cy[1]/2*ygap[t-1]+(.98-cy[1])*ygap[t+1]+cy[2]*(rate[t]-pinf[t+1])+ygap_shk[t] + end + + @reinitialize m + + test_simulation(m, "data/M3_TestData.csv") +end + +@testset "E3 to E6.sim" begin + m = E3.newmodel() + @parameters m begin + p_dlp = 0.0050000000000000 + p_dly = 0.0045000000000000 + end + @variables m begin + @delete pinf rate ygap + dlp; dly; dlyn; lp; ly; lyn + end + @shocks m begin + @delete pinf_shk rate_shk ygap_shk + dlp_shk; dly_shk + end + @autoexogenize m begin + @delete pinf => pinf_shk + @delete rate => rate_shk + @delete ygap => ygap_shk + ly = dly_shk + lp = dlp_shk + end + + @equations m begin + @delete _EQ1; + @delete _EQ2; + @delete _EQ3; + dly[t]=(1-0.2-0.2)*p_dly+0.2*dly[t-1]+0.2*dly[t+1]+dly_shk[t] + dlp[t]=(1-0.5)*p_dlp+0.1*dlp[t-2]+0.1*dlp[t-1]+0.1*dlp[t+1]+0.1*dlp[t+2]+0.1*dlp[t+3]+dlp_shk[t] + dlyn[t]=dly[t]+dlp[t] + ly[t]=ly[t-1]+dly[t] + lp[t]=lp[t-1]+dlp[t] + lyn[t]=lyn[t-1]+dlyn[t] + end + + @reinitialize m + + test_simulation(m, "data/M6_TestData.csv") +end From 05add76f3fcea9b72eb895d01a130c90aa120770 Mon Sep 17 00:00:00 2001 From: Jason Jensen Date: Tue, 6 Jun 2023 06:15:15 -0400 Subject: [PATCH 29/32] add extensive model change tests --- test/dynss.jl | 4 +- test/misc.jl | 2 +- test/modelchanges.jl | 530 +++++++++++++++++++++++++++++++++++++++++++ test/runtests.jl | 16 +- test/sim_fo.jl | 14 +- test/simtests.jl | 128 ++--------- test/sstests.jl | 16 +- test/stochsims.jl | 2 +- 8 files changed, 582 insertions(+), 130 deletions(-) create mode 100644 test/modelchanges.jl diff --git a/test/dynss.jl b/test/dynss.jl index 203d33e..a4cd834 100644 --- a/test/dynss.jl +++ b/test/dynss.jl @@ -6,7 +6,7 @@ ################################################################################## @testset "dynss" begin - m = deepcopy(S1.model) + m = getS1() # reset parameters m.α = 0.5 m.β = 0.6 @@ -48,7 +48,7 @@ end @testset "dynss2" begin - m = deepcopy(S2.model) + m = getS2() m.α = 0.4 diff --git a/test/misc.jl b/test/misc.jl index 0483e03..6a50575 100644 --- a/test/misc.jl +++ b/test/misc.jl @@ -13,7 +13,7 @@ end @testset "fcset" begin - m = deepcopy(E7A.model) + m = getE7A() empty!(m.sstate.constraints) @steadystate m lc = 1 diff --git a/test/modelchanges.jl b/test/modelchanges.jl new file mode 100644 index 0000000..e8789ba --- /dev/null +++ b/test/modelchanges.jl @@ -0,0 +1,530 @@ +# This set of tests reruns the existing test, but uses example models which have been constructed +# from the existing sample models. + +function constructE1!(model::Model) + model.flags.linear = true + model.options.substitutions = false + + @parameters model begin + α = 0.5 + β = 0.5 + end + + @variables model y + + @shocks model y_shk + + @autoexogenize model y = y_shk + + @equations model begin + y[t] = α * y[t - 1] + β * y[t + 1] + y_shk[t] + end + + @reinitialize model + + return model +end + +function constructE2!(model::Model) + model.flags.linear = true + model.options.substitutions = false + + # add parameters + @parameters model begin + cp = [0.5, 0.02] + cr = [0.75, 1.5, 0.5] + cy = [0.5, -0.02] + end + + # add variables: a list of symbols + @variables model begin + pinf + rate + ygap + end + + # add shocks: a list of symbols + @shocks model begin + pinf_shk + rate_shk + ygap_shk + end + + # autoexogenize: define a mapping of variables to shocks + @autoexogenize model begin + pinf = pinf_shk + rate = rate_shk + ygap = ygap_shk + end + + # add equations: a sequence of expressions, such that + # use y[t+1] for expectations/leads + # use y[t] for contemporaneous + # use y[t-1] for lags + # each expression must have exactly one "=" + @equations model begin + pinf[t]=cp[1]*pinf[t-1]+(.98-cp[1])*pinf[t+1]+cp[2]*ygap[t]+pinf_shk[t] + rate[t]=cr[1]*rate[t-1]+(1-cr[1])*(cr[2]*pinf[t]+cr[3]*ygap[t])+rate_shk[t] + ygap[t]=cy[1]*ygap[t-1]+(.98-cy[1])*ygap[t+1]+cy[2]*(rate[t]-pinf[t+1])+ygap_shk[t] + end + + # call initialize! to construct internal structures + @reinitialize model + return model +end + +function constructE3!(model::Model) + model.flags.linear = true + model.options.substitutions = false + + # add parameters + @parameters model begin + cp = [0.5, 0.02] + cr = [0.75, 1.5, 0.5] + cy = [0.5, -0.02] + end + + # add variables: a list of symbols + @variables model begin + pinf + rate + ygap + end + + # add shocks: a list of symbols + @shocks model begin + pinf_shk + rate_shk + ygap_shk + end + + # autoexogenize: define a mapping of variables to shocks + @autoexogenize model begin + pinf = pinf_shk + rate = rate_shk + ygap = ygap_shk + end + + # add equations: a sequence of expressions, such that + # use y[t+1] for expectations/leads + # use y[t] for contemporaneous + # use y[t-1] for lags + # each expression must have exactly one "=" + @equations model begin + pinf[t]=cp[1]*pinf[t-1]+0.3*pinf[t+1]+0.05*pinf[t+2]+0.05*pinf[t+3]+cp[2]*ygap[t]+pinf_shk[t] + rate[t]=cr[1]*rate[t-1]+(1-cr[1])*(cr[2]*pinf[t]+cr[3]*ygap[t])+rate_shk[t] + ygap[t]=cy[1]/2*ygap[t-2]+cy[1]/2*ygap[t-1]+(.98-cy[1])*ygap[t+1]+cy[2]*(rate[t]-pinf[t+1])+ygap_shk[t] + end + + # call initialize! to construct internal structures + @reinitialize model + return model +end + +function constructE3nl!(model::Model) + model.flags.linear = false + model.options.maxiter = 200 + model.options.substitutions = false + + @parameters model begin + cp = [0.5, 0.02] + cr = [0.75, 1.5, 0.5] + cy = [0.5, -0.02] + end + + @variables model begin + pinf + rate + ygap + end + + @shocks model begin + pinf_shk + rate_shk + ygap_shk + end + + @autoexogenize model begin + pinf = pinf_shk + rate = rate_shk + ygap = ygap_shk + end + + @equations model begin + exp(pinf[t])=exp(cp[1]*pinf[t-1]+0.3*pinf[t+1]+0.05*pinf[t+2]+0.05*pinf[t+3]+cp[2]*ygap[t]+pinf_shk[t]) + rate[t]=cr[1]*rate[t-1]+(1-cr[1])*(cr[2]*pinf[t]+cr[3]*ygap[t])+rate_shk[t] + ygap[t]=cy[1]/2*ygap[t-2]+cy[1]/2*ygap[t-1]+(.98-cy[1])*ygap[t+1]+cy[2]*(rate[t]-pinf[t+1])+ygap_shk[t] + end + + @reinitialize model + return model +end + +function constructE6!(model::Model) + model.flags.linear = true + model.options.substitutions = false + + @parameters model begin + p_dlp = 0.0050000000000000 + p_dly = 0.0045000000000000 + end + + @variables model begin + dlp; dly; dlyn; lp; ly; lyn + end + + @shocks model begin + dlp_shk; dly_shk + end + + @autoexogenize model begin + ly = dly_shk + lp = dlp_shk + end + + @equations model begin + dly[t]=(1-0.2-0.2)*p_dly+0.2*dly[t-1]+0.2*dly[t+1]+dly_shk[t] + dlp[t]=(1-0.5)*p_dlp+0.1*dlp[t-2]+0.1*dlp[t-1]+0.1*dlp[t+1]+0.1*dlp[t+2]+0.1*dlp[t+3]+dlp_shk[t] + dlyn[t]=dly[t]+dlp[t] + ly[t]=ly[t-1]+dly[t] + lp[t]=lp[t-1]+dlp[t] + lyn[t]=lyn[t-1]+dlyn[t] + end + + @reinitialize model + return model +end + +function constructE7!(model::Model) + model.flags.linear = false + model.options.substitutions = true + + @parameters model begin + delta = 0.1000000000000000 + p_dlc_ss = 0.0040000000000000 + p_dlinv_ss = 0.0040000000000000 + p_growth = 0.0040000000000000 + end + + @variables model begin + dlc; dlinv; dly; lc; linv; + lk; ly; + end + + @shocks model begin + dlc_shk; dlinv_shk; + end + + @autoexogenize model begin + lc = dlc_shk + linv = dlinv_shk + end + + @equations model begin + dlc[t]=(1-0.2-0.2)*p_dlc_ss+0.2*dlc[t-1]+0.2*dlc[t+1]+dlc_shk[t] + dlinv[t]=(1-0.5)*p_dlinv_ss+0.1*dlinv[t-2]+0.1*dlinv[t-1]+0.1*dlinv[t+1]+0.1*dlinv[t+2]+0.1*dlinv[t+3]+dlinv_shk[t] + lc[t]=lc[t-1]+dlc[t] + linv[t]=linv[t-1]+dlinv[t] + ly[t]=log(exp(lc[t])+exp(linv[t])) + dly[t]=ly[t]-ly[t-1] + lk[t]=log((1-delta)*exp(lk[t-1])+exp(linv[t])) + end + + @reinitialize model + return model +end + +function constructE7A!(model::Model) + model.flags.linear = false + model.options.substitutions = false + + @parameters model begin + delta = 0.1000000000000000 + p_dlc_ss = 0.0040000000000000 + p_dlinv_ss = 0.0040000000000000 + p_growth = 0.0040000000000000 + end + + @variables model begin + dlc; dlinv; dly; lc; linv; + lk; ly; + end + + @shocks model begin + dlc_shk; dlinv_shk; + end + + @autoexogenize model begin + lc = dlc_shk + linv = dlinv_shk + end + + @equations model begin + dlc[t] = (1 - 0.2 - 0.2) * p_dlc_ss + 0.2 * dlc[t - 1] + 0.2 * dlc[t + 1] + dlc_shk[t] + dlinv[t] = (1 - 0.5) * p_dlinv_ss + 0.1 * dlinv[t - 2] + 0.1 * dlinv[t - 1] + 0.1 * dlinv[t + 1] + 0.1 * dlinv[t + 2] + 0.1 * dlinv[t + 3] + dlinv_shk[t] + lc[t] = lc[t - 1] + dlc[t] + linv[t] = linv[t - 1] + dlinv[t] + ly[t] = log(exp(lc[t]) + exp(linv[t])) + dly[t] = ly[t] - ly[t - 1] + lk[t] = log((1 - delta) * exp(lk[t - 1]) + exp(linv[t])) + end + + @reinitialize model + + @steadystate model linv = lc - 7; + @steadystate model lc = 14; + + return model +end + +function constructS1!(model::Model) + model.flags.linear = true + model.options.substitutions = false + + @variables model a b c + + @shocks model b_shk c_shk + + @parameters model begin + a_ss = 1.2 + α = 0.5 + β = 0.8 + q = 2 + end + + @equations model begin + a[t] = b[t] + c[t] + b[t] = @sstate(b) * (1 - α) + α * b[t-1] + b_shk[t] + c[t] = q * @sstate(b) * (1 - β) + β * c[t-1] + c_shk[t] + end + + @reinitialize model + + @steadystate model a = a_ss + + return model +end + +function constructS2!(model::Model) + model.flags.linear = false + model.options.substitutions = false + + @parameters model begin + α = 0.5 + x_ss = 3.1 + end + + @variables model begin + y + @log x + # @shock x_shk + end + + @shocks model begin + x_shk + end + + @autoexogenize model begin + x = x_shk + end + + @equations model begin + y[t] = (1 - α) * 2 * @sstate(x) + (α) * @movav(y[t-1], 4) + log(x[t]) = (1 - α) * log(x_ss) + (α) * @movav(log(x[t-1]), 2) + x_shk[t] + end + + @reinitialize model + + return model +end + +function deconstructE1!(model::Model) + @variables model @delete y + + @shocks model @delete y_shk + + @autoexogenize model begin + @delete y => y_shk + end + + @equations model begin + @delete _EQ1 + end +end + +function deconstructE2!(model::Model) + @variables model @delete pinf rate ygap + + @shocks model @delete pinf_shk rate_shk ygap_shk + + @autoexogenize model begin + @delete pinf = pinf_shk + @delete rate = rate_shk + @delete ygap = ygap_shk + end + + @equations model begin + @delete _EQ1 _EQ2 _EQ3 + end +end + +function deconstructE3!(model::Model) + @variables model @delete pinf rate ygap + + @shocks model @delete pinf_shk rate_shk ygap_shk + + @autoexogenize model begin + @delete pinf = pinf_shk + @delete rate = rate_shk + @delete ygap = ygap_shk + end + + @equations model begin + @delete _EQ1 _EQ2 _EQ3 + end +end + +function deconstructE3nl!(model::Model) + @variables model @delete pinf rate ygap + + @shocks model @delete pinf_shk rate_shk ygap_shk + + @autoexogenize model begin + @delete pinf = pinf_shk + @delete rate = rate_shk + @delete ygap = ygap_shk + end + + @equations model begin + @delete _EQ1 _EQ2 _EQ3 + end +end + +function deconstructE6!(model::Model) + @variables model @delete dlp dly dlyn lp ly lyn + + @shocks model @delete dlp_shk dly_shk + + @autoexogenize model begin + @delete ly => dly_shk + @delete lp => dlp_shk + end + + @equations model begin + @delete _EQ1 _EQ2 _EQ3 _EQ4 _EQ5 _EQ6 + end +end + +function deconstructE7!(model::Model) + @variables model @delete dlc dlinv dly lc linv lk ly + + @shocks model @delete dlc_shk dlinv_shk + + @autoexogenize model begin + @delete lc => dlc_shk + @delete linv => dlinv_shk + end + + @equations model begin + @delete _EQ1 _EQ2 _EQ3 _EQ4 _EQ5 _EQ6 _EQ7 + end +end + +function deconstructE7A!(model::Model) + @variables model @delete dlc dlinv dly lc linv lk ly + + @shocks model @delete dlc_shk dlinv_shk + + @autoexogenize model begin + @delete lc => dlc_shk + @delete linv => dlinv_shk + end + + @equations model begin + @delete _EQ1 _EQ2 _EQ3 _EQ4 _EQ5 _EQ6 _EQ7 + end + + @steadystate model begin + @delete _SSEQ1 _SSEQ2 + end +end + +function deconstructS1!(model::Model) + @variables model @delete a b c + + @shocks model @delete b_shk c_shk + + @equations model begin + @delete _EQ1 _EQ2 _EQ3 + end + + @steadystate model begin + @delete _SSEQ1 + end +end + +function deconstructS2!(model::Model) + @variables model @delete y x + + @shocks model @delete x_shk + + @autoexogenize model begin + @delete x => x_shk + end + + @equations model begin + @delete _EQ1 _EQ2 + end +end + + +function getInitialModel(s::Symbol) + s == :E1 && return E1.newmodel() + s == :E2 && return E2.newmodel() + s == :E3 && return E3.newmodel() + s == :E3nl && return E3nl.newmodel() + s == :E6 && return E6.newmodel() + s == :E7 && return E7.newmodel() + s == :E7A && return E7A.newmodel() + s == :S1 && return S1.newmodel() + s == :S2 && return S2.newmodel() +end + +function deconstructedModel() + symbols = (:E1, :E2, :E3, :E3nl, :E6, :E7, :E7A, :S1, :S2) + sym = symbols[rand(1:length(symbols))] + m = getInitialModel(sym) + + sym == :E1 && deconstructE1!(m) + sym == :E2 && deconstructE2!(m) + sym == :E3 && deconstructE3!(m) + sym == :E3nl && deconstructE3nl!(m) + sym == :E6 && deconstructE6!(m) + sym == :E7 && deconstructE7!(m) + sym == :E7A && deconstructE7A!(m) + sym == :S1 && deconstructS1!(m) + sym == :S2 && deconstructS2!(m) + + return m +end + +for i in 1:3 + @info "Model change variations $i" + global getE1() = constructE1!(deconstructedModel()) + global getE2() = constructE2!(deconstructedModel()) + global getE3() = constructE3!(deconstructedModel()) + global getE3nl() = constructE3nl!(deconstructedModel()) + global getE6() = constructE6!(deconstructedModel()) + global getE7() = constructE7!(deconstructedModel()) + global getE7A() = constructE7A!(deconstructedModel()) + global getS1() = constructS1!(deconstructedModel()) + global getS2() = constructS2!(deconstructedModel()) + + include("simtests.jl") + include("sim_fo.jl") + include("dynss.jl") + include("stochsims.jl") + include("sstests.jl") + include("misc.jl") +end + diff --git a/test/runtests.jl b/test/runtests.jl index 1c4278a..2e33988 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -20,6 +20,16 @@ using ModelBaseEcon @using_example S1 @using_example S2 +getE1() = E1.newmodel() +getE2() = E2.newmodel() +getE3() = E3.newmodel() +getE3nl() = E3nl.newmodel() +getE6() = E6.newmodel() +getE7() = E7.newmodel() +getE7A() = E7A.newmodel() +getS1() = S1.newmodel() +getS2() = S2.newmodel() + using Test using Random using Suppressor @@ -50,7 +60,7 @@ import Pardiso end @testset "Plans" begin - m = deepcopy(E1.model) + m = getE1() p = Plan(m, 1:3) @test first(p.range) == 0U @test last(p.range) == 4U @@ -146,7 +156,7 @@ include("simdatatests.jl") include("sstests.jl") @testset "misc" begin - m = deepcopy(E3.model) + m = getE3() sim = m.maxlag .+ (1:10) p = Plan(m, sim) @@ -250,3 +260,5 @@ for sfdef = QuoteNode.(StateSpaceEcon.StackedTimeSolver.sf_libs) include("stochsims.jl") end + +include("modelchanges.jl") \ No newline at end of file diff --git a/test/sim_fo.jl b/test/sim_fo.jl index 7593785..2cfc6c9 100644 --- a/test/sim_fo.jl +++ b/test/sim_fo.jl @@ -85,17 +85,17 @@ end @using_example E2 @testset "E2.fo.unant" begin - run_fo_unant_tests(deepcopy(E2.model)) + run_fo_unant_tests(getE2()) end @using_example E3 @testset "E3.fo.unant" begin - run_fo_unant_tests(deepcopy(E3.model)) + run_fo_unant_tests(getE3()) end @using_example E6 @testset "E6.fo.unant" begin - let m = deepcopy(E6.model) + let m = getE6() # set slopes to 0, otherwise we're not allowed to linearize m.p_dly = 0 m.p_dlp = 0 @@ -199,12 +199,12 @@ end @steadystate m x = 6 test_initdecomp_firstorder(m) end - for m in (E2.model, E3.model) + for m in (getE2(), getE3()) m = deepcopy(m) empty!(m.sstate.constraints) test_initdecomp_firstorder(m) end - let m = deepcopy(E6.model) + let m = getE6() # set slopes to 0, otherwise we're not allowed to linearize m.p_dly = 0 m.p_dlp = 0 @@ -326,12 +326,12 @@ end @steadystate m x = 6 test_initdecomp_stackedtime(m, fctype=fclevel) end - for m in (E2.model, E3.model) + for m in (getE2(), getE3()) m = deepcopy(m) empty!(m.sstate.constraints) test_initdecomp_stackedtime(m) end - let m = deepcopy(E6.model) + let m = getE6() # set slopes to 0, otherwise we're not allowed to linearize m.p_dly = 0 m.p_dlp = 0 diff --git a/test/simtests.jl b/test/simtests.jl index 3c240ad..dfbc313 100644 --- a/test/simtests.jl +++ b/test/simtests.jl @@ -9,7 +9,7 @@ using LinearAlgebra using DelimitedFiles @testset "E1.simple" begin - m = deepcopy(E1.model) + m = getE1() m.α = 0.5 m.β = 1.0 - m.α for T = 6:10 @@ -96,26 +96,26 @@ function test_simulation(m_in, path; atol=1.0e-9) end @testset "E1.sim" begin - test_simulation(E1.model, "data/M1_TestData.csv") + test_simulation(getE1(), "data/M1_TestData.csv") end @testset "E2.sim" begin - test_simulation(E2.model, "data/M2_TestData.csv") + test_simulation(getE2(), "data/M2_TestData.csv") end @testset "E3.sim" begin - test_simulation(E3.model, "data/M3_TestData.csv") + test_simulation(getE3(), "data/M3_TestData.csv") end @testset "E6.sim" begin - test_simulation(E6.model, "data/M6_TestData.csv") + test_simulation(getE6(), "data/M6_TestData.csv") end ############################################################# # linearization tests @testset "linearize" begin - m3 = deepcopy(E3.model) + m3 = getE3() clear_sstate!(m3) @test_throws ModelBaseEcon.LinearizationError linearize!(m3) @@ -124,7 +124,7 @@ end test_simulation(m3, "data/M3_TestData.csv") - m3nl = deepcopy(E3nl.model) + m3nl = getE3nl() m3nl.sstate.values .= m3.sstate.values m3nl.sstate.mask .= m3.sstate.mask @test issssolved(m3nl) @@ -141,7 +141,7 @@ end @test isa(ModelBaseEcon.getevaldata(m3), ModelBaseEcon.LinearizedModelEvaluationData) test_simulation(m3, "data/M3_TestData.csv") - m7 = deepcopy(E7.model) + m7 = getE7() if !issssolved(m7) empty!(m7.sstate.constraints) @steadystate m7 linv = lc - 7 @@ -152,7 +152,7 @@ end # the aux variables are present @test_throws ModelBaseEcon.LinearizationError linearize!(m7) - m7a = deepcopy(E7A.model) + m7a = getE7A() clear_sstate!(m7a) sssolve!(m7a) # non-zeros linear growth in the steady state @@ -165,7 +165,7 @@ end # Tests of unanticipated shocks @testset "E1.unant" begin - m = deepcopy(E1.model) + m = getE1() m.α = m.β = 0.5 p = Plan(m, 1:3) data = zeroarray(m, p) @@ -178,7 +178,7 @@ end end @testset "E2.unant" begin - m = deepcopy(E2.model) + m = getE2() nvars = length(m.variables) nshks = length(m.shocks) var_inds = 1:nvars @@ -251,7 +251,7 @@ end end @testset "E3.unant" begin - m = deepcopy(E3.model) + m = getE3() nvars = length(m.variables) nshks = length(m.shocks) var_inds = 1:nvars @@ -354,7 +354,7 @@ end ed = zerodata(m, p) @test_throws ArgumentError simulate(m, p, ed, fctype=fcnatural) end - let m = deepcopy(E1.model) + let m = getE1() p = Plan(m, 3:17) ed = zerodata(m, p) ed .= rand(Float64, size(ed)) @@ -367,7 +367,7 @@ end R2 = StateSpaceEcon.StackedTimeSolver.stackedtime_R!(similar(R1), x, ed, sd) @test R1 ≈ R2 end - let m = deepcopy(E3.model) + let m = getE3() p = Plan(m, 3:177) ed = zerodata(m, p) ed[3U, m.shocks] .= 0.2 * rand(1, 3) @@ -393,7 +393,7 @@ end R2 = StateSpaceEcon.StackedTimeSolver.stackedtime_R!(similar(R1), x, ed, sd) @test R1 ≈ R2 end - let m = deepcopy(E6.model) + let m = getE6() empty!(m.sstate.constraints) @steadystate m lp = 1.0 @steadystate m lp + lyn = 1.5 @@ -416,7 +416,7 @@ end end @testset "new.unant" begin - m = deepcopy(E1.model) + m = getE1() m.α = m.β = 0.5 test_rng = 20Q1:22Q1 @@ -476,7 +476,7 @@ end @testset "linesearch, deviation" begin # linesearch should get the same results - let m = deepcopy(E3nl.model) + let m = getE3nl() clear_sstate!(m) nvars = ModelBaseEcon.nvariables(m) nshks = ModelBaseEcon.nshocks(m) @@ -508,7 +508,7 @@ end @test isapprox(res01, res01_line; atol=1.0e-9) end - let m = deepcopy(E1.model) + let m = getE1() clear_sstate!(m) nvars = ModelBaseEcon.nvariables(m) nshks = ModelBaseEcon.nshocks(m) @@ -531,7 +531,7 @@ end @test isapprox(res01_dev, res01; atol=1.0e-9) end - let m = deepcopy(E3.model) + let m = getE3() clear_sstate!(m) nvars = ModelBaseEcon.nvariables(m) nshks = ModelBaseEcon.nshocks(m) @@ -573,93 +573,3 @@ end end end -############################################################# -# Tests of equation changes -@testset "E1 to E2.sim" begin - m = E1.newmodel() - @parameters m begin - cp = [0.5, 0.02] - cr = [0.75, 1.5, 0.5] - cy = [0.5, -0.02] - end - @variables m begin - @delete y - pinf - rate - ygap - end - @shocks m begin - @delete y_shk - pinf_shk - rate_shk - ygap_shk - end - @autoexogenize m begin - @delete y => y_shk - pinf = pinf_shk - rate = rate_shk - ygap = ygap_shk - end - @equations m begin - @delete _EQ1; - pinf[t]=cp[1]*pinf[t-1]+(.98-cp[1])*pinf[t+1]+cp[2]*ygap[t]+pinf_shk[t] - rate[t]=cr[1]*rate[t-1]+(1-cr[1])*(cr[2]*pinf[t]+cr[3]*ygap[t])+rate_shk[t] - ygap[t]=cy[1]*ygap[t-1]+(.98-cy[1])*ygap[t+1]+cy[2]*(rate[t]-pinf[t+1])+ygap_shk[t] - end - - @reinitialize m - - test_simulation(m, "data/M2_TestData.csv") -end - -@testset "E2 to E3.sim" begin - m = E2.newmodel() - @equations m begin - :_EQ1 => pinf[t]=cp[1]*pinf[t-1]+0.3*pinf[t+1]+0.05*pinf[t+2]+0.05*pinf[t+3]+cp[2]*ygap[t]+pinf_shk[t] - :_EQ2 => rate[t]=cr[1]*rate[t-1]+(1-cr[1])*(cr[2]*pinf[t]+cr[3]*ygap[t])+rate_shk[t] - :_EQ3 => ygap[t]=cy[1]/2*ygap[t-2]+cy[1]/2*ygap[t-1]+(.98-cy[1])*ygap[t+1]+cy[2]*(rate[t]-pinf[t+1])+ygap_shk[t] - end - - @reinitialize m - - test_simulation(m, "data/M3_TestData.csv") -end - -@testset "E3 to E6.sim" begin - m = E3.newmodel() - @parameters m begin - p_dlp = 0.0050000000000000 - p_dly = 0.0045000000000000 - end - @variables m begin - @delete pinf rate ygap - dlp; dly; dlyn; lp; ly; lyn - end - @shocks m begin - @delete pinf_shk rate_shk ygap_shk - dlp_shk; dly_shk - end - @autoexogenize m begin - @delete pinf => pinf_shk - @delete rate => rate_shk - @delete ygap => ygap_shk - ly = dly_shk - lp = dlp_shk - end - - @equations m begin - @delete _EQ1; - @delete _EQ2; - @delete _EQ3; - dly[t]=(1-0.2-0.2)*p_dly+0.2*dly[t-1]+0.2*dly[t+1]+dly_shk[t] - dlp[t]=(1-0.5)*p_dlp+0.1*dlp[t-2]+0.1*dlp[t-1]+0.1*dlp[t+1]+0.1*dlp[t+2]+0.1*dlp[t+3]+dlp_shk[t] - dlyn[t]=dly[t]+dlp[t] - ly[t]=ly[t-1]+dly[t] - lp[t]=lp[t-1]+dlp[t] - lyn[t]=lyn[t-1]+dlyn[t] - end - - @reinitialize m - - test_simulation(m, "data/M6_TestData.csv") -end diff --git a/test/sstests.jl b/test/sstests.jl index c93b622..9147d2c 100644 --- a/test/sstests.jl +++ b/test/sstests.jl @@ -7,7 +7,7 @@ @testset "E1.sstate" begin - let m = deepcopy(E1.model) + let m = getE1() empty!(m.sstate.constraints) m.α = 0.5 m.β = 1.0 - m.α @@ -53,7 +53,7 @@ end @testset "E1.sstate, auto" begin - let m = deepcopy(E1.model) + let m = getE1() empty!(m.sstate.constraints) m.α = 0.5 m.β = 1.0 - m.α @@ -93,7 +93,7 @@ end end @testset "E1.sstate, lm" begin - let m = deepcopy(E1.model) + let m = getE1() empty!(m.sstate.constraints) m.α = 0.5 m.β = 1.0 - m.α @@ -133,7 +133,7 @@ end end @testset "E2.sstate" begin - let m = deepcopy(E2.model) + let m = getE2() empty!(m.sstate.constraints) clear_sstate!(m) sssolve!(m) @@ -159,7 +159,7 @@ end end @testset "E6.sstate" begin - let m = deepcopy(E6.model) + let m = getE6() m.options.maxiter = 50 clear_sstate!(m) @@ -181,7 +181,7 @@ end end @testset "E7.sstate" begin - let m = deepcopy(E7.model) + let m = getE7() m.options.maxiter = 100 m.options.tol = 1e-9 @@ -319,7 +319,7 @@ end end @testset "presolve, ssZeroSlope" begin - let m = deepcopy(E1.model) + let m = getE1() empty!(m.sstate.constraints) m.α = 0.5 m.β = 1.0 - m.α @@ -329,7 +329,7 @@ end @test m.sstate.mask == [false, false, true, true] end - let m = deepcopy(E1.model) + let m = getE1() empty!(m.sstate.constraints) m.α = 0.5 m.β = 1.0 - m.α diff --git a/test/stochsims.jl b/test/stochsims.jl index 743549a..0559883 100644 --- a/test/stochsims.jl +++ b/test/stochsims.jl @@ -8,7 +8,7 @@ @using_example E6 @testset "stochsims" begin - m = deepcopy(E6.model) + m = getE6() m.options.fctype = fcnatural @steadystate m lp = 0.6 @steadystate m ly = 0.3 From c087abccefa60a8b0ed10a5d053f7836b45ff0b0 Mon Sep 17 00:00:00 2001 From: Boyan Bejanov Date: Tue, 4 Jul 2023 09:12:38 -0400 Subject: [PATCH 30/32] no need to use dicts in _presolve - internal function --- src/stackedtime/solverdata.jl | 4 ++-- src/steadystate/presolve.jl | 39 +++++++++++++++++------------------ 2 files changed, 21 insertions(+), 22 deletions(-) diff --git a/src/stackedtime/solverdata.jl b/src/stackedtime/solverdata.jl index 12e6c38..9e94267 100644 --- a/src/stackedtime/solverdata.jl +++ b/src/stackedtime/solverdata.jl @@ -368,8 +368,8 @@ StackedTimeSolverData(model::Model, plan::Plan, fctype::FinalCondition, variant: # Prep the Jacobian matrix neq = 0 # running counter of equations added to matrix # Model equations are the same for each sim period, just shifted according to t - Jblock = [ti + NT * (var_to_idx[var] - 1) for eqn_pair in equations for (var, ti) in keys(eqn_pair[2].tsrefs)] - Iblock = [i for (i, eqn_pair) in enumerate(equations) for _ in eqn_pair[2].tsrefs] + Jblock = [ti + NT * (var_to_idx[var] - 1) for (_, eqn) in equations for (var, ti) in keys(eqn.tsrefs)] + Iblock = [i for (i, (_, eqn)) in enumerate(equations) for _ in eqn.tsrefs] Tblock = -model.maxlag:model.maxlead for t in sim diff --git a/src/steadystate/presolve.jl b/src/steadystate/presolve.jl index 604f4ed..7881864 100644 --- a/src/steadystate/presolve.jl +++ b/src/steadystate/presolve.jl @@ -1,7 +1,7 @@ ################################################################################## # This file is part of StateSpaceEcon.jl # BSD 3-Clause License -# Copyright (c) 2020-2022, Bank of Canada +# Copyright (c) 2020-2023, Bank of Canada # All rights reserved. ################################################################################## @@ -20,13 +20,13 @@ loop. ### Arguments * `model` - the Model instance * `mask` - a vector of Bool. Defaults to `model.sstate.mask` - * `vals` - a vector of numbers. Defaults to `model.sstate.values` Caller - must specify either both `mask` and `vals` or neither of them. `mask[i]` + * `values` - a vector of numbers. Defaults to `model.sstate.values` Caller + must specify either both `mask` and `values` or neither of them. `mask[i]` equals `true` if and only if the i-th steady state value has alredy been solved for. -`mask` and `vals` are both input and output data to the algorithm. Any vals -that are successfully presolved are updated in place, and their `mask` enties +`mask` and `values` are both input and output data to the algorithm. Any values +that are successfully pre-solved are updated in place, and their `mask` entries are set to `true`. ### Options @@ -52,41 +52,40 @@ function solve1d(sseqn, vals, ind, ::Val{:newton}, tol = 1e-12, maxiter = 5) return newton1!(sseqn.eval_RJ, vals, ind; tol, maxiter) end -function _presolve_equations!(eqns, mask, vals, method, verbose, tol) - eqns_solved = OrderedDict{Symbol, Bool}(key => false for key in keys(eqns)) - eqns_resid = OrderedDict{Symbol, Float64}(key => 0.0 for key in keys(eqns)) - # eqns_resid = zeros(length(eqns)) +function _presolve_equations!(eqns, mask, values, method, verbose, tol) + eqns_solved = falses(length(eqns)) + eqns_resid = zeros(length(eqns)) retval = false # return value: true if any mask changed, false otherwise updated = true # keep track if any mask changed within the loop - while updated && !all(values(eqns_solved)) + while updated && !all(eqns_solved) updated = false # keep track if anything changes this outer iteration - for (eqn_key, sseqn) in eqns - eqns_solved[eqn_key] && continue + for (eqn_idx, (eqn_key, sseqn)) in enumerate(eqns) + eqns_solved[eqn_idx] && continue # mask is true for variables that are already solved unsolved = .!mask[sseqn.vinds] nunsolved = sum(unsolved) if nunsolved == 0 # all variables are solved, yet equation is not marked solved. # check if equation is satisfied - eqns_resid[eqn_key] = R = sseqn.eval_resid(vals[sseqn.vinds]) + eqns_resid[eqn_idx] = R = sseqn.eval_resid(values[sseqn.vinds]) # mark it solved either way, but issue a warning if residual is not zero - eqns_solved[eqn_key] = true + eqns_solved[eqn_idx] = true if verbose && abs(R) > 100tol @warn "Equation $eqn_key has residual $R:\n $sseqn" end elseif nunsolved == 1 # only one variable left unsolved. call the 1d solver ind = findall(unsolved)[1] - _vals = vals[sseqn.vinds] - success = solve1d(sseqn, _vals, ind, Val(method), 0.1tol) + vals = values[sseqn.vinds] + success = solve1d(sseqn, vals, ind, Val(method), 0.1tol) if success - if abs(_vals[ind]) < 1e-10 - _vals[ind] = 0.0 + if abs(vals[ind]) < 1e-10 + vals[ind] = 0.0 end - vals[sseqn.vinds[ind]] = _vals[ind] + values[sseqn.vinds[ind]] = vals[ind] retval = updated = mask[sseqn.vinds[ind]] = true if verbose - @info "Presolved equation $eqn_key for $(sseqn.vsyms[ind]) = $(_vals[ind])\n $sseqn" + @info "Presolved equation $eqn_key for $(sseqn.vsyms[ind]) = $(vals[ind])\n $sseqn" end else if verbose From adc52d379ce62681e0fa40e1262e666ffe0752b9 Mon Sep 17 00:00:00 2001 From: Boyan Bejanov Date: Tue, 4 Jul 2023 09:54:53 -0400 Subject: [PATCH 31/32] ignoremissing=true when comparing data with different ranges --- test/sim_fo.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/sim_fo.jl b/test/sim_fo.jl index 2cfc6c9..28e83cc 100644 --- a/test/sim_fo.jl +++ b/test/sim_fo.jl @@ -258,7 +258,7 @@ function test_initdecomp_stackedtime(m; nonlin=!m.linear, rng=2001Q1:2024Q4, fct # without initdecomp the resulting range is shorter # and the numbers will be identical for linear models @test rangeof(ref) ≠ rangeof(res1a) - @test compare(ref, res1a; atol, quiet=true) == !nonlin + @test compare(ref, res1a; atol, quiet=true, ignoremissing=true) == !nonlin #test 2 - only 1 shock at a time for shk in shks @@ -291,7 +291,7 @@ function test_initdecomp_stackedtime(m; nonlin=!m.linear, rng=2001Q1:2024Q4, fct # without initdecomp the resulting range is shorter # and the numbers are totally different @test rangeof(ref) ≠ rangeof(res1a) - @test compare(ref, res1a; atol, quiet=true, nans=true) == (m.maxlag == 0) + @test compare(ref, res1a; atol, quiet=true, nans=true, ignoremissing=true) == (m.maxlag == 0) end # test 3 - all shocks and init From de30228c1c677a4856cb2bfbbf5a971cf67437a5 Mon Sep 17 00:00:00 2001 From: Boyan Bejanov Date: Thu, 13 Jul 2023 10:05:24 -0400 Subject: [PATCH 32/32] version = "0.5" + update compat section --- Project.toml | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Project.toml b/Project.toml index 84bcdcc..ada71fa 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "StateSpaceEcon" uuid = "e4c825b0-b65c-11ea-0b5a-6176b64e7b7f" authors = ["Atai Akunov ", "Boyan Bejanov ", "Nicholas Labelle St-Pierre "] -version = "0.4.2" +version = "0.5.0" [deps] DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab" @@ -23,11 +23,11 @@ TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" DiffResults = "1.0" ForwardDiff = "0.10" JLD2 = "0.4" -ModelBaseEcon = "0.5.2" +ModelBaseEcon = "0.6" Pardiso = "0.5" Suppressor = "0.2" TimerOutputs = "0.5" -TimeSeriesEcon = "0.5.1" +TimeSeriesEcon = "0.6" julia = "1.7" [extras]