From 755addb9977c287c3a4eaf6895f81f6d512a2a6a Mon Sep 17 00:00:00 2001 From: MichelJuillard Date: Wed, 2 Oct 2024 09:38:29 +0200 Subject: [PATCH] switch to NonlinearSolver --- src/Dynare.jl | 1 + .../conditional_simulation.jl | 83 ++++------ src/perfectforesight/perfectforesight.jl | 149 +++++++----------- src/steady_state/SteadyState.jl | 148 ++++++++++------- test/models/initialization/neoclassical1.mod | 7 +- test/runtests.jl | 2 + test/test_scenario.jl | 27 ++-- 7 files changed, 193 insertions(+), 224 deletions(-) diff --git a/src/Dynare.jl b/src/Dynare.jl index 60011e00..9a073930 100644 --- a/src/Dynare.jl +++ b/src/Dynare.jl @@ -36,6 +36,7 @@ include("DynareParser.jl") export parser include("DynarePreprocessor.jl") export dynare_preprocess +include("dynare_nonlinear_solvers.jl") include("steady_state/SteadyState.jl") export steadystate! include("dynare_table.jl") diff --git a/src/perfectforesight/conditional_simulation.jl b/src/perfectforesight/conditional_simulation.jl index 6f28ba1c..1f361631 100644 --- a/src/perfectforesight/conditional_simulation.jl +++ b/src/perfectforesight/conditional_simulation.jl @@ -4,7 +4,6 @@ function set_future_information!(y::Vector{Float64}, x::Vector{Float64}, context exogenous_nbr = m.exogenous_nbr flips = context.work.scenario[infoperiod] symboltable = context.symboltable - for (period, f1) in flips p = length(infoperiod:period) for (s, f2) in f1 @@ -27,12 +26,14 @@ function perfectforesight_core_conditional!( y0::AbstractVector{Float64}, initialvalues::Vector{<:Real}, terminalvalues::Vector{<:Real}, - solve_algo, - linear_solve_algo::LinearSolveAlgo, dynamic_ws::DynamicWs, flipinfo::FlipInformation, infoperiod; + linear_solve_algo = nothing, maxit = 50, + method = TrustRegion(), + nonlinear_solve_algo = "NonlinearSolve", + show_trace = true, tolf = 1e-5, tolx = 1e-5 @@ -74,7 +75,7 @@ function perfectforesight_core_conditional!( flipinfo.ix_stack ) - J! = make_pf_jacobian( + j! = make_pf_jacobian( DFunctions.dynamic_derivatives!, initialvalues, terminalvalues, @@ -92,62 +93,31 @@ function perfectforesight_core_conditional!( flipinfo, nzval1 ) - - df = OnceDifferentiable(f!, J!, y0, residuals, JJ) + @debug "$(now()): start nlsolve" set_future_information!(y0, exogenous, context, periods, infoperiod) flip!(y0, exogenous, flipinfo.ix_stack) - f!(residuals, y0) - show_trace = (("JULIA_DEBUG" => "Dynare") in ENV) ? true : false - if linear_solve_algo == pardiso - @show "Pardiso" - if isnothing(Base.get_extension(Dynare, :PardisoSolver)) - error("You must load Pardiso with 'using MKL, Pardiso'") - end - ls1!(x, A, b) = linear_solver!(PardisoLS(), x, A, b) - res = nlsolve(df, y0, method = :robust_trust_region, show_trace = show_trace, ftol = tolf, xtol = tolx, iterations= maxit, linsolve = ls1!) - else - ls2!(x, A, b) = linear_solver!(IluLS(), x, A, b) - res = nlsolve(df, y0, method = :robust_trust_region, show_trace = show_trace, ftol = tolf, xtol = tolx, iterations= maxit, linsolve = ls2!) - end - fj = NonlinearFunction(f!, jac = j!, jac_prototype = JJ) - - prob = NonlinearProblem(fj, y0, params) - @debug "$(now()): start nlsolve" - show_trace = (("JULIA_DEBUG" => "Dynare") in ENV) ? true : false - #= - if linear_solve_algo == pardiso - @show "Pardiso" - if isnothing(Base.get_extension(Dynare, :PardisoSolver)) - error("You must load Pardiso with 'using MKL, Pardiso'") + @debug "$(now()): end nlsolve" + let res + try + residuals = similar(y0) + f!(residuals, y0, work.params) + res = dynare_nonlinear_solvers(f!, j!, JJ, [], [], y0, exogenous, work.params, + linear_solve_algo = linear_solve_algo, + maxit = maxit, + method = method, + solve_algo = nonlinear_solve_algo, + show_trace = show_trace, + tolf = tolf, + tolx = tolx) + catch e + @show e end - ls1!(x, A, b) = linear_solver!(PardisoLS(), x, A, b) - res = nlsolve(df, y0, method = :robust_trust_region, show_trace = show_trace, ftol = tolf, xtol = tolx, iterations= maxit, linsolve = ls1!) - #= - elseif linear_solve_algo == mumps - @show "MUMPS" - #= - if isnothing(Base.get_extension(Dynare, :PardisoSolver)) - error("You must load Pardiso with 'using MKL, Pardiso'") - en=# - MPI.Init() - ls!(x, A, b) = linear_solver!(MumpsLS(), x, A, b) - res = nlsolve(df, y0, method = :robust_trust_region, show_trace = show_trace, ftol = tolf, xtol = tolx, iterations= maxit, linsolve = ls!) - =# - else -# ls2!(x, A, b) = linear_solver!(IluLS(), x, A, b) - # res = nlsolve(df, y0, method = :robust_trust_region, show_trace = show_trace, ftol = tolf, xtol = tolx, iterations= maxit, linsolve = ls2!) - # res = NonlinearSolve.solve(prob, solve_algo, show_trace=Val(false), abstol = tolf) - end - =# - @show solve_algo - res = NonlinearSolve.solve(prob, solve_algo, show_trace = Val(show_trace), abstol = tolf) - print_nlsolver_results(res) - @debug "$(now()): end nlsolve" - flip!(res.zero, exogenous, flipinfo.ix_stack) - make_simulation_results!(context::Context, res.zero, exogenous, terminalvalues, periods) + flip!(res, exogenous, flipinfo.ix_stack) + make_simulation_results!(context::Context, res, exogenous, terminalvalues, periods) + end end function make_pf_residuals( @@ -163,7 +133,7 @@ function make_pf_residuals( permutations::Vector{Tuple{Int64,Int64}}, ix_stack::SparseVector{Int, Int} ) where T <: Real - function f!(residuals::AbstractVector{T}, y::AbstractVector{T}) + function f!(residuals, y, params) flip!(y, exogenous, ix_stack) get_residuals!( residuals, @@ -185,7 +155,7 @@ function make_pf_residuals( return f! end -function flip!(y, x, ix_stack::SparseVector{Int, Int}) +function flip!(y, x, ix_stack) I, J = findnz(ix_stack) for (iy, ix) in zip(I, J) y[iy], x[ix] = x[ix], y[iy] @@ -213,6 +183,7 @@ function make_pf_jacobian( function J!( A::SparseArrays.SparseMatrixCSC{Float64,Int64}, y::AbstractVector{Float64}, + params ) updateJacobian!( A, diff --git a/src/perfectforesight/perfectforesight.jl b/src/perfectforesight/perfectforesight.jl index 7c36d1e8..b4595fc7 100644 --- a/src/perfectforesight/perfectforesight.jl +++ b/src/perfectforesight/perfectforesight.jl @@ -1,13 +1,20 @@ include("makeA.jl") @enum PerfectForesightAlgo trustregionA -#@enum LinearSolveAlgo ilu pardiso mumps @enum LinearSolveAlgo ilu pardiso @enum InitializationAlgo initvalfile steadystate firstorder linearinterpolation +struct DynarePerfectForesightSimulationFailed <: Exception end +function Base.showerror(io::IO, e::DynarePerfectForesightSimulationFailed) + print(""" + Dynare couldn't compute the perfectf foresight simulation. + Either there is no solution or the guess values + are too far from the solution + """) +end + abstract type LinearSolver end struct IluLS <: LinearSolver end -#struct MumpsLS <: LinearSolver end struct PardisoLS <: LinearSolver end function linear_solver!(::IluLS, @@ -17,17 +24,6 @@ function linear_solver!(::IluLS, x .= A\b end -#= -using MUMPS, MPI -function linear_solver!(::MumpsLS, - x::AbstractVector{Float64}, - A::AbstractMatrix{Float64}, - b::AbstractVector{Float64}) - x = MUMPS.solve(A, b) - return x -end -=# - abstract type NonLinearSolver end struct PathNLS <: NonLinearSolver end struct DefaultNLS <: NonLinearSolver end @@ -38,11 +34,13 @@ mutable struct PerfectForesightOptions display::Bool homotopy::Bool initialization_algo::InitializationAlgo - solve_algo + solve_algo::String linear_solve_algo::LinearSolveAlgo maxit::Int mcp::Bool + method periods::Int + show_trace::Bool tolf::Float64 tolx::Float64 end @@ -53,18 +51,18 @@ function PerfectForesightOptions(context::Context, field::Dict{String,Any}) display = true homotopy = false initialization_algo = steadystate - solve_algo = NewtonRaphson(linesearch = LineSearchesJL(method = NonlinearSolve.BackTracking())) + solve_algo = "NonlinearSolve" linear_solve_algo = ilu maxit = 50 mcp = false + method = TrustRegion(radius_update_scheme=RadiusUpdateSchemes.NLsolve) periods = context.work.perfect_foresight_setup["periods"] + show_trace = false tolf = 1e-5 tolx = 1e-5 if haskey(field, "options") for (k, v) in pairs(field["options"]) - if k == "stack_solve_algo" - algo = v::Int64 - elseif k == "noprint" + if k == "noprint" display = false elseif k == "print" display = true @@ -74,8 +72,6 @@ function PerfectForesightOptions(context::Context, field::Dict{String,Any}) mcp = true elseif k == "no_homotopy" homotopy = false - elseif k == "solve_algo" - linear_solve_algo = v elseif k == "maxit" maxit = v elseif k == "tolf" @@ -88,7 +84,7 @@ function PerfectForesightOptions(context::Context, field::Dict{String,Any}) return PerfectForesightOptions(algo, datafile, display, homotopy, initialization_algo, solve_algo, linear_solve_algo, maxit, mcp, - periods, tolf, tolx) + method, periods, show_trace, tolf, tolx) end struct FlipInformation @@ -272,13 +268,15 @@ end - `periods::Int`: number of periods in the simulation [required] - `context::Context=context`: context in which the simulation is computed - `display::Bool=true`: whether to display the results - - `solve_algo = NewtonRaphson(linesearch = LineSearchesJL(method = NonlinearSolve.BackTracking()))`: nonlinear solver algorithm + - `solve_algo = one of "NLsolve", "NonlinearSolve", "PATHSolver" - `linear_solve_algo::LinearSolveAlgo=ilu`: algorithm used for the solution of the linear problem. Either `ilu` or `pardiso`. `ilu` is the sparse linear solver used by default in Julia. To use the Pardiso solver, write `using Pardiso` before running Dynare. - `maxit::Int=50` maximum number of iterations + - `method`: solver method inside `solve_algo` package - `mcp::Bool=false`L whether to solve a mixed complementarity problem with occasionally binding constraints + - `show_trace::Bool=false`: whether to show solver iterations - `tolf::Float64=1e-5`: tolerance for the norm of residualts - `tolx::Float64=1e-5`: tolerance for the norm of the change in the result """ @@ -288,20 +286,20 @@ function perfect_foresight!(;context::Context = context, display::Bool = true, homotopy::Bool = false, initialization_algo::InitializationAlgo = steadystate, - solve_algo = NewtonRaphson(linesearch = LineSearchesJL(method = NonlinearSolve.BackTracking())), + solve_algo = "NonlinearSolve", linear_solve_algo::LinearSolveAlgo = ilu, maxit::Int = 50, mcp::Bool = false, + method = TrustRegion(radius_update_scheme=RadiusUpdateSchemes.NLsolve), periods::Int, + show_trace::Bool = false, tolf::Float64 = 1e-5, tolx::Float64 = 1e-5) options = PerfectForesightOptions(algo, datafile, display, homotopy, initialization_algo, - solve_algo, - linear_solve_algo, maxit, - mcp, periods, tolf, tolx) - + solve_algo, linear_solve_algo, maxit, + mcp, method, periods, show_trace, tolf, tolx) scenario = context.work.scenario check_scenario(scenario) if isempty(scenario) || length(scenario) == 1 @@ -423,11 +421,14 @@ function _perfect_foresight!(context::Context, options::PerfectForesightOptions) guess_values, initial_values, terminal_values, - options.linear_solve_algo, dynamic_ws, perfect_foresight_ws.flipinfo, 1, + linear_solve_algo = nothing, maxit = options.maxit, + method = options.method, + nonlinear_solve_algo = options.solve_algo, + show_trace = options.show_trace, tolf = options.tolf, tolx = options.tolx ) @@ -439,10 +440,12 @@ function _perfect_foresight!(context::Context, options::PerfectForesightOptions) guess_values, initial_values, terminal_values, - options.solve_algo, - options.linear_solve_algo, dynamic_ws, + linear_solve_algo = options.linear_solve_algo, maxit = options.maxit, + method = options.method, + nonlinear_solve_algo = "NonlinearSolve", + show_trace = false, tolf = options.tolf, tolx = options.tolx ) @@ -597,10 +600,12 @@ function perfectforesight_core!( y0::AbstractVector{Float64}, initialvalues::Vector{<:Real}, terminalvalues::Vector{<:Real}, - solve_algo, - linear_solve_algo::LinearSolveAlgo, dynamic_ws::DynamicWs; maxit = 50, + nonlinear_solve_algo = "NonlinearSolve", + linear_solve_algo = nothing, + method = TrustRegion(), + show_trace = true, tolf = 1e-5, tolx = 1e-5 ) @@ -662,62 +667,17 @@ function perfectforesight_core!( j!(JJt, y, p) JJ .= transpose(JJt) end - - f!(residuals, y0, params) - j!(JJ, y0, params) - serialize("test.jls", JJ) - JJt = copy(transpose(JJ)) - jj!(JJt, y0, params) - lp = LinearProblem(JJ, residuals) - lp3 = LinearProblem(JJt, residuals) - dy0 = JJ\residuals - @show residuals[1:20] - @show lp - dy1 = LinearSolve.solve(lp) - @show lp - @show residuals[1:20] - dy2 = LinearSolve.solve(lp, PardisoJL()) - @show lp - @show residuals[1:20] - dy3 = LinearSolve.solve(lp3, PardisoJL()) - @show norm(dy1-dy2) - @show norm(dy1-dy3) - @show norm(dy1-dy0) - - # df = OnceDifferentiable(f!, J!, y0, residuals, JJ) - fj = NonlinearFunction(f!, jac = j!, jac_prototype = JJ) - - prob = NonlinearProblem(fj, y0, params) - @debug "$(now()): start nlsolve" - show_trace = (("JULIA_DEBUG" => "Dynare") in ENV) ? true : false - if linear_solve_algo == pardiso - @show "Pardiso" - if isnothing(Base.get_extension(Dynare, :PardisoSolver)) - error("You must load Pardiso with 'using MKL, Pardiso'") - end - ls1!(x, A, b) = linear_solver!(PardisoLS(), x, A, b) - res = nlsolve(df, y0, method = :robust_trust_region, show_trace = show_trace, ftol = tolf, xtol = tolx, iterations= maxit, linsolve = ls1!) - #= - elseif linear_solve_algo == mumps - @show "MUMPS" - #= - if isnothing(Base.get_extension(Dynare, :PardisoSolver)) - error("You must load Pardiso with 'using MKL, Pardiso'") - en=# - MPI.Init() - ls!(x, A, b) = linear_solver!(MumpsLS(), x, A, b) - res = nlsolve(df, y0, method = :robust_trust_region, show_trace = show_trace, ftol = tolf, xtol = tolx, iterations= maxit, linsolve = ls!) - =# - else -# ls2!(x, A, b) = linear_solver!(IluLS(), x, A, b) - # res = nlsolve(df, y0, method = :robust_trust_region, show_trace = show_trace, ftol = tolf, xtol = tolx, iterations= maxit, linsolve = ls2!) - # res = NonlinearSolve.solve(prob, solve_algo, show_trace=Val(false), abstol = tolf) - end - @show solve_algo - res = NonlinearSolve.solve(prob, solve_algo, show_trace = Val(show_trace), abstol = tolf) - print_nlsolver_results(res) + + results = dynare_nonlinear_solvers(f!, j!, JJ, [], [], y0, exogenous, work.params, + linear_solve_algo = linear_solve_algo, + maxit = maxit, + method = method, + solve_algo = nonlinear_solve_algo, + show_trace = show_trace, + tolf = tolf, + tolx = tolx) @debug "$(now()): end nlsolve" - make_simulation_results!(context::Context, res.u, exogenous, terminalvalues, periods) + make_simulation_results!(context::Context, results, exogenous, terminalvalues, periods) end function make_pf_residuals( @@ -732,7 +692,7 @@ function make_pf_residuals( temp_vec::AbstractVector{T}, permutations::Vector{Tuple{Int64,Int64}}, ) where T <: Real - function f!(residuals::AbstractVector{T}, y::AbstractVector{T}, params) + function f!(residuals::AbstractVector{T}, y::AbstractVector{T}, params::AbstractVector{T}) get_residuals!( residuals, vec(y), @@ -811,8 +771,9 @@ function get_residuals!( temp_vec::AbstractVector{Float64}; permutations::Vector{Tuple{Int64,Int64}} = Tuple{Int64,Int64}[], ) + n = m.endogenous_nbr - + rx = 1:m.exogenous_nbr @views get_residuals_1!( residuals, @@ -972,10 +933,7 @@ function print_nlsolver_results(r) @info @sprintf "Results of Nonlinear Solver Algorithm" @info @sprintf " * Algorithm: %s" r.alg @info @sprintf " * Inf-norm of residuals: %f" norm(r.resid) -# @info @sprintf " * Iterations: %d" r.iterations @info @sprintf " * Convergence: %s" r.retcode -# @info @sprintf " * |x - x'| < %.1e: %s" r.xtol r.x_converged -# @info @sprintf " * |f(x)| < %.1e: %s" r.ftol r.f_converged @info @sprintf " * Function Calls (f): %d" r.stats.nf return end @@ -1083,6 +1041,7 @@ end function _recursive_perfect_foresight!(context::Context, infoperiod, initial_values, options) datafile = options.datafile periods = options.periods - infoperiod + 1 + m = context.models[1] ncol = m.n_bkwrd + m.n_current + m.n_fwrd + 2 * m.n_both tmp_nbr = m.dynamic_tmp_nbr::Vector{Int64} @@ -1148,11 +1107,14 @@ function _recursive_perfect_foresight!(context::Context, infoperiod, initial_val guess_values, initial_values, terminal_values, - options.linear_solve_algo, dynamic_ws, perfect_foresight_ws.flipinfo, infoperiod, + linear_solve_algo = options.linear_solve_algo, maxit = options.maxit, + method = options.method, + nonlinear_solve_algo = options.solve_algo, + show_trace = options.show_trace, tolf = options.tolf, tolx = options.tolx ) @@ -1164,6 +1126,7 @@ function _recursive_perfect_foresight!(context::Context, infoperiod, initial_val guess_values, initial_values, terminal_values, + options.solve_algo, options.linear_solve_algo, dynamic_ws, maxit = options.maxit, diff --git a/src/steady_state/SteadyState.jl b/src/steady_state/SteadyState.jl index 8410dcd4..39e58922 100644 --- a/src/steady_state/SteadyState.jl +++ b/src/steady_state/SteadyState.jl @@ -31,7 +31,7 @@ SteadyOptions type maxit::Int64 - maximum number of iterations [50] tolf::Float64 - tolerance criterium for equations error [1e-8] tolx::Float64 - norm difference in x between two successive iterates under which convergence is declared. [0] - solve_algo::NonLinearSolveAlgos - algorithm for nonlinear equations solver [trustregion] + solve_algo::NonLinearSolveAlgos - algorithm for nonlinear equations solver ["NonlinearSolve"] homotopy_mode::HomotopyModes - homotopy mode [None] homotopy_steps::Int64 - homotopy steps [10] nocheck::Bool - don't check steady state values provided by the user [false] @@ -105,16 +105,19 @@ function steadystate!(; context::Context=context, homotopy_mode = SimultaneousFixedSteps, homotopy_steps = 0, maxit = 50, + linear_solve_algo = nothing, + method = TrustRegion(), nocheck = false, - nonlinear_solve_algo = TrustRegion(), + nonlinear_solve_algo = "NonlinearSolve", + show_trace = false, tolf = cbrt(eps()), - tolx = 0.0 + tolx = cbrt(eps()) ) model = context.models[1] modfileinfo = context.modfileinfo trends = context.results.model_results[1].trends work = context.work - + set_or_zero!(trends.endogenous_steady_state, work.initval_endogenous, model.endogenous_nbr) @@ -130,7 +133,7 @@ function steadystate!(; context::Context=context, model.exogenous_nbr) end if homotopy_steps == 0 - compute_steady_state!(context, maxit = maxit, nocheck = nocheck, nonlinear_solve_algo = nonlinear_solve_algo, tolf = tolf) + compute_steady_state!(context, maxit = maxit, nocheck = nocheck, linear_solve_algo = linear_solve_algo, method = method, nonlinear_solve_algo = nonlinear_solve_algo, show_trace = show_trace, tolf = tolf, tolx = tolx) else homotopy_steady!(context, homotopy_mode, homotopy_steps, maxit = maxit, nonlinear_solve_algo, tolf = tolf) end @@ -179,7 +182,7 @@ function check_steadystate_model(context::Context, nocheck) end end -function compute_steady_state!(context::Context; maxit = 50, nocheck = false, nonlinear_solve_algo = TrustRegion(), tolf = cbrt(eps())) +function compute_steady_state!(context::Context; maxit = 50, linear_solve_algo = nothing, method = TrustRegion(), nocheck = false, nonlinear_solve_algo = "NonlinearSolve", show_trace = false, tolf = cbrt(eps()), tolx = cbrt(eps())) model = context.models[1] modfileinfo = context.modfileinfo trends = context.results.model_results[1].trends @@ -189,6 +192,7 @@ function compute_steady_state!(context::Context; maxit = 50, nocheck = false, no isempty(trends.endogenous_steady_state) && (trends.endogenous_steady_state = Vector{Float64}(undef, endogenous_nbr)) isempty(trends.exogenous_steady_state) && (trends.exogenous_steady_state = Vector{Float64}(undef, exogenous_nbr)) if check_steadystate_model(context, nocheck) + evaluate_steady_state!(trends.endogenous_steady_state, trends.exogenous_steady_state, work.params) @@ -220,7 +224,14 @@ function compute_steady_state!(context::Context; maxit = 50, nocheck = false, no (x0 .= Float64.(trends.endogenous_steady_state)) if !nocheck trends.endogenous_steady_state .= - solve_steady_state!(context, x0, exogenous, maxit = maxit, nonlinear_solve_algo = nonlinear_solve_algo, tolf = tolf) + solve_steady_state!(context, x0, exogenous, + maxit = maxit, + linear_solve_algo = linear_solve_algo, + method = method, + nonlinear_solve_algo = nonlinear_solve_algo, + show_trace = show_trace, + tolf = tolf, + tolx = tolx) end # terminal steady state if modfileinfo.has_endval @@ -230,7 +241,14 @@ function compute_steady_state!(context::Context; maxit = 50, nocheck = false, no !isempty(trends.endogenous_terminal_steady_state) && (x0 .= Float64.(trends.endogenous_terminal_steady_state)) !nocheck && (trends.endogenous_terminal_steady_state .= - solve_steady_state!(context, x0, exogenous, maxit = maxit, nonlinear_solve_algo = nonlinear_solve_algo, tolf = tolf)) + solve_steady_state!(context, x0, exogenous, + maxit = maxit, + linear_solve_algo = linear_solve_algo, + method = method, + nonlinear_solve_algo = nonlinear_solve_algo, + show_trace = show_trace, + tolf = tolf, + tolx = tolx)) end end end @@ -238,7 +256,7 @@ end """ steadystate_display(steady_state, context; title_complement = "") - + Display the steady state of the model """ function steadystate_display(steady_state::AbstractVector{<:Real}, @@ -294,47 +312,72 @@ function solve_steady_state!(context::Context, x0::AbstractVector{<:Real}, exogenous::AbstractVector{<:Real}; maxit = 50, - nonlinear_solve_algo = TrustRegion(), - tolf = cbrt(eps())) - @show "OK1" - try - return solve_steady_state_!(context, x0, exogenous; maxit = maxit, nonlinear_solve_algo = nonlinear_solve_algo, tolf = tolf) - catch e - if !isempty(x0) && isa(e, Dynare.DynareSteadyStateComputationFailed) - i = 1 - while i <= maxit - x00 = rand(0.95:0.01:1.05, length(x0)).*x0 - try - return solve_steady_state_!(context, x00, exogenous, maxit = maxit, nonlinear_solve_algo = nonlinear_solve_algo, tolf = tolf) - catch - end - i += 1 - end - if i > maxit - rethrow(e) - end - else - rethrow(e) - end - end + linear_solve_algo = nothing, + method = TrustRegion(), + nonlinear_solve_algo = "NonlinearSolve", + show_trace = false, + tolf = cbrt(eps()), + tolx = cbrt(eps())) + try + return solve_steady_state_!(context, x0, exogenous, + maxit = maxit, + linear_solve_algo = linear_solve_algo, + method = method, + nonlinear_solve_algo = nonlinear_solve_algo, + show_trace = show_trace, + tolf = tolf, + tolx = tolx) + catch e + @show e + if !isempty(x0) && isa(e, Dynare.DynareSteadyStateComputationFailed) + i = 1 + while i <= maxit + x00 = rand(0.95:0.01:1.05, length(x0)).*x0 + try + return solve_steady_state_!(context, x00, exogenous, + linear_solve_algo = linear_solve_algo, + maxit = maxit, + method = method, + nonlinear_solve_algo = nonlinear_solve_algo, + show_trace = show_trace, + tolf = tolf, + tolx = tolx) + catch + end + i += 1 + end + if i > maxit + rethrow(e) + end + else + rethrow(e) + end + end end using LinearSolve """ solve_steady_state_!(context::Context; x0::Vector{Float64}, exogenous::AbstractVector{<:Real}; + linear_solve_algo = nothing, maxit = 50, - tolf = cbrt(eps())) + method = TrustRegion(), + show_trace = show_trace, + tolf = cbrt(eps()), + tolx = cbrt(eps())) Solve the static model to obtain the steady state """ function solve_steady_state_!(context::Context, x0::AbstractVector{<:Real}, exogenous::AbstractVector{<:Real}; + linear_solve_algo = nothing, maxit = 50, - nonlinear_solve_algo = TrustRegion(), - tolf = cbrt(eps())) - @show "OK2" + method = TrustRegion(), + nonlinear_solve_algo = "NonlinearSolve", + show_trace = false, + tolf = cbrt(eps()), + tolx = tolx) ws = StaticWs(context) model = context.models[1] work = context.work @@ -352,28 +395,17 @@ function solve_steady_state_!(context::Context, work.params) results = context.results.model_results[1] - j!(A, x0, params) - display(Matrix(A)) - lp = LinearProblem(A, residuals) - dy1 = LinearSolve.solve(lp) - dy2 = LinearSolve.solve(lp, PardisoJL()) - @show norm(dy1-dy2) - - # of = OnceDifferentiable(f!, J!, vec(x0), residuals, A) - # result = nlsolve(of, x0; method = :robust_trust_region, show_trace = false, ftol = tolf, iterations = maxit) - fj = NonlinearFunction(f!, jac = j!, jac_prototype = A) - @show x0 - prob = NonlinearProblem(fj, x0, params) - - result = NonlinearSolve.solve(prob, nonlinear_solve_algo, show_trace=Val(true), abstol = tolf) -# result = NonlinearSolve.solve(prob, show_trace=Val(true), abstol = tolf) - - @debug result - @show result.retcode - if result.retcode == ReturnCode.Success - return result.u - else - @debug "Steady state computation failed with\n $result" + try + return dynare_nonlinear_solvers(f!, j!, A, [], [], x0, exogenous, work.params, + linear_solve_algo = linear_solve_algo, + maxit = maxit, + method = method, + solve_algo = nonlinear_solve_algo, + show_trace = show_trace, + tolf = tolf, + tolx = tolx) + catch e + @debug "Steady state computation failed" throw(DynareSteadyStateComputationFailed()) end end diff --git a/test/models/initialization/neoclassical1.mod b/test/models/initialization/neoclassical1.mod index 25b600c4..15ad1887 100644 --- a/test/models/initialization/neoclassical1.mod +++ b/test/models/initialization/neoclassical1.mod @@ -34,14 +34,15 @@ end; // display the steady-state steady; -// initial conditionls +// initial conditions histval; k(0) = 0.5*((1-beta*(1-delta))/(beta*alpha))^(1/(alpha-1)); end; // computing the simulation -perfect_foresight_setup(periods=200); -perfect_foresight_solver; +//perfect_foresight_setup(periods=200); +//perfect_foresight_solver; +perfect_foresight!(periods=200); // ploting results plot(simulation("k")) plot(simulation("c")) diff --git a/test/runtests.jl b/test/runtests.jl index 074efcc9..c6a0268c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,4 +1,6 @@ using Dynare +using Pardiso +using PATHSolver using Test # CSV.write fails with AxisArrays diff --git a/test/test_scenario.jl b/test/test_scenario.jl index 53594067..b99bb702 100644 --- a/test/test_scenario.jl +++ b/test/test_scenario.jl @@ -1,7 +1,7 @@ using AxisArrayTables using Dynare using ExtendedDates -using JLD2 +#using JLD2 using LinearAlgebra using Serialization using Test @@ -142,7 +142,7 @@ context = @dynare "models/example1pf/example1pf_conditional" "stoponerror" periods = 4 f!, J!, guess_values, initial_values, terminal_values, dynamic_variables, residuals, JJ, exogenous, temp_vec, nzval, nzval1 = make_f_J(context, context.work.scenario, periods, 1) - f!(residuals, guess_values) + f!(residuals, guess_values, context.work.params) @test guess_values[1] == context.results.model_results[1].trends.endogenous_steady_state[1] @test exogenous[2] == 0 @@ -150,8 +150,7 @@ context = @dynare "models/example1pf/example1pf_conditional" "stoponerror" @test exogenous[6] == 0 @test guess_values[14] == context.results.model_results[1].trends.endogenous_steady_state[2] - J!(JJ, guess_values) - + J!(JJ, guess_values, context.work.params) J1target = zeros(24) J1target[6] = -1 @test JJ[:, 1] == J1target @@ -169,7 +168,7 @@ context = @dynare "models/example1pf/example1pf_conditional" "stoponerror" guess_values[1] = 0 exogenous[1] = 1.0 - f!(residuals, guess_values) + f!(residuals, guess_values, context.work.params) y, c, k, a, h, b = context.results.model_results[1].trends.endogenous_steady_state beta, rho, alpha, delta, theta, psi, tau = context.work.params @@ -181,7 +180,7 @@ context = @dynare "models/example1pf/example1pf_conditional" "stoponerror" @test residuals[4] ≈ k - (exp(b)*(y_star - c) + (1 - delta)*k) @test all(abs.(residuals[5:24]) .< 2*eps()) - J!(JJ, guess_values) + J!(JJ, guess_values, context.work.params) Dy = -JJ\residuals end @@ -200,7 +199,7 @@ context = @dynare "models/example1pf/example1pf_conditional" "stoponerror" @test guess_values[1] == 0 @test exogenous[1] == 1.1 - f!(residuals, guess_values) + f!(residuals, guess_values, context.work.params) y, c, k, a, h, b = context.results.model_results[1].trends.endogenous_steady_state beta, rho, alpha, delta, theta, psi, tau = context.work.params @@ -212,7 +211,7 @@ context = @dynare "models/example1pf/example1pf_conditional" "stoponerror" @test residuals[4] ≈ k - (exp(b)*(y_star - c) + (1 - delta)*k) @test all(abs.(residuals[5:600]) .< 2*eps()) - J!(JJ, guess_values) + J!(JJ, guess_values, context.work.params) Dy = -JJ\residuals @@ -220,8 +219,8 @@ context = @dynare "models/example1pf/example1pf_conditional" "stoponerror" for i=1:4 y += Dy - f!(residuals, y) - J!(JJ, y) + f!(residuals, y, context.work.params) + J!(JJ, y, context.work.params) Dy = -JJ\residuals end end @@ -236,17 +235,17 @@ context = @dynare "models/example1pf/example1pf_conditional" "stoponerror" Dynare.set_future_information!(guess_values, exogenous, context, periods, 1) Dynare.flip!(guess_values, exogenous, FI.ix_stack) - f!(residuals, guess_values) + f!(residuals, guess_values, context.work.params) - J!(JJ, guess_values) + J!(JJ, guess_values, context.work.params) Dy = -JJ\residuals y = copy(guess_values) for i=1:4 y += Dy - f!(residuals, y) - J!(JJ, y) + f!(residuals, y, context.work.params) + J!(JJ, y, context.work.params) Dy = -JJ\residuals end