From 893c9599aac9d076230a21b78c948b87d9e5305b Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sat, 16 Mar 2024 12:48:36 -0700 Subject: [PATCH 01/23] const SEM=StructuralEquationModels for resolving the scope --- src/StructuralEquationModels.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/StructuralEquationModels.jl b/src/StructuralEquationModels.jl index e9d529ded..ba6e5d58e 100644 --- a/src/StructuralEquationModels.jl +++ b/src/StructuralEquationModels.jl @@ -18,6 +18,8 @@ using LinearAlgebra, import DataFrames: DataFrame export StenoGraphs, @StenoGraph, meld +const SEM = StructuralEquationModels + # type hierarchy include("types.jl") include("objective_gradient_hessian.jl") From 18e57fac4eff855fb6a1e84217f937e69cdf2b17 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sat, 9 Mar 2024 13:10:56 -0800 Subject: [PATCH 02/23] obj_grad_hess: simplify mapreduce * destructure lambda function args * don't use zip(): mapreduce support multi-arg functions --- src/objective_gradient_hessian.jl | 60 ++++++++++++++----------------- 1 file changed, 26 insertions(+), 34 deletions(-) diff --git a/src/objective_gradient_hessian.jl b/src/objective_gradient_hessian.jl index 11bc5669f..53d68ec2c 100644 --- a/src/objective_gradient_hessian.jl +++ b/src/objective_gradient_hessian.jl @@ -87,9 +87,10 @@ end function objective!(loss::SemLoss, par, model) return mapreduce( - fun_weight -> fun_weight[2] * objective!(fun_weight[1], par, model), + (fun, weight) -> weight * objective!(fun, par, model), +, - zip(loss.functions, loss.weights), + loss.functions, + loss.weights, ) end @@ -108,19 +109,19 @@ end function objective_gradient!(gradient, loss::SemLoss, par, model) return mapreduce( - fun_weight -> - objective_gradient_wrap_(gradient, fun_weight[1], par, model, fun_weight[2]), + (fun, weight) -> objective_gradient_wrap_(gradient, fun, par, model, weight), +, - zip(loss.functions, loss.weights), + loss.functions, + loss.weights, ) end function objective_hessian!(hessian, loss::SemLoss, par, model) return mapreduce( - fun_weight -> - objective_hessian_wrap_(hessian, fun_weight[1], par, model, fun_weight[2]), + (fun, weight) -> objective_hessian_wrap_(hessian, fun, par, model, weight), +, - zip(loss.functions, loss.weights), + loss.functions, + loss.weights, ) end @@ -134,16 +135,11 @@ end function objective_gradient_hessian!(gradient, hessian, loss::SemLoss, par, model) return mapreduce( - fun_weight -> objective_gradient_hessian_wrap_( - gradient, - hessian, - fun_weight[1], - par, - model, - fun_weight[2], - ), + (fun, weight) -> + objective_gradient_hessian_wrap_(gradient, hessian, fun, par, model, weight), +, - zip(loss.functions, loss.weights), + loss.functions, + loss.weights, ) end @@ -174,9 +170,10 @@ end function objective!(ensemble::SemEnsemble, par) return mapreduce( - model_weight -> model_weight[2] * objective!(model_weight[1], par), + (model, weight) -> weight * objective!(model, par), +, - zip(ensemble.sems, ensemble.weights), + ensemble.sems, + ensemble.weights, ) end @@ -201,20 +198,20 @@ end function objective_gradient!(gradient, ensemble::SemEnsemble, par) fill!(gradient, zero(eltype(gradient))) return mapreduce( - model_weight -> - objective_gradient_wrap_(gradient, model_weight[1], par, model_weight[2]), + (model, weight) -> objective_gradient_wrap_(gradient, model, par, weight), +, - zip(ensemble.sems, ensemble.weights), + ensemble.sems, + ensemble.weights, ) end function objective_hessian!(hessian, ensemble::SemEnsemble, par) fill!(hessian, zero(eltype(hessian))) return mapreduce( - model_weight -> - objective_hessian_wrap_(hessian, model_weight[1], par, model_weight[2]), + (model, weight) -> objective_hessian_wrap_(hessian, model, par, weight), +, - zip(ensemble.sems, ensemble.weights), + ensemble.sems, + ensemble.weights, ) end @@ -236,16 +233,11 @@ function objective_gradient_hessian!(gradient, hessian, ensemble::SemEnsemble, p fill!(gradient, zero(eltype(gradient))) fill!(hessian, zero(eltype(hessian))) return mapreduce( - model_weight -> objective_gradient_hessian_wrap_( - gradient, - hessian, - model_weight[1], - par, - model, - model_weight[2], - ), + (model, weight) -> + objective_gradient_hessian_wrap_(gradient, hessian, model, par, model, weight), +, - zip(ensemble.sems, ensemble.weights), + ensemble.sems, + ensemble.weights, ) end From 5780be0f97bc2bc7113edd6e8d49c7121beb1204 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sat, 20 Apr 2024 12:41:22 -0700 Subject: [PATCH 03/23] cov_and_mean(): use StatsBase.mean_and_cov() --- src/StructuralEquationModels.jl | 1 + src/additional_functions/helper.jl | 7 ++----- 2 files changed, 3 insertions(+), 5 deletions(-) diff --git a/src/StructuralEquationModels.jl b/src/StructuralEquationModels.jl index ba6e5d58e..547549ccb 100644 --- a/src/StructuralEquationModels.jl +++ b/src/StructuralEquationModels.jl @@ -4,6 +4,7 @@ using LinearAlgebra, Optim, NLSolversBase, Statistics, + StatsBase, SparseArrays, Symbolics, NLopt, diff --git a/src/additional_functions/helper.jl b/src/additional_functions/helper.jl index 62100907d..e976a1e34 100644 --- a/src/additional_functions/helper.jl +++ b/src/additional_functions/helper.jl @@ -108,11 +108,8 @@ function sparse_outer_mul!(C, A, B::Vector, ind) #computes A*S*B -> C, where ind end function cov_and_mean(rows; corrected = false) - data = transpose(reduce(hcat, rows)) - size(rows, 1) > 1 ? obs_cov = Statistics.cov(data; corrected = corrected) : - obs_cov = reshape([0.0], 1, 1) - obs_mean = vec(Statistics.mean(data, dims = 1)) - return obs_cov, obs_mean + obs_mean, obs_cov = StatsBase.mean_and_cov(reduce(hcat, rows), 2, corrected = corrected) + return obs_cov, vec(obs_mean) end function duplication_matrix(nobs) From dae0d83e899c59203cde800b4814f00a1bf83fb0 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sat, 9 Mar 2024 13:38:48 -0800 Subject: [PATCH 04/23] ==: use && avoids full evalulation when the end result is known --- src/frontend/specification/RAMMatrices.jl | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/frontend/specification/RAMMatrices.jl b/src/frontend/specification/RAMMatrices.jl index 8124595ed..c942216a9 100644 --- a/src/frontend/specification/RAMMatrices.jl +++ b/src/frontend/specification/RAMMatrices.jl @@ -50,7 +50,7 @@ end import Base.== function ==(c1::RAMConstant, c2::RAMConstant) - res = ((c1.matrix == c2.matrix) & (c1.index == c2.index) & (c1.value == c2.value)) + res = ((c1.matrix == c2.matrix) && (c1.index == c2.index) && (c1.value == c2.value)) return res end @@ -425,13 +425,13 @@ end function ==(mat1::RAMMatrices, mat2::RAMMatrices) res = ( - (mat1.A_ind == mat2.A_ind) & - (mat1.S_ind == mat2.S_ind) & - (mat1.F_ind == mat2.F_ind) & - (mat1.M_ind == mat2.M_ind) & - (mat1.parameters == mat2.parameters) & - (mat1.colnames == mat2.colnames) & - (mat1.size_F == mat2.size_F) & + (mat1.A_ind == mat2.A_ind) && + (mat1.S_ind == mat2.S_ind) && + (mat1.F_ind == mat2.F_ind) && + (mat1.M_ind == mat2.M_ind) && + (mat1.parameters == mat2.parameters) && + (mat1.colnames == mat2.colnames) && + (mat1.size_F == mat2.size_F) && (mat1.constants == mat2.constants) ) return res From ec5310af6bde4eb7ea57462010e0f3e644db9a80 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sat, 9 Mar 2024 15:30:46 -0800 Subject: [PATCH 05/23] use Fix1 instead of anonymous function --- src/frontend/fit/standard_errors/hessian.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/frontend/fit/standard_errors/hessian.jl b/src/frontend/fit/standard_errors/hessian.jl index bb84efe71..31b14fcb7 100644 --- a/src/frontend/fit/standard_errors/hessian.jl +++ b/src/frontend/fit/standard_errors/hessian.jl @@ -17,7 +17,7 @@ function se_hessian(sem_fit::SemFit; hessian = :finitediff) hessian!(H, sem_fit.model, sem_fit.solution) elseif hessian == :finitediff H = FiniteDiff.finite_difference_hessian( - x -> objective!(sem_fit.model, x), + Base.Fix1(objective!, sem_fit.model), sem_fit.solution, ) elseif hessian == :optimizer From 47ce5a6d4f483121e5a77565835b89ccf7686ff0 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sat, 9 Mar 2024 15:30:57 -0800 Subject: [PATCH 06/23] fix typo --- src/frontend/fit/standard_errors/hessian.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/frontend/fit/standard_errors/hessian.jl b/src/frontend/fit/standard_errors/hessian.jl index 31b14fcb7..78c407d3d 100644 --- a/src/frontend/fit/standard_errors/hessian.jl +++ b/src/frontend/fit/standard_errors/hessian.jl @@ -33,7 +33,7 @@ function se_hessian(sem_fit::SemFit; hessian = :finitediff) ), ) else - throw(ArgumentError("I dont know how to compute `$hessian` standard-errors")) + throw(ArgumentError("I don't know how to compute `$hessian` standard-errors")) end invH = c * inv(H) From dabb3f9999c028e37c029080faecf11196d66aa7 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sat, 9 Mar 2024 15:31:34 -0800 Subject: [PATCH 07/23] rename vars for type stability --- src/frontend/specification/Sem.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/frontend/specification/Sem.jl b/src/frontend/specification/Sem.jl index a627b5f09..208ef3000 100644 --- a/src/frontend/specification/Sem.jl +++ b/src/frontend/specification/Sem.jl @@ -9,11 +9,11 @@ function Sem(; optimizer::D = SemOptimizerOptim, kwargs..., ) where {O, I, L, D} - kwargs = Dict{Symbol, Any}(kwargs...) + kwdict = Dict{Symbol, Any}(kwargs...) - set_field_type_kwargs!(kwargs, observed, imply, loss, optimizer, O, I, D) + set_field_type_kwargs!(kwdict, observed, imply, loss, optimizer, O, I, D) - observed, imply, loss, optimizer = get_fields!(kwargs, observed, imply, loss, optimizer) + observed, imply, loss, optimizer = get_fields!(kwdict, observed, imply, loss, optimizer) sem = Sem(observed, imply, loss, optimizer) @@ -27,11 +27,11 @@ function SemFiniteDiff(; optimizer::D = SemOptimizerOptim, kwargs..., ) where {O, I, L, D} - kwargs = Dict{Symbol, Any}(kwargs...) + kwdict = Dict{Symbol, Any}(kwargs...) - set_field_type_kwargs!(kwargs, observed, imply, loss, optimizer, O, I, D) + set_field_type_kwargs!(kwdict, observed, imply, loss, optimizer, O, I, D) - observed, imply, loss, optimizer = get_fields!(kwargs, observed, imply, loss, optimizer) + observed, imply, loss, optimizer = get_fields!(kwdict, observed, imply, loss, optimizer) sem = SemFiniteDiff(observed, imply, loss, optimizer) From 38349d832ac93fec67b247beed1282c4e1181e62 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sat, 9 Mar 2024 18:47:13 -0800 Subject: [PATCH 08/23] reorder_data(): optimize use dict for faster observed-to-specified matches --- src/observed/data.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/observed/data.jl b/src/observed/data.jl index 0dc4b07fb..974979337 100644 --- a/src/observed/data.jl +++ b/src/observed/data.jl @@ -152,8 +152,8 @@ function reorder_data(data::AbstractArray, spec_colnames, obs_colnames) if spec_colnames == obs_colnames return data else - new_position = [findall(x .== obs_colnames)[1] for x in spec_colnames] - data = data[:, new_position] - return data + obs_positions = Dict(col => i for (i, col) in enumerate(obs_colnames)) + new_positions = [obs_positions[col] for col in spec_colnames] + return data[:, new_positions] end end From 4236bd5511506df32cd04f03722c11a3b0c74887 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sat, 9 Mar 2024 18:47:51 -0800 Subject: [PATCH 09/23] SemObservedData(): cleanup code --- src/observed/data.jl | 30 ++++++++---------------------- 1 file changed, 8 insertions(+), 22 deletions(-) diff --git a/src/observed/data.jl b/src/observed/data.jl index 974979337..77bbc8e8b 100644 --- a/src/observed/data.jl +++ b/src/observed/data.jl @@ -103,28 +103,14 @@ function SemObservedData(; data = Matrix(data) end - n_obs, n_man = Float64.(size(data)) - - if compute_covariance - obs_cov = Statistics.cov(data) - else - obs_cov = nothing - end - - # if a meanstructure is needed, compute observed means - if meanstructure - obs_mean = vcat(Statistics.mean(data, dims = 1)...) - else - obs_mean = nothing - end - - if rowwise - data_rowwise = [data[i, :] for i in 1:convert(Int64, n_obs)] - else - data_rowwise = nothing - end - - return SemObservedData(data, obs_cov, obs_mean, n_man, n_obs, data_rowwise) + return SemObservedData( + data, + compute_covariance ? Statistics.cov(data) : nothing, + meanstructure ? vec(Statistics.mean(data, dims = 1)) : nothing, + Float64.(size(data, 2)), + Float64.(size(data, 1)), + rowwise ? [data[i, :] for i in axes(data, 1)] : nothing, + ) end ############################################################################################ From d3b7326445ac9774e25a5e85d4b56e56059cbe6c Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Mon, 11 Mar 2024 21:58:28 -0700 Subject: [PATCH 10/23] fix typo --- src/types.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/types.jl b/src/types.jl index 641b34c60..2e40e48df 100644 --- a/src/types.jl +++ b/src/types.jl @@ -45,7 +45,7 @@ function SemLoss(functions...; loss_weights = nothing, kwargs...) return SemLoss(functions, loss_weights) end -# weights for loss functions or models. If the weight is nothing, multiplication returs second argument +# weights for loss functions or models. If the weight is nothing, multiplication returns the second argument struct SemWeight{T} w::T end From f9fcef01bca7e0010632468755f3609364095111 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Tue, 19 Mar 2024 20:30:26 -0700 Subject: [PATCH 11/23] tiny simplification --- src/frontend/fit/fitmeasures/n_obs.jl | 2 +- src/types.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/frontend/fit/fitmeasures/n_obs.jl b/src/frontend/fit/fitmeasures/n_obs.jl index 3a403ff3d..cd4bdca30 100644 --- a/src/frontend/fit/fitmeasures/n_obs.jl +++ b/src/frontend/fit/fitmeasures/n_obs.jl @@ -13,4 +13,4 @@ n_obs(sem_fit::SemFit) = n_obs(sem_fit.model) n_obs(model::AbstractSemSingle) = n_obs(model.observed) -n_obs(model::SemEnsemble) = sum(n_obs.(model.sems)) +n_obs(model::SemEnsemble) = sum(n_obs, model.sems) diff --git a/src/types.jl b/src/types.jl index 2e40e48df..ca20424b6 100644 --- a/src/types.jl +++ b/src/types.jl @@ -170,7 +170,7 @@ function SemEnsemble(models...; optimizer = SemOptimizerOptim, weights = nothing # default weights if isnothing(weights) - nobs_total = sum(n_obs.(models)) + nobs_total = sum(n_obs, models) weights = [n_obs(model) / nobs_total for model in models] end From 8f887042290fb6bfffdc019fc1e22a6d4adde466 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Wed, 3 Apr 2024 00:45:54 -0700 Subject: [PATCH 12/23] fix typo --- src/loss/WLS/WLS.jl | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/loss/WLS/WLS.jl b/src/loss/WLS/WLS.jl index 3260c8579..18649ab68 100644 --- a/src/loss/WLS/WLS.jl +++ b/src/loss/WLS/WLS.jl @@ -10,20 +10,20 @@ Weighted least squares estimation. SemWLS(; observed, - meanstructure = false, - wls_weight_matrix = nothing, - wls_weight_matrix_mean = nothing, - approximate_hessian = false, + meanstructure = false, + wls_weight_matrix = nothing, + wls_weight_matrix_mean = nothing, + approximate_hessian = false, kwargs...) # Arguments - `observed`: the `SemObserved` part of the model - `meanstructure::Bool`: does the model have a meanstructure? - `approximate_hessian::Bool`: should the hessian be swapped for an approximation -- `wls_weight_matrix`: the weight matrix for weighted least squares. - Defaults to GLS estimation (``0.5*(D^T*kron(S,S)*D)`` where D is the duplication matrix - and S is the inverse ob the observed covariance matrix) -- `wls_weight_matrix_mean`: the weight matrix for the mean part of weighted least squares. +- `wls_weight_matrix`: the weight matrix for weighted least squares. + Defaults to GLS estimation (``0.5*(D^T*kron(S,S)*D)`` where D is the duplication matrix + and S is the inverse of the observed covariance matrix) +- `wls_weight_matrix_mean`: the weight matrix for the mean part of weighted least squares. Defaults to GLS estimation (the inverse of the observed covariance matrix) # Examples From 87ff2504b324014c67af1510f128ea9ec5c9a064 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Mon, 11 Mar 2024 22:07:56 -0700 Subject: [PATCH 13/23] reorder_obs_cov/mean(): cleanup --- src/observed/covariance.jl | 11 +++-------- 1 file changed, 3 insertions(+), 8 deletions(-) diff --git a/src/observed/covariance.jl b/src/observed/covariance.jl index aba1ef629..1bc1158be 100644 --- a/src/observed/covariance.jl +++ b/src/observed/covariance.jl @@ -107,13 +107,8 @@ function reorder_obs_cov(obs_cov, spec_colnames, obs_colnames) if spec_colnames == obs_colnames return obs_cov else - new_position = [findall(x .== obs_colnames)[1] for x in spec_colnames] - indices = reshape( - [CartesianIndex(i, j) for j in new_position for i in new_position], - size(obs_cov, 1), - size(obs_cov, 1), - ) - obs_cov = obs_cov[indices] + new_position = [findfirst(==(x), obs_colnames) for x in spec_colnames] + obs_cov = obs_cov[new_position, new_position] return obs_cov end end @@ -125,7 +120,7 @@ function reorder_obs_mean(obs_mean, spec_colnames, obs_colnames) if spec_colnames == obs_colnames return obs_mean else - new_position = [findall(x .== obs_colnames)[1] for x in spec_colnames] + new_position = [findfirst(==(x), obs_colnames) for x in spec_colnames] obs_mean = obs_mean[new_position] return obs_mean end From af939be11478c758afd5b73b2b0709e527e83caf Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sun, 17 Mar 2024 00:09:15 -0700 Subject: [PATCH 14/23] remove no-op method --- src/observed/covariance.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/observed/covariance.jl b/src/observed/covariance.jl index 1bc1158be..ea59a576a 100644 --- a/src/observed/covariance.jl +++ b/src/observed/covariance.jl @@ -76,7 +76,8 @@ function SemObservedCovariance(; if !isnothing(spec_colnames) obs_cov = reorder_obs_cov(obs_cov, spec_colnames, obs_colnames) - obs_mean = reorder_obs_mean(obs_mean, spec_colnames, obs_colnames) + isnothing(obs_mean) || + (obs_mean = reorder_obs_mean(obs_mean, spec_colnames, obs_colnames)) end n_man = Float64(size(obs_cov, 1)) @@ -114,7 +115,6 @@ function reorder_obs_cov(obs_cov, spec_colnames, obs_colnames) end # reorder means ---------------------------------------------------------------------------- -reorder_obs_mean(obs_mean::Nothing, spec_colnames, obs_colnames) = nothing function reorder_obs_mean(obs_mean, spec_colnames, obs_colnames) if spec_colnames == obs_colnames From 2803b55f8b6c301497387a9d08e3c77808373c6c Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sat, 16 Mar 2024 12:47:20 -0700 Subject: [PATCH 15/23] remove module spec --- src/diff/Empty.jl | 2 +- src/imply/empty.jl | 3 +-- 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/src/diff/Empty.jl b/src/diff/Empty.jl index 03923a851..57fa9ee98 100644 --- a/src/diff/Empty.jl +++ b/src/diff/Empty.jl @@ -34,5 +34,5 @@ update_observed(optimizer::SemOptimizerEmpty, observed::SemObserved; kwargs...) ############################################################################################ function Base.show(io::IO, struct_inst::SemOptimizerEmpty) - StructuralEquationModels.print_type_name(io, struct_inst) + print_type_name(io, struct_inst) end diff --git a/src/imply/empty.jl b/src/imply/empty.jl index a29fb3cac..ba8580d16 100644 --- a/src/imply/empty.jl +++ b/src/imply/empty.jl @@ -36,11 +36,10 @@ end function ImplyEmpty(; specification, kwargs...) ram_matrices = RAMMatrices(specification) - identifier = StructuralEquationModels.identifier(ram_matrices) n_par = length(ram_matrices.parameters) - return ImplyEmpty(identifier, n_par) + return ImplyEmpty(identifier(ram_matrices), n_par) end ############################################################################################ From caf35e25863389f88bbda9672b39d3e59a27ed9a Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sat, 20 Apr 2024 12:41:35 -0700 Subject: [PATCH 16/23] RAMMatrices: cleanup indexing params in arrays * introduce ArrayParamsMap type for clarity * annotate types of the corresp. RAMMatrices fields * remove cartesian indexing since it is not used * rename get_parameter_indices() into array_parameters_map() * use the single pass over the array for performance --- src/additional_functions/parameters.jl | 23 +++++++++++++---------- src/frontend/specification/RAMMatrices.jl | 18 +++++++++++------- 2 files changed, 24 insertions(+), 17 deletions(-) diff --git a/src/additional_functions/parameters.jl b/src/additional_functions/parameters.jl index 4305d56e7..44fa06ee3 100644 --- a/src/additional_functions/parameters.jl +++ b/src/additional_functions/parameters.jl @@ -18,14 +18,20 @@ function fill_A_S_M(A, S, M, A_indices, S_indices, M_indices, parameters) end end -function get_parameter_indices(parameters, M; linear = true, kwargs...) - M_indices = [findall(x -> (x == par), M) for par in parameters] - - if linear - M_indices = cartesian2linear.(M_indices, [M]) +# build the map from the index of the parameter to the linear indices +# of this parameter occurences in M +# returns ArrayParamsMap object +function array_parameters_map(parameters::AbstractVector, M::AbstractArray) + params_index = Dict(param => i for (i, param) in enumerate(parameters)) + T = Base.eltype(eachindex(M)) + res = [Vector{T}() for _ in eachindex(parameters)] + for (i, val) in enumerate(M) + par_ind = get(params_index, val, nothing) + if !isnothing(par_ind) + push!(res[par_ind], i) + end end - - return M_indices + return res end function eachindex_lower(M; linear_indices = false, kwargs...) @@ -49,9 +55,6 @@ function linear2cartesian(ind_lin, dims) return ind_cart end -cartesian2linear(ind_cart, A::AbstractArray) = cartesian2linear(ind_cart, size(A)) -linear2cartesian(ind_linear, A::AbstractArray) = linear2cartesian(ind_linear, size(A)) - function set_constants!(M, M_pre) for index in eachindex(M) δ = tryparse(Float64, string(M[index])) diff --git a/src/frontend/specification/RAMMatrices.jl b/src/frontend/specification/RAMMatrices.jl index c942216a9..e0fcc575c 100644 --- a/src/frontend/specification/RAMMatrices.jl +++ b/src/frontend/specification/RAMMatrices.jl @@ -2,11 +2,15 @@ ### Type ############################################################################################ +# map from parameter index to linear indices of matrix/vector positions where it occurs +AbstractArrayParamsMap = AbstractVector{<:AbstractVector{<:Integer}} +ArrayParamsMap = Vector{Vector{Int}} + struct RAMMatrices - A_ind::Any - S_ind::Any - F_ind::Any - M_ind::Any + A_ind::ArrayParamsMap + S_ind::ArrayParamsMap + F_ind::Vector{Int} + M_ind::Union{ArrayParamsMap, Nothing} parameters::Any colnames::Any constants::Any @@ -18,9 +22,9 @@ end ############################################################################################ function RAMMatrices(; A, S, F, M = nothing, parameters, colnames) - A_indices = get_parameter_indices(parameters, A) - S_indices = get_parameter_indices(parameters, S) - isnothing(M) ? M_indices = nothing : M_indices = get_parameter_indices(parameters, M) + A_indices = array_parameters_map(parameters, A) + S_indices = array_parameters_map(parameters, S) + M_indices = !isnothing(M) ? array_parameters_map(parameters, M) : nothing F_indices = findall([any(isone.(col)) for col in eachcol(F)]) constants = get_RAMConstants(A, S, M) return RAMMatrices( From 29f7f5bed7b9f67868157790d1524f17141974d3 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sat, 9 Mar 2024 13:14:09 -0800 Subject: [PATCH 17/23] fill_A_S_M(): use `@inbounds` --- src/additional_functions/parameters.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/additional_functions/parameters.jl b/src/additional_functions/parameters.jl index 44fa06ee3..b066126e3 100644 --- a/src/additional_functions/parameters.jl +++ b/src/additional_functions/parameters.jl @@ -1,5 +1,5 @@ function fill_A_S_M(A, S, M, A_indices, S_indices, M_indices, parameters) - for (iA, iS, par) in zip(A_indices, S_indices, parameters) + @inbounds for (iA, iS, par) in zip(A_indices, S_indices, parameters) for index_A in iA A[index_A] = par end @@ -10,7 +10,7 @@ function fill_A_S_M(A, S, M, A_indices, S_indices, M_indices, parameters) end if !isnothing(M) - for (iM, par) in zip(M_indices, parameters) + @inbounds for (iM, par) in zip(M_indices, parameters) for index_M in iM M[index_M] = par end From ef25f8c096434a74b84b0dbd41ddf04c3bc731c3 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sun, 10 Mar 2024 12:23:49 -0700 Subject: [PATCH 18/23] fill_A_S_M!(): add !, specify argtypes to match julia naming convention --- src/additional_functions/parameters.jl | 19 +++++++++++++++++-- src/imply/RAM/generic.jl | 6 +++--- src/imply/RAM/symbolic.jl | 2 +- 3 files changed, 21 insertions(+), 6 deletions(-) diff --git a/src/additional_functions/parameters.jl b/src/additional_functions/parameters.jl index b066126e3..1b7d4643f 100644 --- a/src/additional_functions/parameters.jl +++ b/src/additional_functions/parameters.jl @@ -1,4 +1,13 @@ -function fill_A_S_M(A, S, M, A_indices, S_indices, M_indices, parameters) +# fill A, S, and M matrices with the parameter values according to the parameters map +function fill_A_S_M!( + A::AbstractMatrix, + S::AbstractMatrix, + M::Union{AbstractVector, Nothing}, + A_indices::AbstractArrayParamsMap, + S_indices::AbstractArrayParamsMap, + M_indices::Union{AbstractArrayParamsMap, Nothing}, + parameters::AbstractVector, +) @inbounds for (iA, iS, par) in zip(A_indices, S_indices, parameters) for index_A in iA A[index_A] = par @@ -88,12 +97,18 @@ function get_matrix_derivative(M_indices, parameters, n_long) return ∇M end -function fill_matrix(M, M_indices, parameters) +# fill M with parameters +function fill_matrix!( + M::AbstractMatrix, + M_indices::AbstractArrayParamsMap, + parameters::AbstractVector, +) for (iM, par) in zip(M_indices, parameters) for index_M in iM M[index_M] = par end end + return M end function get_partition(A_indices, S_indices) diff --git a/src/imply/RAM/generic.jl b/src/imply/RAM/generic.jl index 7bba224eb..b9d6fad69 100644 --- a/src/imply/RAM/generic.jl +++ b/src/imply/RAM/generic.jl @@ -222,7 +222,7 @@ gradient!(imply::RAM, par, model::AbstractSemSingle) = # objective and gradient function objective!(imply::RAM, parameters, model, has_meanstructure::Val{T}) where {T} - fill_A_S_M( + fill_A_S_M!( imply.A, imply.S, imply.M, @@ -250,7 +250,7 @@ function gradient!( model::AbstractSemSingle, has_meanstructure::Val{T}, ) where {T} - fill_A_S_M( + fill_A_S_M!( imply.A, imply.S, imply.M, @@ -346,7 +346,7 @@ function check_acyclic(A_pre, n_par, A_indices) A_rand = copy(A_pre) randpar = rand(n_par) - fill_matrix(A_rand, A_indices, randpar) + fill_matrix!(A_rand, A_indices, randpar) # check if the model is acyclic acyclic = isone(det(I - A_rand)) diff --git a/src/imply/RAM/symbolic.jl b/src/imply/RAM/symbolic.jl index 79faf509f..6ce7b6ae5 100644 --- a/src/imply/RAM/symbolic.jl +++ b/src/imply/RAM/symbolic.jl @@ -112,7 +112,7 @@ function RAMSymbolic(; F[CartesianIndex.(1:n_var, ram_matrices.F_ind)] .= 1.0 set_RAMConstants!(A, S, M, ram_matrices.constants) - fill_A_S_M(A, S, M, ram_matrices.A_ind, ram_matrices.S_ind, ram_matrices.M_ind, par) + fill_A_S_M!(A, S, M, ram_matrices.A_ind, ram_matrices.S_ind, ram_matrices.M_ind, par) A, S, F = sparse(A), sparse(S), sparse(F) From 4794301c5c0575deafbed4d823c9f508acff0000 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sat, 20 Apr 2024 12:30:16 -0700 Subject: [PATCH 19/23] fix dangling whitespace --- src/additional_functions/helper.jl | 2 +- src/additional_functions/start_val/start_fabin3.jl | 6 +++--- src/frontend/fit/standard_errors/hessian.jl | 4 ++-- src/frontend/specification/EnsembleParameterTable.jl | 2 +- src/frontend/specification/ParameterTable.jl | 12 ++++++------ src/observed/missing.jl | 2 +- src/types.jl | 2 +- .../recover_parameters/recover_parameters_twofact.jl | 2 +- 8 files changed, 16 insertions(+), 16 deletions(-) diff --git a/src/additional_functions/helper.jl b/src/additional_functions/helper.jl index e976a1e34..05253cd8c 100644 --- a/src/additional_functions/helper.jl +++ b/src/additional_functions/helper.jl @@ -10,7 +10,7 @@ function neumann_series(mat::SparseMatrixCSC) return inverse end -#= +#= function make_onelement_array(A) isa(A, Array) ? nothing : (A = [A]) return A diff --git a/src/additional_functions/start_val/start_fabin3.jl b/src/additional_functions/start_val/start_fabin3.jl index eaf899127..bce6471c7 100644 --- a/src/additional_functions/start_val/start_fabin3.jl +++ b/src/additional_functions/start_val/start_fabin3.jl @@ -1,6 +1,6 @@ """ start_fabin3(model) - + Return a vector of FABIN 3 starting values (see Hägglund 1982). Not available for ensemble models. """ @@ -58,8 +58,8 @@ function start_fabin3(ram_matrices::RAMMatrices, Σ, μ) if !isnothing(M) in_M = length.(M_ind) .!= 0 in_any = in_A .| in_S .| in_M - else - in_any = in_A .| in_S + else + in_any = in_A .| in_S end if !all(in_any) diff --git a/src/frontend/fit/standard_errors/hessian.jl b/src/frontend/fit/standard_errors/hessian.jl index 78c407d3d..396d3b98c 100644 --- a/src/frontend/fit/standard_errors/hessian.jl +++ b/src/frontend/fit/standard_errors/hessian.jl @@ -4,9 +4,9 @@ Return hessian based standard errors. # Arguments -- `hessian`: how to compute the hessian. Options are +- `hessian`: how to compute the hessian. Options are - `:analytic`: (only if an analytic hessian for the model can be computed) - - `:finitediff`: for finite difference approximation + - `:finitediff`: for finite difference approximation """ function se_hessian(sem_fit::SemFit; hessian = :finitediff) c = H_scaling(sem_fit.model) diff --git a/src/frontend/specification/EnsembleParameterTable.jl b/src/frontend/specification/EnsembleParameterTable.jl index fcbf188fa..79283953f 100644 --- a/src/frontend/specification/EnsembleParameterTable.jl +++ b/src/frontend/specification/EnsembleParameterTable.jl @@ -27,7 +27,7 @@ function Dict(partable::EnsembleParameterTable) end #= function DataFrame( - partable::ParameterTable; + partable::ParameterTable; columns = nothing) if isnothing(columns) columns = keys(partable.columns) end out = DataFrame([key => partable.columns[key] for key in columns]) diff --git a/src/frontend/specification/ParameterTable.jl b/src/frontend/specification/ParameterTable.jl index fd3ed71da..1910d666e 100644 --- a/src/frontend/specification/ParameterTable.jl +++ b/src/frontend/specification/ParameterTable.jl @@ -224,9 +224,9 @@ update_partable!(partable::AbstractParameterTable, sem_fit::SemFit, vec, column) # update estimates ------------------------------------------------------------------------- """ update_estimate!( - partable::AbstractParameterTable, + partable::AbstractParameterTable, sem_fit::SemFit) - + Write parameter estimates from `sem_fit` to the `:estimate` column of `partable` """ update_estimate!(partable::AbstractParameterTable, sem_fit::SemFit) = @@ -236,7 +236,7 @@ update_estimate!(partable::AbstractParameterTable, sem_fit::SemFit) = """ update_start!(partable::AbstractParameterTable, sem_fit::SemFit) update_start!(partable::AbstractParameterTable, model::AbstractSem, start_val; kwargs...) - + Write starting values from `sem_fit` or `start_val` to the `:estimate` column of `partable`. # Arguments @@ -262,10 +262,10 @@ end # update partable standard errors ---------------------------------------------------------- """ update_se_hessian!( - partable::AbstractParameterTable, - sem_fit::SemFit; + partable::AbstractParameterTable, + sem_fit::SemFit; hessian = :finitediff) - + Write hessian standard errors computed for `sem_fit` to the `:se` column of `partable` # Arguments diff --git a/src/observed/missing.jl b/src/observed/missing.jl index cf1620ad9..6cfd09391 100644 --- a/src/observed/missing.jl +++ b/src/observed/missing.jl @@ -33,7 +33,7 @@ For observed data with missing values. - `get_data(::SemObservedMissing)` -> observed data - `data_rowwise(::SemObservedMissing)` -> observed data as vector per observation, with missing values deleted -- `patterns(::SemObservedMissing)` -> indices of non-missing variables per missing patterns +- `patterns(::SemObservedMissing)` -> indices of non-missing variables per missing patterns - `patterns_not(::SemObservedMissing)` -> indices of missing variables per missing pattern - `rows(::SemObservedMissing)` -> row indices of observed data points that belong to each pattern - `pattern_n_obs(::SemObservedMissing)` -> number of data points per pattern diff --git a/src/types.jl b/src/types.jl index ca20424b6..803bc733a 100644 --- a/src/types.jl +++ b/src/types.jl @@ -147,7 +147,7 @@ Constructor for ensemble models. - `weights::Vector`: Weights for each model. Defaults to the number of observed data points. All additional kwargs are passed down to the constructor for the optimizer field. - + Returns a SemEnsemble with fields - `n::Int`: Number of models. - `sems::Tuple`: `AbstractSem`s. diff --git a/test/examples/recover_parameters/recover_parameters_twofact.jl b/test/examples/recover_parameters/recover_parameters_twofact.jl index f7e555813..06f67ec92 100644 --- a/test/examples/recover_parameters/recover_parameters_twofact.jl +++ b/test/examples/recover_parameters/recover_parameters_twofact.jl @@ -66,7 +66,7 @@ semobserved = SemObservedData(data = x, specification = nothing) loss_ml = SemLoss(SEM.SemML(; observed = semobserved, n_par = length(start))) optimizer = SemOptimizerOptim( - BFGS(; linesearch = BackTracking(order = 3), alphaguess = InitialHagerZhang()),# m = 100), + BFGS(; linesearch = BackTracking(order = 3), alphaguess = InitialHagerZhang()),# m = 100), Optim.Options(; f_tol = 1e-10, x_tol = 1.5e-8), ) From e084533baa6fc4ca524eca63ce851c229641827c Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Mon, 11 Mar 2024 14:35:43 -0700 Subject: [PATCH 20/23] remove spurious "using SEM" --- src/observed/covariance.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/observed/covariance.jl b/src/observed/covariance.jl index ea59a576a..5785ba404 100644 --- a/src/observed/covariance.jl +++ b/src/observed/covariance.jl @@ -45,7 +45,7 @@ struct SemObservedCovariance{B, C, D, O} <: SemObserved n_man::D n_obs::O end -using StructuralEquationModels + function SemObservedCovariance(; specification, obs_cov, From b901d978240eced71924acd15f99bc015ef45bac Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sat, 20 Apr 2024 12:30:03 -0700 Subject: [PATCH 21/23] tests/examples: import -> using no declarations, so import is not required --- test/examples/multigroup/build_models.jl | 19 +++++++++---------- test/examples/multigroup/multigroup.jl | 13 +++++-------- .../political_democracy/constructor.jl | 2 +- .../political_democracy.jl | 2 +- .../recover_parameters_twofact.jl | 4 ++-- test/unit_tests/data_input_formats.jl | 3 +-- 6 files changed, 19 insertions(+), 24 deletions(-) diff --git a/test/examples/multigroup/build_models.jl b/test/examples/multigroup/build_models.jl index 01200f406..45fc72e6e 100644 --- a/test/examples/multigroup/build_models.jl +++ b/test/examples/multigroup/build_models.jl @@ -134,16 +134,15 @@ struct UserSemML <: SemLossFunction end ### functors ############################################################################################ -import LinearAlgebra: isposdef, logdet, tr, inv -import StructuralEquationModels: Σ, obs_cov, objective! - -function objective!(semml::UserSemML, parameters, model::AbstractSem) - let Σ = Σ(imply(model)), Σₒ = obs_cov(observed(model)) - if !isposdef(Σ) - return Inf - else - return logdet(Σ) + tr(inv(Σ) * Σₒ) - end +using LinearAlgebra: isposdef, logdet, tr, inv + +function SEM.objective!(semml::UserSemML, parameters, model::AbstractSem) + Σ = imply(model).Σ + Σₒ = SEM.obs_cov(observed(model)) + if !isposdef(Σ) + return Inf + else + return logdet(Σ) + tr(inv(Σ) * Σₒ) end end diff --git a/test/examples/multigroup/multigroup.jl b/test/examples/multigroup/multigroup.jl index d5309d40c..759c24eda 100644 --- a/test/examples/multigroup/multigroup.jl +++ b/test/examples/multigroup/multigroup.jl @@ -1,12 +1,9 @@ using StructuralEquationModels, Test, FiniteDiff -import LinearAlgebra: diagind, LowerTriangular -# import StructuralEquationModels as SEM -include( - joinpath( - chop(dirname(pathof(StructuralEquationModels)), tail = 3), - "test/examples/helper.jl", - ), -) +using LinearAlgebra: diagind, LowerTriangular + +const SEM = StructuralEquationModels + +include(joinpath(chop(dirname(pathof(SEM)), tail = 3), "test/examples/helper.jl")) dat = example_data("holzinger_swineford") dat_missing = example_data("holzinger_swineford_missing") diff --git a/test/examples/political_democracy/constructor.jl b/test/examples/political_democracy/constructor.jl index 417e199cd..e13bc3816 100644 --- a/test/examples/political_democracy/constructor.jl +++ b/test/examples/political_democracy/constructor.jl @@ -1,4 +1,4 @@ -import Statistics: cov, mean +using Statistics: cov, mean ############################################################################################ ### models w.o. meanstructure diff --git a/test/examples/political_democracy/political_democracy.jl b/test/examples/political_democracy/political_democracy.jl index 093b07588..389800745 100644 --- a/test/examples/political_democracy/political_democracy.jl +++ b/test/examples/political_democracy/political_democracy.jl @@ -1,5 +1,5 @@ using StructuralEquationModels, Test, FiniteDiff -# import StructuralEquationModels as SEM + include( joinpath( chop(dirname(pathof(StructuralEquationModels)), tail = 3), diff --git a/test/examples/recover_parameters/recover_parameters_twofact.jl b/test/examples/recover_parameters/recover_parameters_twofact.jl index 06f67ec92..dc47a9055 100644 --- a/test/examples/recover_parameters/recover_parameters_twofact.jl +++ b/test/examples/recover_parameters/recover_parameters_twofact.jl @@ -1,5 +1,5 @@ using StructuralEquationModels, Distributions, Random, Optim, LineSearches -import StructuralEquationModels as SEM + include( joinpath( chop(dirname(pathof(StructuralEquationModels)), tail = 3), @@ -63,7 +63,7 @@ Random.seed!(1234) x = transpose(rand(true_dist, 100000)) semobserved = SemObservedData(data = x, specification = nothing) -loss_ml = SemLoss(SEM.SemML(; observed = semobserved, n_par = length(start))) +loss_ml = SemLoss(SemML(; observed = semobserved, n_par = length(start))) optimizer = SemOptimizerOptim( BFGS(; linesearch = BackTracking(order = 3), alphaguess = InitialHagerZhang()),# m = 100), diff --git a/test/unit_tests/data_input_formats.jl b/test/unit_tests/data_input_formats.jl index 540e43c25..bef7ba8df 100644 --- a/test/unit_tests/data_input_formats.jl +++ b/test/unit_tests/data_input_formats.jl @@ -1,6 +1,5 @@ using StructuralEquationModels, Test, Statistics -import StructuralEquationModels: obs_cov, obs_mean, get_data - +using StructuralEquationModels: obs_cov, obs_mean, get_data ### model specification -------------------------------------------------------------------- spec = ParameterTable(nothing) From a37a9a741449a9b4913195c3502880e2bbe65f63 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Wed, 20 Mar 2024 16:41:48 -0700 Subject: [PATCH 22/23] test: use proper partable --- test/examples/multigroup/build_models.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/examples/multigroup/build_models.jl b/test/examples/multigroup/build_models.jl index 45fc72e6e..2f0fc749f 100644 --- a/test/examples/multigroup/build_models.jl +++ b/test/examples/multigroup/build_models.jl @@ -86,7 +86,7 @@ grad_fd = FiniteDiff.finite_difference_gradient( solution = sem_fit(model_ml_multigroup) update_estimate!(partable_s, solution) @test compare_estimates( - partable, + partable_s, solution_lav[:parameter_estimates_ml]; atol = 1e-4, lav_groups = Dict(:Pasteur => 1, :Grant_White => 2), From 47ae50740759f80de0aefff3371ac9ac316d495c Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sat, 20 Apr 2024 12:47:40 -0700 Subject: [PATCH 23/23] tests: turn JuliaFormatter check to warn --- test/runtests.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index a21ff9b87..c3b15475f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,7 +1,9 @@ using Test, SafeTestsets, JuliaFormatter, StructuralEquationModels @testset "JuliaFormatter.jl" begin - @test format(StructuralEquationModels; verbose = false, overwrite = false) + if !format(StructuralEquationModels; verbose = false, overwrite = false) + @warn "Julia code formatting style inconsistencies detected." + end end @time @safetestset "Unit Tests" begin include("unit_tests/unit_tests.jl")