eliminate abs(NSSS) < 1e-12

thorek1 · Sep 23, 2023 · 3a44e78 · 3a44e78
1 parent f1db144
commit 3a44e78
Show file tree

Hide file tree

Showing 3 changed files with 11 additions and 153 deletions.
diff --git a/docs/src/unfinished_docs/todo.md b/docs/src/unfinished_docs/todo.md
@@ -5,7 +5,6 @@
 - [ ] implement occasionally binding constraints with shocks
 - [ ] recheck get function examples and docs
 - [ ] riccati with analytical derivatives (much faster if sparse) instead of implicit diff
-- [ ] set to 0 SS values < 1e-12
 - [ ] write method of moments how to
 - [ ] autocorr and covariance with derivatives. return 3d array
 - [ ] Docs: document outputs and associated functions to work with function
@@ -71,6 +70,7 @@
 - [ ] Find any SS by optimising over both SS guesses and parameter inputs
 - [ ] weed out SS solver and saved objects
 
+- [x] set to 0 SS values < 1e-12
 - [x] sylvester with analytical derivatives (much faster if sparse) instead of implicit diff - yes but there are still way too large matrices being realised. implicitdiff is better here
 - [x] autocorr to statistics output and in general for higher order pruned sols
 - [x] fix product moments and test for cases with more than 2 shocks

diff --git a/src/MacroModelling.jl b/src/MacroModelling.jl
@@ -847,8 +847,7 @@ function levenberg_marquardt(f::Function,
     λ̂̅¹::T   =       0.9815,
     λ̂̅²::T   =       1.0,
     transformation_level::S = 3,
-    backtracking_order::S = 2,
-    ) where {T <: AbstractFloat, S <: Integer}
+    backtracking_order::S = 2) where {T <: AbstractFloat, S <: Integer}
     # issues with optimization: https://www.gurobi.com/documentation/8.1/refman/numerics_gurobi_guidelines.html
 
     @assert size(lower_bounds) == size(upper_bounds) == size(initial_guess)
