Merge pull request #45 from bankofcanada/dev

Dev 0.4.2

.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

Project.toml

@@ -1,7 +1,7 @@
 name = "StateSpaceEcon"
 uuid = "e4c825b0-b65c-11ea-0b5a-6176b64e7b7f"
 authors = ["Atai Akunov <[email protected]>", "Boyan Bejanov <[email protected]>", "Nicholas Labelle St-Pierre <[email protected]>"]
-version = "0.4.1"
+version = "0.5.0"
 
 [deps]
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
@@ -10,19 +10,24 @@ 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"
+SuiteSparse = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"
 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"
-ModelBaseEcon = "0.5.2"
+ModelBaseEcon = "0.6"
+Pardiso = "0.5"
 Suppressor = "0.2"
-TimeSeriesEcon = "0.5.1"
+TimerOutputs = "0.5"
+TimeSeriesEcon = "0.6"
 julia = "1.7"
 
 [extras]

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]...]

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,14 @@ import ModelBaseEcon.hassolverdata
 import ModelBaseEcon.getsolverdata
 import ModelBaseEcon.setsolverdata!
 
+using Pardiso
+
+
+using SuiteSparse
+using SuiteSparse.UMFPACK
+
+include("stackedtime/sparse.jl")
+
 include("stackedtime/fctypes.jl")
 include("stackedtime/misc.jl")
 include("stackedtime/solverdata.jl")
@@ -61,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
 

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)

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
@@ -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) -> A \ b,
     linesearch::Bool=getoption(m, :linesearch, false),
     fctype::FinalCondition=getoption(m, :fctype, fcgiven),
     _debug=false
@@ -75,11 +74,11 @@ 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.
-    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)
@@ -109,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
@@ -129,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
@@ -149,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)
@@ -161,23 +160,23 @@ 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)
                         end
                     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])
@@ -196,29 +195,27 @@ 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)
-
+        sf_solve!(Jac_control, SDMAT)
     end
 
     # now split and decorate the shock-decomposition matrix
     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_data = result.sd[v] 
+            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, :]
@@ -240,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