From 57ef2e13553ee445e7026b505d31c5115f4fb073 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Mon, 22 Apr 2024 08:57:24 -0700 Subject: [PATCH 01/29] obj!()/grad!(): avoid tmp array creation --- src/imply/RAM/generic.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/imply/RAM/generic.jl b/src/imply/RAM/generic.jl index b9d6fad69..907d24103 100644 --- a/src/imply/RAM/generic.jl +++ b/src/imply/RAM/generic.jl @@ -232,7 +232,8 @@ function objective!(imply::RAM, parameters, model, has_meanstructure::Val{T}) wh parameters, ) - imply.I_A .= I - imply.A + @. imply.I_A = -imply.A + @view(imply.I_A[diagind(imply.I_A)]) .+= 1 copyto!(imply.F⨉I_A⁻¹, imply.F) rdiv!(imply.F⨉I_A⁻¹, factorize(imply.I_A)) @@ -260,8 +261,8 @@ function gradient!( parameters, ) - imply.I_A .= I - imply.A - copyto!(imply.I_A⁻¹, imply.I_A) + @. imply.I_A = -imply.A + @view(imply.I_A[diagind(imply.I_A)]) .+= 1 imply.I_A⁻¹ .= LinearAlgebra.inv!(factorize(imply.I_A⁻¹)) imply.F⨉I_A⁻¹ .= imply.F * imply.I_A⁻¹ From c4ecdc49cb211c28362c0b1b597d4ec25957a18f Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sat, 9 Mar 2024 15:38:07 -0800 Subject: [PATCH 02/29] grad!(): avoid extra array copying --- src/imply/RAM/generic.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/imply/RAM/generic.jl b/src/imply/RAM/generic.jl index 907d24103..f240814ab 100644 --- a/src/imply/RAM/generic.jl +++ b/src/imply/RAM/generic.jl @@ -264,8 +264,8 @@ function gradient!( @. imply.I_A = -imply.A @view(imply.I_A[diagind(imply.I_A)]) .+= 1 - imply.I_A⁻¹ .= LinearAlgebra.inv!(factorize(imply.I_A⁻¹)) - imply.F⨉I_A⁻¹ .= imply.F * imply.I_A⁻¹ + imply.I_A⁻¹ = LinearAlgebra.inv!(factorize(imply.I_A)) + mul!(imply.F⨉I_A⁻¹, imply.F, imply.I_A⁻¹) Σ_RAM!(imply.Σ, imply.F⨉I_A⁻¹, imply.S, imply.F⨉I_A⁻¹S) From d877244d5d66aaf6768bf0b5a3456ba700eb42d1 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sat, 9 Mar 2024 15:38:33 -0800 Subject: [PATCH 03/29] check_acyclic(): use istril/u() --- src/imply/RAM/generic.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/imply/RAM/generic.jl b/src/imply/RAM/generic.jl index f240814ab..2089cf3da 100644 --- a/src/imply/RAM/generic.jl +++ b/src/imply/RAM/generic.jl @@ -353,9 +353,9 @@ function check_acyclic(A_pre, n_par, A_indices) acyclic = isone(det(I - A_rand)) # check if A is lower or upper triangular - if iszero(A_rand[.!tril(ones(Bool, size(A_pre)...))]) + if istril(A_rand) A_pre = LowerTriangular(A_pre) - elseif iszero(A_rand[.!tril(ones(Bool, size(A_pre)...))']) + elseif istriu(A_rand) A_pre = UpperTriangular(A_pre) elseif acyclic @info "Your model is acyclic, specifying the A Matrix as either Upper or Lower Triangular can have great performance benefits.\n" maxlog = From 244f61a13eaa4290a67ca372e709f6a8b6e1fb5a Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sat, 9 Mar 2024 15:40:57 -0800 Subject: [PATCH 04/29] grad!(SemML): reduce * ops --- src/loss/ML/ML.jl | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/src/loss/ML/ML.jl b/src/loss/ML/ML.jl index e6616258d..879ac1a7c 100644 --- a/src/loss/ML/ML.jl +++ b/src/loss/ML/ML.jl @@ -531,8 +531,7 @@ function gradient!( μ₋ᵀΣ⁻¹ = μ₋' * Σ⁻¹ k = μ₋ᵀΣ⁻¹ * F⨉I_A⁻¹ - gradient += - -2k * ∇M - 2vec(k'M'I_A⁻¹')'∇A - 2vec(k'k * S * I_A⁻¹')'∇A - vec(k'k)'∇S + gradient .+= -2k * ∇M - 2vec(k' * (M' + k * S) * I_A⁻¹')'∇A - vec(k'k)'∇S end return gradient' @@ -582,8 +581,7 @@ function objective_gradient!( μ₋ᵀΣ⁻¹ = μ₋' * Σ⁻¹ k = μ₋ᵀΣ⁻¹ * F⨉I_A⁻¹ - gradient += - -2k * ∇M - 2vec(k'M'I_A⁻¹')'∇A - 2vec(k'k * S * I_A⁻¹')'∇A - vec(k'k)'∇S + gradient .+= -2k * ∇M - 2vec(k' * (M' + k * S) * I_A⁻¹')'∇A - vec(k'k)'∇S end return objective, gradient' From 580f85d5efa86e6a855ba6d35c5aa4d5031f85ac Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sat, 9 Mar 2024 16:03:22 -0800 Subject: [PATCH 05/29] obj/grad/hess!(SemML): avoid extra arr copying --- src/loss/ML/ML.jl | 60 +++++++++++++---------------------------------- 1 file changed, 16 insertions(+), 44 deletions(-) diff --git a/src/loss/ML/ML.jl b/src/loss/ML/ML.jl index 879ac1a7c..ed9bfad2e 100644 --- a/src/loss/ML/ML.jl +++ b/src/loss/ML/ML.jl @@ -95,13 +95,9 @@ function objective!( copyto!(Σ⁻¹, Σ) Σ_chol = cholesky!(Symmetric(Σ⁻¹); check = false) - - if !isposdef(Σ_chol) - return non_posdef_return(par) - end - + isposdef(Σ_chol) || return non_posdef_return(par) ld = logdet(Σ_chol) - Σ⁻¹ .= LinearAlgebra.inv!(Σ_chol) + Σ⁻¹ = LinearAlgebra.inv!(Σ_chol) #mul!(Σ⁻¹Σₒ, Σ⁻¹, Σₒ) if T @@ -131,12 +127,8 @@ function gradient!( copyto!(Σ⁻¹, Σ) Σ_chol = cholesky!(Symmetric(Σ⁻¹); check = false) - - if !isposdef(Σ_chol) - return ones(eltype(par), size(par)) - end - - Σ⁻¹ .= LinearAlgebra.inv!(Σ_chol) + isposdef(Σ_chol) || return ones(eltype(par), size(par)) + Σ⁻¹ = LinearAlgebra.inv!(Σ_chol) mul!(Σ⁻¹Σₒ, Σ⁻¹, Σₒ) if T @@ -168,12 +160,8 @@ function hessian!( copyto!(Σ⁻¹, Σ) Σ_chol = cholesky!(Symmetric(Σ⁻¹); check = false) - - if !isposdef(Σ_chol) - return diagm(fill(one(eltype(par)), length(par))) - end - - Σ⁻¹ .= LinearAlgebra.inv!(Σ_chol) + isposdef(Σ_chol) || return diagm(fill(one(eltype(par)), length(par))) + Σ⁻¹ = LinearAlgebra.inv!(Σ_chol) if semml.approximate_hessian hessian = 2 * ∇Σ' * kron(Σ⁻¹, Σ⁻¹) * ∇Σ @@ -221,12 +209,11 @@ function objective_gradient!( copyto!(Σ⁻¹, Σ) Σ_chol = cholesky!(Symmetric(Σ⁻¹); check = false) - if !isposdef(Σ_chol) return non_posdef_return(par), ones(eltype(par), size(par)) else ld = logdet(Σ_chol) - Σ⁻¹ .= LinearAlgebra.inv!(Σ_chol) + Σ⁻¹ = LinearAlgebra.inv!(Σ_chol) mul!(Σ⁻¹Σₒ, Σ⁻¹, Σₒ) if T @@ -263,12 +250,11 @@ function objective_hessian!( copyto!(Σ⁻¹, Σ) Σ_chol = cholesky!(Symmetric(Σ⁻¹); check = false) - if !isposdef(Σ_chol) return non_posdef_return(par), diagm(fill(one(eltype(par)), length(par))) else ld = logdet(Σ_chol) - Σ⁻¹ .= LinearAlgebra.inv!(Σ_chol) + Σ⁻¹ = LinearAlgebra.inv!(Σ_chol) mul!(Σ⁻¹Σₒ, Σ⁻¹, Σₒ) objective = ld + tr(Σ⁻¹Σₒ) @@ -318,12 +304,9 @@ function gradient_hessian!( copyto!(Σ⁻¹, Σ) Σ_chol = cholesky!(Symmetric(Σ⁻¹); check = false) - - if !isposdef(Σ_chol) + isposdef(Σ_chol) || return ones(eltype(par), size(par)), diagm(fill(one(eltype(par)), length(par))) - end - - Σ⁻¹ .= LinearAlgebra.inv!(Σ_chol) + Σ⁻¹ = LinearAlgebra.inv!(Σ_chol) mul!(Σ⁻¹Σₒ, Σ⁻¹, Σₒ) Σ⁻¹ΣₒΣ⁻¹ = Σ⁻¹Σₒ * Σ⁻¹ @@ -373,16 +356,14 @@ function objective_gradient_hessian!( copyto!(Σ⁻¹, Σ) Σ_chol = cholesky!(Symmetric(Σ⁻¹); check = false) - if !isposdef(Σ_chol) objective = non_posdef_return(par) gradient = ones(eltype(par), size(par)) hessian = diagm(fill(one(eltype(par)), length(par))) return objective, gradient, hessian end - ld = logdet(Σ_chol) - Σ⁻¹ .= LinearAlgebra.inv!(Σ_chol) + Σ⁻¹ = LinearAlgebra.inv!(Σ_chol) mul!(Σ⁻¹Σₒ, Σ⁻¹, Σₒ) objective = ld + tr(Σ⁻¹Σₒ) @@ -474,13 +455,9 @@ function objective!( copyto!(Σ⁻¹, Σ) Σ_chol = cholesky!(Symmetric(Σ⁻¹); check = false) - - if !isposdef(Σ_chol) - return non_posdef_return(par) - end - + isposdef(Σ_chol) || return non_posdef_return(par) ld = logdet(Σ_chol) - Σ⁻¹ .= LinearAlgebra.inv!(Σ_chol) + Σ⁻¹ = LinearAlgebra.inv!(Σ_chol) mul!(Σ⁻¹Σₒ, Σ⁻¹, Σₒ) if T @@ -515,12 +492,8 @@ function gradient!( copyto!(Σ⁻¹, Σ) Σ_chol = cholesky!(Symmetric(Σ⁻¹); check = false) - - if !isposdef(Σ_chol) - return ones(eltype(par), size(par)) - end - - Σ⁻¹ .= LinearAlgebra.inv!(Σ_chol) + isposdef(Σ_chol) || return ones(eltype(par), size(par)) + Σ⁻¹ = LinearAlgebra.inv!(Σ_chol) #mul!(Σ⁻¹Σₒ, Σ⁻¹, Σₒ) C = F⨉I_A⁻¹' * (I - Σₒ * Σ⁻¹)' * Σ⁻¹ * F⨉I_A⁻¹ @@ -561,14 +534,13 @@ function objective_gradient!( copyto!(Σ⁻¹, Σ) Σ_chol = cholesky!(Symmetric(Σ⁻¹); check = false) - if !isposdef(Σ_chol) objective = non_posdef_return(par) gradient = ones(eltype(par), size(par)) return objective, gradient else ld = logdet(Σ_chol) - Σ⁻¹ .= LinearAlgebra.inv!(Σ_chol) + Σ⁻¹ = LinearAlgebra.inv!(Σ_chol) mul!(Σ⁻¹Σₒ, Σ⁻¹, Σₒ) objective = ld + tr(Σ⁻¹Σₒ) From 9bbd0100e56faac61b5a7b3c58bfcd4ff015c4fc Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sat, 20 Apr 2024 14:38:13 -0700 Subject: [PATCH 06/29] use .+= to reduce allocs --- src/loss/ML/ML.jl | 8 ++++---- src/loss/WLS/WLS.jl | 8 ++++---- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/src/loss/ML/ML.jl b/src/loss/ML/ML.jl index ed9bfad2e..4c9a1f134 100644 --- a/src/loss/ML/ML.jl +++ b/src/loss/ML/ML.jl @@ -174,7 +174,7 @@ function hessian!( # outer H_outer = 2 * kron(Σ⁻¹ΣₒΣ⁻¹, Σ⁻¹) - kron(Σ⁻¹, Σ⁻¹) hessian = ∇Σ' * H_outer * ∇Σ - hessian += ∇²Σ + hessian .+= ∇²Σ end return hessian @@ -268,7 +268,7 @@ function objective_hessian!( # outer H_outer = 2 * kron(Σ⁻¹ΣₒΣ⁻¹, Σ⁻¹) - kron(Σ⁻¹, Σ⁻¹) hessian = ∇Σ' * H_outer * ∇Σ - hessian += ∇²Σ + hessian .+= ∇²Σ end return objective, hessian @@ -322,7 +322,7 @@ function gradient_hessian!( # outer H_outer = 2 * kron(Σ⁻¹ΣₒΣ⁻¹, Σ⁻¹) - kron(Σ⁻¹, Σ⁻¹) hessian = ∇Σ' * H_outer * ∇Σ - hessian += ∇²Σ + hessian .+= ∇²Σ end return gradient', hessian @@ -381,7 +381,7 @@ function objective_gradient_hessian!( # outer H_outer = 2 * kron(Σ⁻¹ΣₒΣ⁻¹, Σ⁻¹) - kron(Σ⁻¹, Σ⁻¹) hessian = ∇Σ' * H_outer * ∇Σ - hessian += ∇²Σ + hessian .+= ∇²Σ end return objective, gradient', hessian diff --git a/src/loss/WLS/WLS.jl b/src/loss/WLS/WLS.jl index 18649ab68..61c89fc85 100644 --- a/src/loss/WLS/WLS.jl +++ b/src/loss/WLS/WLS.jl @@ -180,7 +180,7 @@ function hessian!( if !semwls.approximate_hessian J = -2 * (σ₋' * semwls.V)' ∇²Σ_function!(∇²Σ, J, par) - hessian += ∇²Σ + hessian .+= ∇²Σ end return hessian end @@ -238,7 +238,7 @@ function objective_hessian!( if !semwls.approximate_hessian J = -2 * (σ₋' * semwls.V)' ∇²Σ_function!(∇²Σ, J, par) - hessian += ∇²Σ + hessian .+= ∇²Σ end return objective, hessian @@ -273,7 +273,7 @@ function gradient_hessian!( if !semwls.approximate_hessian J = -2 * (σ₋' * semwls.V)' ∇²Σ_function!(∇²Σ, J, par) - hessian += ∇²Σ + hessian .+= ∇²Σ end return gradient, hessian @@ -308,7 +308,7 @@ function objective_gradient_hessian!( if !semwls.approximate_hessian J = -2 * (σ₋' * semwls.V)' ∇²Σ_function!(∇²Σ, J, par) - hessian += ∇²Σ + hessian .+= ∇²Σ end return objective, gradient, hessian end From f567d865166ac9cb8d51b19d938d910e45c66ffc Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sun, 10 Mar 2024 13:15:53 -0700 Subject: [PATCH 07/29] ML: optimize kron --- src/loss/ML/ML.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/loss/ML/ML.jl b/src/loss/ML/ML.jl index 4c9a1f134..b6104d68a 100644 --- a/src/loss/ML/ML.jl +++ b/src/loss/ML/ML.jl @@ -172,7 +172,7 @@ function hessian!( J = vec(Σ⁻¹ - Σ⁻¹ΣₒΣ⁻¹)' ∇²Σ_function!(∇²Σ, J, par) # outer - H_outer = 2 * kron(Σ⁻¹ΣₒΣ⁻¹, Σ⁻¹) - kron(Σ⁻¹, Σ⁻¹) + H_outer = kron(2Σ⁻¹ΣₒΣ⁻¹ - Σ⁻¹, Σ⁻¹) hessian = ∇Σ' * H_outer * ∇Σ hessian .+= ∇²Σ end @@ -266,7 +266,7 @@ function objective_hessian!( J = vec(Σ⁻¹ - Σ⁻¹ΣₒΣ⁻¹)' ∇²Σ_function!(∇²Σ, J, par) # outer - H_outer = 2 * kron(Σ⁻¹ΣₒΣ⁻¹, Σ⁻¹) - kron(Σ⁻¹, Σ⁻¹) + H_outer = kron(2Σ⁻¹ΣₒΣ⁻¹ - Σ⁻¹, Σ⁻¹) hessian = ∇Σ' * H_outer * ∇Σ hessian .+= ∇²Σ end @@ -320,7 +320,7 @@ function gradient_hessian!( # inner ∇²Σ_function!(∇²Σ, J, par) # outer - H_outer = 2 * kron(Σ⁻¹ΣₒΣ⁻¹, Σ⁻¹) - kron(Σ⁻¹, Σ⁻¹) + H_outer = kron(2Σ⁻¹ΣₒΣ⁻¹ - Σ⁻¹, Σ⁻¹) hessian = ∇Σ' * H_outer * ∇Σ hessian .+= ∇²Σ end @@ -379,7 +379,7 @@ function objective_gradient_hessian!( # inner ∇²Σ_function!(∇²Σ, J, par) # outer - H_outer = 2 * kron(Σ⁻¹ΣₒΣ⁻¹, Σ⁻¹) - kron(Σ⁻¹, Σ⁻¹) + H_outer = kron(2Σ⁻¹ΣₒΣ⁻¹ - Σ⁻¹, Σ⁻¹) hessian = ∇Σ' * H_outer * ∇Σ hessian .+= ∇²Σ end From b37bc767e4f6f26c27674a3deb301afcb6cdf377 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sun, 10 Mar 2024 13:26:41 -0700 Subject: [PATCH 08/29] ML: optimize C reuse sigma^{-1}*sigma_0 product --- src/loss/ML/ML.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/loss/ML/ML.jl b/src/loss/ML/ML.jl index b6104d68a..7811cda7f 100644 --- a/src/loss/ML/ML.jl +++ b/src/loss/ML/ML.jl @@ -134,10 +134,10 @@ function gradient!( if T μ₋ = μₒ - μ μ₋ᵀΣ⁻¹ = μ₋' * Σ⁻¹ - gradient = vec(Σ⁻¹ * (I - Σₒ * Σ⁻¹ - μ₋ * μ₋ᵀΣ⁻¹))' * ∇Σ - 2 * μ₋ᵀΣ⁻¹ * ∇μ + gradient = vec(Σ⁻¹ - Σ⁻¹Σₒ * Σ⁻¹ - μ₋ᵀΣ⁻¹'μ₋ᵀΣ⁻¹)' * ∇Σ - 2 * μ₋ᵀΣ⁻¹ * ∇μ return gradient' else - gradient = (vec(Σ⁻¹) - vec(Σ⁻¹Σₒ * Σ⁻¹))' * ∇Σ + gradient = vec(Σ⁻¹ - Σ⁻¹Σₒ * Σ⁻¹)' * ∇Σ return gradient' end end @@ -494,9 +494,9 @@ function gradient!( Σ_chol = cholesky!(Symmetric(Σ⁻¹); check = false) isposdef(Σ_chol) || return ones(eltype(par), size(par)) Σ⁻¹ = LinearAlgebra.inv!(Σ_chol) - #mul!(Σ⁻¹Σₒ, Σ⁻¹, Σₒ) + mul!(Σ⁻¹Σₒ, Σ⁻¹, Σₒ) - C = F⨉I_A⁻¹' * (I - Σₒ * Σ⁻¹)' * Σ⁻¹ * F⨉I_A⁻¹ + C = F⨉I_A⁻¹' * (I - Σ⁻¹Σₒ) * Σ⁻¹ * F⨉I_A⁻¹ gradient = 2vec(C * S * I_A⁻¹')'∇A + vec(C)'∇S if T @@ -544,7 +544,7 @@ function objective_gradient!( mul!(Σ⁻¹Σₒ, Σ⁻¹, Σₒ) objective = ld + tr(Σ⁻¹Σₒ) - C = F⨉I_A⁻¹' * (I - Σₒ * Σ⁻¹)' * Σ⁻¹ * F⨉I_A⁻¹ + C = F⨉I_A⁻¹' * (I - Σ⁻¹Σₒ) * Σ⁻¹ * F⨉I_A⁻¹ gradient = 2vec(C * S * I_A⁻¹')'∇A + vec(C)'∇S if T From 84b0076995701e88e4429211bc8007143c92a027 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sat, 9 Mar 2024 18:39:48 -0800 Subject: [PATCH 09/29] start_fabin3: optimize indexing * use sorted searches * use dict to map indicator to param index --- .../start_val/start_fabin3.jl | 116 +++++++++--------- 1 file changed, 61 insertions(+), 55 deletions(-) diff --git a/src/additional_functions/start_val/start_fabin3.jl b/src/additional_functions/start_val/start_fabin3.jl index bce6471c7..bb29ce63b 100644 --- a/src/additional_functions/start_val/start_fabin3.jl +++ b/src/additional_functions/start_val/start_fabin3.jl @@ -46,9 +46,7 @@ function start_fabin3(ram_matrices::RAMMatrices, Σ, μ) # check in which matrix each parameter appears - indices = Vector{CartesianIndex}(undef, n_par) - - start_val = zeros(n_par) + indices = Vector{CartesianIndex{2}}(undef, n_par) #= in_S = length.(S_ind) .!= 0 in_A = length.(A_ind) .!= 0 @@ -68,16 +66,12 @@ function start_fabin3(ram_matrices::RAMMatrices, Σ, μ) # set undirected parameters in S for (i, S_ind) in enumerate(S_ind) - if length(S_ind) != 0 - c_ind = C_indices[S_ind][1] - if c_ind[1] == c_ind[2] - if c_ind[1] ∈ F_ind - index_position = findall(c_ind[1] .== F_ind) - start_val[i] = Σ[index_position[1], index_position[1]] / 2 - else - start_val[i] = 0.05 - end - end # covariances stay 0 + for c_ind in C_indices[S_ind] + (c_ind[1] == c_ind[2]) || continue # covariances stay 0 + pos = searchsortedfirst(F_ind, c_ind[1]) + start_val[i] = + (pos <= length(F_ind)) && (F_ind[pos] == c_ind[1]) ? Σ[pos, pos] / 2 : 0.05 + break # i-th parameter initialized end end @@ -86,77 +80,86 @@ function start_fabin3(ram_matrices::RAMMatrices, Σ, μ) A_ind_c = [linear2cartesian(ind, (n_nod, n_nod)) for ind in A_ind] # ind_Λ = findall([is_in_Λ(ind_vec, F_ind) for ind_vec in A_ind_c]) - for i in findall(.!(1:n_nod .∈ [F_ind])) + for i in setdiff(1:n_nod, F_ind) reference = Int64[] indicators = Int64[] - loadings = Symbol[] - - for (j, ind_c_vec) in enumerate(A_ind_c) - for ind_c in ind_c_vec - if (ind_c[2] == i) & (ind_c[1] ∈ F_ind) - push!(indicators, ind_c[1]) - push!(loadings, parameters[j]) + indicator2parampos = Dict{Int, Int}() + + for (j, Aj_ind_c) in enumerate(A_ind_c) + for ind_c in Aj_ind_c + (ind_c[2] == i) || continue + ind_pos = searchsortedfirst(F_ind, ind_c[1]) + if (ind_pos <= length(F_ind)) && (F_ind[ind_pos] == ind_c[1]) + push!(indicators, ind_pos) + indicator2parampos[ind_pos] = j end end end - for ram_constant in constants - if (ram_constant.matrix == :A) && - (ram_constant.index[2] == i) && - (ram_constant.index[1] ∈ F_ind) - push!(loadings, Symbol("")) - if isone(ram_constant.value) - push!(reference, ram_constant.index[1]) - else - push!(indicators, ram_constant.index[1]) + for ram_const in constants + if (ram_const.matrix == :A) && (ram_const.index[2] == i) + ind_pos = searchsortedfirst(F_ind, ram_const.index[1]) + if (ind_pos <= length(F_ind)) && (F_ind[ind_pos] == ram_const.index[1]) + if isone(ram_const.value) + push!(reference, ind_pos) + else + push!(indicators, ind_pos) + # no parameter associated + end end end end - reference = [findfirst(x -> x == ref, F_ind) for ref in reference] - indicators = [findfirst(x -> x == ind, F_ind) for ind in indicators] - # is there at least one reference indicator? - if size(reference, 1) > 0 - if size(reference, 1) > 1 - @warn "You have more than 1 scaling indicator" - reference = reference[1] + if length(reference) > 0 + if length(reference) > 1 + if isempty(indicator2parampos) # don't warn if entire column is fixed + @warn "You have more than 1 scaling indicator for $(ram_matrices.colnames[i])" + end + ref = reference[1] else - reference = reference[1] + ref = reference[1] end for (j, indicator) in enumerate(indicators) - if indicator != reference - instruments = indicators[.!(indicators .∈ [[reference; indicator]])] + if indicator != ref + instruments = filter(i -> (i != ref) && (i != indicator), indicators) s32 = Σ[instruments, indicator] - s13 = Σ[reference, instruments] + s13 = Σ[ref, instruments] S33 = Σ[instruments, instruments] if size(instruments, 1) == 1 temp = S33[1] / s13[1] λ = s32[1] * temp / (s13[1] * temp) - start_val[isequal.(parameters, loadings[j])] .= λ else temp = S33 \ s13 λ = s32' * temp / (s13' * temp) - start_val[isequal.(parameters, loadings[j])] .= λ end + start_val[indicator2parampos[indicator]] = λ end end # no reference indicator: - else - reference = indicators[1] - λ = zeros(size(indicators, 1)) + elseif length(indicators) > 0 + ref = indicators[1] + λ = Vector{Float64}(undef, length(indicators)) λ[1] = 1.0 for (j, indicator) in enumerate(indicators) - if indicator != reference - instruments = indicators[.!(indicators .∈ [[reference; indicator]])] + if indicator != ref + instruments = filter(i -> (i != ref) && (i != indicator), indicators) s32 = Σ[instruments, indicator] - s13 = Σ[reference, instruments] + s13 = Σ[ref, instruments] S33 = Σ[instruments, instruments] + if length(instruments) == 1 + temp = S33[1] / s13[1] + λ[j] = s32[1] * temp / (s13[1] * temp) + else + temp = S33 \ s13 + λ[j] = s32' * temp / (s13' * temp) + end + if size(instruments, 1) == 1 temp = S33[1] / s13[1] λ[j] = s32[1] * temp / (s13[1] * temp) @@ -179,8 +182,12 @@ function start_fabin3(ram_matrices::RAMMatrices, Σ, μ) λ = λ .* sign(Ψ) .* sqrt(abs(Ψ)) for (j, indicator) in enumerate(indicators) - start_val[isequal.(parameters, loadings[j])] .= λ[j] + if (parampos = get(indicator2parampos, indicator, 0)) != 0 + start_val[parampos] = λ[j] + end end + else + @warn "No scaling indicators for $(ram_matrices.colnames[i])" end end @@ -189,9 +196,9 @@ function start_fabin3(ram_matrices::RAMMatrices, Σ, μ) for (i, M_ind) in enumerate(M_ind) if length(M_ind) != 0 ind = M_ind[1] - if ind[1] ∈ F_ind - index_position = findfirst(ind[1] .== F_ind) - start_val[i] = μ[index_position] + pos = searchsortedfirst(F_ind, ind[1]) + if (pos <= length(F_ind)) && (F_ind[pos] == ind[1]) + start_val[i] = μ[pos] end # latent means stay 0 end end @@ -201,6 +208,5 @@ function start_fabin3(ram_matrices::RAMMatrices, Σ, μ) end function is_in_Λ(ind_vec, F_ind) - res = [!(ind[2] ∈ F_ind) & (ind[1] ∈ F_ind) for ind in ind_vec] - return any(res) + return any(ind -> !(ind[2] ∈ F_ind) && (ind[1] ∈ F_ind), ind_vec) end From 7dbd687b5663a2a1161e9f841dc75b17f80b9b68 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Mon, 22 Apr 2024 03:35:53 -0700 Subject: [PATCH 10/29] fixup start_fabin3.jl --- src/additional_functions/start_val/start_fabin3.jl | 10 +++------- 1 file changed, 3 insertions(+), 7 deletions(-) diff --git a/src/additional_functions/start_val/start_fabin3.jl b/src/additional_functions/start_val/start_fabin3.jl index bb29ce63b..23b9ff44d 100644 --- a/src/additional_functions/start_val/start_fabin3.jl +++ b/src/additional_functions/start_val/start_fabin3.jl @@ -112,14 +112,10 @@ function start_fabin3(ram_matrices::RAMMatrices, Σ, μ) # is there at least one reference indicator? if length(reference) > 0 - if length(reference) > 1 - if isempty(indicator2parampos) # don't warn if entire column is fixed - @warn "You have more than 1 scaling indicator for $(ram_matrices.colnames[i])" - end - ref = reference[1] - else - ref = reference[1] + if (length(reference) > 1) && isempty(indicator2parampos) # don't warn if entire column is fixed + @warn "You have more than 1 scaling indicator for $(ram_matrices.colnames[i])" end + ref = reference[1] for (j, indicator) in enumerate(indicators) if indicator != ref From 177a3424830647770080d71bb358c2a0cb05062f Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sat, 9 Mar 2024 18:45:22 -0800 Subject: [PATCH 11/29] start_fabin3(): optimize math * extract common code to calculate_lambda() * use 3-arg dot() * use in-place ops --- .../start_val/start_fabin3.jl | 67 ++++++++----------- 1 file changed, 27 insertions(+), 40 deletions(-) diff --git a/src/additional_functions/start_val/start_fabin3.jl b/src/additional_functions/start_val/start_fabin3.jl index 23b9ff44d..558b4c40f 100644 --- a/src/additional_functions/start_val/start_fabin3.jl +++ b/src/additional_functions/start_val/start_fabin3.jl @@ -80,6 +80,25 @@ function start_fabin3(ram_matrices::RAMMatrices, Σ, μ) A_ind_c = [linear2cartesian(ind, (n_nod, n_nod)) for ind in A_ind] # ind_Λ = findall([is_in_Λ(ind_vec, F_ind) for ind_vec in A_ind_c]) + function calculate_lambda( + ref::Integer, + indicator::Integer, + indicators::AbstractVector{<:Integer}, + ) + instruments = filter(i -> (i != ref) && (i != indicator), indicators) + if length(instruments) == 1 + s13 = Σ[ref, instruments[1]] + s32 = Σ[instruments[1], indicator] + return s32 / s13 + else + s13 = Σ[ref, instruments] + s32 = Σ[instruments, indicator] + S33 = Σ[instruments, instruments] + temp = S33 \ s13 + return dot(s32, temp) / dot(s13, temp) + end + end + for i in setdiff(1:n_nod, F_ind) reference = Int64[] indicators = Int64[] @@ -118,21 +137,9 @@ function start_fabin3(ram_matrices::RAMMatrices, Σ, μ) ref = reference[1] for (j, indicator) in enumerate(indicators) - if indicator != ref - instruments = filter(i -> (i != ref) && (i != indicator), indicators) - - s32 = Σ[instruments, indicator] - s13 = Σ[ref, instruments] - S33 = Σ[instruments, instruments] - - if size(instruments, 1) == 1 - temp = S33[1] / s13[1] - λ = s32[1] * temp / (s13[1] * temp) - else - temp = S33 \ s13 - λ = s32' * temp / (s13' * temp) - end - start_val[indicator2parampos[indicator]] = λ + if (indicator != ref) && + (parampos = get(indicator2parampos, indicator, 0)) != 0 + start_val[parampos] = calculate_lambda(ref, indicator, indicators) end end # no reference indicator: @@ -142,40 +149,20 @@ function start_fabin3(ram_matrices::RAMMatrices, Σ, μ) λ[1] = 1.0 for (j, indicator) in enumerate(indicators) if indicator != ref - instruments = filter(i -> (i != ref) && (i != indicator), indicators) - - s32 = Σ[instruments, indicator] - s13 = Σ[ref, instruments] - S33 = Σ[instruments, instruments] - - if length(instruments) == 1 - temp = S33[1] / s13[1] - λ[j] = s32[1] * temp / (s13[1] * temp) - else - temp = S33 \ s13 - λ[j] = s32' * temp / (s13' * temp) - end - - if size(instruments, 1) == 1 - temp = S33[1] / s13[1] - λ[j] = s32[1] * temp / (s13[1] * temp) - else - temp = S33 \ s13 - λ[j] = s32' * temp / (s13' * temp) - end + λ[j] = calculate_lambda(ref, indicator, indicators) end end Σ_λ = Σ[indicators, indicators] - D = λ * λ' ./ sum(λ .^ 2) + l₂ = sum(abs2, λ) + D = λ * λ' ./ l₂ θ = (I - D .^ 2) \ (diag(Σ_λ - D * Σ_λ * D)) # 3. psi Σ₁ = Σ_λ - Diagonal(θ) - l₂ = sum(λ .^ 2) - Ψ = sum(sum(λ .* Σ₁, dims = 1) .* λ') ./ l₂^2 + Ψ = dot(λ, Σ₁, λ) / l₂^2 - λ = λ .* sign(Ψ) .* sqrt(abs(Ψ)) + λ .*= sign(Ψ) * sqrt(abs(Ψ)) for (j, indicator) in enumerate(indicators) if (parampos = get(indicator2parampos, indicator, 0)) != 0 From 93a54e6aeedf6e7eb38165c1f6c810e7a5368068 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Mon, 11 Mar 2024 22:36:22 -0700 Subject: [PATCH 12/29] start_fabin3(): directly access imply.ram_matrices --- src/additional_functions/start_val/start_fabin3.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/additional_functions/start_val/start_fabin3.jl b/src/additional_functions/start_val/start_fabin3.jl index 558b4c40f..ee7dcb8cf 100644 --- a/src/additional_functions/start_val/start_fabin3.jl +++ b/src/additional_functions/start_val/start_fabin3.jl @@ -18,7 +18,7 @@ function start_fabin3(model::AbstractSemSingle; kwargs...) end function start_fabin3(observed, imply, optimizer, args...; kwargs...) - return start_fabin3(ram_matrices(imply), obs_cov(observed), obs_mean(observed)) + return start_fabin3(imply.ram_matrices, obs_cov(observed), obs_mean(observed)) end # SemObservedMissing @@ -27,7 +27,7 @@ function start_fabin3(observed::SemObservedMissing, imply, optimizer, args...; k em_mvn(observed; kwargs...) end - return start_fabin3(ram_matrices(imply), observed.em_model.Σ, observed.em_model.μ) + return start_fabin3(imply.ram_matrices, observed.em_model.Σ, observed.em_model.μ) end function start_fabin3(ram_matrices::RAMMatrices, Σ, μ) From be249a11ecd8dc15801d6a6fd12caf18588aa860 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sun, 10 Mar 2024 13:32:28 -0700 Subject: [PATCH 13/29] remove_all_missing(): optimize avoid unnecessary allocations --- src/additional_functions/helper.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/additional_functions/helper.jl b/src/additional_functions/helper.jl index 05253cd8c..46598cfbe 100644 --- a/src/additional_functions/helper.jl +++ b/src/additional_functions/helper.jl @@ -52,10 +52,10 @@ function F_one_person(imp_mean, meandiff, inverse, data, logdet) return F end -function remove_all_missing(data) +function remove_all_missing(data::AbstractMatrix) keep = Vector{Int64}() - for i in 1:size(data, 1) - if any(.!ismissing.(data[i, :])) + for (i, coldata) in zip(axes(data, 1), eachrow(data)) + if any(!ismissing, coldata) push!(keep, i) end end From 461393ad52896585892e65eeae58437442d11f02 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sun, 21 Apr 2024 10:28:56 -0700 Subject: [PATCH 14/29] skipmissing_mean(): optimize * use eachcol() within comprehension --- src/additional_functions/helper.jl | 9 ++------- 1 file changed, 2 insertions(+), 7 deletions(-) diff --git a/src/additional_functions/helper.jl b/src/additional_functions/helper.jl index 46598cfbe..3d55fed42 100644 --- a/src/additional_functions/helper.jl +++ b/src/additional_functions/helper.jl @@ -37,13 +37,8 @@ function get_observed(rowind, data, semobserved; args = (), kwargs = NamedTuple( return observed_vec end -function skipmissing_mean(mat) - means = Vector{Float64}(undef, size(mat, 2)) - for i in 1:size(mat, 2) - @views means[i] = mean(skipmissing(mat[:, i])) - end - return means -end +skipmissing_mean(mat::AbstractMatrix) = + [mean(skipmissing(coldata)) for coldata in eachcol(mat)] function F_one_person(imp_mean, meandiff, inverse, data, logdet) F = logdet From 5c96862c097e0f464de3a1bd578801ca676ed5ad Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sun, 10 Mar 2024 13:46:53 -0700 Subject: [PATCH 15/29] use ternary op as intended --- src/imply/RAM/generic.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/imply/RAM/generic.jl b/src/imply/RAM/generic.jl index 2089cf3da..a0c29185d 100644 --- a/src/imply/RAM/generic.jl +++ b/src/imply/RAM/generic.jl @@ -140,8 +140,7 @@ function RAM(; # get indices A_indices = copy(ram_matrices.A_ind) S_indices = copy(ram_matrices.S_ind) - !isnothing(ram_matrices.M_ind) ? M_indices = copy(ram_matrices.M_ind) : - M_indices = nothing + M_indices = !isnothing(ram_matrices.M_ind) ? copy(ram_matrices.M_ind) : nothing #preallocate arrays A_pre = zeros(n_nod, n_nod) From f3c7dd3403badda3bb9492efcbfbce670b98fd42 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sun, 10 Mar 2024 13:48:22 -0700 Subject: [PATCH 16/29] symbolic: constrain to tril before simplifying --- src/imply/RAM/symbolic.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/imply/RAM/symbolic.jl b/src/imply/RAM/symbolic.jl index 6ce7b6ae5..32b06fbd8 100644 --- a/src/imply/RAM/symbolic.jl +++ b/src/imply/RAM/symbolic.jl @@ -277,12 +277,12 @@ function get_Σ_symbolic_RAM(S, A, F; vech = false) Σ_symbolic = F * invia * S * permutedims(invia) * permutedims(F) Σ_symbolic = Array(Σ_symbolic) # Σ_symbolic = Symbolics.simplify.(Σ_symbolic) - Threads.@threads for i in eachindex(Σ_symbolic) - Σ_symbolic[i] = Symbolics.simplify(Σ_symbolic[i]) - end if vech Σ_symbolic = Σ_symbolic[tril(trues(size(F, 1), size(F, 1)))] end + Threads.@threads for i in eachindex(Σ_symbolic) + Σ_symbolic[i] = Symbolics.simplify(Σ_symbolic[i]) + end return Σ_symbolic end From 17e50f2ce517d66092cda2c55742ec89a798ee78 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sun, 21 Apr 2024 10:27:34 -0700 Subject: [PATCH 17/29] neumann_series(): avoid endless loop if the series don't converge --- src/additional_functions/helper.jl | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/additional_functions/helper.jl b/src/additional_functions/helper.jl index 3d55fed42..abc37207c 100644 --- a/src/additional_functions/helper.jl +++ b/src/additional_functions/helper.jl @@ -1,10 +1,14 @@ -function neumann_series(mat::SparseMatrixCSC) +# Neumann series representation of (I - mat)⁻¹ +function neumann_series(mat::SparseMatrixCSC; maxiter::Integer = size(mat, 1)) inverse = I + mat next_term = mat^2 + n = 1 while nnz(next_term) != 0 + (n <= maxiter) || error("Neumann series did not converge in $maxiter steps") inverse += next_term next_term *= mat + n += 1 end return inverse From 447e56e73392943885f82b325b9b0fbc3c4dc6c6 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Mon, 1 Apr 2024 10:01:19 -0700 Subject: [PATCH 18/29] RAMSymbolic: calc (I-A)^{-1} once as symbolic Neumann series is a very comp.intensive operation --- src/imply/RAM/symbolic.jl | 40 +++++++++++++++++++-------------------- 1 file changed, 20 insertions(+), 20 deletions(-) diff --git a/src/imply/RAM/symbolic.jl b/src/imply/RAM/symbolic.jl index 32b06fbd8..4a92e002e 100644 --- a/src/imply/RAM/symbolic.jl +++ b/src/imply/RAM/symbolic.jl @@ -120,8 +120,10 @@ function RAMSymbolic(; any(loss_types .<: SemWLS) ? vech = true : nothing end + I_A⁻¹ = neumann_series(A) + # Σ - Σ_symbolic = get_Σ_symbolic_RAM(S, A, F; vech = vech) + Σ_symbolic = eval_Σ_symbolic(S, I_A⁻¹, F; vech = vech) #print(Symbolics.build_function(Σ_symbolic)[2]) Σ_function = Symbolics.build_function(Σ_symbolic, par, expression = Val{false})[2] Σ = zeros(size(Σ_symbolic)) @@ -163,7 +165,7 @@ function RAMSymbolic(; # μ if meanstructure has_meanstructure = Val(true) - μ_symbolic = get_μ_symbolic_RAM(M, A, F) + μ_symbolic = eval_μ_symbolic(M, I_A⁻¹, F) μ_function = Symbolics.build_function(μ_symbolic, par, expression = Val{false})[2] μ = zeros(size(μ_symbolic)) if gradient @@ -272,26 +274,24 @@ ram_matrices(imply::RAMSymbolic) = imply.ram_matrices ### additional functions ############################################################################################ -function get_Σ_symbolic_RAM(S, A, F; vech = false) - invia = neumann_series(A) - Σ_symbolic = F * invia * S * permutedims(invia) * permutedims(F) - Σ_symbolic = Array(Σ_symbolic) - # Σ_symbolic = Symbolics.simplify.(Σ_symbolic) - if vech - Σ_symbolic = Σ_symbolic[tril(trues(size(F, 1), size(F, 1)))] - end - Threads.@threads for i in eachindex(Σ_symbolic) - Σ_symbolic[i] = Symbolics.simplify(Σ_symbolic[i]) +# expected covariations of observed vars +function eval_Σ_symbolic(S, I_A⁻¹, F; vech = false) + Σ = F * I_A⁻¹ * S * permutedims(I_A⁻¹) * permutedims(F) + Σ = Array(Σ) + vech && (Σ = Σ[tril(trues(size(F, 1), size(F, 1)))]) + # Σ = Symbolics.simplify.(Σ) + Threads.@threads for i in eachindex(Σ) + Σ[i] = Symbolics.simplify(Σ[i]) end - return Σ_symbolic + return Σ end -function get_μ_symbolic_RAM(M, A, F) - invia = neumann_series(A) - μ_symbolic = F * invia * M - μ_symbolic = Array(μ_symbolic) - Threads.@threads for i in eachindex(μ_symbolic) - μ_symbolic[i] = Symbolics.simplify(μ_symbolic[i]) +# expected means of observed vars +function eval_μ_symbolic(M, I_A⁻¹, F) + μ = F * I_A⁻¹ * M + μ = Array(μ) + Threads.@threads for i in eachindex(μ) + μ[i] = Symbolics.simplify(μ[i]) end - return μ_symbolic + return μ end From ac5e7bb7e9727a2298998a5f2d6581c712efee34 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Thu, 18 Apr 2024 21:46:58 -0700 Subject: [PATCH 19/29] n_obs/man(data): restrict to integer don't allow missing n_obs --- src/imply/RAM/symbolic.jl | 2 +- src/observed/covariance.jl | 10 +++--- src/observed/data.jl | 12 +++---- .../political_democracy/constructor.jl | 4 +-- test/unit_tests/data_input_formats.jl | 32 +++++++++++++------ 5 files changed, 36 insertions(+), 24 deletions(-) diff --git a/src/imply/RAM/symbolic.jl b/src/imply/RAM/symbolic.jl index 4a92e002e..5c5e52112 100644 --- a/src/imply/RAM/symbolic.jl +++ b/src/imply/RAM/symbolic.jl @@ -244,7 +244,7 @@ identifier(imply::RAMSymbolic) = imply.identifier n_par(imply::RAMSymbolic) = imply.n_par function update_observed(imply::RAMSymbolic, observed::SemObserved; kwargs...) - if Int(n_man(observed)) == size(imply.Σ, 1) + if n_man(observed) == size(imply.Σ, 1) return imply else return RAMSymbolic(; observed = observed, kwargs...) diff --git a/src/observed/covariance.jl b/src/observed/covariance.jl index 5785ba404..1b5de9fc2 100644 --- a/src/observed/covariance.jl +++ b/src/observed/covariance.jl @@ -39,11 +39,11 @@ use this if you are sure your covariance matrix is in the right format. ## Additional keyword arguments: - `spec_colnames::Vector{Symbol} = nothing`: overwrites column names of the specification object """ -struct SemObservedCovariance{B, C, D, O} <: SemObserved +struct SemObservedCovariance{B, C} <: SemObserved obs_cov::B obs_mean::C - n_man::D - n_obs::O + n_man::Int + n_obs::Int end function SemObservedCovariance(; @@ -80,9 +80,7 @@ function SemObservedCovariance(; (obs_mean = reorder_obs_mean(obs_mean, spec_colnames, obs_colnames)) end - n_man = Float64(size(obs_cov, 1)) - - return SemObservedCovariance(obs_cov, obs_mean, n_man, n_obs) + return SemObservedCovariance(obs_cov, obs_mean, size(obs_cov, 1), n_obs) end ############################################################################################ diff --git a/src/observed/data.jl b/src/observed/data.jl index 77bbc8e8b..0d9ad3a04 100644 --- a/src/observed/data.jl +++ b/src/observed/data.jl @@ -39,18 +39,18 @@ use this if you are sure your observed data is in the right format. - `compute_covariance::Bool ) = true`: should the covariance of `data` be computed and stored? - `rowwise::Bool = false`: should the data be stored also as vectors per observation """ -struct SemObservedData{A, B, C, D, O, R} <: SemObserved +struct SemObservedData{A, B, C, R} <: SemObserved data::A obs_cov::B obs_mean::C - n_man::D - n_obs::O + n_man::Int + n_obs::Int data_rowwise::R end # error checks function check_arguments_SemObservedData(kwargs...) - # data is a data frame, + # data is a data frame, end @@ -107,8 +107,8 @@ function SemObservedData(; data, compute_covariance ? Statistics.cov(data) : nothing, meanstructure ? vec(Statistics.mean(data, dims = 1)) : nothing, - Float64.(size(data, 2)), - Float64.(size(data, 1)), + size(data, 2), + size(data, 1), rowwise ? [data[i, :] for i in axes(data, 1)] : nothing, ) end diff --git a/test/examples/political_democracy/constructor.jl b/test/examples/political_democracy/constructor.jl index e13bc3816..97f5e42d8 100644 --- a/test/examples/political_democracy/constructor.jl +++ b/test/examples/political_democracy/constructor.jl @@ -12,7 +12,7 @@ model_ml_cov = Sem( obs_cov = cov(Matrix(dat)), obs_colnames = Symbol.(names(dat)), optimizer = semoptimizer, - n_obs = 75.0, + n_obs = 75, ) model_ls_sym = Sem( @@ -255,7 +255,7 @@ model_ml_cov = Sem( obs_colnames = Symbol.(names(dat)), meanstructure = true, optimizer = semoptimizer, - n_obs = 75.0, + n_obs = 75, ) model_ml_sym = Sem( diff --git a/test/unit_tests/data_input_formats.jl b/test/unit_tests/data_input_formats.jl index bef7ba8df..485cf82d2 100644 --- a/test/unit_tests/data_input_formats.jl +++ b/test/unit_tests/data_input_formats.jl @@ -200,13 +200,18 @@ end # errors @test_throws ArgumentError("observed means were passed, but `meanstructure = false`") begin - SemObservedCovariance(specification = nothing, obs_cov = dat_cov, obs_mean = dat_mean) + SemObservedCovariance( + specification = nothing, + obs_cov = dat_cov, + obs_mean = dat_mean, + n_obs = 75, + ) end @test_throws UndefKeywordError(:specification) SemObservedCovariance(obs_cov = dat_cov) @test_throws ArgumentError("no `obs_colnames` were specified") begin - SemObservedCovariance(specification = spec, obs_cov = dat_cov) + SemObservedCovariance(specification = spec, obs_cov = dat_cov, n_obs = 75) end @test_throws ArgumentError("please specify `obs_colnames` as a vector of Symbols") begin @@ -214,6 +219,7 @@ end specification = spec, obs_cov = dat_cov, obs_colnames = names(dat), + n_obs = 75, ) end @@ -222,17 +228,18 @@ observed = SemObservedCovariance( specification = spec, obs_cov = dat_cov, obs_colnames = obs_colnames = Symbol.(names(dat)), + n_obs = 75, ) observed_nospec = - SemObservedCovariance(specification = nothing, obs_cov = dat_cov, n_obs = 75.0) + SemObservedCovariance(specification = nothing, obs_cov = dat_cov, n_obs = 75) all_equal_cov = (obs_cov(observed) == obs_cov(observed_nospec)) @testset "unit tests | SemObservedCovariance | input formats" begin @test all_equal_cov - @test isnothing(n_obs(observed)) - @test n_obs(observed_nospec) == 75.0 + @test n_obs(observed) == 75 + @test n_obs(observed_nospec) == 75 end # shuffle variables @@ -248,6 +255,7 @@ observed_shuffle = SemObservedCovariance( specification = spec, obs_cov = shuffle_dat_cov, obs_colnames = shuffle_names, + n_obs = 75, ) all_equal_cov_suffled = (obs_cov(observed) ≈ obs_cov(observed_shuffle)) @@ -260,7 +268,12 @@ end # errors @test_throws ArgumentError("`meanstructure = true`, but no observed means were passed") begin - SemObservedCovariance(specification = spec, obs_cov = dat_cov, meanstructure = true) + SemObservedCovariance( + specification = spec, + obs_cov = dat_cov, + meanstructure = true, + n_obs = 75, + ) end @test_throws UndefKeywordError SemObservedCovariance( @@ -279,6 +292,7 @@ end obs_cov = dat_cov, obs_colnames = Symbol.(names(dat)), meanstructure = true, + n_obs = 75, ) end @@ -288,7 +302,7 @@ observed = SemObservedCovariance( obs_cov = dat_cov, obs_mean = dat_mean, obs_colnames = Symbol.(names(dat)), - n_obs = 75.0, + n_obs = 75, meanstructure = true, ) @@ -297,7 +311,7 @@ observed_nospec = SemObservedCovariance( obs_cov = dat_cov, obs_mean = dat_mean, meanstructure = true, - n_obs = 75.0, + n_obs = 75, ) all_equal_mean = (obs_mean(observed) == obs_mean(observed_nospec)) @@ -323,7 +337,7 @@ observed_shuffle = SemObservedCovariance( obs_cov = shuffle_dat_cov, obs_mean = shuffle_dat_mean, obs_colnames = shuffle_names, - n_obs = 75.0, + n_obs = 75, meanstructure = true, ) From 46ca0975bc4bf7df7765b04f6d30c007a923ced5 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sun, 10 Mar 2024 13:30:23 -0700 Subject: [PATCH 20/29] refactor get_partition() * convert to param_range() which gets the range for a single matrix * use findfirst/findlast() instead of manual loop --- src/additional_functions/parameters.jl | 106 ++++--------------------- 1 file changed, 16 insertions(+), 90 deletions(-) diff --git a/src/additional_functions/parameters.jl b/src/additional_functions/parameters.jl index 1b7d4643f..144590bd1 100644 --- a/src/additional_functions/parameters.jl +++ b/src/additional_functions/parameters.jl @@ -111,97 +111,23 @@ function fill_matrix!( return M end -function get_partition(A_indices, S_indices) - n_par = length(A_indices) - - first_A = "a" - first_S = "a" - last_A = "a" - last_S = "a" - - for i in 1:n_par - if length(A_indices[i]) != 0 - first_A = i - break - end - end - - for i in 1:n_par - if length(S_indices[i]) != 0 - first_S = i - break - end - end - - for i in n_par + 1 .- (1:n_par) - if length(A_indices[i]) != 0 - last_A = i - break - end - end - - for i in n_par + 1 .- (1:n_par) - if length(S_indices[i]) != 0 - last_S = i - break - end - end - - for i in first_A:last_A - if length(A_indices[i]) == 0 - throw( - ErrorException( - "Your parameter vector is not partitioned into directed and undirected effects", - ), - ) - return nothing - end - end - - for i in first_S:last_S - if length(S_indices[i]) == 0 - throw( - ErrorException( - "Your parameter vector is not partitioned into directed and undirected effects", - ), - ) - return nothing - end - end - - return first_A:last_A, first_S:last_S -end - -function get_partition(M_indices) - n_par = length(M_indices) - - first_M = "a" - last_M = "a" - - for i in 1:n_par - if length(M_indices[i]) != 0 - first_M = i - break - end - end - - for i in n_par + 1 .- (1:n_par) - if length(M_indices[i]) != 0 - last_M = i - break - end - end - - for i in first_M:last_M - if length(M_indices[i]) == 0 - throw( - ErrorException( - "Your parameter vector is not partitioned into directed, undirected and mean effects", - ), - ) - return nothing +# range of parameters that are referenced in the matrix +function param_range(mtx_indices::AbstractArrayParamsMap) + first_i = findfirst(!isempty, mtx_indices) + last_i = findlast(!isempty, mtx_indices) + + if !isnothing(first_i) && !isnothing(last_i) + for i in first_i:last_i + if isempty(mtx_indices[i]) + # TODO show which parameter is missing in which matrix + throw( + ErrorException( + "Your parameter vector is not partitioned into directed and undirected effects", + ), + ) + end end end - return first_M:last_M + return first_i:last_i end From 483efd0af62c6f2e22f39bdad610550b9c240b1a Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sat, 20 Apr 2024 16:16:56 -0700 Subject: [PATCH 21/29] test_gradient(): do tests inside --- test/examples/helper.jl | 12 +++++------- test/examples/multigroup/build_models.jl | 10 +++++----- test/examples/political_democracy/by_parts.jl | 8 ++++---- test/examples/political_democracy/constructor.jl | 8 ++++---- 4 files changed, 18 insertions(+), 20 deletions(-) diff --git a/test/examples/helper.jl b/test/examples/helper.jl index 6c212eadd..dca43695e 100644 --- a/test/examples/helper.jl +++ b/test/examples/helper.jl @@ -1,19 +1,17 @@ function test_gradient(model, parameters; rtol = 1e-10, atol = 0) true_grad = - FiniteDiff.finite_difference_gradient(x -> objective!(model, x)[1], parameters) + FiniteDiff.finite_difference_gradient(Base.Fix1(objective!, model), parameters) gradient = similar(parameters) - gradient .= 1.0 # F and G + fill!(gradient, NaN) gradient!(gradient, model, parameters) - correct1 = isapprox(gradient, true_grad; rtol = rtol, atol = atol) + @test gradient ≈ true_grad rtol = rtol atol = atol # only G - gradient .= 1.0 + fill!(gradient, NaN) objective_gradient!(gradient, model, parameters) - correct2 = isapprox(gradient, true_grad; rtol = rtol, atol = atol) - - return correct1 & correct2 + @test gradient ≈ true_grad rtol = rtol atol = atol end function test_hessian(model, parameters; rtol = 1e-4, atol = 0) diff --git a/test/examples/multigroup/build_models.jl b/test/examples/multigroup/build_models.jl index 2f0fc749f..c0038ccae 100644 --- a/test/examples/multigroup/build_models.jl +++ b/test/examples/multigroup/build_models.jl @@ -10,7 +10,7 @@ model_ml_multigroup = SemEnsemble(model_g1, model_g2; optimizer = semoptimizer) # gradients @testset "ml_gradients_multigroup" begin - @test test_gradient(model_ml_multigroup, start_test; atol = 1e-9) + test_gradient(model_ml_multigroup, start_test; atol = 1e-9) end # fit @@ -71,7 +71,7 @@ model_ml_multigroup = SemEnsemble(model_g1, model_g2; optimizer = semoptimizer) # gradients @testset "ml_gradients_multigroup | sorted" begin - @test test_gradient(model_ml_multigroup, start_test; atol = 1e-2) + test_gradient(model_ml_multigroup, start_test; atol = 1e-2) end grad = similar(start_test) @@ -159,7 +159,7 @@ model_g2 = SemFiniteDiff( model_ml_multigroup = SemEnsemble(model_g1, model_g2; optimizer = semoptimizer) @testset "gradients_user_defined_loss" begin - @test test_gradient(model_ml_multigroup, start_test; atol = 1e-9) + test_gradient(model_ml_multigroup, start_test; atol = 1e-9) end # fit @@ -187,7 +187,7 @@ model_ls_g2 = model_ls_multigroup = SemEnsemble(model_ls_g1, model_ls_g2; optimizer = semoptimizer) @testset "ls_gradients_multigroup" begin - @test test_gradient(model_ls_multigroup, start_test; atol = 1e-9) + test_gradient(model_ls_multigroup, start_test; atol = 1e-9) end @testset "ls_solution_multigroup" begin @@ -281,7 +281,7 @@ if !isnothing(specification_miss_g1) ] @testset "fiml_gradients_multigroup" begin - @test test_gradient(model_ml_multigroup, start_test; atol = 1e-7) + test_gradient(model_ml_multigroup, start_test; atol = 1e-7) end @testset "fiml_solution_multigroup" begin diff --git a/test/examples/political_democracy/by_parts.jl b/test/examples/political_democracy/by_parts.jl index 99d7a6715..4c3567572 100644 --- a/test/examples/political_democracy/by_parts.jl +++ b/test/examples/political_democracy/by_parts.jl @@ -54,7 +54,7 @@ model_names = ["ml", "ls_sym", "ridge", "constant", "ml_sym", "ml_weighted"] for (model, name) in zip(models, model_names) try @testset "$(name)_gradient" begin - @test test_gradient(model, start_test; rtol = 1e-9) + test_gradient(model, start_test; rtol = 1e-9) end catch end @@ -254,7 +254,7 @@ model_names = ["ml", "ls_sym", "ml_sym"] for (model, name) in zip(models, model_names) try @testset "$(name)_gradient_mean" begin - @test test_gradient(model, start_test_mean; rtol = 1e-9) + test_gradient(model, start_test_mean; rtol = 1e-9) end catch end @@ -343,11 +343,11 @@ model_ml_sym = Sem(observed, imply_ram_sym, loss_fiml, optimizer_obj) ############################################################################################ @testset "fiml_gradient" begin - @test test_gradient(model_ml, start_test_mean; atol = 1e-6) + test_gradient(model_ml, start_test_mean; atol = 1e-6) end @testset "fiml_gradient_symbolic" begin - @test test_gradient(model_ml_sym, start_test_mean; atol = 1e-6) + test_gradient(model_ml_sym, start_test_mean; atol = 1e-6) end ############################################################################################ diff --git a/test/examples/political_democracy/constructor.jl b/test/examples/political_democracy/constructor.jl index 97f5e42d8..757b26e60 100644 --- a/test/examples/political_democracy/constructor.jl +++ b/test/examples/political_democracy/constructor.jl @@ -68,7 +68,7 @@ model_names = ["ml", "ml_cov", "ls_sym", "ridge", "constant", "ml_sym", "ml_weig for (model, name) in zip(models, model_names) try @testset "$(name)_gradient" begin - @test test_gradient(model, start_test; rtol = 1e-9) + test_gradient(model, start_test; rtol = 1e-9) end catch end @@ -277,7 +277,7 @@ model_names = ["ml", "ml_cov", "ls_sym", "ml_sym"] for (model, name) in zip(models, model_names) try @testset "$(name)_gradient_mean" begin - @test test_gradient(model, start_test_mean; rtol = 1e-9) + test_gradient(model, start_test_mean; rtol = 1e-9) end catch end @@ -377,11 +377,11 @@ model_ml_sym = Sem( ############################################################################################ @testset "fiml_gradient" begin - @test test_gradient(model_ml, start_test_mean; atol = 1e-6) + test_gradient(model_ml, start_test_mean; atol = 1e-6) end @testset "fiml_gradient_symbolic" begin - @test test_gradient(model_ml_sym, start_test_mean; atol = 1e-6) + test_gradient(model_ml_sym, start_test_mean; atol = 1e-6) end ############################################################################################ From 0c4299df984af5fd2de2158cd5e47c6c67dda944 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sat, 16 Mar 2024 21:57:11 -0700 Subject: [PATCH 22/29] test_hessian(): do tests inside --- test/examples/helper.jl | 22 +++++++++---------- test/examples/political_democracy/by_parts.jl | 4 ++-- .../political_democracy/constructor.jl | 4 ++-- 3 files changed, 14 insertions(+), 16 deletions(-) diff --git a/test/examples/helper.jl b/test/examples/helper.jl index dca43695e..736df37d3 100644 --- a/test/examples/helper.jl +++ b/test/examples/helper.jl @@ -16,31 +16,29 @@ end function test_hessian(model, parameters; rtol = 1e-4, atol = 0) true_hessian = - FiniteDiff.finite_difference_hessian(x -> objective!(model, x)[1], parameters) - hessian = zeros(size(true_hessian)) - hessian .= 1.0 + FiniteDiff.finite_difference_hessian(Base.Fix1(objective!, model), parameters) + hessian = similar(parameters, size(true_hessian)) gradient = similar(parameters) # H + fill!(hessian, NaN) hessian!(hessian, model, parameters) - correct1 = isapprox(hessian, true_hessian; rtol = rtol, atol = atol) + @test hessian ≈ true_hessian rtol = rtol atol = atol # F and H - hessian .= 1.0 + fill!(hessian, NaN) objective_hessian!(hessian, model, parameters) - correct2 = isapprox(hessian, true_hessian; rtol = rtol, atol = atol) + @test hessian ≈ true_hessian rtol = rtol atol = atol # G and H - hessian .= 1.0 + fill!(hessian, NaN) gradient_hessian!(gradient, hessian, model, parameters) - correct3 = isapprox(hessian, true_hessian; rtol = rtol, atol = atol) + @test hessian ≈ true_hessian rtol = rtol atol = atol # F, G and H - hessian .= 1.0 + fill!(hessian, NaN) objective_gradient_hessian!(gradient, hessian, model, parameters) - correct4 = isapprox(hessian, true_hessian; rtol = rtol, atol = atol) - - return correct1 & correct2 & correct3 & correct4 + @test hessian ≈ true_hessian rtol = rtol atol = atol end fitmeasure_names_ml = Dict( diff --git a/test/examples/political_democracy/by_parts.jl b/test/examples/political_democracy/by_parts.jl index 4c3567572..8ca6df041 100644 --- a/test/examples/political_democracy/by_parts.jl +++ b/test/examples/political_democracy/by_parts.jl @@ -179,11 +179,11 @@ if semoptimizer == SemOptimizerOptim Sem(observed, imply_sym_hessian, loss_ml, SemOptimizerOptim(algorithm = Newton())) @testset "ml_hessians" begin - @test test_hessian(model_ml, start_test; atol = 1e-4) + test_hessian(model_ml, start_test; atol = 1e-4) end @testset "ls_hessians" begin - @test test_hessian(model_ls, start_test; atol = 1e-4) + test_hessian(model_ls, start_test; atol = 1e-4) end @testset "ml_solution_hessian" begin diff --git a/test/examples/political_democracy/constructor.jl b/test/examples/political_democracy/constructor.jl index 757b26e60..f53682c95 100644 --- a/test/examples/political_democracy/constructor.jl +++ b/test/examples/political_democracy/constructor.jl @@ -197,11 +197,11 @@ if semoptimizer == SemOptimizerOptim ) @testset "ml_hessians" begin - @test test_hessian(model_ml, start_test; atol = 1e-4) + test_hessian(model_ml, start_test; atol = 1e-4) end @testset "ls_hessians" begin - @test test_hessian(model_ls, start_test; atol = 1e-4) + test_hessian(model_ls, start_test; atol = 1e-4) end @testset "ml_solution_hessian" begin From cac6607d7f323809e0ddd30a453b01763dbf51cc Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sat, 16 Mar 2024 21:54:46 -0700 Subject: [PATCH 23/29] comp_fitmeasures() -> test_fitmeasures() --- test/examples/helper.jl | 9 ++-- test/examples/multigroup/build_models.jl | 50 ++++++++---------- test/examples/political_democracy/by_parts.jl | 52 +++++++------------ .../political_democracy/constructor.jl | 52 +++++++------------ 4 files changed, 62 insertions(+), 101 deletions(-) diff --git a/test/examples/helper.jl b/test/examples/helper.jl index 736df37d3..3bb4e217a 100644 --- a/test/examples/helper.jl +++ b/test/examples/helper.jl @@ -66,13 +66,10 @@ function test_fitmeasures( atol = 0, fitmeasure_names = fitmeasure_names_ml, ) - correct = [] - for key in keys(fitmeasure_names) - measure = measures[key] - measure_lav = measures_lav.x[measures_lav[:, 1].==fitmeasure_names[key]][1] - push!(correct, isapprox(measure, measure_lav; rtol = rtol, atol = atol)) + @testset "$name" for (key, name) in pairs(fitmeasure_names) + measure_lav = measures_lav.x[findfirst(==(name), measures_lav[!, 1])] + @test measures[key] ≈ measure_lav rtol = rtol atol = atol end - return correct end function compare_estimates( diff --git a/test/examples/multigroup/build_models.jl b/test/examples/multigroup/build_models.jl index c0038ccae..e563025a7 100644 --- a/test/examples/multigroup/build_models.jl +++ b/test/examples/multigroup/build_models.jl @@ -27,13 +27,11 @@ end @testset "fitmeasures/se_ml" begin solution_ml = sem_fit(model_ml_multigroup) - @test all( - test_fitmeasures( - fit_measures(solution_ml), - solution_lav[:fitmeasures_ml]; - rtol = 1e-2, - atol = 1e-7, - ), + test_fitmeasures( + fit_measures(solution_ml), + solution_lav[:fitmeasures_ml]; + rtol = 1e-2, + atol = 1e-7, ) update_partable!( @@ -95,13 +93,11 @@ end @testset "fitmeasures/se_ml | sorted" begin solution_ml = sem_fit(model_ml_multigroup) - @test all( - test_fitmeasures( - fit_measures(solution_ml), - solution_lav[:fitmeasures_ml]; - rtol = 1e-2, - atol = 1e-7, - ), + test_fitmeasures( + fit_measures(solution_ml), + solution_lav[:fitmeasures_ml]; + rtol = 1e-2, + atol = 1e-7, ) update_partable!( @@ -203,14 +199,12 @@ end @testset "fitmeasures/se_ls" begin solution_ls = sem_fit(model_ls_multigroup) - @test all( - test_fitmeasures( - fit_measures(solution_ls), - solution_lav[:fitmeasures_ls]; - fitmeasure_names = fitmeasure_names_ls, - rtol = 1e-2, - atol = 1e-5, - ), + test_fitmeasures( + fit_measures(solution_ls), + solution_lav[:fitmeasures_ls]; + fitmeasure_names = fitmeasure_names_ls, + rtol = 1e-2, + atol = 1e-5, ) update_partable!( @@ -297,13 +291,11 @@ if !isnothing(specification_miss_g1) @testset "fitmeasures/se_fiml" begin solution = sem_fit(model_ml_multigroup) - @test all( - test_fitmeasures( - fit_measures(solution), - solution_lav[:fitmeasures_fiml]; - rtol = 1e-3, - atol = 0, - ), + test_fitmeasures( + fit_measures(solution), + solution_lav[:fitmeasures_fiml]; + rtol = 1e-3, + atol = 0, ) update_partable!( diff --git a/test/examples/political_democracy/by_parts.jl b/test/examples/political_democracy/by_parts.jl index 8ca6df041..e260e2442 100644 --- a/test/examples/political_democracy/by_parts.jl +++ b/test/examples/political_democracy/by_parts.jl @@ -114,13 +114,7 @@ end @testset "fitmeasures/se_ml" begin solution_ml = sem_fit(model_ml) - @test all( - test_fitmeasures( - fit_measures(solution_ml), - solution_lav[:fitmeasures_ml]; - atol = 1e-3, - ), - ) + test_fitmeasures(fit_measures(solution_ml), solution_lav[:fitmeasures_ml]; atol = 1e-3) update_partable!(partable, identifier(model_ml), se_hessian(solution_ml), :se) @test compare_estimates( @@ -135,13 +129,11 @@ end @testset "fitmeasures/se_ls" begin solution_ls = sem_fit(model_ls_sym) fm = fit_measures(solution_ls) - @test all( - test_fitmeasures( - fm, - solution_lav[:fitmeasures_ls]; - atol = 1e-3, - fitmeasure_names = fitmeasure_names_ls, - ), + test_fitmeasures( + fm, + solution_lav[:fitmeasures_ls]; + atol = 1e-3, + fitmeasure_names = fitmeasure_names_ls, ) @test (fm[:AIC] === missing) & (fm[:BIC] === missing) & (fm[:minus2ll] === missing) @@ -283,12 +275,10 @@ end @testset "fitmeasures/se_ml_mean" begin solution_ml = sem_fit(model_ml) - @test all( - test_fitmeasures( - fit_measures(solution_ml), - solution_lav[:fitmeasures_ml_mean]; - atol = 1e-3, - ), + test_fitmeasures( + fit_measures(solution_ml), + solution_lav[:fitmeasures_ml_mean]; + atol = 1e-3, ) update_partable!(partable_mean, identifier(model_ml), se_hessian(solution_ml), :se) @@ -304,13 +294,11 @@ end @testset "fitmeasures/se_ls_mean" begin solution_ls = sem_fit(model_ls) fm = fit_measures(solution_ls) - @test all( - test_fitmeasures( - fm, - solution_lav[:fitmeasures_ls_mean]; - atol = 1e-3, - fitmeasure_names = fitmeasure_names_ls, - ), + test_fitmeasures( + fm, + solution_lav[:fitmeasures_ls_mean]; + atol = 1e-3, + fitmeasure_names = fitmeasure_names_ls, ) @test (fm[:AIC] === missing) & (fm[:BIC] === missing) & (fm[:minus2ll] === missing) @@ -380,12 +368,10 @@ end @testset "fitmeasures/se_fiml" begin solution_ml = sem_fit(model_ml) - @test all( - test_fitmeasures( - fit_measures(solution_ml), - solution_lav[:fitmeasures_fiml]; - atol = 1e-3, - ), + test_fitmeasures( + fit_measures(solution_ml), + solution_lav[:fitmeasures_fiml]; + atol = 1e-3, ) update_partable!(partable_mean, identifier(model_ml), se_hessian(solution_ml), :se) diff --git a/test/examples/political_democracy/constructor.jl b/test/examples/political_democracy/constructor.jl index f53682c95..b4eacf76d 100644 --- a/test/examples/political_democracy/constructor.jl +++ b/test/examples/political_democracy/constructor.jl @@ -128,13 +128,7 @@ end @testset "fitmeasures/se_ml" begin solution_ml = sem_fit(model_ml) - @test all( - test_fitmeasures( - fit_measures(solution_ml), - solution_lav[:fitmeasures_ml]; - atol = 1e-3, - ), - ) + test_fitmeasures(fit_measures(solution_ml), solution_lav[:fitmeasures_ml]; atol = 1e-3) update_partable!(partable, identifier(model_ml), se_hessian(solution_ml), :se) @test compare_estimates( @@ -149,13 +143,11 @@ end @testset "fitmeasures/se_ls" begin solution_ls = sem_fit(model_ls_sym) fm = fit_measures(solution_ls) - @test all( - test_fitmeasures( - fm, - solution_lav[:fitmeasures_ls]; - atol = 1e-3, - fitmeasure_names = fitmeasure_names_ls, - ), + test_fitmeasures( + fm, + solution_lav[:fitmeasures_ls]; + atol = 1e-3, + fitmeasure_names = fitmeasure_names_ls, ) @test (fm[:AIC] === missing) & (fm[:BIC] === missing) & (fm[:minus2ll] === missing) @@ -306,12 +298,10 @@ end @testset "fitmeasures/se_ml_mean" begin solution_ml = sem_fit(model_ml) - @test all( - test_fitmeasures( - fit_measures(solution_ml), - solution_lav[:fitmeasures_ml_mean]; - atol = 0.002, - ), + test_fitmeasures( + fit_measures(solution_ml), + solution_lav[:fitmeasures_ml_mean]; + atol = 0.002, ) update_partable!(partable_mean, identifier(model_ml), se_hessian(solution_ml), :se) @@ -327,13 +317,11 @@ end @testset "fitmeasures/se_ls_mean" begin solution_ls = sem_fit(model_ls) fm = fit_measures(solution_ls) - @test all( - test_fitmeasures( - fm, - solution_lav[:fitmeasures_ls_mean]; - atol = 1e-3, - fitmeasure_names = fitmeasure_names_ls, - ), + test_fitmeasures( + fm, + solution_lav[:fitmeasures_ls_mean]; + atol = 1e-3, + fitmeasure_names = fitmeasure_names_ls, ) @test (fm[:AIC] === missing) & (fm[:BIC] === missing) & (fm[:minus2ll] === missing) @@ -414,12 +402,10 @@ end @testset "fitmeasures/se_fiml" begin solution_ml = sem_fit(model_ml) - @test all( - test_fitmeasures( - fit_measures(solution_ml), - solution_lav[:fitmeasures_fiml]; - atol = 1e-3, - ), + test_fitmeasures( + fit_measures(solution_ml), + solution_lav[:fitmeasures_fiml]; + atol = 1e-3, ) update_partable!(partable_mean, identifier(model_ml), se_hessian(solution_ml), :se) From 19ceaa7e97fe6585bfc378fb57d54351ed45d503 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sat, 16 Mar 2024 21:57:11 -0700 Subject: [PATCH 24/29] tests: tiny improvements --- test/examples/political_democracy/by_parts.jl | 11 ++++------- 1 file changed, 4 insertions(+), 7 deletions(-) diff --git a/test/examples/political_democracy/by_parts.jl b/test/examples/political_democracy/by_parts.jl index e260e2442..2f6adfe7a 100644 --- a/test/examples/political_democracy/by_parts.jl +++ b/test/examples/political_democracy/by_parts.jl @@ -83,7 +83,7 @@ end solution_ridge = sem_fit(model_ridge) solution_ml = sem_fit(model_ml) # solution_ridge_id = sem_fit(model_ridge_id) - @test abs(solution_ridge.minimum - solution_ml.minimum) < 1 + @test solution_ridge.minimum < solution_ml.minimum + 1 end # test constant objective value @@ -100,12 +100,9 @@ end @testset "ml_solution_weighted" begin solution_ml = sem_fit(model_ml) solution_ml_weighted = sem_fit(model_ml_weighted) - @test isapprox(solution(solution_ml), solution(solution_ml_weighted), rtol = 1e-3) - @test isapprox( - n_obs(model_ml) * StructuralEquationModels.minimum(solution_ml), - StructuralEquationModels.minimum(solution_ml_weighted), - rtol = 1e-6, - ) + @test solution(solution_ml) ≈ solution(solution_ml_weighted) rtol = 1e-3 + @test n_obs(model_ml) * StructuralEquationModels.minimum(solution_ml) ≈ + StructuralEquationModels.minimum(solution_ml_weighted) rtol = 1e-6 end ############################################################################################ From 7563e2091fd00de7bc37bb592974b94e57b87b06 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Fri, 15 Mar 2024 08:36:35 -0700 Subject: [PATCH 25/29] tests: use approx op --- test/examples/recover_parameters/recover_parameters_twofact.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/examples/recover_parameters/recover_parameters_twofact.jl b/test/examples/recover_parameters/recover_parameters_twofact.jl index dc47a9055..68e44ce20 100644 --- a/test/examples/recover_parameters/recover_parameters_twofact.jl +++ b/test/examples/recover_parameters/recover_parameters_twofact.jl @@ -74,4 +74,4 @@ model_ml = Sem(semobserved, imply_ml, loss_ml, optimizer) objective!(model_ml, true_val) solution_ml = sem_fit(model_ml) -@test isapprox(true_val, solution(solution_ml); atol = 0.05) +@test true_val ≈ solution(solution_ml) atol = 0.05 From 964037a0af3459471e5cf2b8a6ea6fd502f34316 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Fri, 15 Mar 2024 12:15:47 -0700 Subject: [PATCH 26/29] tests: relax multithreading check --- test/unit_tests/multithreading.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/unit_tests/multithreading.jl b/test/unit_tests/multithreading.jl index e18e40abd..a3d8ebe9c 100644 --- a/test/unit_tests/multithreading.jl +++ b/test/unit_tests/multithreading.jl @@ -2,6 +2,6 @@ using Test if haskey(ENV, "JULIA_ON_CI") @testset "multithreading_enabled" begin - @test Threads.nthreads() == 8 + @test Threads.nthreads() >= 8 end end From 9e95a9eca90eb48ff4ed73d0d3d92a5d8b67ab00 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sun, 17 Mar 2024 00:42:47 -0700 Subject: [PATCH 27/29] tests: use ismissing() --- test/examples/political_democracy/constructor.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/examples/political_democracy/constructor.jl b/test/examples/political_democracy/constructor.jl index b4eacf76d..0e58940cb 100644 --- a/test/examples/political_democracy/constructor.jl +++ b/test/examples/political_democracy/constructor.jl @@ -149,7 +149,7 @@ end atol = 1e-3, fitmeasure_names = fitmeasure_names_ls, ) - @test (fm[:AIC] === missing) & (fm[:BIC] === missing) & (fm[:minus2ll] === missing) + @test ismissing(fm[:AIC]) && ismissing(fm[:BIC]) && ismissing(fm[:minus2ll]) update_partable!(partable, identifier(model_ls_sym), se_hessian(solution_ls), :se) @test compare_estimates( @@ -323,7 +323,7 @@ end atol = 1e-3, fitmeasure_names = fitmeasure_names_ls, ) - @test (fm[:AIC] === missing) & (fm[:BIC] === missing) & (fm[:minus2ll] === missing) + @test ismissing(fm[:AIC]) && ismissing(fm[:BIC]) && ismissing(fm[:minus2ll]) update_partable!(partable_mean, identifier(model_ls), se_hessian(solution_ls), :se) @test compare_estimates( From 115ccfbace57f5b099aed9f033665e17ad067cf5 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Tue, 19 Mar 2024 20:37:36 -0700 Subject: [PATCH 28/29] tests: use update_se_hessian!() --- src/StructuralEquationModels.jl | 1 + test/examples/multigroup/build_models.jl | 28 +++---------------- test/examples/political_democracy/by_parts.jl | 10 +++---- .../political_democracy/constructor.jl | 10 +++---- test/unit_tests/bootstrap.jl | 3 +- 5 files changed, 16 insertions(+), 36 deletions(-) diff --git a/src/StructuralEquationModels.jl b/src/StructuralEquationModels.jl index 547549ccb..048b7181c 100644 --- a/src/StructuralEquationModels.jl +++ b/src/StructuralEquationModels.jl @@ -142,6 +142,7 @@ export AbstractSem, update_partable!, update_estimate!, update_start!, + update_se_hessian!, Fixed, fixed, Start, diff --git a/test/examples/multigroup/build_models.jl b/test/examples/multigroup/build_models.jl index e563025a7..9b97300df 100644 --- a/test/examples/multigroup/build_models.jl +++ b/test/examples/multigroup/build_models.jl @@ -34,12 +34,7 @@ end atol = 1e-7, ) - update_partable!( - partable, - identifier(model_ml_multigroup), - se_hessian(solution_ml), - :se, - ) + update_se_hessian!(partable, solution_ml) @test compare_estimates( partable, solution_lav[:parameter_estimates_ml]; @@ -100,12 +95,7 @@ end atol = 1e-7, ) - update_partable!( - partable_s, - identifier(model_ml_multigroup), - se_hessian(solution_ml), - :se, - ) + update_se_hessian!(partable_s, solution_ml) @test compare_estimates( partable_s, solution_lav[:parameter_estimates_ml]; @@ -207,12 +197,7 @@ end atol = 1e-5, ) - update_partable!( - partable, - identifier(model_ls_multigroup), - se_hessian(solution_ls), - :se, - ) + update_se_hessian!(partable, solution_ls) @test compare_estimates( partable, solution_lav[:parameter_estimates_ls]; @@ -298,12 +283,7 @@ if !isnothing(specification_miss_g1) atol = 0, ) - update_partable!( - partable_miss, - identifier(model_ml_multigroup), - se_hessian(solution), - :se, - ) + update_se_hessian!(partable_miss, solution) @test compare_estimates( partable_miss, solution_lav[:parameter_estimates_fiml]; diff --git a/test/examples/political_democracy/by_parts.jl b/test/examples/political_democracy/by_parts.jl index 2f6adfe7a..c071d9e00 100644 --- a/test/examples/political_democracy/by_parts.jl +++ b/test/examples/political_democracy/by_parts.jl @@ -113,7 +113,7 @@ end solution_ml = sem_fit(model_ml) test_fitmeasures(fit_measures(solution_ml), solution_lav[:fitmeasures_ml]; atol = 1e-3) - update_partable!(partable, identifier(model_ml), se_hessian(solution_ml), :se) + update_se_hessian!(partable, solution_ml) @test compare_estimates( partable, solution_lav[:parameter_estimates_ml]; @@ -134,7 +134,7 @@ end ) @test (fm[:AIC] === missing) & (fm[:BIC] === missing) & (fm[:minus2ll] === missing) - update_partable!(partable, identifier(model_ls_sym), se_hessian(solution_ls), :se) + update_se_hessian!(partable, solution_ls) @test compare_estimates( partable, solution_lav[:parameter_estimates_ls]; @@ -278,7 +278,7 @@ end atol = 1e-3, ) - update_partable!(partable_mean, identifier(model_ml), se_hessian(solution_ml), :se) + update_se_hessian!(partable_mean, solution_ml) @test compare_estimates( partable_mean, solution_lav[:parameter_estimates_ml_mean]; @@ -299,7 +299,7 @@ end ) @test (fm[:AIC] === missing) & (fm[:BIC] === missing) & (fm[:minus2ll] === missing) - update_partable!(partable_mean, identifier(model_ls), se_hessian(solution_ls), :se) + update_se_hessian!(partable_mean, solution_ls) @test compare_estimates( partable_mean, solution_lav[:parameter_estimates_ls_mean]; @@ -371,7 +371,7 @@ end atol = 1e-3, ) - update_partable!(partable_mean, identifier(model_ml), se_hessian(solution_ml), :se) + update_se_hessian!(partable_mean, solution_ml) @test compare_estimates( partable_mean, solution_lav[:parameter_estimates_fiml]; diff --git a/test/examples/political_democracy/constructor.jl b/test/examples/political_democracy/constructor.jl index 0e58940cb..4ca1994bd 100644 --- a/test/examples/political_democracy/constructor.jl +++ b/test/examples/political_democracy/constructor.jl @@ -130,7 +130,7 @@ end solution_ml = sem_fit(model_ml) test_fitmeasures(fit_measures(solution_ml), solution_lav[:fitmeasures_ml]; atol = 1e-3) - update_partable!(partable, identifier(model_ml), se_hessian(solution_ml), :se) + update_se_hessian!(partable, solution_ml) @test compare_estimates( partable, solution_lav[:parameter_estimates_ml]; @@ -151,7 +151,7 @@ end ) @test ismissing(fm[:AIC]) && ismissing(fm[:BIC]) && ismissing(fm[:minus2ll]) - update_partable!(partable, identifier(model_ls_sym), se_hessian(solution_ls), :se) + update_se_hessian!(partable, solution_ls) @test compare_estimates( partable, solution_lav[:parameter_estimates_ls]; @@ -304,7 +304,7 @@ end atol = 0.002, ) - update_partable!(partable_mean, identifier(model_ml), se_hessian(solution_ml), :se) + update_se_hessian!(partable_mean, solution_ml) @test compare_estimates( partable_mean, solution_lav[:parameter_estimates_ml_mean]; @@ -325,7 +325,7 @@ end ) @test ismissing(fm[:AIC]) && ismissing(fm[:BIC]) && ismissing(fm[:minus2ll]) - update_partable!(partable_mean, identifier(model_ls), se_hessian(solution_ls), :se) + update_se_hessian!(partable_mean, solution_ls) @test compare_estimates( partable_mean, solution_lav[:parameter_estimates_ls_mean]; @@ -408,7 +408,7 @@ end atol = 1e-3, ) - update_partable!(partable_mean, identifier(model_ml), se_hessian(solution_ml), :se) + update_se_hessian!(partable_mean, solution_ml) @test compare_estimates( partable_mean, solution_lav[:parameter_estimates_fiml]; diff --git a/test/unit_tests/bootstrap.jl b/test/unit_tests/bootstrap.jl index deaafddce..f30092865 100644 --- a/test/unit_tests/bootstrap.jl +++ b/test/unit_tests/bootstrap.jl @@ -1,6 +1,5 @@ solution_ml = sem_fit(model_ml) bs = se_bootstrap(solution_ml; n_boot = 20) -se = se_hessian(solution_ml) -update_partable!(partable, solution_ml, se, :se) +update_se_hessian!(partable, solution_ml) update_partable!(partable, solution_ml, bs, :se_boot) From d9884ae64f31a5caecdf191c24738af6d5b304bb Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sun, 10 Mar 2024 13:14:05 -0700 Subject: [PATCH 29/29] matrix_gradient(): refactor * rename from get_matrix_derivative(): it's not a getter, and gradient is a better term * construct sparse matrix directly, which is much more efficient * parameters arg is not needed --- src/additional_functions/parameters.jl | 22 +++++++++++++--------- src/imply/RAM/generic.jl | 13 +++---------- 2 files changed, 16 insertions(+), 19 deletions(-) diff --git a/src/additional_functions/parameters.jl b/src/additional_functions/parameters.jl index 144590bd1..8d01b3747 100644 --- a/src/additional_functions/parameters.jl +++ b/src/additional_functions/parameters.jl @@ -86,15 +86,19 @@ function check_constants(M) return false end -function get_matrix_derivative(M_indices, parameters, n_long) - ∇M = [ - sparsevec(M_indices[i], ones(length(M_indices[i])), n_long) for - i in 1:length(parameters) - ] - - ∇M = reduce(hcat, ∇M) - - return ∇M +# construct length(M)×length(parameters) sparse matrix of 1s at the positions, +# where the corresponding parameter occurs in the M matrix +function matrix_gradient(M_indices::ArrayParamsMap, M_length::Integer) + rowval = reduce(vcat, M_indices) + colptr = + pushfirst!(accumulate((ptr, M_ind) -> ptr + length(M_ind), M_indices, init = 1), 1) + return SparseMatrixCSC( + M_length, + length(M_indices), + colptr, + rowval, + ones(length(rowval)), + ) end # fill M with parameters diff --git a/src/imply/RAM/generic.jl b/src/imply/RAM/generic.jl index a0c29185d..00c0d0ef9 100644 --- a/src/imply/RAM/generic.jl +++ b/src/imply/RAM/generic.jl @@ -158,8 +158,8 @@ function RAM(; I_A = similar(A_pre) if gradient - ∇A = get_matrix_derivative(A_indices, parameters, n_nod^2) - ∇S = get_matrix_derivative(S_indices, parameters, n_nod^2) + ∇A = matrix_gradient(A_indices, n_nod^2) + ∇S = matrix_gradient(S_indices, n_nod^2) else ∇A = nothing ∇S = nothing @@ -168,15 +168,8 @@ function RAM(; # μ if meanstructure has_meanstructure = Val(true) - - if gradient - ∇M = get_matrix_derivative(M_indices, parameters, n_nod) - else - ∇M = nothing - end - + ∇M = gradient ? matrix_gradient(M_indices, n_nod) : nothing μ = zeros(n_var) - else has_meanstructure = Val(false) M_indices = nothing