From d315cb2291fa182d25953f5b67d5228a9e192dd8 Mon Sep 17 00:00:00 2001 From: thorek1 Date: Wed, 18 Dec 2024 20:30:28 +0100 Subject: [PATCH] timeroutputs commented out --- Project.toml | 2 - src/MacroModelling.jl | 222 +++++----- src/algorithms/lyapunov.jl | 52 +-- src/algorithms/quadratic_matrix_equation.jl | 94 ++-- src/algorithms/sylvester.jl | 186 ++++---- src/filter/inversion.jl | 216 +++++----- src/filter/kalman.jl | 74 ++-- src/get_functions.jl | 24 +- src/perturbation.jl | 448 ++++++++++---------- 9 files changed, 660 insertions(+), 658 deletions(-) diff --git a/Project.toml b/Project.toml index ea6813ee..4c370652 100644 --- a/Project.toml +++ b/Project.toml @@ -40,7 +40,6 @@ Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb" SymPyPythonCall = "bc8888f7-b21e-4b7c-a06a-5d9c9496438c" Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" ThreadedSparseArrays = "59d54670-b8ac-4d81-ab7a-bb56233e17ab" -TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" Unicode = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5" [weakdeps] @@ -94,7 +93,6 @@ SymPyPythonCall = "0.2, 0.3, 0.4" Symbolics = "5, 6" Test = "1" ThreadedSparseArrays = "0.2.3, 0.3" -TimerOutputs = "0.5" Turing = "0.30, 0.31, 0.32, 0.33, 0.34, 0.35" Unicode = "1" Zygote = "0.6" diff --git a/src/MacroModelling.jl b/src/MacroModelling.jl index 60ea2727..cb90041d 100644 --- a/src/MacroModelling.jl +++ b/src/MacroModelling.jl @@ -10,8 +10,8 @@ import SpecialFunctions import SymPyPythonCall as SPyPyC import Symbolics import Accessors -import TimerOutputs -import TimerOutputs: TimerOutput, @timeit, @timeit_debug +# import TimerOutputs +# import TimerOutputs: TimerOutput, @timeit, @timeit_debug # import NaNMath # import Memoization: @memoize # import LRUCache: LRU @@ -1059,11 +1059,11 @@ end function compressed_kron³(a::AbstractMatrix{T}; rowmask::Vector{Int} = Int[], colmask::Vector{Int} = Int[], - timer::TimerOutput = TimerOutput(), + # timer::TimerOutput = TimerOutput(), tol::AbstractFloat= eps()) where T <: Real - @timeit_debug timer "Compressed 3rd kronecker power" begin + # @timeit_debug timer "Compressed 3rd kronecker power" begin - @timeit_debug timer "Preallocation" begin + # @timeit_debug timer "Preallocation" begin a_is_adjoint = typeof(a) <: ℒ.Adjoint{T,Matrix{T}} @@ -1113,18 +1113,18 @@ function compressed_kron³(a::AbstractMatrix{T}; # k̄ = Threads.Atomic{Int}(0) # effectively slower than the non-threaded version k = 0 - end # timeit_debug + # end # timeit_debug - @timeit_debug timer "findnz" begin + # @timeit_debug timer "findnz" begin # Find unique non-zero row and column indices rowinds, colinds, _ = findnz(a) ui = unique(rowinds) uj = unique(colinds) - end # timeit_debug + # end # timeit_debug - @timeit_debug timer "Loop" begin + # @timeit_debug timer "Loop" begin # Triple nested loops for (i1 ≤ j1 ≤ k1) and (i2 ≤ j2 ≤ k2) # Polyester.@batch threadlocal=(Vector{Int}(), Vector{Int}(), Vector{T}()) for i1 in ui # Polyester.@batch minbatch = 10 for i1 in ui @@ -1214,9 +1214,9 @@ function compressed_kron³(a::AbstractMatrix{T}; end end - end # timeit_debug + # end # timeit_debug - @timeit_debug timer "Resize" begin + # @timeit_debug timer "Resize" begin # out = map(fetch, threadlocal) @@ -1232,8 +1232,8 @@ function compressed_kron³(a::AbstractMatrix{T}; resize!(J, k) resize!(V, k) - end # timeit_debug - end # timeit_debug + # end # timeit_debug + # end # timeit_debug # Create the sparse matrix from the collected indices and values if VERSION >= v"1.10" @@ -4082,7 +4082,7 @@ function calculate_second_order_stochastic_steady_state(parameters::Vector{M}, pruning::Bool = false, quadratic_matrix_equation_solver::Symbol = :schur, sylvester_algorithm::Symbol = :doubling, - timer::TimerOutput = TimerOutput(), + # timer::TimerOutput = TimerOutput(), tol::AbstractFloat = 1e-12)::Tuple{Vector{M}, Bool, Vector{M}, M, AbstractMatrix{M}, SparseMatrixCSC{M}, AbstractMatrix{M}, SparseMatrixCSC{M}} where M # @timeit_debug timer "Calculate NSSS" begin @@ -4211,9 +4211,9 @@ function calculate_second_order_stochastic_steady_state(::Val{:Newton}, 𝐒₂::AbstractSparseMatrix{Float64}, x::Vector{Float64}, 𝓂::ℳ; - timer::TimerOutput = TimerOutput(), + # timer::TimerOutput = TimerOutput(), tol::AbstractFloat = 1e-14) - @timeit_debug timer "Setup matrices" begin + # @timeit_debug timer "Setup matrices" begin nᵉ = 𝓂.timings.nExo @@ -4231,9 +4231,9 @@ function calculate_second_order_stochastic_steady_state(::Val{:Newton}, max_iters = 100 # SSS .= 𝐒₁ * aug_state + 𝐒₂ * ℒ.kron(aug_state, aug_state) / 2 + 𝐒₃ * ℒ.kron(ℒ.kron(aug_state,aug_state),aug_state) / 6 - end # timeit_debug + # end # timeit_debug - @timeit_debug timer "Iterations" begin + # @timeit_debug timer "Iterations" begin for i in 1:max_iters ∂x = (A + B * ℒ.kron(vcat(x,1), ℒ.I(𝓂.timings.nPast_not_future_and_mixed)) - ℒ.I(𝓂.timings.nPast_not_future_and_mixed)) @@ -4256,7 +4256,7 @@ function calculate_second_order_stochastic_steady_state(::Val{:Newton}, ℒ.axpy!(-1, Δx, x) end - end # timeit_debug + # end # timeit_debug return x, isapprox(A * x + B̂ * ℒ.kron(vcat(x,1), vcat(x,1)) / 2, x, rtol = tol) end @@ -4269,7 +4269,7 @@ function calculate_second_order_stochastic_steady_state(::Val{:Newton}, 𝐒₂::AbstractSparseMatrix{ℱ.Dual{Z,S,N}}, x::Vector{ℱ.Dual{Z,S,N}}, 𝓂::ℳ; - timer::TimerOutput = TimerOutput(), + # timer::TimerOutput = TimerOutput(), tol::AbstractFloat = 1e-14) where {Z,S,N} 𝐒₁̂ = ℱ.value.(𝐒₁) @@ -4339,10 +4339,10 @@ function rrule(::typeof(calculate_second_order_stochastic_steady_state), 𝐒₂::AbstractSparseMatrix{Float64}, x::Vector{Float64}, 𝓂::ℳ; - timer::TimerOutput = TimerOutput(), + # timer::TimerOutput = TimerOutput(), tol::AbstractFloat = 1e-14) - @timeit_debug timer "Calculate SSS - forward" begin - @timeit_debug timer "Setup indices" begin + # @timeit_debug timer "Calculate SSS - forward" begin + # @timeit_debug timer "Setup indices" begin nᵉ = 𝓂.timings.nExo @@ -4357,9 +4357,9 @@ function rrule(::typeof(calculate_second_order_stochastic_steady_state), B = 𝐒₂[𝓂.timings.past_not_future_and_mixed_idx,kron_s⁺_s] B̂ = 𝐒₂[𝓂.timings.past_not_future_and_mixed_idx,kron_s⁺_s⁺] - end # timeit_debug + # end # timeit_debug - @timeit_debug timer "Iterations" begin + # @timeit_debug timer "Iterations" begin max_iters = 100 # SSS .= 𝐒₁ * aug_state + 𝐒₂ * ℒ.kron(aug_state, aug_state) / 2 + 𝐒₃ * ℒ.kron(ℒ.kron(aug_state,aug_state),aug_state) / 6 @@ -4389,11 +4389,11 @@ function rrule(::typeof(calculate_second_order_stochastic_steady_state), ∂𝐒₁ = zero(𝐒₁) ∂𝐒₂ = zero(𝐒₂) - end # timeit_debug - end # timeit_debug + # end # timeit_debug + # end # timeit_debug function second_order_stochastic_steady_state_pullback(∂x) - @timeit_debug timer "Calculate SSS - pullback" begin + # @timeit_debug timer "Calculate SSS - pullback" begin S = -∂x[1]' / (A + B * ℒ.kron(vcat(x,1), ℒ.I(𝓂.timings.nPast_not_future_and_mixed)) - ℒ.I(𝓂.timings.nPast_not_future_and_mixed)) @@ -4401,7 +4401,7 @@ function rrule(::typeof(calculate_second_order_stochastic_steady_state), ∂𝐒₂[𝓂.timings.past_not_future_and_mixed_idx,kron_s⁺_s⁺] = S' * ℒ.kron(vcat(x,1), vcat(x,1))' / 2 - end # timeit_debug + # end # timeit_debug return NoTangent(), NoTangent(), ∂𝐒₁, ∂𝐒₂, NoTangent(), NoTangent(), NoTangent() end @@ -4417,7 +4417,7 @@ function calculate_third_order_stochastic_steady_state( parameters::Vector{M}, pruning::Bool = false, quadratic_matrix_equation_solver::Symbol = :schur, sylvester_algorithm::Symbol = :bicgstab, - timer::TimerOutput = TimerOutput(), + # timer::TimerOutput = TimerOutput(), tol::AbstractFloat = 1e-12)::Tuple{Vector{M}, Bool, Vector{M}, M, AbstractMatrix{M}, SparseMatrixCSC{M}, SparseMatrixCSC{M}, AbstractMatrix{M}, SparseMatrixCSC{M}, SparseMatrixCSC{M}} where M SS_and_pars, (solution_error, iters) = get_NSSS_and_parameters(𝓂, parameters, verbose = verbose, timer = timer) @@ -4550,7 +4550,7 @@ function calculate_third_order_stochastic_steady_state(::Val{:Newton}, 𝐒₃::AbstractSparseMatrix{Float64}, x::Vector{Float64}, 𝓂::ℳ; - timer::TimerOutput = TimerOutput(), + # timer::TimerOutput = TimerOutput(), tol::AbstractFloat = 1e-14) nᵉ = 𝓂.timings.nExo @@ -4746,22 +4746,22 @@ end function solve!(𝓂::ℳ; - parameters::ParameterType = nothing, - dynamics::Bool = false, - algorithm::Symbol = :first_order, - obc::Bool = false, - verbose::Bool = false, - silent::Bool = false, - timer::TimerOutput = TimerOutput(), - tol::AbstractFloat = 1e-12) + parameters::ParameterType = nothing, + dynamics::Bool = false, + algorithm::Symbol = :first_order, + obc::Bool = false, + verbose::Bool = false, + silent::Bool = false, + # timer::TimerOutput = TimerOutput(), + tol::AbstractFloat = 1e-12) @assert algorithm ∈ all_available_algorithms - @timeit_debug timer "Write parameter inputs" begin + # @timeit_debug timer "Write parameter inputs" begin write_parameters_input!(𝓂, parameters, verbose = verbose) - end # timeit_debug + # end # timeit_debug if 𝓂.solution.perturbation.second_order_auxilliary_matrices.𝛔 == SparseMatrixCSC{Int, Int64}(ℒ.I,0,0) && algorithm ∈ [:second_order, :pruned_second_order] @@ -4785,21 +4785,21 @@ function solve!(𝓂::ℳ; ((:third_order == algorithm) && ((:third_order ∈ 𝓂.solution.outdated_algorithms) || (obc && obc_not_solved))) || ((:pruned_third_order == algorithm) && ((:pruned_third_order ∈ 𝓂.solution.outdated_algorithms) || (obc && obc_not_solved))) - @timeit_debug timer "Solve for NSSS (if necessary)" begin + # @timeit_debug timer "Solve for NSSS (if necessary)" begin SS_and_pars, (solution_error, iters) = 𝓂.solution.outdated_NSSS ? get_NSSS_and_parameters(𝓂, 𝓂.parameter_values, verbose = verbose) : (𝓂.solution.non_stochastic_steady_state, (eps(), 0)) - end # timeit_debug + # end # timeit_debug @assert solution_error < tol "Could not find non stochastic steady steady." - @timeit_debug timer "Calculate Jacobian" begin + # @timeit_debug timer "Calculate Jacobian" begin ∇₁ = calculate_jacobian(𝓂.parameter_values, SS_and_pars, 𝓂)# |> Matrix - end # timeit_debug + # end # timeit_debug - @timeit_debug timer "Calculate first order solution" begin + # @timeit_debug timer "Calculate first order solution" begin if algorithm == :first_order_doubling qme_solver = :doubling @@ -4811,7 +4811,7 @@ function solve!(𝓂::ℳ; if solved 𝓂.solution.perturbation.qme_solution = qme_sol end - end # timeit_debug + # end # timeit_debug @assert solved "Could not find stable first order solution." @@ -6065,9 +6065,9 @@ end function calculate_jacobian(parameters::Vector{M}, SS_and_pars::Vector{N}, - 𝓂::ℳ; - timer::TimerOutput = TimerOutput())::Matrix{M} where {M,N} - @timeit_debug timer "Calculate jacobian" begin + 𝓂::ℳ)::Matrix{M} where {M,N} + # timer::TimerOutput = TimerOutput())::Matrix{M} where {M,N} + # @timeit_debug timer "Calculate jacobian" begin SS = SS_and_pars[1:end - length(𝓂.calibration_equations)] calibrated_parameters = SS_and_pars[(end - length(𝓂.calibration_equations)+1):end] @@ -6092,7 +6092,7 @@ function calculate_jacobian(parameters::Vector{M}, # lk = ReentrantLock() - @timeit_debug timer "Loop" begin + # @timeit_debug timer "Loop" begin Polyester.@batch minbatch = 200 for f in 𝓂.model_jacobian[1] # for f in 𝓂.model_jacobian[1] @@ -6116,8 +6116,8 @@ function calculate_jacobian(parameters::Vector{M}, 𝓂.model_jacobian[3][𝓂.model_jacobian[2]] .= vals - end # timeit_debug - end # timeit_debug + # end # timeit_debug + # end # timeit_debug return 𝓂.model_jacobian[3] end @@ -6126,16 +6126,16 @@ end function rrule(::typeof(calculate_jacobian), parameters, SS_and_pars, - 𝓂; - timer::TimerOutput = TimerOutput()) - @timeit_debug timer "Calculate jacobian - forward" begin + 𝓂)#; + # timer::TimerOutput = TimerOutput()) + # @timeit_debug timer "Calculate jacobian - forward" begin jacobian = calculate_jacobian(parameters, SS_and_pars, 𝓂) - end # timeit_debug + # end # timeit_debug function calculate_jacobian_pullback(∂∇₁) - @timeit_debug timer "Calculate jacobian - reverse" begin + # @timeit_debug timer "Calculate jacobian - reverse" begin X = [parameters; SS_and_pars] # vals = Float64[] @@ -6148,7 +6148,7 @@ function rrule(::typeof(calculate_jacobian), # lk = ReentrantLock() - @timeit_debug timer "Loop" begin + # @timeit_debug timer "Loop" begin Polyester.@batch minbatch = 200 for f in 𝓂.model_jacobian_SS_and_pars_vars[1] out = f(X) @@ -6173,8 +6173,8 @@ function rrule(::typeof(calculate_jacobian), ∂parameters_and_SS_and_pars = analytical_jacobian_SS_and_pars_vars[:,cols_unique] * v∂∇₁ - end # timeit_debug - end # timeit_debug + # end # timeit_debug + # end # timeit_debug return NoTangent(), ∂parameters_and_SS_and_pars[1:length(parameters)], ∂parameters_and_SS_and_pars[length(parameters)+1:end], NoTangent() end @@ -6289,8 +6289,11 @@ function rrule(::typeof(calculate_hessian), parameters, SS_and_pars, 𝓂) end -function calculate_third_order_derivatives(parameters::Vector{M}, SS_and_pars::Vector{N}, 𝓂::ℳ; timer::TimerOutput = TimerOutput()) where {M,N} - @timeit_debug timer "3rd order derivatives" begin +function calculate_third_order_derivatives(parameters::Vector{M}, + SS_and_pars::Vector{N}, + 𝓂::ℳ) where {M,N} #; + # timer::TimerOutput = TimerOutput()) where {M,N} + # @timeit_debug timer "3rd order derivatives" begin SS = SS_and_pars[1:end - length(𝓂.calibration_equations)] calibrated_parameters = SS_and_pars[(end - length(𝓂.calibration_equations)+1):end] @@ -6327,7 +6330,7 @@ function calculate_third_order_derivatives(parameters::Vector{M}, SS_and_pars::V # lk = ReentrantLock() - @timeit_debug timer "Loop" begin + # @timeit_debug timer "Loop" begin Polyester.@batch minbatch = 200 for f in 𝓂.model_third_order_derivatives[1] out = f(X) @@ -6342,14 +6345,14 @@ function calculate_third_order_derivatives(parameters::Vector{M}, SS_and_pars::V # end end - end # timeit_debug + # end # timeit_debug - @timeit_debug timer "Allocation" begin + # @timeit_debug timer "Allocation" begin Accessors.@reset 𝓂.model_third_order_derivatives[2].nzval = vals - end # timeit_debug - end # timeit_debug + # end # timeit_debug + # end # timeit_debug return 𝓂.model_third_order_derivatives[2]# * 𝓂.solution.perturbation.third_order_auxilliary_matrices.𝐔∇₃ @@ -6371,18 +6374,19 @@ function calculate_third_order_derivatives(parameters::Vector{M}, SS_and_pars::V end -function rrule(::typeof(calculate_third_order_derivatives), parameters, SS_and_pars, 𝓂; timer::TimerOutput = TimerOutput()) - @timeit_debug timer "3rd order derivatives - forward" begin +function rrule(::typeof(calculate_third_order_derivatives), parameters, SS_and_pars, 𝓂) # ; + # timer::TimerOutput = TimerOutput()) + # @timeit_debug timer "3rd order derivatives - forward" begin third_order_derivatives = calculate_third_order_derivatives(parameters, SS_and_pars, 𝓂, timer = timer) - end # timeit_debug + # end # timeit_debug function calculate_third_order_derivatives_pullback(∂∇₁) - @timeit_debug timer "3rd order derivatives - pullback" begin + # @timeit_debug timer "3rd order derivatives - pullback" begin X = [parameters; SS_and_pars] vals = zeros(Float64, length(𝓂.model_third_order_derivatives_SS_and_pars_vars[1])) - @timeit_debug timer "Loop" begin + # @timeit_debug timer "Loop" begin Polyester.@batch minbatch = 200 for f in 𝓂.model_third_order_derivatives_SS_and_pars_vars[1] out = f(X) @@ -6390,13 +6394,13 @@ function rrule(::typeof(calculate_third_order_derivatives), parameters, SS_and_p @inbounds vals[out[2]] = out[1] end - end # timeit_debug - @timeit_debug timer "Allocation" begin + # end # timeit_debug + # @timeit_debug timer "Allocation" begin Accessors.@reset 𝓂.model_third_order_derivatives_SS_and_pars_vars[2].nzval = vals - end # timeit_debug - @timeit_debug timer "Post process" begin + # end # timeit_debug + # @timeit_debug timer "Post process" begin analytical_third_order_derivatives_SS_and_pars_vars = 𝓂.model_third_order_derivatives_SS_and_pars_vars[2] |> ThreadedSparseArrays.ThreadedSparseMatrixCSC @@ -6406,8 +6410,8 @@ function rrule(::typeof(calculate_third_order_derivatives), parameters, SS_and_p ∂parameters_and_SS_and_pars = analytical_third_order_derivatives_SS_and_pars_vars[:,cols_unique] * v∂∇₁ - end # timeit_debug - end # timeit_debug + # end # timeit_debug + # end # timeit_debug return NoTangent(), ∂parameters_and_SS_and_pars[1:length(parameters)], ∂parameters_and_SS_and_pars[length(parameters)+1:end], NoTangent() end @@ -6991,9 +6995,9 @@ end function get_NSSS_and_parameters(𝓂::ℳ, parameter_values::Vector{S}; verbose::Bool = false, - timer::TimerOutput = TimerOutput(), + # timer::TimerOutput = TimerOutput(), tol::AbstractFloat = 1e-12) where S <: Float64 - @timeit_debug timer "Calculate NSSS" begin + # @timeit_debug timer "Calculate NSSS" begin SS_and_pars, (solution_error, iters) = 𝓂.SS_solve_func(parameter_values, 𝓂, verbose, false, 𝓂.solver_parameters) if solution_error > tol || isnan(solution_error) @@ -7004,7 +7008,7 @@ function get_NSSS_and_parameters(𝓂::ℳ, return (SS_and_pars, (10, iters))#, x -> (NoTangent(), NoTangent(), NoTangent(), NoTangent()) end - end # timeit_debug + # end # timeit_debug return SS_and_pars, (solution_error, iters) end @@ -7013,19 +7017,19 @@ function rrule(::typeof(get_NSSS_and_parameters), 𝓂, parameter_values; verbose = false, - timer::TimerOutput = TimerOutput(), + # timer::TimerOutput = TimerOutput(), tol::AbstractFloat = 1e-12) - @timeit_debug timer "Calculate NSSS - forward" begin + # @timeit_debug timer "Calculate NSSS - forward" begin SS_and_pars, (solution_error, iters) = 𝓂.SS_solve_func(parameter_values, 𝓂, verbose, false, 𝓂.solver_parameters) - end # timeit_debug + # end # timeit_debug if solution_error > tol || isnan(solution_error) return (SS_and_pars, (solution_error, iters)), x -> (NoTangent(), NoTangent(), NoTangent(), NoTangent()) end - @timeit_debug timer "Calculate NSSS - pullback" begin + # @timeit_debug timer "Calculate NSSS - pullback" begin SS_and_pars_names_lead_lag = vcat(Symbol.(string.(sort(union(𝓂.var,𝓂.exo_past,𝓂.exo_future)))), 𝓂.calibration_equations_parameters) @@ -7050,7 +7054,7 @@ function rrule(::typeof(get_NSSS_and_parameters), # lk = ReentrantLock() - @timeit_debug timer "Loop - parameter derivatives" begin + # @timeit_debug timer "Loop - parameter derivatives" begin Polyester.@batch minbatch = 200 for f in 𝓂.∂SS_equations_∂parameters[1] out = f(X) @@ -7069,7 +7073,7 @@ function rrule(::typeof(get_NSSS_and_parameters), ∂SS_equations_∂parameters = 𝓂.∂SS_equations_∂parameters[2] - end # timeit_debug + # end # timeit_debug # vals = Float64[] @@ -7081,7 +7085,7 @@ function rrule(::typeof(get_NSSS_and_parameters), # lk = ReentrantLock() - @timeit_debug timer "Loop - NSSS derivatives" begin + # @timeit_debug timer "Loop - NSSS derivatives" begin Polyester.@batch minbatch = 200 for f in 𝓂.∂SS_equations_∂SS_and_pars[1] out = f(X) @@ -7101,12 +7105,12 @@ function rrule(::typeof(get_NSSS_and_parameters), ∂SS_equations_∂SS_and_pars = 𝓂.∂SS_equations_∂SS_and_pars[3] - end # timeit_debug + # end # timeit_debug # ∂SS_equations_∂parameters = 𝓂.∂SS_equations_∂parameters(parameter_values, SS_and_pars[indexin(unknowns, SS_and_pars_names_lead_lag)]) |> Matrix # ∂SS_equations_∂SS_and_pars = 𝓂.∂SS_equations_∂SS_and_pars(parameter_values, SS_and_pars[indexin(unknowns, SS_and_pars_names_lead_lag)]) |> Matrix - @timeit_debug timer "Implicit diff - mat inv" begin + # @timeit_debug timer "Implicit diff - mat inv" begin ∂SS_equations_∂SS_and_pars_lu = RF.lu!(∂SS_equations_∂SS_and_pars, check = false) @@ -7124,8 +7128,8 @@ function rrule(::typeof(get_NSSS_and_parameters), end end - end # timeit_debug - end # timeit_debug + # end # timeit_debug + # end # timeit_debug # try block-gmres here function get_non_stochastic_steady_state_pullback(∂SS_and_pars) @@ -7144,7 +7148,7 @@ end function get_NSSS_and_parameters(𝓂::ℳ, parameter_values_dual::Vector{ℱ.Dual{Z,S,N}}; verbose::Bool = false, - timer::TimerOutput = TimerOutput(), + # timer::TimerOutput = TimerOutput(), tol::AbstractFloat = 1e-12) where {Z,S,N} parameter_values = ℱ.value.(parameter_values_dual) @@ -7280,9 +7284,9 @@ function get_relevant_steady_state_and_state_update(::Val{:second_order}, 𝓂::ℳ, tol::AbstractFloat; quadratic_matrix_equation_solver::Symbol = :schur, - sylvester_algorithm::Symbol = :doubling, - verbose::Bool = false, - timer::TimerOutput = TimerOutput()) where S <: Real + sylvester_algorithm::Symbol = :doubling, + # timer::TimerOutput = TimerOutput(), + verbose::Bool = false) where S <: Real sss, converged, SS_and_pars, solution_error, ∇₁, ∇₂, 𝐒₁, 𝐒₂ = calculate_second_order_stochastic_steady_state(parameter_values, 𝓂, timer = timer, @@ -7311,9 +7315,9 @@ function get_relevant_steady_state_and_state_update(::Val{:pruned_second_order}, 𝓂::ℳ, tol::AbstractFloat; quadratic_matrix_equation_solver::Symbol = :schur, - sylvester_algorithm::Symbol = :doubling, - verbose::Bool = false, - timer::TimerOutput = TimerOutput())::Tuple{timings, Vector{S}, Union{Matrix{S},Vector{AbstractMatrix{S}}}, Vector{Vector{S}}, Bool} where S <: Real + sylvester_algorithm::Symbol = :doubling, + # timer::TimerOutput = TimerOutput(), + verbose::Bool = false)::Tuple{timings, Vector{S}, Union{Matrix{S},Vector{AbstractMatrix{S}}}, Vector{Vector{S}}, Bool} where S <: Real sss, converged, SS_and_pars, solution_error, ∇₁, ∇₂, 𝐒₁, 𝐒₂ = calculate_second_order_stochastic_steady_state(parameter_values, 𝓂, pruning = true, @@ -7344,9 +7348,9 @@ function get_relevant_steady_state_and_state_update(::Val{:third_order}, 𝓂::ℳ, tol::AbstractFloat; quadratic_matrix_equation_solver::Symbol = :schur, - sylvester_algorithm::Symbol = :bicgstab, - verbose::Bool = false, - timer::TimerOutput = TimerOutput()) where S <: Real + sylvester_algorithm::Symbol = :bicgstab, + # timer::TimerOutput = TimerOutput(), + verbose::Bool = false) where S <: Real sss, converged, SS_and_pars, solution_error, ∇₁, ∇₂, ∇₃, 𝐒₁, 𝐒₂, 𝐒₃ = calculate_third_order_stochastic_steady_state(parameter_values, 𝓂, timer = timer, @@ -7375,9 +7379,9 @@ function get_relevant_steady_state_and_state_update(::Val{:pruned_third_order}, 𝓂::ℳ, tol::AbstractFloat; quadratic_matrix_equation_solver::Symbol = :schur, - sylvester_algorithm::Symbol = :bicgstab, - verbose::Bool = false, - timer::TimerOutput = TimerOutput())::Tuple{timings, Vector{S}, Union{Matrix{S},Vector{AbstractMatrix{S}}}, Vector{Vector{S}}, Bool} where S <: Real + sylvester_algorithm::Symbol = :bicgstab, + # timer::TimerOutput = TimerOutput(), + verbose::Bool = false)::Tuple{timings, Vector{S}, Union{Matrix{S},Vector{AbstractMatrix{S}}}, Vector{Vector{S}}, Bool} where S <: Real sss, converged, SS_and_pars, solution_error, ∇₁, ∇₂, ∇₃, 𝐒₁, 𝐒₂, 𝐒₃ = calculate_third_order_stochastic_steady_state(parameter_values, 𝓂, pruning = true, @@ -7406,9 +7410,9 @@ function get_relevant_steady_state_and_state_update(::Val{:first_order}, 𝓂::ℳ, tol::AbstractFloat; quadratic_matrix_equation_solver::Symbol = :schur, - sylvester_algorithm::Symbol = :bicgstab, - verbose::Bool = false, - timer::TimerOutput = TimerOutput())::Tuple{timings, Vector{S}, Union{Matrix{S},Vector{AbstractMatrix{S}}}, Vector{Vector{Float64}}, Bool} where S <: Real + sylvester_algorithm::Symbol = :bicgstab, + # timer::TimerOutput = TimerOutput(), + verbose::Bool = false)::Tuple{timings, Vector{S}, Union{Matrix{S},Vector{AbstractMatrix{S}}}, Vector{Vector{Float64}}, Bool} where S <: Real SS_and_pars, (solution_error, iters) = get_NSSS_and_parameters(𝓂, parameter_values, tol = tol, timer = timer, verbose = verbose) state = zeros(𝓂.timings.nVars) diff --git a/src/algorithms/lyapunov.jl b/src/algorithms/lyapunov.jl index 09ab1225..6d88e60c 100644 --- a/src/algorithms/lyapunov.jl +++ b/src/algorithms/lyapunov.jl @@ -13,10 +13,10 @@ function solve_lyapunov_equation(A::AbstractMatrix{Float64}, C::AbstractMatrix{Float64}; lyapunov_algorithm::Symbol = :doubling, tol::AbstractFloat = 1e-12, - verbose::Bool = false, - timer::TimerOutput = TimerOutput()) - @timeit_debug timer "Solve lyapunov equation" begin - @timeit_debug timer "Choose matrix formats" begin + # timer::TimerOutput = TimerOutput(), + verbose::Bool = false) + # @timeit_debug timer "Solve lyapunov equation" begin + # @timeit_debug timer "Choose matrix formats" begin if lyapunov_algorithm ≠ :bartels_stewart A = choose_matrix_format(A) @@ -25,8 +25,8 @@ function solve_lyapunov_equation(A::AbstractMatrix{Float64}, C = choose_matrix_format(C, density_threshold = 0.0) # C = collect(C) # C is always dense because the output will be dense in all of these cases as we use this function to compute dense covariance matrices - end # timeit_debug - @timeit_debug timer "Solve" begin + # end # timeit_debug + # @timeit_debug timer "Solve" begin X, i, reached_tol = solve_lyapunov_equation(A, C, Val(lyapunov_algorithm), @@ -65,8 +65,8 @@ function solve_lyapunov_equation(A::AbstractMatrix{Float64}, end end - end # timeit_debug - end # timeit_debug + # end # timeit_debug + # end # timeit_debug # if (reached_tol > tol) println("Lyapunov failed: $reached_tol") end @@ -78,8 +78,8 @@ function rrule(::typeof(solve_lyapunov_equation), C::AbstractMatrix{Float64}; lyapunov_algorithm::Symbol = :doubling, tol::AbstractFloat = 1e-12, - verbose::Bool = false, - timer::TimerOutput = TimerOutput()) + # timer::TimerOutput = TimerOutput(), + verbose::Bool = false) P, solved = solve_lyapunov_equation(A, C, lyapunov_algorithm = lyapunov_algorithm, # tol = tol, @@ -106,8 +106,8 @@ function solve_lyapunov_equation( A::AbstractMatrix{ℱ.Dual{Z,S,N}}, C::AbstractMatrix{ℱ.Dual{Z,S,N}}; lyapunov_algorithm::Symbol = :doubling, tol::AbstractFloat = 1e-12, - verbose::Bool = false, - timer::TimerOutput = TimerOutput()) where {Z,S,N} + # timer::TimerOutput = TimerOutput(), + verbose::Bool = false) where {Z,S,N} # unpack: AoS -> SoA Â = ℱ.value.(A) Ĉ = ℱ.value.(C) @@ -147,8 +147,8 @@ end function solve_lyapunov_equation(A::Union{ℒ.Adjoint{Float64,Matrix{Float64}},DenseMatrix{Float64}}, C::Union{ℒ.Adjoint{Float64,Matrix{Float64}},DenseMatrix{Float64}}, ::Val{:bartels_stewart}; - tol::AbstractFloat = 1e-12, - timer::TimerOutput = TimerOutput()) + # timer::TimerOutput = TimerOutput(), + tol::AbstractFloat = 1e-12) 𝐂 = try MatrixEquations.lyapd(A, C) catch @@ -175,8 +175,8 @@ end function solve_lyapunov_equation( A::AbstractSparseMatrix{Float64}, C::AbstractSparseMatrix{Float64}, ::Val{:doubling}; - tol::Float64 = 1e-14, - timer::TimerOutput = TimerOutput()) + # timer::TimerOutput = TimerOutput(), + tol::Float64 = 1e-14) 𝐂 = copy(C) 𝐀 = copy(A) @@ -222,8 +222,8 @@ end function solve_lyapunov_equation( A::Union{ℒ.Adjoint{Float64,Matrix{Float64}},DenseMatrix{Float64}}, C::AbstractSparseMatrix{Float64}, ::Val{:doubling}; - tol::Float64 = 1e-14, - timer::TimerOutput = TimerOutput()) + # timer::TimerOutput = TimerOutput(), + tol::Float64 = 1e-14) 𝐂 = copy(C) 𝐀 = copy(A) @@ -272,8 +272,8 @@ end function solve_lyapunov_equation( A::AbstractSparseMatrix{Float64}, C::Union{ℒ.Adjoint{Float64,Matrix{Float64}},DenseMatrix{Float64}}, ::Val{:doubling}; - tol::Float64 = 1e-14, - timer::TimerOutput = TimerOutput()) + # timer::TimerOutput = TimerOutput(), + tol::Float64 = 1e-14) 𝐂 = copy(C) 𝐀 = copy(A) 𝐂A = collect(𝐀) @@ -334,8 +334,8 @@ end function solve_lyapunov_equation( A::Union{ℒ.Adjoint{Float64,Matrix{Float64}},DenseMatrix{Float64}}, C::Union{ℒ.Adjoint{Float64,Matrix{Float64}},DenseMatrix{Float64}}, ::Val{:doubling}; - tol::Float64 = 1e-14, - timer::TimerOutput = TimerOutput()) + # timer::TimerOutput = TimerOutput(), + tol::Float64 = 1e-14) 𝐂 = copy(C) 𝐂¹ = copy(C) 𝐀 = copy(A) @@ -391,8 +391,8 @@ end function solve_lyapunov_equation(A::AbstractMatrix{Float64}, C::Union{ℒ.Adjoint{Float64,Matrix{Float64}},DenseMatrix{Float64}}, ::Val{:bicgstab}; - tol::Float64 = 1e-14, - timer::TimerOutput = TimerOutput()) + # timer::TimerOutput = TimerOutput(), + tol::Float64 = 1e-14) tmp̄ = similar(C) 𝐗 = similar(C) @@ -431,8 +431,8 @@ end function solve_lyapunov_equation(A::AbstractMatrix{Float64}, C::Union{ℒ.Adjoint{Float64,Matrix{Float64}},DenseMatrix{Float64}}, ::Val{:gmres}; - tol::Float64 = 1e-14, - timer::TimerOutput = TimerOutput()) + # timer::TimerOutput = TimerOutput(), + tol::Float64 = 1e-14) tmp̄ = similar(C) 𝐗 = similar(C) diff --git a/src/algorithms/quadratic_matrix_equation.jl b/src/algorithms/quadratic_matrix_equation.jl index f7f3643d..4d3fd580 100644 --- a/src/algorithms/quadratic_matrix_equation.jl +++ b/src/algorithms/quadratic_matrix_equation.jl @@ -12,7 +12,7 @@ function solve_quadratic_matrix_equation(A::AbstractMatrix{R}, T::timings; initial_guess::AbstractMatrix{R} = zeros(0,0), quadratic_matrix_equation_solver::Symbol = :schur, - timer::TimerOutput = TimerOutput(), + # timer::TimerOutput = TimerOutput(), tol::AbstractFloat = 1e-8, # 1e-14 is too tight verbose::Bool = false) where R <: Real if length(initial_guess) > 0 @@ -81,9 +81,9 @@ function solve_quadratic_matrix_equation(A::AbstractMatrix{R}, T::timings; initial_guess::AbstractMatrix{R} = zeros(0,0), tol::AbstractFloat = 1e-14, - timer::TimerOutput = TimerOutput(), + # timer::TimerOutput = TimerOutput(), verbose::Bool = false) where R <: Real - @timeit_debug timer "Prepare indice" begin + # @timeit_debug timer "Prepare indice" begin comb = union(T.future_not_past_and_mixed_idx, T.past_not_future_idx) sort!(comb) @@ -92,8 +92,8 @@ function solve_quadratic_matrix_equation(A::AbstractMatrix{R}, past_not_future_and_mixed_in_comb = indexin(T.past_not_future_and_mixed_idx, comb) indices_past_not_future_in_comb = indexin(T.past_not_future_idx, comb) - end # timeit_debug - @timeit_debug timer "Assemble matrices" begin + # end # timeit_debug + # @timeit_debug timer "Assemble matrices" begin Ã₊ = A[:,future_not_past_and_mixed_in_comb] @@ -115,8 +115,8 @@ function solve_quadratic_matrix_equation(A::AbstractMatrix{R}, ℒ.rmul!(Ã₀₊,-1) E = vcat(hcat(Ã₋,Ã₀₊), hcat(Z₋, I₊)) - end # timeit_debug - @timeit_debug timer "Schur decomposition" begin + # end # timeit_debug + # @timeit_debug timer "Schur decomposition" begin # this is the companion form and by itself the linearisation of the matrix polynomial used in the linear time iteration method. see: https://opus4.kobv.de/opus4-matheon/files/209/240.pdf schdcmp = try @@ -128,8 +128,8 @@ function solve_quadratic_matrix_equation(A::AbstractMatrix{R}, eigenselect = abs.(schdcmp.β ./ schdcmp.α) .< 1 - end # timeit_debug - @timeit_debug timer "Reorder Schur decomposition" begin + # end # timeit_debug + # @timeit_debug timer "Reorder Schur decomposition" begin try ℒ.ordschur!(schdcmp, eigenselect) @@ -138,8 +138,8 @@ function solve_quadratic_matrix_equation(A::AbstractMatrix{R}, return A, false, 0, 1.0 end - end # timeit_debug - @timeit_debug timer "Postprocess" begin + # end # timeit_debug + # @timeit_debug timer "Postprocess" begin Z₂₁ = schdcmp.Z[T.nPast_not_future_and_mixed+1:end, 1:T.nPast_not_future_and_mixed] Z₁₁ = schdcmp.Z[1:T.nPast_not_future_and_mixed, 1:T.nPast_not_future_and_mixed] @@ -147,7 +147,7 @@ function solve_quadratic_matrix_equation(A::AbstractMatrix{R}, S₁₁ = schdcmp.S[1:T.nPast_not_future_and_mixed, 1:T.nPast_not_future_and_mixed] T₁₁ = schdcmp.T[1:T.nPast_not_future_and_mixed, 1:T.nPast_not_future_and_mixed] - @timeit_debug timer "Matrix inversions" begin + # @timeit_debug timer "Matrix inversions" begin Ẑ₁₁ = ℒ.lu(Z₁₁, check = false) @@ -163,8 +163,8 @@ function solve_quadratic_matrix_equation(A::AbstractMatrix{R}, return A, false, 0, 1.0 end - end # timeit_debug - @timeit_debug timer "Matrix divisions" begin + # end # timeit_debug + # @timeit_debug timer "Matrix divisions" begin # D = Z₂₁ / Ẑ₁₁ ℒ.rdiv!(Z₂₁, Ẑ₁₁) @@ -178,8 +178,8 @@ function solve_quadratic_matrix_equation(A::AbstractMatrix{R}, sol = vcat(L[T.not_mixed_in_past_idx,:], D) - end # timeit_debug - end # timeit_debug + # end # timeit_debug + # end # timeit_debug X = sol[T.dynamic_order,:] * ℒ.I(length(comb))[past_not_future_and_mixed_in_comb,:] @@ -210,13 +210,13 @@ function solve_quadratic_matrix_equation(A::AbstractMatrix{R}, T::timings; initial_guess::AbstractMatrix{R} = zeros(0,0), tol::AbstractFloat = 1e-14, - timer::TimerOutput = TimerOutput(), + # timer::TimerOutput = TimerOutput(), verbose::Bool = false, max_iter::Int = 100) where R <: Real # Johannes Huber, Alexander Meyer-Gohde, Johanna Saecker (2024). Solving Linear DSGE Models with Structure Preserving Doubling Methods. # https://www.imfs-frankfurt.de/forschung/imfs-working-papers/details.html?tx_mmpublications_publicationsdetail%5Bcontroller%5D=Publication&tx_mmpublications_publicationsdetail%5Bpublication%5D=461&cHash=f53244e0345a27419a9d40a3af98c02f # https://arxiv.org/abs/2212.09491 - @timeit_debug timer "Invert B" begin + # @timeit_debug timer "Invert B" begin guess_provided = true @@ -244,8 +244,8 @@ function solve_quadratic_matrix_equation(A::AbstractMatrix{R}, X = -E - initial_guess Y = -F - end # timeit_debug - @timeit_debug timer "Prellocate" begin + # end # timeit_debug + # @timeit_debug timer "Prellocate" begin X_new = similar(X) Y_new = similar(Y) @@ -265,19 +265,19 @@ function solve_quadratic_matrix_equation(A::AbstractMatrix{R}, solved = false iter = max_iter - end # timeit_debug - @timeit_debug timer "Loop" begin + # end # timeit_debug + # @timeit_debug timer "Loop" begin for i in 1:max_iter - @timeit_debug timer "Compute EI" begin + # @timeit_debug timer "Compute EI" begin # Compute EI = I - Y * X ℒ.mul!(temp1, Y, X) ℒ.axpby!(1, II, -1, temp1) # temp1 = II - Y * X - end # timeit_debug - @timeit_debug timer "Invert EI" begin + # end # timeit_debug + # @timeit_debug timer "Invert EI" begin fEI = ℒ.lu!(temp1, check = false) @@ -285,16 +285,16 @@ function solve_quadratic_matrix_equation(A::AbstractMatrix{R}, return A, false, iter, 1.0 end - end # timeit_debug - @timeit_debug timer "Compute E" begin + # end # timeit_debug + # @timeit_debug timer "Compute E" begin # Compute E = E * EI * E ℒ.ldiv!(temp3, fEI, E) ℒ.mul!(E_new, E, temp3) # E_new = E / fEI * E - end # timeit_debug - @timeit_debug timer "Compute FI" begin + # end # timeit_debug + # @timeit_debug timer "Compute FI" begin # Compute FI = I - X * Y ℒ.mul!(temp2, X, Y) @@ -303,8 +303,8 @@ function solve_quadratic_matrix_equation(A::AbstractMatrix{R}, # ℒ.mul!(E_new, X, Y, -1, 1) # temp1 .= II - X * Y - end # timeit_debug - @timeit_debug timer "Invert FI" begin + # end # timeit_debug + # @timeit_debug timer "Invert FI" begin fFI = ℒ.lu!(temp2, check = false) @@ -312,16 +312,16 @@ function solve_quadratic_matrix_equation(A::AbstractMatrix{R}, return A, false, iter, 1.0 end - end # timeit_debug - @timeit_debug timer "Compute F" begin + # end # timeit_debug + # @timeit_debug timer "Compute F" begin # Compute F = F * FI * F ℒ.ldiv!(temp3, fFI, F) ℒ.mul!(F_new, F, temp3) # F_new = F / fFI * F - end # timeit_debug - @timeit_debug timer "Compute X_new" begin + # end # timeit_debug + # @timeit_debug timer "Compute X_new" begin # Compute X_new = X + F * FI * X * E ℒ.mul!(temp3, X, E) @@ -335,8 +335,8 @@ function solve_quadratic_matrix_equation(A::AbstractMatrix{R}, ℒ.axpy!(1, X, X_new) # X_new += X - end # timeit_debug - @timeit_debug timer "Compute Y_new" begin + # end # timeit_debug + # @timeit_debug timer "Compute Y_new" begin # Compute Y_new = Y + E * EI * Y * F ℒ.mul!(X, Y, F) # use X as temporary storage @@ -359,17 +359,17 @@ function solve_quadratic_matrix_equation(A::AbstractMatrix{R}, break end - end # timeit_debug - @timeit_debug timer "Copy" begin + # end # timeit_debug + # @timeit_debug timer "Copy" begin # Update values for the next iteration copy!(X, X_new) copy!(Y, Y_new) copy!(E, E_new) copy!(F, F_new) - end # timeit_debug + # end # timeit_debug end - end # timeit_debug + # end # timeit_debug ℒ.axpy!(1, initial_guess, X_new) @@ -406,15 +406,15 @@ end # max_iter::Int = 5000) where R <: Real # # Iterate over: A * X * X + C + B * X = (A * X + B) * X + C = X + (A * X + B) \ C -# @timeit_debug timer "Preallocate" begin +# # @timeit_debug timer "Preallocate" begin # F = similar(C) # t = similar(C) # initial_guess = length(initial_guess) > 0 ? initial_guess : zero(A) -# end # timeit_debug -# @timeit_debug timer "Loop" begin +# # end # timeit_debug +# # @timeit_debug timer "Loop" begin # sol = @suppress begin # speedmapping(initial_guess; m! = (F̄, F) -> begin @@ -438,7 +438,7 @@ end # reached_tol = ℒ.norm(AXX) / AXXnorm -# end # timeit_debug +# # end # timeit_debug # # if reached_tol > tol # # println("QME: linear time iteration $reached_tol") @@ -510,7 +510,7 @@ function solve_quadratic_matrix_equation(A::AbstractMatrix{ℱ.Dual{Z,S,N}}, initial_guess::AbstractMatrix{<:Real} = zeros(0,0), tol::AbstractFloat = 1e-8, quadratic_matrix_equation_solver::Symbol = :schur, - timer::TimerOutput = TimerOutput(), + # timer::TimerOutput = TimerOutput(), verbose::Bool = false) where {Z,S,N} # unpack: AoS -> SoA Â = ℱ.value.(A) @@ -551,7 +551,7 @@ function solve_quadratic_matrix_equation(A::AbstractMatrix{ℱ.Dual{Z,S,N}}, if ℒ.norm(CC) < eps() continue end - dX, solved = solve_sylvester_equation(AA, -X, -CC, sylvester_algorithm = :sylvester) + dX, solved = solve_sylvester_equation(AA, -X, -CC, sylvester_algorithm = :doubling) X̃[:,i] = vec(dX) end diff --git a/src/algorithms/sylvester.jl b/src/algorithms/sylvester.jl index a385610b..9280b21c 100644 --- a/src/algorithms/sylvester.jl +++ b/src/algorithms/sylvester.jl @@ -17,9 +17,9 @@ function solve_sylvester_equation(A::M, initial_guess::AbstractMatrix{<:AbstractFloat} = zeros(0,0), sylvester_algorithm::Symbol = :doubling, tol::AbstractFloat = 1e-10, - timer::TimerOutput = TimerOutput(), + # timer::TimerOutput = TimerOutput(), verbose::Bool = false) where {M <: AbstractMatrix{Float64}, N <: AbstractMatrix{Float64}, O <: AbstractMatrix{Float64}} - @timeit_debug timer "Choose matrix formats" begin + # @timeit_debug timer "Choose matrix formats" begin if sylvester_algorithm == :bartels_stewart b = collect(B) @@ -37,8 +37,8 @@ function solve_sylvester_equation(A::M, c = choose_matrix_format(C)# |> sparse end - end # timeit_debug - @timeit_debug timer "Check if guess solves it already" begin + # end # timeit_debug + # @timeit_debug timer "Check if guess solves it already" begin if length(initial_guess) > 0 𝐂 = a * initial_guess * b + c - initial_guess @@ -54,8 +54,8 @@ function solve_sylvester_equation(A::M, end end - end # timeit_debug - @timeit_debug timer "Solve sylvester equation" begin + # end # timeit_debug + # @timeit_debug timer "Solve sylvester equation" begin x, i, reached_tol = solve_sylvester_equation(a, b, c, Val(sylvester_algorithm), initial_guess = initial_guess, @@ -200,7 +200,7 @@ function solve_sylvester_equation(A::M, # end # end - end # timeit_debug + # end # timeit_debug X = choose_matrix_format(x)# |> sparse @@ -217,7 +217,7 @@ function rrule(::typeof(solve_sylvester_equation), initial_guess::AbstractMatrix{<:AbstractFloat} = zeros(0,0), sylvester_algorithm::Symbol = :doubling, tol::AbstractFloat = 1e-12, - timer::TimerOutput = TimerOutput(), + # timer::TimerOutput = TimerOutput(), verbose::Bool = false) where {M <: AbstractMatrix{Float64}, N <: AbstractMatrix{Float64}, O <: AbstractMatrix{Float64}} P, solved = solve_sylvester_equation(A, B, C, @@ -248,7 +248,7 @@ function solve_sylvester_equation( A::AbstractMatrix{ℱ.Dual{Z,S,N}}, initial_guess::AbstractMatrix{<:AbstractFloat} = zeros(0,0), sylvester_algorithm::Symbol = :doubling, tol::AbstractFloat = 1e-12, - timer::TimerOutput = TimerOutput(), + # timer::TimerOutput = TimerOutput(), verbose::Bool = false) where {Z,S,N} # unpack: AoS -> SoA Â = ℱ.value.(A) @@ -293,7 +293,7 @@ function solve_sylvester_equation( A::AbstractSparseMatrix{Float64}, C::AbstractSparseMatrix{Float64}, ::Val{:doubling}; initial_guess::AbstractMatrix{<:AbstractFloat} = zeros(0,0), - timer::TimerOutput = TimerOutput(), + # timer::TimerOutput = TimerOutput(), verbose::Bool = false, tol::Float64 = 1e-14) # see doi:10.1016/j.aml.2009.01.012 @@ -360,7 +360,7 @@ function solve_sylvester_equation( A::AbstractSparseMatrix{Float64}, C::Matrix{Float64}, ::Val{:doubling}; initial_guess::AbstractMatrix{<:AbstractFloat} = zeros(0,0), - timer::TimerOutput = TimerOutput(), + # timer::TimerOutput = TimerOutput(), verbose::Bool = false, tol::Float64 = 1e-14) # see doi:10.1016/j.aml.2009.01.012 @@ -444,12 +444,12 @@ function solve_sylvester_equation( A::Matrix{Float64}, C::Matrix{Float64}, ::Val{:doubling}; initial_guess::AbstractMatrix{<:AbstractFloat} = zeros(0,0), - timer::TimerOutput = TimerOutput(), + # timer::TimerOutput = TimerOutput(), verbose::Bool = false, tol::Float64 = 1e-14) # see doi:10.1016/j.aml.2009.01.012 - @timeit_debug timer "Doubling solve" begin - @timeit_debug timer "Setup buffers" begin + # @timeit_debug timer "Doubling solve" begin + # @timeit_debug timer "Setup buffers" begin # guess_provided = true @@ -473,33 +473,33 @@ function solve_sylvester_equation( A::Matrix{Float64}, iters = max_iter - end # timeit_debug + # end # timeit_debug for i in 1:max_iter - @timeit_debug timer "Update C" begin + # @timeit_debug timer "Update C" begin ℒ.mul!(𝐂B, 𝐂, 𝐁) ℒ.mul!(𝐂¹, 𝐀, 𝐂B) ℒ.axpy!(1, 𝐂, 𝐂¹) # 𝐂¹ = 𝐀 * 𝐂 * 𝐁 + 𝐂 - end # timeit_debug + # end # timeit_debug - @timeit_debug timer "Square A" begin + # @timeit_debug timer "Square A" begin ℒ.mul!(𝐀¹,𝐀,𝐀) copy!(𝐀,𝐀¹) - end # timeit_debug + # end # timeit_debug # 𝐀 = 𝐀^2 - @timeit_debug timer "Square B" begin + # @timeit_debug timer "Square B" begin 𝐁 = 𝐁^2 # ℒ.mul!(𝐁¹,𝐁,𝐁) # copy!(𝐁,𝐁¹) - end # timeit_debug + # end # timeit_debug # droptol!(𝐀, eps()) - @timeit_debug timer "droptol B" begin + # @timeit_debug timer "droptol B" begin droptol!(𝐁, eps()) - end # timeit_debug + # end # timeit_debug if i % 2 == 0 normdiff = ℒ.norm(𝐂¹ - 𝐂) @@ -510,12 +510,12 @@ function solve_sylvester_equation( A::Matrix{Float64}, end end - @timeit_debug timer "Copy C" begin + # @timeit_debug timer "Copy C" begin copy!(𝐂,𝐂¹) - end # timeit_debug + # end # timeit_debug end - @timeit_debug timer "Finalise" begin + # @timeit_debug timer "Finalise" begin # ℒ.mul!(𝐂B, 𝐂, 𝐁) # ℒ.mul!(𝐂¹, 𝐀, 𝐂B) # ℒ.axpy!(1, 𝐂, 𝐂¹) @@ -531,8 +531,8 @@ function solve_sylvester_equation( A::Matrix{Float64}, reached_tol = ℒ.norm(A * 𝐂 * B + C - 𝐂) / max(ℒ.norm(𝐂), ℒ.norm(C)) - end # timeit_debug - end # timeit_debug + # end # timeit_debug + # end # timeit_debug # if reached_tol > tol # println("Sylvester: doubling $reached_tol") @@ -547,7 +547,7 @@ function solve_sylvester_equation( A::AbstractSparseMatrix{Float64}, C::Matrix{Float64}, ::Val{:doubling}; initial_guess::AbstractMatrix{<:AbstractFloat} = zeros(0,0), - timer::TimerOutput = TimerOutput(), + # timer::TimerOutput = TimerOutput(), verbose::Bool = false, tol::Float64 = 1e-14) # see doi:10.1016/j.aml.2009.01.012 On Smith-type iterative algorithms for the Stein matrix equation @@ -628,7 +628,7 @@ function solve_sylvester_equation( A::Matrix{Float64}, C::AbstractSparseMatrix{Float64}, ::Val{:doubling}; initial_guess::AbstractMatrix{<:AbstractFloat} = zeros(0,0), - timer::TimerOutput = TimerOutput(), + # timer::TimerOutput = TimerOutput(), verbose::Bool = false, tol::Float64 = 1e-14) # see doi:10.1016/j.aml.2009.01.012 @@ -710,7 +710,7 @@ function solve_sylvester_equation( A::AbstractSparseMatrix{Float64}, C::AbstractSparseMatrix{Float64}, ::Val{:doubling}; initial_guess::AbstractMatrix{<:AbstractFloat} = zeros(0,0), - timer::TimerOutput = TimerOutput(), + # timer::TimerOutput = TimerOutput(), verbose::Bool = false, tol::Float64 = 1e-14) # see doi:10.1016/j.aml.2009.01.012 @@ -791,7 +791,7 @@ function solve_sylvester_equation( A::Matrix{Float64}, C::AbstractSparseMatrix{Float64}, ::Val{:doubling}; initial_guess::AbstractMatrix{<:AbstractFloat} = zeros(0,0), - timer::TimerOutput = TimerOutput(), + # timer::TimerOutput = TimerOutput(), verbose::Bool = false, tol::Float64 = 1e-14) # see doi:10.1016/j.aml.2009.01.012 @@ -872,11 +872,11 @@ function solve_sylvester_equation( A::Union{ℒ.Adjoint{Float64,Matrix{Float64} C::Union{ℒ.Adjoint{Float64,Matrix{Float64}},DenseMatrix{Float64}}, ::Val{:doubling}; initial_guess::AbstractMatrix{<:AbstractFloat} = zeros(0,0), - timer::TimerOutput = TimerOutput(), + # timer::TimerOutput = TimerOutput(), verbose::Bool = false, tol::Float64 = 1e-14) # see doi:10.1016/j.aml.2009.01.012 - @timeit_debug timer "Setup buffers" begin + # @timeit_debug timer "Setup buffers" begin # guess_provided = true @@ -900,24 +900,24 @@ function solve_sylvester_equation( A::Union{ℒ.Adjoint{Float64,Matrix{Float64} iters = max_iter - end # timeit_debug + # end # timeit_debug for i in 1:max_iter - @timeit_debug timer "Update C" begin + # @timeit_debug timer "Update C" begin ℒ.mul!(𝐂B, 𝐂, 𝐁) ℒ.mul!(𝐂¹, 𝐀, 𝐂B) ℒ.axpy!(1, 𝐂, 𝐂¹) # 𝐂¹ = 𝐀 * 𝐂 * 𝐁 + 𝐂 - end # timeit_debug + # end # timeit_debug - @timeit_debug timer "Square A" begin + # @timeit_debug timer "Square A" begin ℒ.mul!(𝐀¹,𝐀,𝐀) copy!(𝐀,𝐀¹) - end # timeit_debug - @timeit_debug timer "Square B" begin + # end # timeit_debug + # @timeit_debug timer "Square B" begin ℒ.mul!(𝐁¹,𝐁,𝐁) copy!(𝐁,𝐁¹) - end # timeit_debug + # end # timeit_debug # 𝐀 = 𝐀^2 # 𝐁 = 𝐁^2 @@ -933,12 +933,12 @@ function solve_sylvester_equation( A::Union{ℒ.Adjoint{Float64,Matrix{Float64} end end - @timeit_debug timer "Copy C" begin + # @timeit_debug timer "Copy C" begin copy!(𝐂,𝐂¹) - end # timeit_debug + # end # timeit_debug end - @timeit_debug timer "Finalise" begin + # @timeit_debug timer "Finalise" begin # ℒ.mul!(𝐂B, 𝐂, 𝐁) # ℒ.mul!(𝐂¹, 𝐀, 𝐂B) @@ -955,7 +955,7 @@ function solve_sylvester_equation( A::Union{ℒ.Adjoint{Float64,Matrix{Float64} reached_tol = ℒ.norm(A * 𝐂 * B + C - 𝐂) / max(ℒ.norm(𝐂), ℒ.norm(C)) - end # timeit_debug + # end # timeit_debug # if reached_tol > tol # println("Sylvester: doubling $reached_tol") @@ -970,7 +970,7 @@ function solve_sylvester_equation(A::DenseMatrix{Float64}, C::DenseMatrix{Float64}, ::Val{:bartels_stewart}; initial_guess::AbstractMatrix{<:AbstractFloat} = zeros(0,0), - timer::TimerOutput = TimerOutput(), + # timer::TimerOutput = TimerOutput(), verbose::Bool = false, tol::AbstractFloat = 1e-12) # guess_provided = true @@ -1011,10 +1011,10 @@ function solve_sylvester_equation(A::DenseMatrix{Float64}, C::DenseMatrix{Float64}, ::Val{:bicgstab}; initial_guess::AbstractMatrix{<:AbstractFloat} = zeros(0,0), - timer::TimerOutput = TimerOutput(), + # timer::TimerOutput = TimerOutput(), verbose::Bool = false, tol::Float64 = 1e-14) - @timeit_debug timer "Preallocate matrices" begin + # @timeit_debug timer "Preallocate matrices" begin # guess_provided = true @@ -1028,28 +1028,28 @@ function solve_sylvester_equation(A::DenseMatrix{Float64}, tmp̄ = zero(C) 𝐗 = zero(C) - end # timeit_debug + # end # timeit_debug # idxs = findall(abs.(C) .> eps()) # this does not work since the problem is incomplete function sylvester!(sol,𝐱) - @timeit_debug timer "Copy1" begin + # @timeit_debug timer "Copy1" begin # @inbounds 𝐗[idxs] = 𝐱 copyto!(𝐗, 𝐱) - end # timeit_debug + # end # timeit_debug # 𝐗 = @view reshape(𝐱, size(𝐗)) - @timeit_debug timer "Mul1" begin + # @timeit_debug timer "Mul1" begin # tmp̄ = A * 𝐗 * B ℒ.mul!(tmp̄, A, 𝐗) - end # timeit_debug - @timeit_debug timer "Mul2" begin + # end # timeit_debug + # @timeit_debug timer "Mul2" begin ℒ.mul!(𝐗, tmp̄, B, -1, 1) # ℒ.axpby!(-1, tmp̄, 1, 𝐗) - end # timeit_debug - @timeit_debug timer "Copy2" begin + # end # timeit_debug + # @timeit_debug timer "Copy2" begin # @inbounds @views copyto!(sol, 𝐗[idxs]) copyto!(sol, 𝐗) - end # timeit_debug + # end # timeit_debug # sol = @view reshape(𝐗, size(sol)) end @@ -1081,7 +1081,7 @@ function solve_sylvester_equation(A::DenseMatrix{Float64}, # precond = LinearOperators.LinearOperator(Float64, length(C), length(C), true, true, preconditioner!) - @timeit_debug timer "BICGSTAB solve" begin + # @timeit_debug timer "BICGSTAB solve" begin # if length(init) == 0 # 𝐂, info = Krylov.bicgstab(sylvester, C[idxs], rtol = tol / 10, atol = tol / 10)#, M = precond) 𝐂, info = Krylov.bicgstab(sylvester, [vec(𝐂¹);], @@ -1093,9 +1093,9 @@ function solve_sylvester_equation(A::DenseMatrix{Float64}, # else # 𝐂, info = Krylov.bicgstab(sylvester, [vec(C);], [vec(init);], rtol = tol / 10) # end - end # timeit_debug + # end # timeit_debug - @timeit_debug timer "Postprocess" begin + # @timeit_debug timer "Postprocess" begin # @inbounds 𝐗[idxs] = 𝐂 copyto!(𝐗, 𝐂) @@ -1112,7 +1112,7 @@ function solve_sylvester_equation(A::DenseMatrix{Float64}, reached_tol = ℒ.norm(A * 𝐗 * B + C - 𝐗) / max(ℒ.norm(𝐗), ℒ.norm(C)) - end # timeit_debug + # end # timeit_debug if !(typeof(C) <: DenseMatrix) 𝐗 = choose_matrix_format(𝐗, density_threshold = 1.0) @@ -1131,10 +1131,10 @@ function solve_sylvester_equation(A::DenseMatrix{Float64}, C::DenseMatrix{Float64}, ::Val{:dqgmres}; initial_guess::AbstractMatrix{<:AbstractFloat} = zeros(0,0), - timer::TimerOutput = TimerOutput(), + # timer::TimerOutput = TimerOutput(), verbose::Bool = false, tol::Float64 = 1e-14) - @timeit_debug timer "Preallocate matrices" begin + # @timeit_debug timer "Preallocate matrices" begin # guess_provided = true @@ -1148,24 +1148,24 @@ function solve_sylvester_equation(A::DenseMatrix{Float64}, tmp̄ = similar(C) 𝐗 = similar(C) - end # timeit_debug + # end # timeit_debug function sylvester!(sol,𝐱) - @timeit_debug timer "Copy1" begin + # @timeit_debug timer "Copy1" begin copyto!(𝐗, 𝐱) - end # timeit_debug + # end # timeit_debug # 𝐗 = @view reshape(𝐱, size(𝐗)) - @timeit_debug timer "Mul1" begin + # @timeit_debug timer "Mul1" begin # tmp̄ = A * 𝐗 * B ℒ.mul!(tmp̄, A, 𝐗) - end # timeit_debug - @timeit_debug timer "Mul2" begin + # end # timeit_debug + # @timeit_debug timer "Mul2" begin ℒ.mul!(𝐗, tmp̄, B, -1, 1) # ℒ.axpby!(-1, tmp̄, 1, 𝐗) - end # timeit_debug - @timeit_debug timer "Copy2" begin + # end # timeit_debug + # @timeit_debug timer "Copy2" begin copyto!(sol, 𝐗) - end # timeit_debug + # end # timeit_debug # sol = @view reshape(𝐗, size(sol)) end @@ -1192,7 +1192,7 @@ function solve_sylvester_equation(A::DenseMatrix{Float64}, # precond = LinearOperators.LinearOperator(Float64, length(C), length(C), true, true, preconditioner!) - @timeit_debug timer "DQGMRES solve" begin + # @timeit_debug timer "DQGMRES solve" begin # if length(init) == 0 𝐂, info = Krylov.dqgmres(sylvester, [vec(𝐂¹);], # [vec(initial_guess);], # start value helps @@ -1203,9 +1203,9 @@ function solve_sylvester_equation(A::DenseMatrix{Float64}, # else # 𝐂, info = Krylov.gmres(sylvester, [vec(C);], [vec(init);], rtol = tol / 10)#, restart = true, M = precond) # end - end # timeit_debug + # end # timeit_debug - @timeit_debug timer "Postprocess" begin + # @timeit_debug timer "Postprocess" begin copyto!(𝐗, 𝐂) @@ -1222,7 +1222,7 @@ function solve_sylvester_equation(A::DenseMatrix{Float64}, reached_tol = ℒ.norm(A * 𝐗 * B + C - 𝐗) / max(ℒ.norm(𝐗), ℒ.norm(C)) - end # timeit_debug + # end # timeit_debug # if reached_tol > tol # println("Sylvester: dqgmres $reached_tol") @@ -1237,10 +1237,10 @@ function solve_sylvester_equation(A::DenseMatrix{Float64}, C::DenseMatrix{Float64}, ::Val{:gmres}; initial_guess::AbstractMatrix{<:AbstractFloat} = zeros(0,0), - timer::TimerOutput = TimerOutput(), + # timer::TimerOutput = TimerOutput(), verbose::Bool = false, tol::Float64 = 1e-14) - @timeit_debug timer "Preallocate matrices" begin + # @timeit_debug timer "Preallocate matrices" begin # guess_provided = true @@ -1254,24 +1254,24 @@ function solve_sylvester_equation(A::DenseMatrix{Float64}, tmp̄ = similar(C) 𝐗 = similar(C) - end # timeit_debug + # end # timeit_debug function sylvester!(sol,𝐱) - @timeit_debug timer "Copy1" begin + # @timeit_debug timer "Copy1" begin copyto!(𝐗, 𝐱) - end # timeit_debug + # end # timeit_debug # 𝐗 = @view reshape(𝐱, size(𝐗)) - @timeit_debug timer "Mul1" begin + # @timeit_debug timer "Mul1" begin # tmp̄ = A * 𝐗 * B ℒ.mul!(tmp̄, A, 𝐗) - end # timeit_debug - @timeit_debug timer "Mul2" begin + # end # timeit_debug + # @timeit_debug timer "Mul2" begin ℒ.mul!(𝐗, tmp̄, B, -1, 1) # ℒ.axpby!(-1, tmp̄, 1, 𝐗) - end # timeit_debug - @timeit_debug timer "Copy2" begin + # end # timeit_debug + # @timeit_debug timer "Copy2" begin copyto!(sol, 𝐗) - end # timeit_debug + # end # timeit_debug # sol = @view reshape(𝐗, size(sol)) end @@ -1298,7 +1298,7 @@ function solve_sylvester_equation(A::DenseMatrix{Float64}, # precond = LinearOperators.LinearOperator(Float64, length(C), length(C), true, true, preconditioner!) - @timeit_debug timer "GMRES solve" begin + # @timeit_debug timer "GMRES solve" begin # if length(init) == 0 𝐂, info = Krylov.gmres(sylvester, [vec(𝐂¹);], # [vec(initial_guess);], # start value helps @@ -1309,9 +1309,9 @@ function solve_sylvester_equation(A::DenseMatrix{Float64}, # else # 𝐂, info = Krylov.gmres(sylvester, [vec(C);], [vec(init);], rtol = tol / 10)#, restart = true, M = precond) # end - end # timeit_debug + # end # timeit_debug - @timeit_debug timer "Postprocess" begin + # @timeit_debug timer "Postprocess" begin copyto!(𝐗, 𝐂) @@ -1328,7 +1328,7 @@ function solve_sylvester_equation(A::DenseMatrix{Float64}, reached_tol = ℒ.norm(A * 𝐗 * B + C - 𝐗) / max(ℒ.norm(𝐗), ℒ.norm(C)) - end # timeit_debug + # end # timeit_debug # if reached_tol > tol # println("Sylvester: gmres $reached_tol") @@ -1365,7 +1365,7 @@ end # iters = max_iter # for i in 1:max_iter -# @timeit_debug timer "Update" begin +# # @timeit_debug timer "Update" begin # ℒ.mul!(𝐂B, 𝐂, B) # ℒ.mul!(𝐂¹, A, 𝐂B) # ℒ.axpy!(1, 𝐂⁰, 𝐂¹) @@ -1380,7 +1380,7 @@ end # end # copyto!(𝐂, 𝐂¹) -# end # timeit_debug +# # end # timeit_debug # end # # ℒ.mul!(𝐂B, 𝐂, B) diff --git a/src/filter/inversion.jl b/src/filter/inversion.jl index a5ef95d3..61e3e830 100644 --- a/src/filter/inversion.jl +++ b/src/filter/inversion.jl @@ -9,8 +9,8 @@ function calculate_loglikelihood(::Val{:inversion}, state, warmup_iterations, filter_algorithm, - verbose; - timer::TimerOutput = TimerOutput()) + verbose) #; + # timer::TimerOutput = TimerOutput()) return calculate_inversion_filter_loglikelihood(Val(algorithm), state, 𝐒, @@ -31,12 +31,12 @@ function calculate_inversion_filter_loglikelihood(::Val{:first_order}, data_in_deviations::Matrix{Float64}, observables::Union{Vector{String}, Vector{Symbol}}, T::timings; - timer::TimerOutput = TimerOutput(), + # timer::TimerOutput = TimerOutput(), warmup_iterations::Int = 0, presample_periods::Int = 0, verbose::Bool = false, filter_algorithm::Symbol = :LagrangeNewton) - @timeit_debug timer "Inversion filter" begin + # @timeit_debug timer "Inversion filter" begin # first order state = copy(state[1]) @@ -121,7 +121,7 @@ function calculate_inversion_filter_loglikelihood(::Val{:first_order}, 𝐒obs = 𝐒[cond_var_idx,1:end-T.nExo] - @timeit_debug timer "Loop" begin + # @timeit_debug timer "Loop" begin for i in axes(data_in_deviations,2) @views ℒ.mul!(y, 𝐒obs, state[T.past_not_future_and_mixed_idx]) @views ℒ.axpby!(1, data_in_deviations[:,i], -1, y) @@ -139,8 +139,8 @@ function calculate_inversion_filter_loglikelihood(::Val{:first_order}, end # TODO: use subset of observables and states when propagating states (see kalman filter) - end # timeit_debug - end # timeit_debug + # end # timeit_debug + # end # timeit_debug return -(logabsdets + shocks² + (length(observables) * (warmup_iterations + n_obs - presample_periods)) * log(2 * 3.141592653589793)) / 2 # return -(logabsdets + (length(observables) * (warmup_iterations + n_obs - presample_periods)) * log(2 * 3.141592653589793)) / 2 @@ -155,12 +155,12 @@ function rrule(::typeof(calculate_inversion_filter_loglikelihood), data_in_deviations::Matrix{Float64}, observables::Union{Vector{String}, Vector{Symbol}}, T::timings; - timer::TimerOutput = TimerOutput(), + # timer::TimerOutput = TimerOutput(), warmup_iterations::Int = 0, presample_periods::Int = 0, verbose::Bool = false, filter_algorithm::Symbol = :LagrangeNewton) - @timeit_debug timer "Inversion filter - forward" begin + # @timeit_debug timer "Inversion filter - forward" begin # first order state = copy(state[1]) @@ -266,11 +266,11 @@ function rrule(::typeof(calculate_inversion_filter_loglikelihood), ∂𝐒t⁻ = copy(tmp2) # ∂𝐒obs_idx = copy(tmp1) - end # timeit_debug + # end # timeit_debug # TODO: optimize allocations # pullback function inversion_pullback(∂llh) - @timeit_debug timer "Inversion filter - pullback" begin + # @timeit_debug timer "Inversion filter - pullback" begin for t in reverse(axes(data_in_deviations,2)) ∂state[t⁻] .= M² * ∂state[t⁻] @@ -311,7 +311,7 @@ function rrule(::typeof(calculate_inversion_filter_loglikelihood), ∂𝐒[obs_idx,end-T.nExo+1:end] -= (size(data_in_deviations,2) - presample_periods) * invjac' / 2 - end # timeit_debug + # end # timeit_debug return NoTangent(), NoTangent(), [∂state * ∂llh], ∂𝐒 * ∂llh, ∂data_in_deviations * ∂llh, NoTangent(), NoTangent(), NoTangent(), NoTangent() end @@ -328,13 +328,13 @@ function calculate_inversion_filter_loglikelihood(::Val{:pruned_second_order}, data_in_deviations::Matrix{Float64}, observables::Union{Vector{String}, Vector{Symbol}}, T::timings; - timer::TimerOutput = TimerOutput(), + # timer::TimerOutput = TimerOutput(), warmup_iterations::Int = 0, presample_periods::Int = 0, verbose::Bool = false, filter_algorithm::Symbol = :LagrangeNewton) - @timeit_debug timer "Pruned 2nd - Inversion filter" begin - @timeit_debug timer "Preallocation" begin + # @timeit_debug timer "Pruned 2nd - Inversion filter" begin + # @timeit_debug timer "Preallocation" begin precision_factor = 1.0 @@ -410,8 +410,8 @@ function calculate_inversion_filter_loglikelihood(::Val{:pruned_second_order}, init_guess = zeros(size(𝐒ⁱ, 2)) - end # timeit_debug - @timeit_debug timer "Loop" begin + # end # timeit_debug + # @timeit_debug timer "Loop" begin for i in axes(data_in_deviations, 2) # state¹⁻ = state₁ @@ -441,7 +441,7 @@ function calculate_inversion_filter_loglikelihood(::Val{:pruned_second_order}, init_guess *= 0 - @timeit_debug timer "Find shocks" begin + # @timeit_debug timer "Find shocks" begin x, matched = find_shocks(Val(filter_algorithm), init_guess, kron_buffer, @@ -452,7 +452,7 @@ function calculate_inversion_filter_loglikelihood(::Val{:pruned_second_order}, shock_independent, # max_iter = 100 ) - end # timeit_debug + # end # timeit_debug # if matched println("$filter_algorithm: $matched; current x: $x") end # if !matched @@ -546,8 +546,8 @@ function calculate_inversion_filter_loglikelihood(::Val{:pruned_second_order}, ℒ.mul!(state₂, 𝐒⁻², kronaug_state₁, 1/2, 1) end - end # timeit_debug - end # timeit_debug + # end # timeit_debug + # end # timeit_debug # See: https://pcubaborda.net/documents/CGIZ-final.pdf and Fair and Taylor (1983) return -(logabsdets + shocks² + (length(observables) * (warmup_iterations + n_obs - presample_periods)) * log(2 * 3.141592653589793)) / 2 @@ -562,13 +562,13 @@ function rrule(::typeof(calculate_inversion_filter_loglikelihood), data_in_deviations::Matrix{Float64}, observables::Union{Vector{String}, Vector{Symbol}}, T::timings; - timer::TimerOutput = TimerOutput(), + # timer::TimerOutput = TimerOutput(), warmup_iterations::Int = 0, presample_periods::Int = 0, verbose::Bool = false, filter_algorithm::Symbol = :LagrangeNewton)# where S <: Real - @timeit_debug timer "Inversion filter pruned 2nd - forward" begin - @timeit_debug timer "Preallocation" begin + # @timeit_debug timer "Inversion filter pruned 2nd - forward" begin + # @timeit_debug timer "Preallocation" begin precision_factor = 1.0 @@ -672,8 +672,8 @@ function rrule(::typeof(calculate_inversion_filter_loglikelihood), lI = -2 * vec(ℒ.I(size(𝐒ⁱ, 2))) - end # timeit_debug - @timeit_debug timer "Main loop" begin + # end # timeit_debug + # @timeit_debug timer "Main loop" begin for i in axes(data_in_deviations,2) # state¹⁻ = state₁ @@ -703,7 +703,7 @@ function rrule(::typeof(calculate_inversion_filter_loglikelihood), init_guess *= 0 - @timeit_debug timer "Find shocks" begin + # @timeit_debug timer "Find shocks" begin x[i], matched = find_shocks(Val(filter_algorithm), init_guess, kronxx[i], @@ -714,7 +714,7 @@ function rrule(::typeof(calculate_inversion_filter_loglikelihood), shock_independent, # max_iter = 100 ) - end # timeit_debug + # end # timeit_debug if !matched if verbose println("Inversion filter failed at step $i") end @@ -779,8 +779,8 @@ function rrule(::typeof(calculate_inversion_filter_loglikelihood), ℒ.mul!(state₂, 𝐒⁻², kronaug_state₁, 1/2, 1) end - end # timeit_debug - end # timeit_debug + # end # timeit_debug + # end # timeit_debug ∂data_in_deviations = similar(data_in_deviations) @@ -797,8 +797,8 @@ function rrule(::typeof(calculate_inversion_filter_loglikelihood), ∂kronstate¹⁻_vol = zero(ℒ.kron(state¹⁻_vol, state¹⁻_vol)) function inversion_filter_loglikelihood_pullback(∂llh) - @timeit_debug timer "Inversion filter pruned 2nd - pullback" begin - @timeit_debug timer "Preallocation" begin + # @timeit_debug timer "Inversion filter pruned 2nd - pullback" begin + # @timeit_debug timer "Preallocation" begin ∂𝐒ⁱ = zero(𝐒ⁱ) ∂𝐒ⁱ²ᵉ = zero(𝐒ⁱ²ᵉ) @@ -821,8 +821,8 @@ function rrule(::typeof(calculate_inversion_filter_loglikelihood), kronSλ = zeros(length(cond_var_idx) * T.nExo) kronxS = zeros(T.nExo * length(cond_var_idx)) - end # timeit_debug - @timeit_debug timer "Main loop" begin + # end # timeit_debug + # @timeit_debug timer "Main loop" begin for i in reverse(axes(data_in_deviations,2)) # state₁, state₂ = [𝐒⁻¹ * aug_state₁[i], 𝐒⁻¹ * aug_state₂[i] + 𝐒⁻² * ℒ.kron(aug_state₁[i], aug_state₁[i]) / 2] @@ -973,8 +973,8 @@ function rrule(::typeof(calculate_inversion_filter_loglikelihood), ℒ.axpy!(1, ∂state¹⁻_vol[1:end-1], ∂state[1]) end - end # timeit_debug - @timeit_debug timer "Post allocation" begin + # end # timeit_debug + # @timeit_debug timer "Post allocation" begin ∂𝐒 = [zero(𝐒[1]), zeros(size(𝐒[2]))] @@ -1001,8 +1001,8 @@ function rrule(::typeof(calculate_inversion_filter_loglikelihood), ∂state[1] = ℒ.I(T.nVars)[:,T.past_not_future_and_mixed_idx] * ∂state[1] * ∂llh ∂state[2] = ℒ.I(T.nVars)[:,T.past_not_future_and_mixed_idx] * ∂state[2] * ∂llh - end # timeit_debug - end # timeit_debug + # end # timeit_debug + # end # timeit_debug return NoTangent(), NoTangent(), ∂state, ∂𝐒, ∂data_in_deviations, NoTangent(), NoTangent(), NoTangent(), NoTangent(), NoTangent() end @@ -1022,13 +1022,13 @@ function calculate_inversion_filter_loglikelihood(::Val{:second_order}, data_in_deviations::Matrix{Float64}, observables::Union{Vector{String}, Vector{Symbol}}, T::timings; - timer::TimerOutput = TimerOutput(), + # timer::TimerOutput = TimerOutput(), warmup_iterations::Int = 0, presample_periods::Int = 0, verbose::Bool = false, filter_algorithm::Symbol = :LagrangeNewton)# where S <: Real - @timeit_debug timer "2nd - Inversion filter" begin - @timeit_debug timer "Preallocation" begin + # @timeit_debug timer "2nd - Inversion filter" begin + # @timeit_debug timer "Preallocation" begin precision_factor = 1.0 @@ -1102,8 +1102,8 @@ function calculate_inversion_filter_loglikelihood(::Val{:second_order}, init_guess = zeros(size(𝐒ⁱ, 2)) - end # timeit_debug - @timeit_debug timer "Loop" begin + # end # timeit_debug + # @timeit_debug timer "Loop" begin for i in axes(data_in_deviations,2) # state¹⁻ = state#[T.past_not_future_and_mixed_idx] @@ -1128,7 +1128,7 @@ function calculate_inversion_filter_loglikelihood(::Val{:second_order}, init_guess *= 0 - @timeit_debug timer "Find shocks" begin + # @timeit_debug timer "Find shocks" begin x, matched = find_shocks(Val(filter_algorithm), init_guess, kron_buffer, @@ -1139,7 +1139,7 @@ function calculate_inversion_filter_loglikelihood(::Val{:second_order}, shock_independent, # max_iter = 100 ) - end # timeit_debug + # end # timeit_debug # if !matched # x, matched = find_shocks(Val(:COBYLA), @@ -1232,8 +1232,8 @@ function calculate_inversion_filter_loglikelihood(::Val{:second_order}, ℒ.mul!(state, 𝐒⁻², kronaug_state, 1/2 ,1) end - end # timeit_debug - end # timeit_debug + # end # timeit_debug + # end # timeit_debug # See: https://pcubaborda.net/documents/CGIZ-final.pdf return -(logabsdets + shocks² + (length(observables) * (warmup_iterations + n_obs - presample_periods)) * log(2 * 3.141592653589793)) / 2 @@ -1248,14 +1248,14 @@ function rrule(::typeof(calculate_inversion_filter_loglikelihood), data_in_deviations::Matrix{Float64}, observables::Union{Vector{String}, Vector{Symbol}}, T::timings; - timer::TimerOutput = TimerOutput(), + # timer::TimerOutput = TimerOutput(), warmup_iterations::Int = 0, presample_periods::Int = 0, verbose::Bool = false, filter_algorithm::Symbol = :LagrangeNewton)# where S <: Real - @timeit_debug timer "Inversion filter 2nd - forward" begin + # @timeit_debug timer "Inversion filter 2nd - forward" begin - @timeit_debug timer "Preallocation" begin + # @timeit_debug timer "Preallocation" begin precision_factor = 1.0 @@ -1357,8 +1357,8 @@ function rrule(::typeof(calculate_inversion_filter_loglikelihood), init_guess = zeros(size(𝐒ⁱ, 2)) - end # timeit_debug - @timeit_debug timer "Main loop" begin + # end # timeit_debug + # @timeit_debug timer "Main loop" begin @inbounds for i in axes(data_in_deviations,2) # aug_state[i][1:T.nPast_not_future_and_mixed] = state¹⁻ @@ -1384,7 +1384,7 @@ function rrule(::typeof(calculate_inversion_filter_loglikelihood), init_guess *= 0 - @timeit_debug timer "Find shocks" begin + # @timeit_debug timer "Find shocks" begin x[i], matched = find_shocks(Val(filter_algorithm), init_guess, kronxx[i], @@ -1395,7 +1395,7 @@ function rrule(::typeof(calculate_inversion_filter_loglikelihood), shock_independent, # max_iter = 100 ) - end # timeit_debug + # end # timeit_debug if !matched if verbose println("Inversion filter failed at step $i") end @@ -1459,8 +1459,8 @@ function rrule(::typeof(calculate_inversion_filter_loglikelihood), ℒ.mul!(state¹⁻, 𝐒⁻², kronaug_state[i], 1/2 ,1) end - end # timeit_debug - end # timeit_debug + # end # timeit_debug + # end # timeit_debug ∂aug_state = zero(aug_state[1]) @@ -1477,9 +1477,9 @@ function rrule(::typeof(calculate_inversion_filter_loglikelihood), ∂kronIx = zero(ℒ.kron(ℒ.I(length(x[1])), x[1])) function inversion_filter_loglikelihood_pullback(∂llh) - @timeit_debug timer "Inversion filter 2nd - pullback" begin + # @timeit_debug timer "Inversion filter 2nd - pullback" begin - @timeit_debug timer "Preallocation" begin + # @timeit_debug timer "Preallocation" begin ∂𝐒ⁱ = zero(𝐒ⁱ) ∂𝐒ⁱ²ᵉ = zero(𝐒ⁱ²ᵉ) @@ -1504,8 +1504,8 @@ function rrule(::typeof(calculate_inversion_filter_loglikelihood), kronSλ = zeros(length(cond_var_idx) * T.nExo) kronxS = zeros(T.nExo * length(cond_var_idx)) - end # timeit_debug - @timeit_debug timer "Main loop" begin + # end # timeit_debug + # @timeit_debug timer "Main loop" begin for i in reverse(axes(data_in_deviations,2)) # stt = 𝐒⁻¹ * aug_state + 𝐒⁻² * ℒ.kron(aug_state, aug_state) / 2 @@ -1640,8 +1640,8 @@ function rrule(::typeof(calculate_inversion_filter_loglikelihood), ∂state += ∂state¹⁻_vol[1:end-1] end - end # timeit_debug - @timeit_debug timer "Post allocation" begin + # end # timeit_debug + # @timeit_debug timer "Post allocation" begin ∂𝐒 = [copy(𝐒[1]) * 0, copy(𝐒[2]) * 0] @@ -1660,8 +1660,8 @@ function rrule(::typeof(calculate_inversion_filter_loglikelihood), return NoTangent(), NoTangent(), ℒ.I(T.nVars)[:,T.past_not_future_and_mixed_idx] * ∂state * ∂llh, ∂𝐒, ∂data_in_deviations * ∂llh, NoTangent(), NoTangent(), NoTangent(), NoTangent(), NoTangent() end - end # timeit_debug - end # timeit_debug + # end # timeit_debug + # end # timeit_debug # See: https://pcubaborda.net/documents/CGIZ-final.pdf llh = -(logabsdets + shocks² + (length(observables) * (warmup_iterations + n_obs - presample_periods)) * log(2 * 3.141592653589793)) / 2 @@ -1678,12 +1678,12 @@ function calculate_inversion_filter_loglikelihood(::Val{:pruned_third_order}, data_in_deviations::Matrix{Float64}, observables::Union{Vector{String}, Vector{Symbol}}, T::timings; - timer::TimerOutput = TimerOutput(), + # timer::TimerOutput = TimerOutput(), warmup_iterations::Int = 0, presample_periods::Int = 0, verbose::Bool = false, filter_algorithm::Symbol = :LagrangeNewton) - @timeit_debug timer "Inversion filter" begin + # @timeit_debug timer "Inversion filter" begin precision_factor = 1.0 @@ -1791,7 +1791,7 @@ function calculate_inversion_filter_loglikelihood(::Val{:pruned_third_order}, kron_buffer4 = ℒ.kron(ℒ.kron(J, J), zeros(T.nExo)) - @timeit_debug timer "Loop" begin + # @timeit_debug timer "Loop" begin for i in axes(data_in_deviations,2) state¹⁻ = state[1] @@ -1839,7 +1839,7 @@ function calculate_inversion_filter_loglikelihood(::Val{:pruned_third_order}, # ) # println(x²) - @timeit_debug timer "Find shocks" begin + # @timeit_debug timer "Find shocks" begin x, matched = find_shocks(Val(filter_algorithm), init_guess, kron_buffer, @@ -1854,7 +1854,7 @@ function calculate_inversion_filter_loglikelihood(::Val{:pruned_third_order}, shock_independent, # max_iter = 200 ) - end # timeit_debug + # end # timeit_debug # println(x) # println("$filter_algorithm: $matched; current x: $x, $(ℒ.norm(x))") @@ -2016,8 +2016,8 @@ function calculate_inversion_filter_loglikelihood(::Val{:pruned_third_order}, state = [𝐒⁻¹ * aug_state₁, 𝐒⁻¹ * aug_state₂ + 𝐒⁻² * kron_aug_state₁ / 2, 𝐒⁻¹ * aug_state₃ + 𝐒⁻² * ℒ.kron(aug_state₁̂, aug_state₂) + 𝐒⁻³ * ℒ.kron(kron_aug_state₁,aug_state₁) / 6] end - end # timeit_debug - end # timeit_debug + # end # timeit_debug + # end # timeit_debug # See: https://pcubaborda.net/documents/CGIZ-final.pdf return -(logabsdets + shocks² + (length(observables) * (warmup_iterations + n_obs - presample_periods)) * log(2 * 3.141592653589793)) / 2 @@ -2033,12 +2033,12 @@ function rrule(::typeof(calculate_inversion_filter_loglikelihood), data_in_deviations::Matrix{Float64}, observables::Union{Vector{String}, Vector{Symbol}}, T::timings; - timer::TimerOutput = TimerOutput(), + # timer::TimerOutput = TimerOutput(), warmup_iterations::Int = 0, presample_periods::Int = 0, verbose::Bool = false, filter_algorithm::Symbol = :LagrangeNewton) - @timeit_debug timer "Inversion filter - forward" begin + # @timeit_debug timer "Inversion filter - forward" begin precision_factor = 1.0 n_obs = size(data_in_deviations,2) @@ -2193,7 +2193,7 @@ function rrule(::typeof(calculate_inversion_filter_loglikelihood), 𝐒ⁱ³ᵉ = 𝐒³ᵉ / 6 - @timeit_debug timer "Loop" begin + # @timeit_debug timer "Loop" begin for i in axes(data_in_deviations,2) state¹⁻ = state₁ @@ -2225,7 +2225,7 @@ function rrule(::typeof(calculate_inversion_filter_loglikelihood), init_guess = zeros(size(𝐒ⁱ, 2)) - @timeit_debug timer "Find shocks" begin + # @timeit_debug timer "Find shocks" begin x[i], matched = find_shocks(Val(filter_algorithm), init_guess, kronxx[i], @@ -2240,7 +2240,7 @@ function rrule(::typeof(calculate_inversion_filter_loglikelihood), shock_independent, # max_iter = 100 ) - end # timeit_debug + # end # timeit_debug if !matched if verbose println("Inversion filter failed at step $i") end @@ -2287,7 +2287,7 @@ function rrule(::typeof(calculate_inversion_filter_loglikelihood), state₁, state₂, state₃ = [𝐒⁻¹ * aug_state₁[i], 𝐒⁻¹ * aug_state₂[i] + 𝐒⁻² * kron_aug_state₁[i] / 2, 𝐒⁻¹ * aug_state₃[i] + 𝐒⁻² * ℒ.kron(aug_state₁̂[i], aug_state₂[i]) + 𝐒⁻³ * ℒ.kron(kron_aug_state₁[i], aug_state₁[i]) / 6] end - end # timeit_debug + # end # timeit_debug # See: https://pcubaborda.net/documents/CGIZ-final.pdf llh = -(logabsdets + shocks² + (length(observables) * (warmup_iterations + n_obs - presample_periods)) * log(2 * 3.141592653589793)) / 2 @@ -2298,10 +2298,10 @@ function rrule(::typeof(calculate_inversion_filter_loglikelihood), ∂data_in_deviations = similar(data_in_deviations) - end # timeit_debug + # end # timeit_debug function inversion_filter_loglikelihood_pullback(∂llh) - @timeit_debug timer "Inversion filter - pullback" begin + # @timeit_debug timer "Inversion filter - pullback" begin ∂𝐒ⁱ = zero(𝐒ⁱ) ∂𝐒²ᵉ = zero(𝐒²ᵉ) ∂𝐒ⁱ³ᵉ = zero(𝐒ⁱ³ᵉ) @@ -2329,7 +2329,7 @@ function rrule(::typeof(calculate_inversion_filter_loglikelihood), ∂kronstate¹⁻_vol = zeros(length(state¹⁻_vol)^2) ∂state = [zeros(T.nPast_not_future_and_mixed), zeros(T.nPast_not_future_and_mixed), zeros(T.nPast_not_future_and_mixed)] - @timeit_debug timer "Loop" begin + # @timeit_debug timer "Loop" begin for i in reverse(axes(data_in_deviations,2)) # state₁ = 𝐒⁻¹ * aug_state₁[i] ∂𝐒⁻¹ += ∂state[1] * aug_state₁[i]' @@ -2535,7 +2535,7 @@ function rrule(::typeof(calculate_inversion_filter_loglikelihood), # state¹⁻_vol = vcat(state¹⁻, 1) ∂state[1] += ∂state¹⁻_vol[1:end-1] end - end # timeit_debug + # end # timeit_debug ∂𝐒 = [copy(𝐒[1]) * 0, copy(𝐒[2]) * 0, copy(𝐒[3]) * 0] @@ -2565,7 +2565,7 @@ function rrule(::typeof(calculate_inversion_filter_loglikelihood), ∂state[2] = ℒ.I(T.nVars)[:,T.past_not_future_and_mixed_idx] * ∂state[2] * ∂llh ∂state[3] = ℒ.I(T.nVars)[:,T.past_not_future_and_mixed_idx] * ∂state[3] * ∂llh - end # timeit_debug + # end # timeit_debug return NoTangent(), NoTangent(), ∂state, ∂𝐒, ∂data_in_deviations * ∂llh, NoTangent(), NoTangent(), NoTangent(), NoTangent(), NoTangent() end @@ -2583,13 +2583,13 @@ function calculate_inversion_filter_loglikelihood(::Val{:third_order}, data_in_deviations::Matrix{Float64}, observables::Union{Vector{String}, Vector{Symbol}}, T::timings; - timer::TimerOutput = TimerOutput(), + # timer::TimerOutput = TimerOutput(), warmup_iterations::Int = 0, presample_periods::Int = 0, verbose::Bool = false, filter_algorithm::Symbol = :LagrangeNewton) - @timeit_debug timer "3rd - Inversion filter" begin - @timeit_debug timer "Preallocation" begin + # @timeit_debug timer "3rd - Inversion filter" begin + # @timeit_debug timer "Preallocation" begin precision_factor = 1.0 @@ -2690,8 +2690,8 @@ function calculate_inversion_filter_loglikelihood(::Val{:third_order}, II = sparse(ℒ.I(T.nExo^2)) - end # timeit_debug - @timeit_debug timer "Loop" begin + # end # timeit_debug + # @timeit_debug timer "Loop" begin for i in axes(data_in_deviations,2) state¹⁻ = state @@ -2716,7 +2716,7 @@ function calculate_inversion_filter_loglikelihood(::Val{:third_order}, init_guess = zeros(size(𝐒ⁱ, 2)) - @timeit_debug timer "Find shocks" begin + # @timeit_debug timer "Find shocks" begin x, matched = find_shocks(Val(filter_algorithm), init_guess, kron_buffer, @@ -2731,7 +2731,7 @@ function calculate_inversion_filter_loglikelihood(::Val{:third_order}, shock_independent, # max_iter = 200 ) - end # timeit_debug + # end # timeit_debug # println("$filter_algorithm: $matched; current x: $x, $(ℒ.norm(x))") # if !matched @@ -2888,8 +2888,8 @@ function calculate_inversion_filter_loglikelihood(::Val{:third_order}, # state = state_update(state, x) end - end # timeit_debug - end # timeit_debug + # end # timeit_debug + # end # timeit_debug # See: https://pcubaborda.net/documents/CGIZ-final.pdf return -(logabsdets + shocks² + (length(observables) * (warmup_iterations + n_obs - presample_periods)) * log(2 * 3.141592653589793)) / 2 @@ -2905,13 +2905,13 @@ function rrule(::typeof(calculate_inversion_filter_loglikelihood), data_in_deviations::Matrix{Float64}, observables::Union{Vector{String}, Vector{Symbol}}, T::timings; - timer::TimerOutput = TimerOutput(), + # timer::TimerOutput = TimerOutput(), warmup_iterations::Int = 0, presample_periods::Int = 0, verbose::Bool = false, filter_algorithm::Symbol = :LagrangeNewton) - @timeit_debug timer "Inversion filter pruned 2nd - forward" begin - @timeit_debug timer "Preallocation" begin + # @timeit_debug timer "Inversion filter pruned 2nd - forward" begin + # @timeit_debug timer "Preallocation" begin precision_factor = 1.0 @@ -3047,8 +3047,8 @@ function rrule(::typeof(calculate_inversion_filter_loglikelihood), 𝐒ⁱ³ᵉ = 𝐒³ᵉ / 6 - end # timeit_debug - @timeit_debug timer "Main loop" begin + # end # timeit_debug + # @timeit_debug timer "Main loop" begin for i in axes(data_in_deviations,2) state¹⁻ = stt @@ -3069,7 +3069,7 @@ function rrule(::typeof(calculate_inversion_filter_loglikelihood), init_guess = zeros(size(𝐒ⁱ, 2)) - @timeit_debug timer "Find shocks" begin + # @timeit_debug timer "Find shocks" begin x[i], matched = find_shocks(Val(filter_algorithm), init_guess, kronxx[i], @@ -3084,7 +3084,7 @@ function rrule(::typeof(calculate_inversion_filter_loglikelihood), shock_independent, # max_iter = 100 ) - end # timeit_debug + # end # timeit_debug if !matched if verbose println("Inversion filter failed at step $i") end @@ -3130,8 +3130,8 @@ function rrule(::typeof(calculate_inversion_filter_loglikelihood), # See: https://pcubaborda.net/documents/CGIZ-final.pdf llh = -(logabsdets + shocks² + (length(observables) * (warmup_iterations + n_obs - presample_periods)) * log(2 * 3.141592653589793)) / 2 - end # timeit_debug - end # timeit_debug + # end # timeit_debug + # end # timeit_debug ∂state = similar(state) @@ -3140,8 +3140,8 @@ function rrule(::typeof(calculate_inversion_filter_loglikelihood), ∂data_in_deviations = similar(data_in_deviations) function inversion_filter_loglikelihood_pullback(∂llh) - @timeit_debug timer "Inversion filter pruned 2nd - pullback" begin - @timeit_debug timer "Preallocation" begin + # @timeit_debug timer "Inversion filter pruned 2nd - pullback" begin + # @timeit_debug timer "Preallocation" begin ∂𝐒ⁱ = zero(𝐒ⁱ) ∂𝐒²ᵉ = zero(𝐒²ᵉ) @@ -3166,8 +3166,8 @@ function rrule(::typeof(calculate_inversion_filter_loglikelihood), ∂kronstate¹⁻_vol = zeros(length(state¹⁻_vol)^2) ∂state = zeros(T.nPast_not_future_and_mixed) - end # timeit_debug - @timeit_debug timer "Main loop" begin + # end # timeit_debug + # @timeit_debug timer "Main loop" begin for i in reverse(axes(data_in_deviations,2)) # stt = 𝐒⁻¹ * aug_state[i] + 𝐒⁻² * ℒ.kron(aug_state[i], aug_state[i]) / 2 + 𝐒⁻³ * ℒ.kron(ℒ.kron(aug_state[i],aug_state[i]),aug_state[i]) / 6 @@ -3318,8 +3318,8 @@ function rrule(::typeof(calculate_inversion_filter_loglikelihood), ∂state += ∂state¹⁻_vol[1:end-1] end - end # timeit_debug - @timeit_debug timer "Post allocation" begin + # end # timeit_debug + # @timeit_debug timer "Post allocation" begin ∂𝐒 = [copy(𝐒[1]) * 0, copy(𝐒[2]) * 0, copy(𝐒[3]) * 0] @@ -3345,8 +3345,8 @@ function rrule(::typeof(calculate_inversion_filter_loglikelihood), return NoTangent(), NoTangent(), ℒ.I(T.nVars)[:,T.past_not_future_and_mixed_idx] * ∂state * ∂llh, ∂𝐒, ∂data_in_deviations * ∂llh, NoTangent(), NoTangent(), NoTangent(), NoTangent(), NoTangent() end - end # timeit_debug - end # timeit_debug + # end # timeit_debug + # end # timeit_debug return llh, inversion_filter_loglikelihood_pullback end diff --git a/src/filter/kalman.jl b/src/filter/kalman.jl index 7cb612eb..c0b3fafc 100644 --- a/src/filter/kalman.jl +++ b/src/filter/kalman.jl @@ -10,8 +10,8 @@ function calculate_loglikelihood(::Val{:kalman}, state, warmup_iterations, filter_algorithm, - verbose; - timer::TimerOutput = TimerOutput()) + verbose) #; + # timer::TimerOutput = TimerOutput()) return calculate_kalman_filter_loglikelihood(observables, 𝐒, data_in_deviations, @@ -26,7 +26,7 @@ function calculate_kalman_filter_loglikelihood(observables::Vector{Symbol}, 𝐒::Union{Matrix{S},Vector{AbstractMatrix{S}}}, data_in_deviations::Matrix{S}, T::timings; - timer::TimerOutput = TimerOutput(), + # timer::TimerOutput = TimerOutput(), presample_periods::Int = 0, initial_covariance::Symbol = :theoretical, verbose::Bool = false)::S where S <: Real @@ -39,7 +39,7 @@ function calculate_kalman_filter_loglikelihood(observables::Vector{String}, 𝐒::Union{Matrix{S},Vector{AbstractMatrix{S}}}, data_in_deviations::Matrix{S}, T::timings; - timer::TimerOutput = TimerOutput(), + # timer::TimerOutput = TimerOutput(), presample_periods::Int = 0, initial_covariance::Symbol = :theoretical, verbose::Bool = false)::S where S <: Real @@ -52,7 +52,7 @@ function calculate_kalman_filter_loglikelihood(observables_index::Vector{Int}, 𝐒::Union{Matrix{S},Vector{AbstractMatrix{S}}}, data_in_deviations::Matrix{S}, T::timings; - timer::TimerOutput = TimerOutput(), + # timer::TimerOutput = TimerOutput(), presample_periods::Int = 0, initial_covariance::Symbol = :theoretical, lyapunov_algorithm::Symbol = :doubling, @@ -77,9 +77,9 @@ end function get_initial_covariance(::Val{:theoretical}, A::AbstractMatrix{S}, B::AbstractMatrix{S}; - lyapunov_algorithm::Symbol = :doubling, - verbose::Bool = false, - timer::TimerOutput = TimerOutput())::Matrix{S} where S <: Real + lyapunov_algorithm::Symbol = :doubling, + # timer::TimerOutput = TimerOutput(), + verbose::Bool = false)::Matrix{S} where S <: Real P, _ = solve_lyapunov_equation(A, B, lyapunov_algorithm = lyapunov_algorithm, verbose = verbose, timer = timer) return P end @@ -89,9 +89,9 @@ end function get_initial_covariance(::Val{:diagonal}, A::AbstractMatrix{S}, B::AbstractMatrix{S}; - lyapunov_algorithm::Symbol = :doubling, - verbose::Bool = false, - timer::TimerOutput = TimerOutput())::Matrix{S} where S <: Real + lyapunov_algorithm::Symbol = :doubling, + # timer::TimerOutput = TimerOutput(), + verbose::Bool = false)::Matrix{S} where S <: Real P = @ignore_derivatives collect(ℒ.I(size(A, 1)) * 10.0) return P end @@ -103,9 +103,9 @@ function run_kalman_iterations(A::Matrix{S}, P::Matrix{S}, data_in_deviations::Matrix{S}; presample_periods::Int = 0, - verbose::Bool = false, - timer::TimerOutput = TimerOutput())::S where S <: Float64 - @timeit_debug timer "Calculate Kalman filter" begin + # timer::TimerOutput = TimerOutput(), + verbose::Bool = false)::S where S <: Float64 + # @timeit_debug timer "Calculate Kalman filter" begin u = zeros(S, size(C,2)) @@ -127,7 +127,7 @@ function run_kalman_iterations(A::Matrix{S}, tmp = similar(P) Ptmp = similar(P) - @timeit_debug timer "Loop" begin + # @timeit_debug timer "Loop" begin for t in 1:size(data_in_deviations, 2) if !all(isfinite.(z)) if verbose println("KF not finite at step $t") end @@ -141,9 +141,9 @@ function run_kalman_iterations(A::Matrix{S}, mul!(F, Ctmp, C') # F = C * P * C' - @timeit_debug timer "LU factorisation" begin + # @timeit_debug timer "LU factorisation" begin luF = RF.lu!(F, check = false) ### has to be LU since F will always be symmetric and positive semi-definite but not positive definite (due to linear dependencies) - end # timeit_debug + # end # timeit_debug if !ℒ.issuccess(luF) if verbose println("KF factorisation failed step $t") end @@ -160,7 +160,7 @@ function run_kalman_iterations(A::Matrix{S}, # invF = inv(luF) ### - @timeit_debug timer "LU div" begin + # @timeit_debug timer "LU div" begin if t > presample_periods ℒ.ldiv!(ztmp, luF, z) loglik += log(Fdet) + ℒ.dot(z', ztmp) ### @@ -175,8 +175,8 @@ function run_kalman_iterations(A::Matrix{S}, # K = P * Ct / luF # K = P * C' * invF - end # timeit_debug - @timeit_debug timer "Matmul" begin + # end # timeit_debug + # @timeit_debug timer "Matmul" begin mul!(tmp, K, C) mul!(Ptmp, tmp, P) @@ -195,11 +195,11 @@ function run_kalman_iterations(A::Matrix{S}, mul!(z, C, u) # z = C * u - end # timeit_debug + # end # timeit_debug end - end # timeit_debug - end # timeit_debug + # end # timeit_debug + # end # timeit_debug return -(loglik + ((size(data_in_deviations, 2) - presample_periods) * size(data_in_deviations, 1)) * log(2 * 3.141592653589793)) / 2 end @@ -212,9 +212,9 @@ function run_kalman_iterations(A::Matrix{S}, P::Matrix{S}, data_in_deviations::Matrix{S}; presample_periods::Int = 0, - verbose::Bool = false, - timer::TimerOutput = TimerOutput())::S where S <: ℱ.Dual - @timeit_debug timer "Calculate Kalman filter - forward mode AD" begin + # timer::TimerOutput = TimerOutput(), + verbose::Bool = false)::S where S <: ℱ.Dual + # @timeit_debug timer "Calculate Kalman filter - forward mode AD" begin u = zeros(S, size(C,2)) z = C * u @@ -265,7 +265,7 @@ function run_kalman_iterations(A::Matrix{S}, z = C * u end - end # timeit_debug + # end # timeit_debug return -(loglik + ((size(data_in_deviations, 2) - presample_periods) * size(data_in_deviations, 1)) * log(2 * 3.141592653589793)) / 2 end @@ -278,9 +278,9 @@ function rrule(::typeof(run_kalman_iterations), P, data_in_deviations; presample_periods = 0, - verbose::Bool = false, - timer::TimerOutput = TimerOutput()) - @timeit_debug timer "Calculate Kalman filter - forward" begin + # timer::TimerOutput = TimerOutput(), + verbose::Bool = false) + # @timeit_debug timer "Calculate Kalman filter - forward" begin T = size(data_in_deviations, 2) + 1 z = zeros(size(data_in_deviations, 1)) @@ -309,7 +309,7 @@ function rrule(::typeof(run_kalman_iterations), loglik = 0.0 - @timeit_debug timer "Loop" begin + # @timeit_debug timer "Loop" begin for t in 2:T if !all(isfinite.(z)) @@ -385,19 +385,19 @@ function rrule(::typeof(run_kalman_iterations), vtmp = zero(v[1]) Ptmp = zero(P[1]) - end # timeit_debug - end # timeit_debug + # end # timeit_debug + # end # timeit_debug # pullback function kalman_pullback(∂llh) - @timeit_debug timer "Calculate Kalman filter - reverse" begin + # @timeit_debug timer "Calculate Kalman filter - reverse" begin ℒ.rmul!(∂A, 0) ℒ.rmul!(∂Faccum, 0) ℒ.rmul!(∂P, 0) ℒ.rmul!(∂ū, 0) ℒ.rmul!(∂𝐁, 0) - @timeit_debug timer "Loop" begin + # @timeit_debug timer "Loop" begin for t in T:-1:2 if t > presample_periods + 1 # ∂llh∂F @@ -534,8 +534,8 @@ function rrule(::typeof(run_kalman_iterations), ℒ.rmul!(∂𝐁, -∂llh/2) ℒ.rmul!(∂data_in_deviations, -∂llh/2) - end # timeit_debug - end # timeit_debug + # end # timeit_debug + # end # timeit_debug return NoTangent(), ∂A, ∂𝐁, NoTangent(), ∂P, ∂data_in_deviations, NoTangent() end diff --git a/src/get_functions.jl b/src/get_functions.jl index 1176028c..7d126fba 100644 --- a/src/get_functions.jl +++ b/src/get_functions.jl @@ -994,7 +994,7 @@ function get_irf(𝓂::ℳ; timer::TimerOutput = TimerOutput(), verbose::Bool = false) - @timeit_debug timer "Wrangling inputs" begin + # @timeit_debug timer "Wrangling inputs" begin shocks = shocks isa KeyedArray ? axiskeys(shocks,1) isa Vector{String} ? rekey(shocks, 1 => axiskeys(shocks,1) .|> Meta.parse .|> replace_indices) : shocks : shocks @@ -1046,9 +1046,9 @@ function get_irf(𝓂::ℳ; occasionally_binding_constraints = length(𝓂.obc_violation_equations) > 0 end - end # timeit_debug + # end # timeit_debug - @timeit_debug timer "Solve model" begin + # @timeit_debug timer "Solve model" begin solve!(𝓂, parameters = parameters, @@ -1058,13 +1058,13 @@ function get_irf(𝓂::ℳ; obc = occasionally_binding_constraints || obc_shocks_included, timer = timer) - end # timeit_debug + # end # timeit_debug - @timeit_debug timer "Get relevant steady state" begin + # @timeit_debug timer "Get relevant steady state" begin reference_steady_state, NSSS, SSS_delta = get_relevant_steady_states(𝓂, algorithm, verbose = verbose) - end # timeit_debug + # end # timeit_debug unspecified_initial_state = initial_state == [0.0] @@ -1103,7 +1103,7 @@ function get_irf(𝓂::ℳ; end if generalised_irf - @timeit_debug timer "Calculate IRFs" begin + # @timeit_debug timer "Calculate IRFs" begin girfs = girf(state_update, initial_state, levels ? reference_steady_state + SSS_delta : SSS_delta, @@ -1113,7 +1113,7 @@ function get_irf(𝓂::ℳ; variables = variables, shock_size = shock_size, negative_shock = negative_shock)#, warmup_periods::Int = 100, draws::Int = 50, iterations_to_steady_state::Int = 500) - end # timeit_debug + # end # timeit_debug return girfs else @@ -1184,7 +1184,7 @@ function get_irf(𝓂::ℳ; return present_states, present_shocks, solved end - @timeit_debug timer "Calculate IRFs" begin + # @timeit_debug timer "Calculate IRFs" begin irfs = irf(state_update, obc_state_update, initial_state, @@ -1195,9 +1195,9 @@ function get_irf(𝓂::ℳ; variables = variables, shock_size = shock_size, negative_shock = negative_shock) - end # timeit_debug + # end # timeit_debug else - @timeit_debug timer "Calculate IRFs" begin + # @timeit_debug timer "Calculate IRFs" begin irfs = irf(state_update, initial_state, levels ? reference_steady_state + SSS_delta : SSS_delta, @@ -1207,7 +1207,7 @@ function get_irf(𝓂::ℳ; variables = variables, shock_size = shock_size, negative_shock = negative_shock) - end # timeit_debug + # end # timeit_debug end return irfs diff --git a/src/perturbation.jl b/src/perturbation.jl index 7ce85fdd..a48ed91f 100644 --- a/src/perturbation.jl +++ b/src/perturbation.jl @@ -2,10 +2,10 @@ function calculate_first_order_solution(∇₁::Matrix{Float64}; T::timings, quadratic_matrix_equation_solver::Symbol = :schur, verbose::Bool = false, - initial_guess::AbstractMatrix{<:AbstractFloat} = zeros(0,0), - timer::TimerOutput = TimerOutput())::Tuple{Matrix{Float64}, Matrix{Float64}, Bool} - @timeit_debug timer "Calculate 1st order solution" begin - @timeit_debug timer "Preprocessing" begin + # timer::TimerOutput = TimerOutput(), + initial_guess::AbstractMatrix{<:AbstractFloat} = zeros(0,0))::Tuple{Matrix{Float64}, Matrix{Float64}, Bool} + # @timeit_debug timer "Calculate 1st order solution" begin + # @timeit_debug timer "Preprocessing" begin dynIndex = T.nPresent_only+1:T.nVars @@ -24,8 +24,8 @@ function calculate_first_order_solution(∇₁::Matrix{Float64}; ∇₋ = ∇₁[:,T.nFuture_not_past_and_mixed + T.nVars .+ range(1, T.nPast_not_future_and_mixed)] ∇ₑ = ∇₁[:,(T.nFuture_not_past_and_mixed + T.nVars + T.nPast_not_future_and_mixed + 1):end] - end # timeit_debug - @timeit_debug timer "Invert ∇₀" begin + # end # timeit_debug + # @timeit_debug timer "Invert ∇₀" begin Q = ℒ.qr!(∇₀[:,T.present_only_idx]) @@ -33,15 +33,15 @@ function calculate_first_order_solution(∇₁::Matrix{Float64}; A₀ = Q.Q' * ∇₀ A₋ = Q.Q' * ∇₋ - end # timeit_debug - @timeit_debug timer "Sort matrices" begin + # end # timeit_debug + # @timeit_debug timer "Sort matrices" begin Ã₊ = A₊[dynIndex,:] * Ir[future_not_past_and_mixed_in_comb,:] Ã₀ = A₀[dynIndex, comb] Ã₋ = A₋[dynIndex,:] * Ir[past_not_future_and_mixed_in_comb,:] - end # timeit_debug - @timeit_debug timer "Quadratic matrix equation solve" begin + # end # timeit_debug + # @timeit_debug timer "Quadratic matrix equation solve" begin sol, solved = solve_quadratic_matrix_equation(Ã₊, Ã₀, Ã₋, T, @@ -55,9 +55,9 @@ function calculate_first_order_solution(∇₁::Matrix{Float64}; return zeros(T.nVars,T.nPast_not_future_and_mixed + T.nExo), sol, false end - end # timeit_debug - @timeit_debug timer "Postprocessing" begin - @timeit_debug timer "Setup matrices" begin + # end # timeit_debug + # @timeit_debug timer "Postprocessing" begin + # @timeit_debug timer "Setup matrices" begin sol_compact = sol[reverse_dynamic_order, past_not_future_and_mixed_in_comb] @@ -70,8 +70,8 @@ function calculate_first_order_solution(∇₁::Matrix{Float64}; Ã₀ᵤ = A₀[1:T.nPresent_only, T.present_but_not_only_idx] A₋ᵤ = A₋[1:T.nPresent_only,:] - end # timeit_debug - @timeit_debug timer "Invert Ā₀ᵤ" begin + # end # timeit_debug + # @timeit_debug timer "Invert Ā₀ᵤ" begin Ā̂₀ᵤ = ℒ.lu!(Ā₀ᵤ, check = false) @@ -91,9 +91,9 @@ function calculate_first_order_solution(∇₁::Matrix{Float64}; A = vcat(A₋ᵤ, sol_compact)[T.reorder,:] - end # timeit_debug - end # timeit_debug - @timeit_debug timer "Exogenous part solution" begin + # end # timeit_debug + # end # timeit_debug + # @timeit_debug timer "Exogenous part solution" begin M = A[T.future_not_past_and_mixed_idx,:] * ℒ.I(T.nVars)[T.past_not_future_and_mixed_idx,:] @@ -109,8 +109,8 @@ function calculate_first_order_solution(∇₁::Matrix{Float64}; ℒ.ldiv!(C, ∇ₑ) ℒ.rmul!(∇ₑ, -1) - end # timeit_debug - end # timeit_debug + # end # timeit_debug + # end # timeit_debug return hcat(A, ∇ₑ), sol, true end @@ -121,11 +121,11 @@ function rrule(::typeof(calculate_first_order_solution), T::timings, quadratic_matrix_equation_solver::Symbol = :schur, verbose::Bool = false, - initial_guess::AbstractMatrix{<:AbstractFloat} = zeros(0,0), - timer::TimerOutput = TimerOutput()) + # timer::TimerOutput = TimerOutput(), + initial_guess::AbstractMatrix{<:AbstractFloat} = zeros(0,0)) # Forward pass to compute the output and intermediate values needed for the backward pass - @timeit_debug timer "Calculate 1st order solution" begin - @timeit_debug timer "Preprocessing" begin + # @timeit_debug timer "Calculate 1st order solution" begin + # @timeit_debug timer "Preprocessing" begin dynIndex = T.nPresent_only+1:T.nVars @@ -144,8 +144,8 @@ function rrule(::typeof(calculate_first_order_solution), ∇₋ = ∇₁[:,T.nFuture_not_past_and_mixed + T.nVars .+ range(1, T.nPast_not_future_and_mixed)] ∇̂ₑ = ∇₁[:,(T.nFuture_not_past_and_mixed + T.nVars + T.nPast_not_future_and_mixed + 1):end] - end # timeit_debug - @timeit_debug timer "Invert ∇₀" begin + # end # timeit_debug + # @timeit_debug timer "Invert ∇₀" begin Q = ℒ.qr!(∇₀[:,T.present_only_idx]) @@ -153,15 +153,15 @@ function rrule(::typeof(calculate_first_order_solution), A₀ = Q.Q' * ∇₀ A₋ = Q.Q' * ∇₋ - end # timeit_debug - @timeit_debug timer "Sort matrices" begin + # end # timeit_debug + # @timeit_debug timer "Sort matrices" begin Ã₊ = A₊[dynIndex,:] * Ir[future_not_past_and_mixed_in_comb,:] Ã₀ = A₀[dynIndex, comb] Ã₋ = A₋[dynIndex,:] * Ir[past_not_future_and_mixed_in_comb,:] - end # timeit_debug - @timeit_debug timer "Quadratic matrix equation solve" begin + # end # timeit_debug + # @timeit_debug timer "Quadratic matrix equation solve" begin sol, solved = solve_quadratic_matrix_equation(Ã₊, Ã₀, Ã₋, T, @@ -174,9 +174,9 @@ function rrule(::typeof(calculate_first_order_solution), return (zeros(T.nVars,T.nPast_not_future_and_mixed + T.nExo), sol, false), x -> NoTangent(), NoTangent(), NoTangent() end - end # timeit_debug - @timeit_debug timer "Postprocessing" begin - @timeit_debug timer "Setup matrices" begin + # end # timeit_debug + # @timeit_debug timer "Postprocessing" begin + # @timeit_debug timer "Setup matrices" begin sol_compact = sol[reverse_dynamic_order, past_not_future_and_mixed_in_comb] @@ -189,8 +189,8 @@ function rrule(::typeof(calculate_first_order_solution), Ã₀ᵤ = A₀[1:T.nPresent_only, T.present_but_not_only_idx] A₋ᵤ = A₋[1:T.nPresent_only,:] - end # timeit_debug - @timeit_debug timer "Invert Ā₀ᵤ" begin + # end # timeit_debug + # @timeit_debug timer "Invert Ā₀ᵤ" begin Ā̂₀ᵤ = ℒ.lu!(Ā₀ᵤ, check = false) @@ -207,9 +207,9 @@ function rrule(::typeof(calculate_first_order_solution), ℒ.rmul!(A₋ᵤ, -1) end - end # timeit_debug - end # timeit_debug - @timeit_debug timer "Exogenous part solution" begin + # end # timeit_debug + # end # timeit_debug + # @timeit_debug timer "Exogenous part solution" begin expand = [ℒ.I(T.nVars)[T.future_not_past_and_mixed_idx,:], ℒ.I(T.nVars)[T.past_not_future_and_mixed_idx,:]] @@ -229,8 +229,8 @@ function rrule(::typeof(calculate_first_order_solution), ℒ.ldiv!(C, ∇̂ₑ) ℒ.rmul!(∇̂ₑ, -1) - end # timeit_debug - end # timeit_debug + # end # timeit_debug + # end # timeit_debug M = inv(C) @@ -276,8 +276,8 @@ function calculate_first_order_solution(∇₁::Matrix{ℱ.Dual{Z,S,N}}; T::timings, quadratic_matrix_equation_solver::Symbol = :schur, verbose::Bool = false, - initial_guess::AbstractMatrix{<:AbstractFloat} = zeros(0,0), - timer::TimerOutput = TimerOutput())::Tuple{Matrix{ℱ.Dual{Z,S,N}}, Matrix{Float64}, Bool} where {Z,S,N} + # timer::TimerOutput = TimerOutput(), + initial_guess::AbstractMatrix{<:AbstractFloat} = zeros(0,0))::Tuple{Matrix{ℱ.Dual{Z,S,N}}, Matrix{Float64}, Bool} where {Z,S,N} ∇̂₁ = ℱ.value.(∇₁) expand = [ℒ.I(T.nVars)[T.future_not_past_and_mixed_idx,:], ℒ.I(T.nVars)[T.past_not_future_and_mixed_idx,:]] @@ -378,9 +378,9 @@ function calculate_second_order_solution(∇₁::AbstractMatrix{S}, #first order initial_guess::AbstractMatrix{Float64} = zeros(0,0), sylvester_algorithm::Symbol = :doubling, tol::AbstractFloat = eps(), - timer::TimerOutput = TimerOutput(), + # timer::TimerOutput = TimerOutput(), verbose::Bool = false)::Tuple{<: AbstractMatrix{S}, Bool} where S <: Real - @timeit_debug timer "Calculate second order solution" begin + # @timeit_debug timer "Calculate second order solution" begin # inspired by Levintal @@ -394,7 +394,7 @@ function calculate_second_order_solution(∇₁::AbstractMatrix{S}, #first order n = T.nVars nₑ₋ = n₋ + 1 + nₑ - @timeit_debug timer "Setup matrices" begin + # @timeit_debug timer "Setup matrices" begin # 1st order solution 𝐒₁ = @views [𝑺₁[:,1:n₋] zeros(n) 𝑺₁[:,n₋+1:end]]# |> sparse @@ -415,9 +415,9 @@ function calculate_second_order_solution(∇₁::AbstractMatrix{S}, #first order ∇₁₊𝐒₁➕∇₁₀ = @views -∇₁[:,1:n₊] * 𝐒₁[i₊,1:n₋] * ℒ.I(n)[i₋,:] - ∇₁[:,range(1,n) .+ n₊] - end # timeit_debug + # end # timeit_debug - @timeit_debug timer "Invert matrix" begin + # @timeit_debug timer "Invert matrix" begin ∇₁₊𝐒₁➕∇₁₀lu = ℒ.lu(∇₁₊𝐒₁➕∇₁₀, check = false) @@ -429,34 +429,34 @@ function calculate_second_order_solution(∇₁::AbstractMatrix{S}, #first order # spinv = inv(∇₁₊𝐒₁➕∇₁₀) # spinv = choose_matrix_format(spinv) - end # timeit_debug + # end # timeit_debug - @timeit_debug timer "Setup second order matrices" begin - @timeit_debug timer "A" begin + # @timeit_debug timer "Setup second order matrices" begin + # @timeit_debug timer "A" begin ∇₁₊ = @views ∇₁[:,1:n₊] * ℒ.I(n)[i₊,:] A = ∇₁₊𝐒₁➕∇₁₀lu \ ∇₁₊ - end # timeit_debug - @timeit_debug timer "C" begin + # end # timeit_debug + # @timeit_debug timer "C" begin # ∇₂⎸k⎸𝐒₁𝐒₁₋╱𝟏ₑ⎹╱𝐒₁╱𝟏ₑ₋➕𝛔k𝐒₁₊╱𝟎⎹ = ∇₂ * (ℒ.kron(⎸𝐒₁𝐒₁₋╱𝟏ₑ⎹╱𝐒₁╱𝟏ₑ₋, ⎸𝐒₁𝐒₁₋╱𝟏ₑ⎹╱𝐒₁╱𝟏ₑ₋) + ℒ.kron(𝐒₁₊╱𝟎, 𝐒₁₊╱𝟎) * M₂.𝛔) * M₂.𝐂₂ ∇₂⎸k⎸𝐒₁𝐒₁₋╱𝟏ₑ⎹╱𝐒₁╱𝟏ₑ₋➕𝛔k𝐒₁₊╱𝟎⎹ = mat_mult_kron(∇₂, ⎸𝐒₁𝐒₁₋╱𝟏ₑ⎹╱𝐒₁╱𝟏ₑ₋, ⎸𝐒₁𝐒₁₋╱𝟏ₑ⎹╱𝐒₁╱𝟏ₑ₋, M₂.𝐂₂) + mat_mult_kron(∇₂, 𝐒₁₊╱𝟎, 𝐒₁₊╱𝟎, M₂.𝛔 * M₂.𝐂₂) C = ∇₁₊𝐒₁➕∇₁₀lu \ ∇₂⎸k⎸𝐒₁𝐒₁₋╱𝟏ₑ⎹╱𝐒₁╱𝟏ₑ₋➕𝛔k𝐒₁₊╱𝟎⎹ - end # timeit_debug - @timeit_debug timer "B" begin + # end # timeit_debug + # @timeit_debug timer "B" begin # 𝐒₁₋╱𝟏ₑ = choose_matrix_format(𝐒₁₋╱𝟏ₑ, density_threshold = 0.0) 𝐒₁₋╱𝟏ₑ = choose_matrix_format(𝐒₁₋╱𝟏ₑ, density_threshold = 0.0) B = mat_mult_kron(M₂.𝐔₂, 𝐒₁₋╱𝟏ₑ, 𝐒₁₋╱𝟏ₑ, M₂.𝐂₂) + M₂.𝐔₂ * M₂.𝛔 * M₂.𝐂₂ - end # timeit_debug - end # timeit_debug - @timeit_debug timer "Solve sylvester equation" begin + # end # timeit_debug + # end # timeit_debug + # @timeit_debug timer "Solve sylvester equation" begin 𝐒₂, solved = solve_sylvester_equation(A, B, C, sylvester_algorithm = sylvester_algorithm, @@ -465,8 +465,8 @@ function calculate_second_order_solution(∇₁::AbstractMatrix{S}, #first order # tol = tol, timer = timer) - end # timeit_debug - # @timeit_debug timer "Refine sylvester equation" begin + # end # timeit_debug + # # @timeit_debug timer "Refine sylvester equation" begin # # if !solved && !(sylvester_algorithm == :doubling) # # 𝐒₂, solved = solve_sylvester_equation(A, B, C, @@ -479,8 +479,8 @@ function calculate_second_order_solution(∇₁::AbstractMatrix{S}, #first order # # timer = timer) # # end - # end # timeit_debug - @timeit_debug timer "Post-process" begin + # # end # timeit_debug + # @timeit_debug timer "Post-process" begin # 𝐒₂ *= M₂.𝐔₂ @@ -490,8 +490,8 @@ function calculate_second_order_solution(∇₁::AbstractMatrix{S}, #first order return 𝐒₂, false end - end # timeit_debug - end # timeit_debug + # end # timeit_debug + # end # timeit_debug return 𝐒₂, true end @@ -508,9 +508,9 @@ function rrule(::typeof(calculate_second_order_solution), initial_guess::AbstractMatrix{Float64} = zeros(0,0), sylvester_algorithm::Symbol = :doubling, tol::AbstractFloat = eps(), - timer::TimerOutput = TimerOutput(), + # timer::TimerOutput = TimerOutput(), verbose::Bool = false) - @timeit_debug timer "Second order solution - forward" begin + # @timeit_debug timer "Second order solution - forward" begin # inspired by Levintal # Indices and number of variables @@ -523,7 +523,7 @@ function rrule(::typeof(calculate_second_order_solution), n = T.nVars nₑ₋ = n₋ + 1 + nₑ - @timeit_debug timer "Setup matrices" begin + # @timeit_debug timer "Setup matrices" begin # 1st order solution 𝐒₁ = @views [𝑺₁[:,1:n₋] zeros(n) 𝑺₁[:,n₋+1:end]]# |> sparse @@ -541,8 +541,8 @@ function rrule(::typeof(calculate_second_order_solution), ∇₁₊𝐒₁➕∇₁₀ = @views -∇₁[:,1:n₊] * 𝐒₁[i₊,1:n₋] * ℒ.I(n)[i₋,:] - ∇₁[:,range(1,n) .+ n₊] - end # timeit_debug - @timeit_debug timer "Invert matrix" begin + # end # timeit_debug + # @timeit_debug timer "Invert matrix" begin ∇₁₊𝐒₁➕∇₁₀lu = ℒ.lu(∇₁₊𝐒₁➕∇₁₀, check = false) @@ -554,33 +554,33 @@ function rrule(::typeof(calculate_second_order_solution), spinv = inv(∇₁₊𝐒₁➕∇₁₀lu) spinv = choose_matrix_format(spinv) - end # timeit_debug - @timeit_debug timer "Setup second order matrices" begin - @timeit_debug timer "A" begin + # end # timeit_debug + # @timeit_debug timer "Setup second order matrices" begin + # @timeit_debug timer "A" begin ∇₁₊ = @views ∇₁[:,1:n₊] * ℒ.I(n)[i₊,:] A = spinv * ∇₁₊ - end # timeit_debug - @timeit_debug timer "C" begin + # end # timeit_debug + # @timeit_debug timer "C" begin # ∇₂⎸k⎸𝐒₁𝐒₁₋╱𝟏ₑ⎹╱𝐒₁╱𝟏ₑ₋➕𝛔k𝐒₁₊╱𝟎⎹ = ∇₂ * (ℒ.kron(⎸𝐒₁𝐒₁₋╱𝟏ₑ⎹╱𝐒₁╱𝟏ₑ₋, ⎸𝐒₁𝐒₁₋╱𝟏ₑ⎹╱𝐒₁╱𝟏ₑ₋) + ℒ.kron(𝐒₁₊╱𝟎, 𝐒₁₊╱𝟎) * M₂.𝛔) * M₂.𝐂₂ ∇₂⎸k⎸𝐒₁𝐒₁₋╱𝟏ₑ⎹╱𝐒₁╱𝟏ₑ₋➕𝛔k𝐒₁₊╱𝟎⎹ = mat_mult_kron(∇₂, ⎸𝐒₁𝐒₁₋╱𝟏ₑ⎹╱𝐒₁╱𝟏ₑ₋, ⎸𝐒₁𝐒₁₋╱𝟏ₑ⎹╱𝐒₁╱𝟏ₑ₋, M₂.𝐂₂) + mat_mult_kron(∇₂, 𝐒₁₊╱𝟎, 𝐒₁₊╱𝟎, M₂.𝛔 * M₂.𝐂₂) C = spinv * ∇₂⎸k⎸𝐒₁𝐒₁₋╱𝟏ₑ⎹╱𝐒₁╱𝟏ₑ₋➕𝛔k𝐒₁₊╱𝟎⎹ - end # timeit_debug - @timeit_debug timer "B" begin + # end # timeit_debug + # @timeit_debug timer "B" begin # 𝐒₁₋╱𝟏ₑ = choose_matrix_format(𝐒₁₋╱𝟏ₑ, density_threshold = 0.0) 𝐒₁₋╱𝟏ₑ = choose_matrix_format(𝐒₁₋╱𝟏ₑ, density_threshold = 0.0) B = mat_mult_kron(M₂.𝐔₂, 𝐒₁₋╱𝟏ₑ, 𝐒₁₋╱𝟏ₑ, M₂.𝐂₂) + M₂.𝐔₂ * M₂.𝛔 * M₂.𝐂₂ - end # timeit_debug - end # timeit_debug - @timeit_debug timer "Solve sylvester equation" begin + # end # timeit_debug + # end # timeit_debug + # @timeit_debug timer "Solve sylvester equation" begin 𝐒₂, solved = solve_sylvester_equation(A, B, C, sylvester_algorithm = sylvester_algorithm, @@ -589,14 +589,14 @@ function rrule(::typeof(calculate_second_order_solution), # tol = tol, timer = timer) - end # timeit_debug - @timeit_debug timer "Post-process" begin + # end # timeit_debug + # @timeit_debug timer "Post-process" begin if !solved return (𝐒₂, solved), x -> NoTangent(), NoTangent(), NoTangent(), NoTangent(), NoTangent(), NoTangent(), NoTangent(), NoTangent(), NoTangent() end - end # timeit_debug + # end # timeit_debug # sp⎸𝐒₁𝐒₁₋╱𝟏ₑ⎹╱𝐒₁╱𝟏ₑ₋t = choose_matrix_format(⎸𝐒₁𝐒₁₋╱𝟏ₑ⎹╱𝐒₁╱𝟏ₑ₋', density_threshold = 1.0) @@ -610,12 +610,12 @@ function rrule(::typeof(calculate_second_order_solution), ∇₂t = choose_matrix_format(∇₂', density_threshold = 1.0) - end # timeit_debug + # end # timeit_debug function second_order_solution_pullback(∂𝐒₂_solved) - @timeit_debug timer "Second order solution - pullback" begin + # @timeit_debug timer "Second order solution - pullback" begin - @timeit_debug timer "Preallocate" begin + # @timeit_debug timer "Preallocate" begin ∂∇₂ = zeros(size(∇₂)) ∂∇₁ = zero(∇₁) ∂𝐒₁ = zero(𝐒₁) @@ -624,13 +624,13 @@ function rrule(::typeof(calculate_second_order_solution), ∂𝐒₁₊╱𝟎 = zeros(size(𝐒₁₊╱𝟎)) ∂⎸𝐒₁𝐒₁₋╱𝟏ₑ⎹╱𝐒₁╱𝟏ₑ₋ = zeros(size(⎸𝐒₁𝐒₁₋╱𝟏ₑ⎹╱𝐒₁╱𝟏ₑ₋)) - end # timeit_debug + # end # timeit_debug ∂𝐒₂ = ∂𝐒₂_solved[1] # ∂𝐒₂ *= 𝐔₂t - @timeit_debug timer "Sylvester" begin + # @timeit_debug timer "Sylvester" begin ∂C, solved = solve_sylvester_equation(A', B', ∂𝐒₂, sylvester_algorithm = sylvester_algorithm, @@ -642,9 +642,9 @@ function rrule(::typeof(calculate_second_order_solution), return (𝐒₂, solved), x -> NoTangent(), NoTangent(), NoTangent(), NoTangent(), NoTangent(), NoTangent(), NoTangent(), NoTangent(), NoTangent() end - end # timeit_debug + # end # timeit_debug - @timeit_debug timer "Matmul" begin + # @timeit_debug timer "Matmul" begin ∂C = choose_matrix_format(∂C) # Dense @@ -655,15 +655,15 @@ function rrule(::typeof(calculate_second_order_solution), # B = (M₂.𝐔₂ * ℒ.kron(𝐒₁₋╱𝟏ₑ, 𝐒₁₋╱𝟏ₑ) + M₂.𝐔₂ * M₂.𝛔) * M₂.𝐂₂ ∂kron𝐒₁₋╱𝟏ₑ = 𝐔₂t * ∂B * 𝐂₂t - end # timeit_debug + # end # timeit_debug - @timeit_debug timer "Kron adjoint" begin + # @timeit_debug timer "Kron adjoint" begin fill_kron_adjoint!(∂𝐒₁₋╱𝟏ₑ, ∂𝐒₁₋╱𝟏ₑ, ∂kron𝐒₁₋╱𝟏ₑ, 𝐒₁₋╱𝟏ₑ, 𝐒₁₋╱𝟏ₑ) - end # timeit_debug + # end # timeit_debug - @timeit_debug timer "Matmul2" begin + # @timeit_debug timer "Matmul2" begin # A = spinv * ∇₁₊ ∂∇₁₊ = spinv' * ∂A @@ -677,9 +677,9 @@ function rrule(::typeof(calculate_second_order_solution), ∂spinv += ∂C * ∇₂⎸k⎸𝐒₁𝐒₁₋╱𝟏ₑ⎹╱𝐒₁╱𝟏ₑ₋➕𝛔k𝐒₁₊╱𝟎⎹' - end # timeit_debug + # end # timeit_debug - @timeit_debug timer "Matmul3" begin + # @timeit_debug timer "Matmul3" begin # ∇₂⎸k⎸𝐒₁𝐒₁₋╱𝟏ₑ⎹╱𝐒₁╱𝟏ₑ₋➕𝛔k𝐒₁₊╱𝟎⎹ = ∇₂ * ℒ.kron(⎸𝐒₁𝐒₁₋╱𝟏ₑ⎹╱𝐒₁╱𝟏ₑ₋, ⎸𝐒₁𝐒₁₋╱𝟏ₑ⎹╱𝐒₁╱𝟏ₑ₋) * M₂.𝐂₂ + ∇₂ * ℒ.kron(𝐒₁₊╱𝟎, 𝐒₁₊╱𝟎) * M₂.𝛔 * M₂.𝐂₂ # kron⎸𝐒₁𝐒₁₋╱𝟏ₑ⎹╱𝐒₁╱𝟏ₑ₋ = choose_matrix_format(ℒ.kron(sp⎸𝐒₁𝐒₁₋╱𝟏ₑ⎹╱𝐒₁╱𝟏ₑ₋t, sp⎸𝐒₁𝐒₁₋╱𝟏ₑ⎹╱𝐒₁╱𝟏ₑ₋t), density_threshold = 1.0) @@ -696,29 +696,29 @@ function rrule(::typeof(calculate_second_order_solution), ∂∇₂ += mat_mult_kron(∂∇₂⎸k⎸𝐒₁𝐒₁₋╱𝟏ₑ⎹╱𝐒₁╱𝟏ₑ₋➕𝛔k𝐒₁₊╱𝟎⎹𝐂₂, ⎸𝐒₁𝐒₁₋╱𝟏ₑ⎹╱𝐒₁╱𝟏ₑ₋', ⎸𝐒₁𝐒₁₋╱𝟏ₑ⎹╱𝐒₁╱𝟏ₑ₋') - end # timeit_debug + # end # timeit_debug - @timeit_debug timer "Matmul4" begin + # @timeit_debug timer "Matmul4" begin ∂kron𝐒₁₊╱𝟎 = ∇₂t * ∂∇₂⎸k⎸𝐒₁𝐒₁₋╱𝟏ₑ⎹╱𝐒₁╱𝟏ₑ₋➕𝛔k𝐒₁₊╱𝟎⎹𝐂₂ * 𝛔t - end # timeit_debug + # end # timeit_debug - @timeit_debug timer "Kron adjoint 2" begin + # @timeit_debug timer "Kron adjoint 2" begin fill_kron_adjoint!(∂𝐒₁₊╱𝟎, ∂𝐒₁₊╱𝟎, ∂kron𝐒₁₊╱𝟎, 𝐒₁₊╱𝟎, 𝐒₁₊╱𝟎) - end # timeit_debug + # end # timeit_debug ∂kron⎸𝐒₁𝐒₁₋╱𝟏ₑ⎹╱𝐒₁╱𝟏ₑ₋ = ∇₂t * ∂∇₂⎸k⎸𝐒₁𝐒₁₋╱𝟏ₑ⎹╱𝐒₁╱𝟏ₑ₋➕𝛔k𝐒₁₊╱𝟎⎹𝐂₂ - @timeit_debug timer "Kron adjoint 3" begin + # @timeit_debug timer "Kron adjoint 3" begin fill_kron_adjoint!(∂⎸𝐒₁𝐒₁₋╱𝟏ₑ⎹╱𝐒₁╱𝟏ₑ₋, ∂⎸𝐒₁𝐒₁₋╱𝟏ₑ⎹╱𝐒₁╱𝟏ₑ₋, ∂kron⎸𝐒₁𝐒₁₋╱𝟏ₑ⎹╱𝐒₁╱𝟏ₑ₋, ⎸𝐒₁𝐒₁₋╱𝟏ₑ⎹╱𝐒₁╱𝟏ₑ₋, ⎸𝐒₁𝐒₁₋╱𝟏ₑ⎹╱𝐒₁╱𝟏ₑ₋) # filling dense is much faster - end # timeit_debug + # end # timeit_debug - @timeit_debug timer "Matmul5" begin + # @timeit_debug timer "Matmul5" begin # spinv = sparse(inv(∇₁₊𝐒₁➕∇₁₀)) ∂∇₁₊𝐒₁➕∇₁₀ = -spinv' * ∂spinv * spinv' @@ -748,9 +748,9 @@ function rrule(::typeof(calculate_second_order_solution), # 𝐒₁ = [𝑺₁[:,1:n₋] zeros(n) 𝑺₁[:,n₋+1:end]] ∂𝑺₁ = [∂𝐒₁[:,1:n₋] ∂𝐒₁[:,n₋+2:end]] - end # timeit_debug + # end # timeit_debug - end # timeit_debug + # end # timeit_debug return NoTangent(), ∂∇₁, ∂∇₂, ∂𝑺₁, NoTangent(), NoTangent(), NoTangent(), NoTangent(), NoTangent() end @@ -772,10 +772,10 @@ function calculate_third_order_solution(∇₁::AbstractMatrix{<: Real}, #first T::timings, initial_guess::AbstractMatrix{Float64} = zeros(0,0), sylvester_algorithm::Symbol = :bicgstab, - timer::TimerOutput = TimerOutput(), + # timer::TimerOutput = TimerOutput(), tol::AbstractFloat = 1e-12, # sylvester tol verbose::Bool = false) - @timeit_debug timer "Calculate third order solution" begin + # @timeit_debug timer "Calculate third order solution" begin # inspired by Levintal # Indices and number of variables @@ -788,7 +788,7 @@ function calculate_third_order_solution(∇₁::AbstractMatrix{<: Real}, #first n = T.nVars nₑ₋ = n₋ + 1 + nₑ - @timeit_debug timer "Setup matrices" begin + # @timeit_debug timer "Setup matrices" begin # 1st order solution 𝐒₁ = @views [𝑺₁[:,1:n₋] zeros(n) 𝑺₁[:,n₋+1:end]]# |> sparse @@ -807,8 +807,8 @@ function calculate_third_order_solution(∇₁::AbstractMatrix{<: Real}, #first ∇₁₊𝐒₁➕∇₁₀ = @views -∇₁[:,1:n₊] * 𝐒₁[i₊,1:n₋] * ℒ.I(n)[i₋,:] - ∇₁[:,range(1,n) .+ n₊] - end # timeit_debug - @timeit_debug timer "Invert matrix" begin + # end # timeit_debug + # @timeit_debug timer "Invert matrix" begin ∇₁₊𝐒₁➕∇₁₀lu = ℒ.lu(∇₁₊𝐒₁➕∇₁₀, check = false) @@ -820,45 +820,45 @@ function calculate_third_order_solution(∇₁::AbstractMatrix{<: Real}, #first # spinv = inv(∇₁₊𝐒₁➕∇₁₀) # spinv = choose_matrix_format(spinv) - end # timeit_debug + # end # timeit_debug ∇₁₊ = @views ∇₁[:,1:n₊] * ℒ.I(n)[i₊,:] A = ∇₁₊𝐒₁➕∇₁₀lu \ ∇₁₊ - @timeit_debug timer "Setup B" begin - @timeit_debug timer "Add tmpkron" begin + # @timeit_debug timer "Setup B" begin + # @timeit_debug timer "Add tmpkron" begin tmpkron = ℒ.kron(𝐒₁₋╱𝟏ₑ, M₂.𝛔) kron𝐒₁₋╱𝟏ₑ = ℒ.kron(𝐒₁₋╱𝟏ₑ, 𝐒₁₋╱𝟏ₑ) B = tmpkron - end # timeit_debug - @timeit_debug timer "Step 1" begin + # end # timeit_debug + # @timeit_debug timer "Step 1" begin B += M₃.𝐏₁ₗ̄ * tmpkron * M₃.𝐏₁ᵣ̃ - end # timeit_debug - @timeit_debug timer "Step 2" begin + # end # timeit_debug + # @timeit_debug timer "Step 2" begin B += M₃.𝐏₂ₗ̄ * tmpkron * M₃.𝐏₂ᵣ̃ - end # timeit_debug - @timeit_debug timer "Mult" begin + # end # timeit_debug + # @timeit_debug timer "Mult" begin B *= M₃.𝐂₃ B = choose_matrix_format(M₃.𝐔₃ * B) - end # timeit_debug - @timeit_debug timer "3rd Kronecker power" begin + # end # timeit_debug + # @timeit_debug timer "3rd Kronecker power" begin # B += mat_mult_kron(M₃.𝐔₃, collect(𝐒₁₋╱𝟏ₑ), collect(ℒ.kron(𝐒₁₋╱𝟏ₑ, 𝐒₁₋╱𝟏ₑ)), M₃.𝐂₃) # slower than direct compression B += compressed_kron³(𝐒₁₋╱𝟏ₑ, timer = timer) - end # timeit_debug - end # timeit_debug - @timeit_debug timer "Setup C" begin - @timeit_debug timer "Initialise smaller matrices" begin + # end # timeit_debug + # end # timeit_debug + # @timeit_debug timer "Setup C" begin + # @timeit_debug timer "Initialise smaller matrices" begin ⎸𝐒₂k𝐒₁₋╱𝟏ₑ➕𝐒₁𝐒₂₋⎹╱𝐒₂╱𝟎 = @views [(𝐒₂ * kron𝐒₁₋╱𝟏ₑ + 𝐒₁ * [𝐒₂[i₋,:] ; zeros(nₑ + 1, nₑ₋^2)])[i₊,:] 𝐒₂ @@ -872,8 +872,8 @@ function calculate_third_order_solution(∇₁::AbstractMatrix{<: Real}, #first aux = M₃.𝐒𝐏 * ⎸𝐒₁𝐒₁₋╱𝟏ₑ⎹╱𝐒₁╱𝟏ₑ₋ # aux = choose_matrix_format(aux, density_threshold = 1.0, min_length = 10) - end # timeit_debug - @timeit_debug timer "∇₃" begin + # end # timeit_debug + # @timeit_debug timer "∇₃" begin tmpkron = ℒ.kron(⎸𝐒₁𝐒₁₋╱𝟏ₑ⎹╱𝐒₁╱𝟏ₑ₋, ℒ.kron(𝐒₁₊╱𝟎, 𝐒₁₊╱𝟎) * M₂.𝛔) @@ -881,8 +881,8 @@ function calculate_third_order_solution(∇₁::AbstractMatrix{<: Real}, #first 𝐗₃ = 𝐔∇₃ * tmpkron + 𝐔∇₃ * M₃.𝐏₁ₗ̂ * tmpkron * M₃.𝐏₁ᵣ̃ + 𝐔∇₃ * M₃.𝐏₂ₗ̂ * tmpkron * M₃.𝐏₂ᵣ̃ - end # timeit_debug - @timeit_debug timer "∇₂ & ∇₁₊" begin + # end # timeit_debug + # @timeit_debug timer "∇₂ & ∇₁₊" begin 𝐒₂₊╱𝟎 = choose_matrix_format(𝐒₂₊╱𝟎, density_threshold = 1.0, min_length = 10) @@ -893,55 +893,55 @@ function calculate_third_order_solution(∇₁::AbstractMatrix{<: Real}, #first 𝐒₂₋╱𝟎 = [𝐒₂[i₋,:] ; zeros(size(𝐒₁)[2] - n₋, nₑ₋^2)] - @timeit_debug timer "Step 1" begin + # @timeit_debug timer "Step 1" begin out2 = mat_mult_kron(∇₂, ⎸𝐒₁𝐒₁₋╱𝟏ₑ⎹╱𝐒₁╱𝟏ₑ₋, ⎸𝐒₂k𝐒₁₋╱𝟏ₑ➕𝐒₁𝐒₂₋⎹╱𝐒₂╱𝟎) # this help - end # timeit_debug - @timeit_debug timer "Step 2" begin + # end # timeit_debug + # @timeit_debug timer "Step 2" begin out2 += ∇₂ * tmpkron1 * tmpkron2# |> findnz - end # timeit_debug - @timeit_debug timer "Step 3" begin + # end # timeit_debug + # @timeit_debug timer "Step 3" begin out2 += ∇₂ * tmpkron1 * M₃.𝐏₁ₗ * tmpkron2 * M₃.𝐏₁ᵣ# |> findnz - end # timeit_debug - @timeit_debug timer "Step 4" begin + # end # timeit_debug + # @timeit_debug timer "Step 4" begin # out2 += ∇₂ * ℒ.kron(⎸𝐒₁𝐒₁₋╱𝟏ₑ⎹╱𝐒₁╱𝟏ₑ₋, 𝐒₂₊╱𝟎 * M₂.𝛔)# |> findnz out2 += mat_mult_kron(∇₂, ⎸𝐒₁𝐒₁₋╱𝟏ₑ⎹╱𝐒₁╱𝟏ₑ₋, collect(𝐒₂₊╱𝟎 * M₂.𝛔))# |> findnz - end # timeit_debug - @timeit_debug timer "Step 5" begin + # end # timeit_debug + # @timeit_debug timer "Step 5" begin # out2 += ∇₁₊ * mat_mult_kron(𝐒₂, collect(𝐒₁₋╱𝟏ₑ), collect(𝐒₂₋╱𝟎)) # out2 += mat_mult_kron(∇₁₊ * 𝐒₂, collect(𝐒₁₋╱𝟏ₑ), collect(𝐒₂₋╱𝟎)) # out2 += ∇₁₊ * 𝐒₂ * ℒ.kron(𝐒₁₋╱𝟏ₑ, 𝐒₂₋╱𝟎) 𝐒₁₋╱𝟏ₑ = choose_matrix_format(𝐒₁₋╱𝟏ₑ, density_threshold = 0.0) out2 += ∇₁₊ * mat_mult_kron(𝐒₂, 𝐒₁₋╱𝟏ₑ, 𝐒₂₋╱𝟎) - end # timeit_debug - @timeit_debug timer "Mult" begin + # end # timeit_debug + # @timeit_debug timer "Mult" begin 𝐗₃ += out2 * M₃.𝐏 𝐗₃ *= M₃.𝐂₃ - end # timeit_debug - end # timeit_debug - @timeit_debug timer "3rd Kronecker power" begin + # end # timeit_debug + # end # timeit_debug + # @timeit_debug timer "3rd Kronecker power" begin # 𝐗₃ += mat_mult_kron(∇₃, collect(aux), collect(ℒ.kron(aux, aux)), M₃.𝐂₃) # slower than direct compression 𝐗₃ += ∇₃ * compressed_kron³(aux, rowmask = unique(findnz(∇₃)[2]), timer = timer) - end # timeit_debug - @timeit_debug timer "Mult 2" begin + # end # timeit_debug + # @timeit_debug timer "Mult 2" begin C = ∇₁₊𝐒₁➕∇₁₀lu \ 𝐗₃# * M₃.𝐂₃ - end # timeit_debug - end # timeit_debug - @timeit_debug timer "Solve sylvester equation" begin + # end # timeit_debug + # end # timeit_debug + # @timeit_debug timer "Solve sylvester equation" begin 𝐒₃, solved = solve_sylvester_equation(A, B, C, sylvester_algorithm = sylvester_algorithm, @@ -950,8 +950,8 @@ function calculate_third_order_solution(∇₁::AbstractMatrix{<: Real}, #first # tol = tol, timer = timer) - end # timeit_debug - # @timeit_debug timer "Refine sylvester equation" begin + # end # timeit_debug + # # @timeit_debug timer "Refine sylvester equation" begin # if !solved # 𝐒₃, solved = solve_sylvester_equation(A, B, C, @@ -966,15 +966,15 @@ function calculate_third_order_solution(∇₁::AbstractMatrix{<: Real}, #first return 𝐒₃, solved end - # end # timeit_debug - @timeit_debug timer "Post-process" begin + # # end # timeit_debug + # @timeit_debug timer "Post-process" begin # 𝐒₃ *= M₃.𝐔₃ # 𝐒₃ = sparse(𝐒₃) - end # timeit_debug - end # timeit_debug + # end # timeit_debug + # end # timeit_debug return 𝐒₃, solved end @@ -993,10 +993,10 @@ function rrule(::typeof(calculate_third_order_solution), T::timings, initial_guess::AbstractMatrix{Float64} = zeros(0,0), sylvester_algorithm::Symbol = :bicgstab, - timer::TimerOutput = TimerOutput(), + # timer::TimerOutput = TimerOutput(), tol::AbstractFloat = eps(), verbose::Bool = false) - @timeit_debug timer "Third order solution - forward" begin + # @timeit_debug timer "Third order solution - forward" begin # inspired by Levintal # Indices and number of variables @@ -1009,7 +1009,7 @@ function rrule(::typeof(calculate_third_order_solution), n = T.nVars nₑ₋ = n₋ + 1 + nₑ - @timeit_debug timer "Setup matrices" begin + # @timeit_debug timer "Setup matrices" begin # 1st order solution 𝐒₁ = @views [𝑺₁[:,1:n₋] zeros(n) 𝑺₁[:,n₋+1:end]]# |> sparse @@ -1028,8 +1028,8 @@ function rrule(::typeof(calculate_third_order_solution), ∇₁₊𝐒₁➕∇₁₀ = @views -∇₁[:,1:n₊] * 𝐒₁[i₊,1:n₋] * ℒ.I(n)[i₋,:] - ∇₁[:,range(1,n) .+ n₊] - end # timeit_debug - @timeit_debug timer "Invert matrix" begin + # end # timeit_debug + # @timeit_debug timer "Invert matrix" begin ∇₁₊𝐒₁➕∇₁₀lu = ℒ.lu(∇₁₊𝐒₁➕∇₁₀, check = false) @@ -1041,7 +1041,7 @@ function rrule(::typeof(calculate_third_order_solution), spinv = inv(∇₁₊𝐒₁➕∇₁₀lu) spinv = choose_matrix_format(spinv) - end # timeit_debug + # end # timeit_debug ∇₁₊ = @views ∇₁[:,1:n₊] * ℒ.I(n)[i₊,:] @@ -1051,36 +1051,36 @@ function rrule(::typeof(calculate_third_order_solution), tmpkron = choose_matrix_format(ℒ.kron(𝐒₁₋╱𝟏ₑ,M₂.𝛔), density_threshold = 1.0) kron𝐒₁₋╱𝟏ₑ = ℒ.kron(𝐒₁₋╱𝟏ₑ,𝐒₁₋╱𝟏ₑ) - @timeit_debug timer "Setup B" begin - @timeit_debug timer "Add tmpkron" begin + # @timeit_debug timer "Setup B" begin + # @timeit_debug timer "Add tmpkron" begin B = tmpkron - end # timeit_debug - @timeit_debug timer "Step 1" begin + # end # timeit_debug + # @timeit_debug timer "Step 1" begin B += M₃.𝐏₁ₗ̄ * tmpkron * M₃.𝐏₁ᵣ̃ - end # timeit_debug - @timeit_debug timer "Step 2" begin + # end # timeit_debug + # @timeit_debug timer "Step 2" begin B += M₃.𝐏₂ₗ̄ * tmpkron * M₃.𝐏₂ᵣ̃ - end # timeit_debug - @timeit_debug timer "Mult" begin + # end # timeit_debug + # @timeit_debug timer "Mult" begin B *= M₃.𝐂₃ B = choose_matrix_format(M₃.𝐔₃ * B) - end # timeit_debug - @timeit_debug timer "3rd Kronecker power" begin + # end # timeit_debug + # @timeit_debug timer "3rd Kronecker power" begin B += compressed_kron³(𝐒₁₋╱𝟏ₑ, timer = timer) - end # timeit_debug - end # timeit_debug - @timeit_debug timer "Setup C" begin - @timeit_debug timer "Initialise smaller matrices" begin + # end # timeit_debug + # end # timeit_debug + # @timeit_debug timer "Setup C" begin + # @timeit_debug timer "Initialise smaller matrices" begin ⎸𝐒₂k𝐒₁₋╱𝟏ₑ➕𝐒₁𝐒₂₋⎹╱𝐒₂╱𝟎 = @views [(𝐒₂ * kron𝐒₁₋╱𝟏ₑ + 𝐒₁ * [𝐒₂[i₋,:] ; zeros(nₑ + 1, nₑ₋^2)])[i₊,:] 𝐒₂ @@ -1093,8 +1093,8 @@ function rrule(::typeof(calculate_third_order_solution), aux = M₃.𝐒𝐏 * ⎸𝐒₁𝐒₁₋╱𝟏ₑ⎹╱𝐒₁╱𝟏ₑ₋ - end # timeit_debug - @timeit_debug timer "∇₃" begin + # end # timeit_debug + # @timeit_debug timer "∇₃" begin tmpkron0 = ℒ.kron(𝐒₁₊╱𝟎, 𝐒₁₊╱𝟎) tmpkron22 = ℒ.kron(⎸𝐒₁𝐒₁₋╱𝟏ₑ⎹╱𝐒₁╱𝟏ₑ₋, tmpkron0 * M₂.𝛔) @@ -1103,8 +1103,8 @@ function rrule(::typeof(calculate_third_order_solution), 𝐗₃ = 𝐔∇₃ * tmpkron22 + 𝐔∇₃ * M₃.𝐏₁ₗ̂ * tmpkron22 * M₃.𝐏₁ᵣ̃ + 𝐔∇₃ * M₃.𝐏₂ₗ̂ * tmpkron22 * M₃.𝐏₂ᵣ̃ - end # timeit_debug - @timeit_debug timer "∇₂ & ∇₁₊" begin + # end # timeit_debug + # @timeit_debug timer "∇₂ & ∇₁₊" begin 𝐒₂₊╱𝟎 = choose_matrix_format(𝐒₂₊╱𝟎, density_threshold = 1.0, min_length = 10) @@ -1117,57 +1117,57 @@ function rrule(::typeof(calculate_third_order_solution), 𝐒₂₋╱𝟎 = choose_matrix_format(𝐒₂₋╱𝟎, density_threshold = 1.0, min_length = 10) - @timeit_debug timer "Step 1" begin + # @timeit_debug timer "Step 1" begin # tmpkron10 = ℒ.kron(⎸𝐒₁𝐒₁₋╱𝟏ₑ⎹╱𝐒₁╱𝟏ₑ₋, ⎸𝐒₂k𝐒₁₋╱𝟏ₑ➕𝐒₁𝐒₂₋⎹╱𝐒₂╱𝟎) # out2 = ∇₂ * tmpkron10 out2 = mat_mult_kron(∇₂, ⎸𝐒₁𝐒₁₋╱𝟏ₑ⎹╱𝐒₁╱𝟏ₑ₋, ⎸𝐒₂k𝐒₁₋╱𝟏ₑ➕𝐒₁𝐒₂₋⎹╱𝐒₂╱𝟎) # this help - end # timeit_debug - @timeit_debug timer "Step 2" begin + # end # timeit_debug + # @timeit_debug timer "Step 2" begin out2 += ∇₂ * tmpkron1 * tmpkron2# |> findnz - end # timeit_debug - @timeit_debug timer "Step 3" begin + # end # timeit_debug + # @timeit_debug timer "Step 3" begin out2 += ∇₂ * tmpkron1 * M₃.𝐏₁ₗ * tmpkron2 * M₃.𝐏₁ᵣ# |> findnz - end # timeit_debug - @timeit_debug timer "Step 4" begin + # end # timeit_debug + # @timeit_debug timer "Step 4" begin 𝐒₂₊╱𝟎𝛔 = 𝐒₂₊╱𝟎 * M₂.𝛔 tmpkron11 = ℒ.kron(⎸𝐒₁𝐒₁₋╱𝟏ₑ⎹╱𝐒₁╱𝟏ₑ₋, 𝐒₂₊╱𝟎𝛔) out2 += ∇₂ * tmpkron11# |> findnz - end # timeit_debug - @timeit_debug timer "Step 5" begin + # end # timeit_debug + # @timeit_debug timer "Step 5" begin tmpkron12 = ℒ.kron(𝐒₁₋╱𝟏ₑ, 𝐒₂₋╱𝟎) out2 += ∇₁₊ * 𝐒₂ * tmpkron12 - end # timeit_debug - @timeit_debug timer "Mult" begin + # end # timeit_debug + # @timeit_debug timer "Mult" begin 𝐗₃ += out2 * M₃.𝐏 𝐗₃ *= M₃.𝐂₃ - end # timeit_debug - end # timeit_debug - @timeit_debug timer "3rd Kronecker power aux" begin + # end # timeit_debug + # end # timeit_debug + # @timeit_debug timer "3rd Kronecker power aux" begin 𝐗₃ += ∇₃ * compressed_kron³(aux, rowmask = unique(findnz(∇₃)[2]), timer = timer) 𝐗₃ = choose_matrix_format(𝐗₃, density_threshold = 1.0, min_length = 10) - end # timeit_debug - @timeit_debug timer "Mult 2" begin + # end # timeit_debug + # @timeit_debug timer "Mult 2" begin C = spinv * 𝐗₃ - end # timeit_debug - end # timeit_debug - @timeit_debug timer "Solve sylvester equation" begin + # end # timeit_debug + # end # timeit_debug + # @timeit_debug timer "Solve sylvester equation" begin 𝐒₃, solved = solve_sylvester_equation(A, B, C, sylvester_algorithm = sylvester_algorithm, @@ -1176,8 +1176,8 @@ function rrule(::typeof(calculate_third_order_solution), # tol = tol, timer = timer) - end # timeit_debug - # @timeit_debug timer "Refine sylvester equation" begin + # end # timeit_debug + # # @timeit_debug timer "Refine sylvester equation" begin # if !solved # 𝐒₃, solved = solve_sylvester_equation(A, B, C, @@ -1194,9 +1194,9 @@ function rrule(::typeof(calculate_third_order_solution), 𝐒₃ = sparse(𝐒₃) - # end # timeit_debug + # # end # timeit_debug - @timeit_debug timer "Preallocate for pullback" begin + # @timeit_debug timer "Preallocate for pullback" begin # At = choose_matrix_format(A')# , density_threshold = 1.0) @@ -1246,8 +1246,8 @@ function rrule(::typeof(calculate_third_order_solution), tmpkron10t = ℒ.kron(⎸𝐒₁𝐒₁₋╱𝟏ₑ⎹╱𝐒₁╱𝟏ₑ₋t, ⎸𝐒₂k𝐒₁₋╱𝟏ₑ➕𝐒₁𝐒₂₋⎹╱𝐒₂╱𝟎t) - end # timeit_debug - end # timeit_debug + # end # timeit_debug + # end # timeit_debug function third_order_solution_pullback(∂𝐒₃_solved) ∂∇₁ = zero(∇₁) @@ -1272,9 +1272,9 @@ function rrule(::typeof(calculate_third_order_solution), ∂∇₁₊ = zero(∇₁₊) ∂𝐒₂₋╱𝟎 = zero(𝐒₂₋╱𝟎) - @timeit_debug timer "Third order solution - pullback" begin + # @timeit_debug timer "Third order solution - pullback" begin - @timeit_debug timer "Solve sylvester equation" begin + # @timeit_debug timer "Solve sylvester equation" begin ∂𝐒₃ = ∂𝐒₃_solved[1] @@ -1292,16 +1292,16 @@ function rrule(::typeof(calculate_third_order_solution), ∂C = choose_matrix_format(∂C, density_threshold = 1.0) - end # timeit_debug - @timeit_debug timer "Step 0" begin + # end # timeit_debug + # @timeit_debug timer "Step 0" begin ∂A = ∂C * B' * 𝐒₃' # ∂B = 𝐒₃' * A' * ∂C ∂B = choose_matrix_format(𝐒₃' * A' * ∂C, density_threshold = 1.0) - end # timeit_debug - @timeit_debug timer "Step 1" begin + # end # timeit_debug + # @timeit_debug timer "Step 1" begin # C = spinv * 𝐗₃ # ∂𝐗₃ = spinv' * ∂C * M₃.𝐂₃' @@ -1324,8 +1324,8 @@ function rrule(::typeof(calculate_third_order_solution), # tmpkron12 = ℒ.kron(𝐒₁₋╱𝟏ₑ, 𝐒₂₋╱𝟎) fill_kron_adjoint!(∂𝐒₁₋╱𝟏ₑ, ∂𝐒₂₋╱𝟎, ∂tmpkron12, 𝐒₁₋╱𝟏ₑ, 𝐒₂₋╱𝟎) - end # timeit_debug - @timeit_debug timer "Step 2" begin + # end # timeit_debug + # @timeit_debug timer "Step 2" begin # ∇₂ * (tmpkron10 + tmpkron1 * tmpkron2 + tmpkron1 * M₃.𝐏₁ₗ * tmpkron2 * M₃.𝐏₁ᵣ + tmpkron11) * M₃.𝐏 * M₃.𝐂₃ #improve this @@ -1348,8 +1348,8 @@ function rrule(::typeof(calculate_third_order_solution), ∂tmpkron10 = ∇₂t * ∂𝐗₃ * 𝐂₃t * 𝐏t - end # timeit_debug - @timeit_debug timer "Step 3" begin + # end # timeit_debug + # @timeit_debug timer "Step 3" begin # tmpkron10 = ℒ.kron(⎸𝐒₁𝐒₁₋╱𝟏ₑ⎹╱𝐒₁╱𝟏ₑ₋, ⎸𝐒₂k𝐒₁₋╱𝟏ₑ➕𝐒₁𝐒₂₋⎹╱𝐒₂╱𝟎) fill_kron_adjoint!(∂⎸𝐒₁𝐒₁₋╱𝟏ₑ⎹╱𝐒₁╱𝟏ₑ₋, ∂⎸𝐒₂k𝐒₁₋╱𝟏ₑ➕𝐒₁𝐒₂₋⎹╱𝐒₂╱𝟎, ∂tmpkron10, ⎸𝐒₁𝐒₁₋╱𝟏ₑ⎹╱𝐒₁╱𝟏ₑ₋, ⎸𝐒₂k𝐒₁₋╱𝟏ₑ➕𝐒₁𝐒₂₋⎹╱𝐒₂╱𝟎) @@ -1379,8 +1379,8 @@ function rrule(::typeof(calculate_third_order_solution), ∂𝐒₂₊╱𝟎 += ∂𝐒₂₊╱𝟎𝛔 * 𝛔t - end # timeit_debug - @timeit_debug timer "Step 4" begin + # end # timeit_debug + # @timeit_debug timer "Step 4" begin # out = (𝐔∇₃ * tmpkron22 # + 𝐔∇₃ * M₃.𝐏₁ₗ̂ * tmpkron22 * M₃.𝐏₁ᵣ̃ @@ -1402,15 +1402,15 @@ function rrule(::typeof(calculate_third_order_solution), # A_mult_kron_power_3_B!(∂∇₃, ∂𝐗₃, aux') # not a good idea because filling an existing matrix one by one is slow # ∂∇₃ += A_mult_kron_power_3_B(∂𝐗₃, aux') # this is slower somehow - end # timeit_debug - @timeit_debug timer "Step 5" begin + # end # timeit_debug + # @timeit_debug timer "Step 5" begin # this is very slow ∂∇₃ += ∂𝐗₃ * compressed_kron³(aux', rowmask = unique(findnz(∂𝐗₃)[2]), timer = timer) # ∂∇₃ += ∂𝐗₃ * ℒ.kron(aux', aux', aux') - end # timeit_debug - @timeit_debug timer "Step 6" begin + # end # timeit_debug + # @timeit_debug timer "Step 6" begin ∂kronkronaux = 𝐔∇₃t * ∂𝐗₃ * 𝐂₃t @@ -1418,8 +1418,8 @@ function rrule(::typeof(calculate_third_order_solution), fill_kron_adjoint!(∂aux, ∂aux, ∂kronaux, aux, aux) - end # timeit_debug - @timeit_debug timer "Step 7" begin + # end # timeit_debug + # @timeit_debug timer "Step 7" begin # aux = M₃.𝐒𝐏 * ⎸𝐒₁𝐒₁₋╱𝟏ₑ⎹╱𝐒₁╱𝟏ₑ₋ ∂⎸𝐒₁𝐒₁₋╱𝟏ₑ⎹╱𝐒₁╱𝟏ₑ₋ += M₃.𝐒𝐏' * ∂aux @@ -1456,8 +1456,8 @@ function rrule(::typeof(calculate_third_order_solution), ∂𝐒₂╱𝟎 = 𝐒₁' * ∂𝐒₂k𝐒₁₋╱𝟏ₑ ∂𝐒₂[i₋,:] += ∂𝐒₂╱𝟎[1:length(i₋),:] - end # timeit_debug - @timeit_debug timer "Step 8" begin + # end # timeit_debug + # @timeit_debug timer "Step 8" begin ### # B = M₃.𝐔₃ * (tmpkron + M₃.𝐏₁ₗ̄ * tmpkron * M₃.𝐏₁ᵣ̃ + M₃.𝐏₂ₗ̄ * tmpkron * M₃.𝐏₂ᵣ̃ + ℒ.kron(𝐒₁₋╱𝟏ₑ, kron𝐒₁₋╱𝟏ₑ)) * M₃.𝐂₃ @@ -1508,8 +1508,8 @@ function rrule(::typeof(calculate_third_order_solution), # 𝐒₁ = [𝑺₁[:,1:n₋] zeros(n) 𝑺₁[:,n₋+1:end]] ∂𝑺₁ = [∂𝐒₁[:,1:n₋] ∂𝐒₁[:,n₋+2:end]] - end # timeit_debug - end # timeit_debug + # end # timeit_debug + # end # timeit_debug return NoTangent(), ∂∇₁, ∂∇₂, ∂∇₃, ∂𝑺₁, ∂𝐒₂, NoTangent(), NoTangent(), NoTangent(), NoTangent(), NoTangent(), NoTangent() end