From 8942f3b879b4f05153e8ae2b33365abe9213d1e8 Mon Sep 17 00:00:00 2001 From: FOJ-0 Date: Mon, 10 Apr 2023 12:24:56 -0400 Subject: [PATCH 1/4] update analytical derivatives --- test/test_analytical_derivatives.jl | 37 +++++++++++++++++++++++++---- 1 file changed, 32 insertions(+), 5 deletions(-) diff --git a/test/test_analytical_derivatives.jl b/test/test_analytical_derivatives.jl index 05ac7773..4ab112b7 100644 --- a/test/test_analytical_derivatives.jl +++ b/test/test_analytical_derivatives.jl @@ -1,9 +1,36 @@ -using Dynare +module TestDerivatives +using Dynare #Use analytical_derivatives branch +using LinearAlgebra +#Model context = @dynare "models/analytical_derivatives/fs2000_sa.mod" "params_derivs_order=1" "notmpterms"; - +model = context.models[1] +results = context.results.model_results[1] +wss = Dynare.StaticWs(context) params = context.work.params -steady_state = context.results.model_results[1].trends.endogenous_steady_state -exogenous = context.results.model_results[1].trends.exogenous_steady_state +steadystate = results.trends.endogenous_steady_state +endogenous = repeat(steadystate, 3) +exogenous = results.trends.exogenous_steady_state + +#Get StaticJacobian +df_dx = Dynare.get_static_jacobian!(wss, params, steadystate, exogenous, model) + +#Get StaticJacobianParams +(df_dp, gp) = Dynare.DFunctions.SparseStaticParametersDerivatives!(steadystate, exogenous, params) + +#Get Jacobian matrix +function Derivatives(df_dx, df_dp) + dense_df_dx = Matrix(df_dx) + + if det(dense_df_dx) == 0 + error("The matrix is not invertible") + end + + #Get dx_dp + df_dx_inv = inv(dense_df_dx) + return -df_dx_inv * df_dp +end + +@show Derivatives(df_dx, df_dp) -(rp, gp) = Dynare.DFunctions.SparseStaticParametersDerivatives!(steady_state, exogenous, params) +end # end module From b632fc60d853bf592aa539a37b9fe9dc4d40ce68 Mon Sep 17 00:00:00 2001 From: FOJ-0 Date: Thu, 13 Apr 2023 09:14:42 -0400 Subject: [PATCH 2/4] add new methods --- test/test_analytical_derivatives.jl | 65 +++++++++++++++++++++++++++-- 1 file changed, 62 insertions(+), 3 deletions(-) diff --git a/test/test_analytical_derivatives.jl b/test/test_analytical_derivatives.jl index 4ab112b7..5ca60cd2 100644 --- a/test/test_analytical_derivatives.jl +++ b/test/test_analytical_derivatives.jl @@ -1,6 +1,8 @@ module TestDerivatives using Dynare #Use analytical_derivatives branch using LinearAlgebra +using FastLapackInterface +using Test #Model context = @dynare "models/analytical_derivatives/fs2000_sa.mod" "params_derivs_order=1" "notmpterms"; @@ -18,7 +20,7 @@ df_dx = Dynare.get_static_jacobian!(wss, params, steadystate, exogenous, model) #Get StaticJacobianParams (df_dp, gp) = Dynare.DFunctions.SparseStaticParametersDerivatives!(steadystate, exogenous, params) -#Get Jacobian matrix +#Get Jacobian matrix using the inverse of df_dx function Derivatives(df_dx, df_dp) dense_df_dx = Matrix(df_dx) @@ -27,10 +29,67 @@ function Derivatives(df_dx, df_dp) end #Get dx_dp - df_dx_inv = inv(dense_df_dx) + df_dx_inv = dense_df_dx \ I#inv(dense_df_dx) return -df_dx_inv * df_dp end -@show Derivatives(df_dx, df_dp) +#Get Jacobian matrix using FastLapackInterface (dense matrix) +function Derivatives2(df_dx, df_dp) + dense_df_dx = Matrix(df_dx) + LU_df_dx = LU(LAPACK.getrf!(LUWs(dense_df_dx), dense_df_dx)...) + if any(diag(LU_df_dx.U) .== 0) + error("The matrix is not invertible") + end + + #Get dx_dp + dx_dp = -(LU_df_dx \ df_dp) + return dx_dp +end + +#Get Jacobian matrix using lu(df_dx) (sparse matrix) +function Derivatives3(df_dx, df_dp) + LU_df_dx = lu(df_dx) + + if any(diag(LU_df_dx.U) .== 0) + error("The matrix is not invertible") + end + + # Get dx_dp + dx_dp = -(LU_df_dx \ df_dp) + return dx_dp +end + +#Get Jacobian matrix using ldiv! (dense matrix) +function Derivatives4(df_dx, df_dp) + dense_df_dx = Matrix(df_dx) + LU_df_dx = lu(dense_df_dx) + + if any(diag(LU_df_dx.U) .== 0) + error("The matrix is not invertible") + end + + # Preallocate dx_dp + n, m = size(df_dp) + dx_dp = Matrix{Float64}(undef, n, m) + + # Get dx_dp + for i in 1:m + dp_col = view(df_dp, :, i) + sol_col = view(dx_dp, :, i) + ldiv!(sol_col, LU_df_dx, dp_col) + sol_col .*= -1 + end + + return dx_dp +end + +# @show Derivatives(df_dx, df_dp); +@time sol1 = Derivatives(df_dx, df_dp); +@time sol2 = Derivatives2(df_dx, df_dp); +@time sol3 = Derivatives3(df_dx, df_dp); +@time sol4 = Derivatives4(df_dx, df_dp); +@test sol1 ≈ sol2 +@test sol1 ≈ sol3 +@test sol1 ≈ sol4 end # end module From be3a9156e8085c5042f675c0f6adcefd8bbe75d5 Mon Sep 17 00:00:00 2001 From: FOJ-0 Date: Thu, 20 Apr 2023 12:23:23 -0400 Subject: [PATCH 3/4] structure --- test/test_analytical_derivatives.jl | 136 +++++++++++++++++++++++----- 1 file changed, 112 insertions(+), 24 deletions(-) diff --git a/test/test_analytical_derivatives.jl b/test/test_analytical_derivatives.jl index 5ca60cd2..4e35e08e 100644 --- a/test/test_analytical_derivatives.jl +++ b/test/test_analytical_derivatives.jl @@ -1,11 +1,15 @@ module TestDerivatives +cd("/home/felix/Dropbox/Dynare/2_Dynare_Fork/Dynare.jl") using Dynare #Use analytical_derivatives branch using LinearAlgebra using FastLapackInterface using Test +using SparseArrays +using SuiteSparse #Model context = @dynare "models/analytical_derivatives/fs2000_sa.mod" "params_derivs_order=1" "notmpterms"; +context = @dynare "/home/felix/Dropbox/Dynare/2_Dynare_Fork/Dynare.jl/test/models/analytical_derivatives/fs2000_sa.mod" "params_derivs_order=1" "notmpterms"; model = context.models[1] results = context.results.model_results[1] wss = Dynare.StaticWs(context) @@ -20,59 +24,133 @@ df_dx = Dynare.get_static_jacobian!(wss, params, steadystate, exogenous, model) #Get StaticJacobianParams (df_dp, gp) = Dynare.DFunctions.SparseStaticParametersDerivatives!(steadystate, exogenous, params) -#Get Jacobian matrix using the inverse of df_dx -function Derivatives(df_dx, df_dp) - dense_df_dx = Matrix(df_dx) + +#1. Get Jacobian matrix using the inverse of df_dx +struct DerivativesWorkspace + dense_df_dx::Matrix{Float64} + df_dx_inv::Matrix{Float64} + result::Matrix{Float64} +end + +function init_derivatives_workspace(n::Int, m::Int) + dense_df_dx = Matrix{Float64}(undef, n, n) + df_dx_inv = Matrix{Float64}(undef, n, n) + result = Matrix{Float64}(undef, n, m) + return DerivativesWorkspace(dense_df_dx, df_dx_inv, result) +end + +function Derivatives(workspace::DerivativesWorkspace, df_dx, df_dp) + dense_df_dx = workspace.dense_df_dx + copyto!(dense_df_dx, df_dx) if det(dense_df_dx) == 0 error("The matrix is not invertible") end - #Get dx_dp - df_dx_inv = dense_df_dx \ I#inv(dense_df_dx) - return -df_dx_inv * df_dp + df_dx_inv = workspace.df_dx_inv + df_dx_inv .= dense_df_dx \ I + + result = workspace.result + return mul!(result, -df_dx_inv, df_dp) +end + + +#2. Get Jacobian matrix using FastLapackInterface (dense matrix) +struct Derivatives2Workspace + dense_df_dx::Matrix{Float64} + LU_df_dx::LU{Float64, Matrix{Float64}} #not used yet + dx_dp::Matrix{Float64} end -#Get Jacobian matrix using FastLapackInterface (dense matrix) -function Derivatives2(df_dx, df_dp) - dense_df_dx = Matrix(df_dx) + +function init_derivatives2_workspace(n::Int, m::Int) + dense_df_dx = Matrix{Float64}(undef, n, n) + LU_df_dx = lu(Matrix{Float64}(I, n, n)) #not used yet + dx_dp = Matrix{Float64}(undef, n, m) + return Derivatives2Workspace(dense_df_dx, LU_df_dx, dx_dp) +end + + +function Derivatives2(workspace::Derivatives2Workspace, df_dx, df_dp) + dense_df_dx = workspace.dense_df_dx + copyto!(dense_df_dx, df_dx) + LU_df_dx = LU(LAPACK.getrf!(LUWs(dense_df_dx), dense_df_dx)...) + if any(diag(LU_df_dx.U) .== 0) error("The matrix is not invertible") end - #Get dx_dp - dx_dp = -(LU_df_dx \ df_dp) + dx_dp = workspace.dx_dp + copyto!(dx_dp, -(LU_df_dx \ df_dp)) + # ldiv!(-1.0, LU_df_dx, df_dp, dx_dp) + return dx_dp end -#Get Jacobian matrix using lu(df_dx) (sparse matrix) -function Derivatives3(df_dx, df_dp) +#3. Get Jacobian matrix using lu(df_dx) (sparse matrix) +struct Derivatives3Workspace + # df_dx::SparseMatrixCSC{Float64, Int64} + LU_df_dx::SuiteSparse.UMFPACK.UmfpackLU{Float64, Int64} #not used yet + dx_dp::Matrix{Float64} +end + +function init_derivatives3_workspace(n::Int, m::Int) + LU_df_dx = SuiteSparse.UMFPACK.lu(sparse(Matrix{Float64}(I, n, n))) #not used yet + dx_dp = Matrix{Float64}(undef, n, m) + return Derivatives3Workspace(LU_df_dx, dx_dp) +end + +function Derivatives3(workspace::Derivatives3Workspace, df_dx, df_dp) + # Perform Sparse LU decomposition + # LU_df_dx = workspace.LU_df_dx LU_df_dx = lu(df_dx) + # copyto!(LU_df_dx.L, lu_df_dx.L) + # copyto!(LU_df_dx.U, lu_df_dx.U) if any(diag(LU_df_dx.U) .== 0) error("The matrix is not invertible") end # Get dx_dp - dx_dp = -(LU_df_dx \ df_dp) + dx_dp = workspace.dx_dp + copyto!(dx_dp, -(LU_df_dx \ df_dp)) + # ldiv!(-1.0, LU_df_dx, df_dp, dx_dp) + return dx_dp end -#Get Jacobian matrix using ldiv! (dense matrix) -function Derivatives4(df_dx, df_dp) - dense_df_dx = Matrix(df_dx) + +#4. Get Jacobian matrix using ldiv! (dense matrix) +struct Derivatives4Workspace + dense_df_dx::Matrix{Float64} + LU_df_dx::LU{Float64, Matrix{Float64}} + dx_dp::Matrix{Float64} +end + +function init_derivatives4_workspace(n::Int, m::Int) + dense_df_dx = Matrix{Float64}(undef, n, n) + LU_df_dx = lu(Matrix{Float64}(I, n, n)) + dx_dp = Matrix{Float64}(undef, n, m) + return Derivatives4Workspace(dense_df_dx, LU_df_dx, dx_dp) +end + +function Derivatives4(workspace::Derivatives4Workspace, df_dx, df_dp) + dense_df_dx = workspace.dense_df_dx + copyto!(dense_df_dx, Matrix(df_dx)) + + # LU_df_dx = workspace.LU_df_dx LU_df_dx = lu(dense_df_dx) + # copyto!(LU_df_dx.L, lu_df_dx.L) + # copyto!(LU_df_dx.U, lu_df_dx.U) if any(diag(LU_df_dx.U) .== 0) error("The matrix is not invertible") end - # Preallocate dx_dp + dx_dp = workspace.dx_dp n, m = size(df_dp) - dx_dp = Matrix{Float64}(undef, n, m) - # Get dx_dp for i in 1:m dp_col = view(df_dp, :, i) sol_col = view(dx_dp, :, i) @@ -84,11 +162,21 @@ function Derivatives4(df_dx, df_dp) end +# Initialize +n = model.endogenous_nbr +m = model.parameter_nbr +workspace1 = init_derivatives_workspace(n, m) +workspace2 = init_derivatives2_workspace(n, m) +workspace3 = init_derivatives3_workspace(n, m) +workspace4 = init_derivatives4_workspace(n, m) + +# Luego puedes utilizar workspace para llamar a la función Derivatives con tus matrices df_dx y df_dp +@time sol1 = Derivatives(workspace1, df_dx, df_dp); +@time sol2 = Derivatives2(workspace2, df_dx, df_dp); +@time sol3 = Derivatives3(workspace3, df_dx, df_dp); +@time sol4 = Derivatives4(workspace4, df_dx, df_dp); + # @show Derivatives(df_dx, df_dp); -@time sol1 = Derivatives(df_dx, df_dp); -@time sol2 = Derivatives2(df_dx, df_dp); -@time sol3 = Derivatives3(df_dx, df_dp); -@time sol4 = Derivatives4(df_dx, df_dp); @test sol1 ≈ sol2 @test sol1 ≈ sol3 @test sol1 ≈ sol4 From 76b219c3daa0b44f8f83a295b782b3542e136888 Mon Sep 17 00:00:00 2001 From: FOJ-0 Date: Fri, 21 Apr 2023 09:02:52 -0400 Subject: [PATCH 4/4] structure --- test/test_analytical_derivatives.jl | 2 -- 1 file changed, 2 deletions(-) diff --git a/test/test_analytical_derivatives.jl b/test/test_analytical_derivatives.jl index 4e35e08e..0b6c8ae7 100644 --- a/test/test_analytical_derivatives.jl +++ b/test/test_analytical_derivatives.jl @@ -1,5 +1,4 @@ module TestDerivatives -cd("/home/felix/Dropbox/Dynare/2_Dynare_Fork/Dynare.jl") using Dynare #Use analytical_derivatives branch using LinearAlgebra using FastLapackInterface @@ -9,7 +8,6 @@ using SuiteSparse #Model context = @dynare "models/analytical_derivatives/fs2000_sa.mod" "params_derivs_order=1" "notmpterms"; -context = @dynare "/home/felix/Dropbox/Dynare/2_Dynare_Fork/Dynare.jl/test/models/analytical_derivatives/fs2000_sa.mod" "params_derivs_order=1" "notmpterms"; model = context.models[1] results = context.results.model_results[1] wss = Dynare.StaticWs(context)