@@ -976,7 +975,9 @@ function levenberg_marquardt(f::Function,
         end
     end
 
-    return undo_transform(current_guess,transformation_level), (iterations, largest_step, largest_residual, f(undo_transform(current_guess,transformation_level)))
+    best_guess = undo_transform(current_guess,transformation_level)
+
+    return best_guess, (iterations, largest_step, largest_residual, f(best_guess))
 end
 
 
@@ -1667,7 +1668,9 @@ function solve_steady_state!(𝓂::ℳ, symbolic_SS, Symbolics::symbolics; verbo
                             $(SS_solve_func...)
                             if scale == 1
                                 # return ComponentVector([$(sort(union(𝓂.var,𝓂.exo_past,𝓂.exo_future))...), $(𝓂.calibration_equations_parameters...)], Axis([sort(union(𝓂.exo_present,𝓂.var))...,𝓂.calibration_equations_parameters...])), solution_error
-                                return [$(Symbol.(replace.(string.(sort(union(𝓂.var,𝓂.exo_past,𝓂.exo_future))), r"ᴸ⁽⁻?[⁰¹²³⁴⁵⁶⁷⁸⁹]+⁾" => ""))...), $(𝓂.calibration_equations_parameters...)] , solution_error
+                                NSSS_solution = [$(Symbol.(replace.(string.(sort(union(𝓂.var,𝓂.exo_past,𝓂.exo_future))), r"ᴸ⁽⁻?[⁰¹²³⁴⁵⁶⁷⁸⁹]+⁾" => ""))...), $(𝓂.calibration_equations_parameters...)]
+                                NSSS_solution[abs.(NSSS_solution) .< 1e-12] .= 0
+                                return NSSS_solution , solution_error
                             end
                         end
                     end
@@ -4281,7 +4284,7 @@ function calculate_covariance(parameters::Vector{<: Real}, 𝓂::ℳ; verbose::B
 
     values = vcat(vec(A), vec(collect(-CC)))
 
-    covar_raw, _ = solve_sylvester_equation_forward(values, coords = coordinates, dims = dimensions, solver = :doubling)
+    covar_raw, _ = solve_sylvester_equation_AD(values, coords = coordinates, dims = dimensions, solver = :doubling)
     # covar_raw, _ = solve_sylvester_equation_AD_direct(values, coords = coordinates, dims = dimensions, solver = :doubling)
     # covar_raw, _ = solve_sylvester_equation_AD_direct([vec(A); vec(-CC)], dims = [size(A), size(CC)], solver = :bicgstab)
     # covar_raw, _ = solve_sylvester_equation_forward([vec(A); vec(-CC)], dims = [size(A), size(CC)])
@@ -4478,7 +4481,7 @@ function solve_sylvester_equation_conditions(ABC::Vector{<: Real},
         A = sparse(coords[1]...,vA,dims[1]...) |> ThreadedSparseArrays.ThreadedSparseMatrixCSC
         C = reshape(ABC[lengthA+1:end],dims[2]...)
         if solver != :doubling
-            B = A'
+            B = A' |> sparse |> ThreadedSparseArrays.ThreadedSparseMatrixCSC
         end
     elseif length(coords) == 3
         lengthA = length(coords[1][1])
@@ -4623,151 +4626,6 @@ function solve_sylvester_equation_forward(abc::Vector{ℱ.Dual{Z,S,N}};
 end
 
 
-# function solve_sylvester_equation_forward(ABC::AbstractVector{ℱ.Dual{Z,S,N}}; 
-#     dims::Vector{Tuple{Int,Int}}, 
-#     sparse_output::Bool = false,
-#     solver::Symbol = :gmres) where {Z,S,N}
-
-#     # unpack: AoS -> SoA
-#     ABCv = ℱ.value.(ABC)
-
-#     # you can play with the dimension here, sometimes it makes sense to transpose
-#     partials = mapreduce(ℱ.partials, hcat, ABC)'
-
-#     val, solved = solve_sylvester_equation_forward(ABCv, dims = dims, sparse_output = sparse_output, solver = solver)
-
-#     # get J(f, vs) * ps (cheating). Write your custom rule here
-#     BB = ℱ.jacobian(x -> solve_sylvester_equation_conditions(x, val, solved, dims = dims), ABCv)
-#     AA = ℱ.jacobian(x -> solve_sylvester_equation_conditions(ABCv, x, solved, dims = dims), val)
-
-#     Â = RF.lu(AA, check = false)
-
-#     if !ℒ.issuccess(Â)
-#         Â = ℒ.svd(AA)
-#     end
-
-#     jvp = -(Â \ BB) * partials
-
-#     # pack: SoA -> AoS
-#     return reshape(map(val, eachrow(jvp)) do v, p
-#         ℱ.Dual{Z}(v, p...) # Z is the tag
-#     end,size(val)), solved
-# end
-
-
-
-# function solve_sylvester_equation_forward(abc::SparseVector{ℱ.Dual{Z,S,N}}; 
-#     dims::Vector{Tuple{Int,Int}}, 
-#     sparse_output::Bool = false,
-#     solver::Symbol = :gmres) where {Z,S,N}
-
-#     # unpack: AoS -> SoA
-#     ABC, partials = separate_values_and_partials_from_sparsevec_dual(abc)
-
-#     # get f(vs)
-#     val, solved = solve_sylvester_equation_forward(ABC, dims = dims, sparse_output = sparse_output, solver = solver)
-
-#     lenA = dims[1][1] * dims[1][2]
-
-#     A = reconstruct_sparse_matrix(ABC[1 : lenA], dims[1])
-#     A¹ = sparse((ABC[1 : lenA]).nzind, (ABC[1 : lenA]).nzind, 1, lenA, lenA)
-
-#     if length(dims) == 3
-#         lenB = dims[2][1] * dims[2][2]
-#         B = reconstruct_sparse_matrix(ABC[lenA .+ (1 : lenB)], dims[2])
-#         B¹ = sparse((ABC[lenA .+ (1 : lenB)]).nzind, (ABC[lenA .+ (1 : lenB)]).nzind, 1, lenB, lenB)
-
-#         jacobian_A = A¹ * ℒ.kron(-val * B, ℒ.I(size(A,1)))
-#         jacobian_B = ℒ.kron(ℒ.I(size(B,1)), -A * val) * B¹
-
-#         b = hcat(jacobian_A', jacobian_B, ℒ.I(length(val)))
-#     elseif length(dims) == 2
-#         B = A'
-#         lenB = lenA
-
-#         jacobian_A = A¹ * ℒ.kron(-val * B, ℒ.I(size(A,1)))
-
-#         b = hcat(jacobian_A', ℒ.I(length(val)))
-#     end
-
-#     # get J(f, vs) * ps (cheating). Write your custom rule here. This used to be the conditions but here they are analytically derived.
-#     # a = reshape(permutedims(reshape(ℒ.I - ℒ.kron(A, B) ,size(B,1), size(A,1), size(A,1), size(B,1)), [2, 3, 4, 1]), size(A,1) * size(B,1), size(A,1) * size(B,1))
-
-#     reshape_matmul = LinearOperators.LinearOperator(Float64, size(b,1) * size(partials,2), size(b,1) * size(partials,2), false, false, 
-#         (sol,𝐱) -> begin 
-#         𝐗 = reshape(𝐱, (size(b,1),size(partials,2)))
-#         sol .= vec(reshape(permutedims(reshape(ℒ.I - ℒ.kron(A, B) ,size(B,1), size(A,1), size(A,1), size(B,1)), [2, 3, 4, 1]), size(A,1) * size(B,1), size(A,1) * size(B,1)) * 𝐗)
-#         return sol
-#     end)
-
-#     X, info = Krylov.gmres(reshape_matmul, -vec(b * partials))#, atol = tol)
-
-#     jvp = reshape(X, (size(b,1),size(partials,2)))
-
-
-#     # pack: SoA -> AoS
-#     return sparse(reshape(map(val, eachrow(jvp)) do v, p
-#         ℱ.Dual{Z}(v, p...) # Z is the tag
-#         end,size(val))), solved
-# end
-
-
-
-
-# function solve_sylvester_equation_forward(abc::DenseVector{ℱ.Dual{Z,S,N}};
-#     dims::Vector{Tuple{Int,Int}}, 
-#     sparse_output::Bool = false,
-#     solver::Symbol = :gmres) where {Z,S,N}
-
-#     # unpack: AoS -> SoA
-#     ABC = ℱ.value.(abc)
-
-#     # you can play with the dimension here, sometimes it makes sense to transpose
-#     partials = mapreduce(ℱ.partials, hcat, abc)'
-
-#     # get f(vs)
-#     val, solved = solve_sylvester_equation_forward(ABC, dims = dims, sparse_output = sparse_output, solver = solver)
-
-#     lenA = dims[1][1] * dims[1][2]
-
-#     A = reshape(ABC[1 : lenA], dims[1])
-
-#     if length(dims) == 3
-#         lenB = dims[2][1] * dims[2][2]
-#         B = reshape(ABC[lenA .+ (1 : lenB)], dims[2])
-
-#         jacobian_A = ℒ.kron(-val * B, ℒ.I(size(A,1)))
-#         jacobian_B = ℒ.kron(ℒ.I(size(B,1)), -A * val)
-
-#         b = hcat(jacobian_A', jacobian_B, ℒ.I(length(val)))
-#     elseif length(dims) == 2
-#         B = A'
-#         jacobian_A = ℒ.kron(-val * B, ℒ.I(size(A,1)))
-
-#         b = hcat(jacobian_A', ℒ.I(length(val)))
-#     end
-
-#     # get J(f, vs) * ps (cheating). Write your custom rule here. This used to be the conditions but here they are analytically derived.
-#     # a = reshape(permutedims(reshape(ℒ.I - ℒ.kron(A, B) ,size(B,1), size(A,1), size(A,1), size(B,1)), [2, 3, 4, 1]), size(A,1) * size(B,1), size(A,1) * size(B,1))
-
-#     reshape_matmul = LinearOperators.LinearOperator(Float64, size(b,1) * size(partials,2), size(b,1) * size(partials,2), false, false, 
-#         (sol,𝐱) -> begin 
-#         𝐗 = reshape(𝐱, (size(b,1),size(partials,2)))
-#         sol .= vec(reshape(permutedims(reshape(ℒ.I - ℒ.kron(A, B) ,size(B,1), size(A,1), size(A,1), size(B,1)), [2, 3, 4, 1]), size(A,1) * size(B,1), size(A,1) * size(B,1)) * 𝐗)
-#         return sol
-#     end)
-
-#     X, info = Krylov.gmres(reshape_matmul, -vec(b * partials))#, atol = tol)
-
-#     jvp = reshape(X, (size(b,1),size(partials,2)))
-
-#     # pack: SoA -> AoS
-#     return reshape(map(val, eachrow(jvp)) do v, p
-#         ℱ.Dual{Z}(v, p...) # Z is the tag
-#         end,size(val)), solved
-# end
-
-
 solve_sylvester_equation_AD = ID.ImplicitFunction(solve_sylvester_equation_forward, 
                                                 solve_sylvester_equation_conditions)
 

diff --git a/src/get_functions.jl b/src/get_functions.jl
@@ -1653,7 +1653,7 @@ function get_variance_decomposition(𝓂::ℳ;
 
         values = vcat(vec(A), vec(collect(-CC)))
 
-        covar_raw, _ = solve_sylvester_equation_forward(values, coords = coordinates, dims = dimensions, solver = :doubling)
+        covar_raw, _ = solve_sylvester_equation_AD(values, coords = coordinates, dims = dimensions, solver = :doubling)
         # covar_raw, _ = solve_sylvester_equation_AD_direct([vec(A); vec(-CC)], dims = [size(A), size(CC)], solver = :bicgstab)
         # covar_raw, _ = solve_sylvester_equation_forward([vec(A); vec(-CC)], dims = [size(A), size(CC)])