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/src/additional_functions/helper.jl b/src/additional_functions/helper.jl index 05253cd8c..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 @@ -37,13 +41,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 @@ -52,10 +51,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 diff --git a/src/additional_functions/parameters.jl b/src/additional_functions/parameters.jl index 1b7d4643f..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 @@ -111,97 +115,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 diff --git a/src/additional_functions/start_val/start_fabin3.jl b/src/additional_functions/start_val/start_fabin3.jl index bce6471c7..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, Σ, μ) @@ -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,101 +80,97 @@ 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])) + 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[] - 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] - else - reference = reference[1] + if length(reference) > 0 + 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 != reference - instruments = indicators[.!(indicators .∈ [[reference; indicator]])] - - s32 = Σ[instruments, indicator] - s13 = Σ[reference, 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 + if (indicator != ref) && + (parampos = get(indicator2parampos, indicator, 0)) != 0 + start_val[parampos] = calculate_lambda(ref, indicator, indicators) 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]])] - - s32 = Σ[instruments, indicator] - s13 = Σ[reference, instruments] - S33 = Σ[instruments, instruments] - - 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 + if indicator != ref + λ[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) - 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 +179,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 +191,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 diff --git a/src/imply/RAM/generic.jl b/src/imply/RAM/generic.jl index b9d6fad69..00c0d0ef9 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) @@ -159,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 @@ -169,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 @@ -232,7 +224,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,11 +253,11 @@ 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⁻¹ + 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) @@ -352,9 +345,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 = diff --git a/src/imply/RAM/symbolic.jl b/src/imply/RAM/symbolic.jl index 6ce7b6ae5..5c5e52112 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 @@ -242,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...) @@ -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) - 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)))] +# 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 diff --git a/src/loss/ML/ML.jl b/src/loss/ML/ML.jl index e6616258d..7811cda7f 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,21 +127,17 @@ 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 μ₋ = μₒ - μ μ₋ᵀΣ⁻¹ = μ₋' * Σ⁻¹ - gradient = vec(Σ⁻¹ * (I - Σₒ * Σ⁻¹ - μ₋ * μ₋ᵀΣ⁻¹))' * ∇Σ - 2 * μ₋ᵀΣ⁻¹ * ∇μ + gradient = vec(Σ⁻¹ - Σ⁻¹Σₒ * Σ⁻¹ - μ₋ᵀΣ⁻¹'μ₋ᵀΣ⁻¹)' * ∇Σ - 2 * μ₋ᵀΣ⁻¹ * ∇μ return gradient' else - gradient = (vec(Σ⁻¹) - vec(Σ⁻¹Σₒ * Σ⁻¹))' * ∇Σ + gradient = vec(Σ⁻¹ - Σ⁻¹Σₒ * Σ⁻¹)' * ∇Σ return gradient' end end @@ -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(Σ⁻¹, Σ⁻¹) * ∇Σ @@ -184,9 +172,9 @@ function hessian!( J = vec(Σ⁻¹ - Σ⁻¹ΣₒΣ⁻¹)' ∇²Σ_function!(∇²Σ, J, par) # outer - H_outer = 2 * kron(Σ⁻¹ΣₒΣ⁻¹, Σ⁻¹) - kron(Σ⁻¹, Σ⁻¹) + H_outer = kron(2Σ⁻¹ΣₒΣ⁻¹ - Σ⁻¹, Σ⁻¹) hessian = ∇Σ' * H_outer * ∇Σ - hessian += ∇²Σ + hessian .+= ∇²Σ end return hessian @@ -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(Σ⁻¹Σₒ) @@ -280,9 +266,9 @@ function objective_hessian!( J = vec(Σ⁻¹ - Σ⁻¹ΣₒΣ⁻¹)' ∇²Σ_function!(∇²Σ, J, par) # outer - H_outer = 2 * kron(Σ⁻¹ΣₒΣ⁻¹, Σ⁻¹) - kron(Σ⁻¹, Σ⁻¹) + H_outer = kron(2Σ⁻¹ΣₒΣ⁻¹ - Σ⁻¹, Σ⁻¹) hessian = ∇Σ' * H_outer * ∇Σ - hessian += ∇²Σ + hessian .+= ∇²Σ end return objective, hessian @@ -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!(Σ⁻¹Σₒ, Σ⁻¹, Σₒ) Σ⁻¹ΣₒΣ⁻¹ = Σ⁻¹Σₒ * Σ⁻¹ @@ -337,9 +320,9 @@ function gradient_hessian!( # inner ∇²Σ_function!(∇²Σ, J, par) # outer - H_outer = 2 * kron(Σ⁻¹ΣₒΣ⁻¹, Σ⁻¹) - kron(Σ⁻¹, Σ⁻¹) + H_outer = kron(2Σ⁻¹ΣₒΣ⁻¹ - Σ⁻¹, Σ⁻¹) hessian = ∇Σ' * H_outer * ∇Σ - hessian += ∇²Σ + hessian .+= ∇²Σ end return gradient', hessian @@ -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(Σ⁻¹Σₒ) @@ -398,9 +379,9 @@ function objective_gradient_hessian!( # inner ∇²Σ_function!(∇²Σ, J, par) # outer - H_outer = 2 * kron(Σ⁻¹ΣₒΣ⁻¹, Σ⁻¹) - kron(Σ⁻¹, Σ⁻¹) + H_outer = kron(2Σ⁻¹ΣₒΣ⁻¹ - Σ⁻¹, Σ⁻¹) hessian = ∇Σ' * H_outer * ∇Σ - hessian += ∇²Σ + hessian .+= ∇²Σ end return objective, gradient', hessian @@ -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,15 +492,11 @@ function gradient!( copyto!(Σ⁻¹, Σ) Σ_chol = cholesky!(Symmetric(Σ⁻¹); check = false) + isposdef(Σ_chol) || return ones(eltype(par), size(par)) + Σ⁻¹ = LinearAlgebra.inv!(Σ_chol) + mul!(Σ⁻¹Σₒ, Σ⁻¹, Σₒ) - if !isposdef(Σ_chol) - return ones(eltype(par), size(par)) - end - - Σ⁻¹ .= LinearAlgebra.inv!(Σ_chol) - #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 @@ -531,8 +504,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' @@ -562,18 +534,17 @@ 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(Σ⁻¹Σₒ) - 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 @@ -582,8 +553,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' 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 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/helper.jl b/test/examples/helper.jl index 6c212eadd..3bb4e217a 100644 --- a/test/examples/helper.jl +++ b/test/examples/helper.jl @@ -1,48 +1,44 @@ 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) 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( @@ -70,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 2f0fc749f..9b97300df 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 @@ -27,21 +27,14 @@ 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!( - 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]; @@ -71,7 +64,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) @@ -95,21 +88,14 @@ 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!( - 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]; @@ -159,7 +145,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 +173,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 @@ -203,22 +189,15 @@ 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!( - 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]; @@ -281,7 +260,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 @@ -297,21 +276,14 @@ 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!( - 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 99d7a6715..c071d9e00 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 @@ -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 ############################################################################################ @@ -114,15 +111,9 @@ 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) + update_se_hessian!(partable, solution_ml) @test compare_estimates( partable, solution_lav[:parameter_estimates_ml]; @@ -135,17 +126,15 @@ 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) - 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]; @@ -179,11 +168,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 @@ -254,7 +243,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 @@ -283,15 +272,13 @@ 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) + update_se_hessian!(partable_mean, solution_ml) @test compare_estimates( partable_mean, solution_lav[:parameter_estimates_ml_mean]; @@ -304,17 +291,15 @@ 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) - 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]; @@ -343,11 +328,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 ############################################################################################ @@ -380,15 +365,13 @@ 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) + 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 e13bc3816..4ca1994bd 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( @@ -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 @@ -128,15 +128,9 @@ 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) + update_se_hessian!(partable, solution_ml) @test compare_estimates( partable, solution_lav[:parameter_estimates_ml]; @@ -149,17 +143,15 @@ 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) + @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]; @@ -197,11 +189,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 @@ -255,7 +247,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( @@ -277,7 +269,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 @@ -306,15 +298,13 @@ 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) + update_se_hessian!(partable_mean, solution_ml) @test compare_estimates( partable_mean, solution_lav[:parameter_estimates_ml_mean]; @@ -327,17 +317,15 @@ 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) + @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]; @@ -377,11 +365,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 ############################################################################################ @@ -414,15 +402,13 @@ 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) + update_se_hessian!(partable_mean, solution_ml) @test compare_estimates( partable_mean, solution_lav[:parameter_estimates_fiml]; 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 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) 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, ) 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