From 22542eaa14262013bb70a95c3cd651533d237ca8 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sun, 28 Apr 2024 01:50:15 -0700 Subject: [PATCH 01/56] SemSpecification base type --- src/frontend/specification/ParameterTable.jl | 2 -- src/frontend/specification/RAMMatrices.jl | 2 +- src/types.jl | 7 +++++++ 3 files changed, 8 insertions(+), 3 deletions(-) diff --git a/src/frontend/specification/ParameterTable.jl b/src/frontend/specification/ParameterTable.jl index 1910d666e..4fc9d1513 100644 --- a/src/frontend/specification/ParameterTable.jl +++ b/src/frontend/specification/ParameterTable.jl @@ -1,5 +1,3 @@ -abstract type AbstractParameterTable end - ############################################################################################ ### Types ############################################################################################ diff --git a/src/frontend/specification/RAMMatrices.jl b/src/frontend/specification/RAMMatrices.jl index e0fcc575c..2faf84a8f 100644 --- a/src/frontend/specification/RAMMatrices.jl +++ b/src/frontend/specification/RAMMatrices.jl @@ -6,7 +6,7 @@ AbstractArrayParamsMap = AbstractVector{<:AbstractVector{<:Integer}} ArrayParamsMap = Vector{Vector{Int}} -struct RAMMatrices +struct RAMMatrices <: SemSpecification A_ind::ArrayParamsMap S_ind::ArrayParamsMap F_ind::Vector{Int} diff --git a/src/types.jl b/src/types.jl index 803bc733a..46b2781fb 100644 --- a/src/types.jl +++ b/src/types.jl @@ -247,3 +247,10 @@ loss(model::AbstractSemSingle) = model.loss Returns the optimizer part of a model. """ optimizer(model::AbstractSemSingle) = model.optimizer + +""" +Base type for all SEM specifications. +""" +abstract type SemSpecification end + +abstract type AbstractParameterTable <: SemSpecification end From 8d4187477554e4972e713006c5b8d1aea804846d Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sun, 28 Apr 2024 01:50:32 -0700 Subject: [PATCH 02/56] SemSpecification: use in methods --- src/imply/RAM/generic.jl | 2 +- src/imply/RAM/symbolic.jl | 2 +- src/observed/covariance.jl | 2 +- src/observed/data.jl | 2 +- src/observed/missing.jl | 2 +- 5 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/imply/RAM/generic.jl b/src/imply/RAM/generic.jl index 00c0d0ef9..c14121bb4 100644 --- a/src/imply/RAM/generic.jl +++ b/src/imply/RAM/generic.jl @@ -121,7 +121,7 @@ using StructuralEquationModels ############################################################################################ function RAM(; - specification, + specification::SemSpecification, #vech = false, gradient = true, meanstructure = false, diff --git a/src/imply/RAM/symbolic.jl b/src/imply/RAM/symbolic.jl index 5c5e52112..fae687e1d 100644 --- a/src/imply/RAM/symbolic.jl +++ b/src/imply/RAM/symbolic.jl @@ -88,7 +88,7 @@ end ############################################################################################ function RAMSymbolic(; - specification, + specification::SemSpecification, loss_types = nothing, vech = false, gradient = true, diff --git a/src/observed/covariance.jl b/src/observed/covariance.jl index 1b5de9fc2..8d73b1a99 100644 --- a/src/observed/covariance.jl +++ b/src/observed/covariance.jl @@ -47,7 +47,7 @@ struct SemObservedCovariance{B, C} <: SemObserved end function SemObservedCovariance(; - specification, + specification::Union{SemSpecification, Nothing}, obs_cov, obs_colnames = nothing, spec_colnames = nothing, diff --git a/src/observed/data.jl b/src/observed/data.jl index 0d9ad3a04..89deefd04 100644 --- a/src/observed/data.jl +++ b/src/observed/data.jl @@ -55,7 +55,7 @@ function check_arguments_SemObservedData(kwargs...) end function SemObservedData(; - specification, + specification::Union{SemSpecification, Nothing}, data, obs_colnames = nothing, spec_colnames = nothing, diff --git a/src/observed/missing.jl b/src/observed/missing.jl index 6cfd09391..439e3d837 100644 --- a/src/observed/missing.jl +++ b/src/observed/missing.jl @@ -86,7 +86,7 @@ end ############################################################################################ function SemObservedMissing(; - specification, + specification::Union{SemSpecification, Nothing}, data, obs_colnames = nothing, spec_colnames = nothing, From 239a9426314d3e2257e84f2fe36b65bfe85372aa Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sun, 28 Apr 2024 13:46:46 -0700 Subject: [PATCH 03/56] rename identifier -> param * identifier() -> param_indices() (Dict{Symbol, Int}) * get_identifier_indices() -> param_to_indices() (Vector{Int}) * parameters -> params (Vector{Symbol}) --- src/StructuralEquationModels.jl | 4 +- src/additional_functions/identifier.jl | 48 +++++----- src/additional_functions/parameters.jl | 16 ++-- .../start_val/start_fabin3.jl | 6 +- .../start_val/start_partable.jl | 8 +- .../start_val/start_simple.jl | 8 +- src/frontend/fit/fitmeasures/n_par.jl | 6 +- .../specification/EnsembleParameterTable.jl | 4 +- src/frontend/specification/ParameterTable.jl | 20 ++-- src/frontend/specification/RAMMatrices.jl | 69 +++++++------- src/frontend/specification/StenoGraphs.jl | 20 ++-- src/frontend/specification/documentation.jl | 14 +-- src/imply/RAM/generic.jl | 21 ++--- src/imply/RAM/symbolic.jl | 12 +-- src/imply/empty.jl | 10 +- src/loss/ML/FIML.jl | 20 ++-- src/loss/regularization/ridge.jl | 2 +- src/objective_gradient_hessian.jl | 92 +++++++++---------- src/types.jl | 14 +-- test/examples/helper.jl | 26 +++--- test/examples/multigroup/build_models.jl | 2 +- test/examples/multigroup/multigroup.jl | 4 +- .../political_democracy.jl | 4 +- .../recover_parameters_twofact.jl | 2 +- test/unit_tests/specification.jl | 14 +-- 25 files changed, 222 insertions(+), 224 deletions(-) diff --git a/src/StructuralEquationModels.jl b/src/StructuralEquationModels.jl index 113022960..3319f049b 100644 --- a/src/StructuralEquationModels.jl +++ b/src/StructuralEquationModels.jl @@ -153,9 +153,9 @@ export AbstractSem, start, Label, label, - get_identifier_indices, + params_to_indices, RAMMatrices, - identifier, + param_indices, fit_measures, AIC, BIC, diff --git a/src/additional_functions/identifier.jl b/src/additional_functions/identifier.jl index fefcc1be5..1b10357d6 100644 --- a/src/additional_functions/identifier.jl +++ b/src/additional_functions/identifier.jl @@ -1,59 +1,59 @@ ############################################################################################ -# get parameter identifier +# get a map from parameters to their indices ############################################################################################ -identifier(sem_fit::SemFit) = identifier(sem_fit.model) -identifier(model::AbstractSemSingle) = identifier(model.imply) -identifier(model::SemEnsemble) = model.identifier +param_indices(sem_fit::SemFit) = param_indices(sem_fit.model) +param_indices(model::AbstractSemSingle) = param_indices(model.imply) +param_indices(model::SemEnsemble) = model.param_indices ############################################################################################ -# construct identifier +# construct a map from parameters to indices ############################################################################################ -identifier(ram_matrices::RAMMatrices) = - Dict{Symbol, Int64}(ram_matrices.parameters .=> 1:length(ram_matrices.parameters)) -function identifier(partable::ParameterTable) - _, _, identifier = get_par_npar_identifier(partable) - return identifier +param_indices(ram_matrices::RAMMatrices) = + Dict(par => i for (i, par) in enumerate(ram_matrices.params)) +function param_indices(partable::ParameterTable) + _, _, param_indices = get_par_npar_indices(partable) + return param_indices end ############################################################################################ # get indices of a Vector of parameter labels ############################################################################################ -get_identifier_indices(parameters, identifier::Dict{Symbol, Int}) = - [identifier[par] for par in parameters] +params_to_indices(params, param_indices::Dict{Symbol, Int}) = + [param_indices[par] for par in params] -get_identifier_indices( - parameters, +params_to_indices( + params, obj::Union{SemFit, AbstractSemSingle, SemEnsemble, SemImply}, -) = get_identifier_indices(parameters, identifier(obj)) +) = params_to_indices(params, params(obj)) -function get_identifier_indices(parameters, obj::Union{ParameterTable, RAMMatrices}) +function params_to_indices(params, obj::Union{ParameterTable, RAMMatrices}) @warn "You are trying to find parameter indices from a ParameterTable or RAMMatrices object. \n If your model contains user-defined types, this may lead to wrong results. \n - To be on the safe side, try to reference parameters by labels or query the indices from - the constructed model (`get_identifier_indices(parameters, model)`)." maxlog = 1 - return get_identifier_indices(parameters, identifier(obj)) + To be on the safe side, try to reference parameters by labels or query the indices from + the constructed model (`params_to_indices(params, model)`)." maxlog = 1 + return params_to_indices(params, params(obj)) end ############################################################################################ # documentation ############################################################################################ """ - get_identifier_indices(parameters, model) + params_to_indices(params, model) -Returns the indices of `parameters`. +Returns the indices of `params`. # Arguments -- `parameters::Vector{Symbol}`: parameter labels +- `params::Vector{Symbol}`: parameter labels - `model`: either a SEM or a fitted SEM # Examples ```julia -parameter_indices = get_identifier_indices([:λ₁, λ₂], my_fitted_sem) +parameter_indices = params_to_indices([:λ₁, λ₂], my_fitted_sem) values = solution(my_fitted_sem)[parameter_indices] ``` """ -function get_identifier_indices end +function params_to_indices end diff --git a/src/additional_functions/parameters.jl b/src/additional_functions/parameters.jl index 8d01b3747..d6e8eb535 100644 --- a/src/additional_functions/parameters.jl +++ b/src/additional_functions/parameters.jl @@ -6,9 +6,9 @@ function fill_A_S_M!( A_indices::AbstractArrayParamsMap, S_indices::AbstractArrayParamsMap, M_indices::Union{AbstractArrayParamsMap, Nothing}, - parameters::AbstractVector, + params::AbstractVector, ) - @inbounds for (iA, iS, par) in zip(A_indices, S_indices, parameters) + @inbounds for (iA, iS, par) in zip(A_indices, S_indices, params) for index_A in iA A[index_A] = par end @@ -19,7 +19,7 @@ function fill_A_S_M!( end if !isnothing(M) - @inbounds for (iM, par) in zip(M_indices, parameters) + @inbounds for (iM, par) in zip(M_indices, params) for index_M in iM M[index_M] = par end @@ -30,10 +30,10 @@ end # 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)) +function array_params_map(params::AbstractVector, M::AbstractArray) + params_index = Dict(param => i for (i, param) in enumerate(params)) T = Base.eltype(eachindex(M)) - res = [Vector{T}() for _ in eachindex(parameters)] + res = [Vector{T}() for _ in eachindex(params)] for (i, val) in enumerate(M) par_ind = get(params_index, val, nothing) if !isnothing(par_ind) @@ -105,9 +105,9 @@ end function fill_matrix!( M::AbstractMatrix, M_indices::AbstractArrayParamsMap, - parameters::AbstractVector, + params::AbstractVector, ) - for (iM, par) in zip(M_indices, parameters) + for (iM, par) in zip(M_indices, params) for index_M in iM M[index_M] = par end diff --git a/src/additional_functions/start_val/start_fabin3.jl b/src/additional_functions/start_val/start_fabin3.jl index ee7dcb8cf..b56ee60a1 100644 --- a/src/additional_functions/start_val/start_fabin3.jl +++ b/src/additional_functions/start_val/start_fabin3.jl @@ -31,13 +31,13 @@ function start_fabin3(observed::SemObservedMissing, imply, optimizer, args...; k end function start_fabin3(ram_matrices::RAMMatrices, Σ, μ) - A_ind, S_ind, F_ind, M_ind, parameters = ram_matrices.A_ind, + A_ind, S_ind, F_ind, M_ind, params = ram_matrices.A_ind, ram_matrices.S_ind, ram_matrices.F_ind, ram_matrices.M_ind, - ram_matrices.parameters + ram_matrices.params - n_par = length(parameters) + n_par = length(params) start_val = zeros(n_par) n_var, n_nod = ram_matrices.size_F n_latent = n_nod - n_var diff --git a/src/additional_functions/start_val/start_partable.jl b/src/additional_functions/start_val/start_partable.jl index 01d06ac71..6fb15e365 100644 --- a/src/additional_functions/start_val/start_partable.jl +++ b/src/additional_functions/start_val/start_partable.jl @@ -1,6 +1,6 @@ """ start_parameter_table(model; parameter_table) - + Return a vector of starting values taken from `parameter_table`. """ function start_parameter_table end @@ -28,10 +28,10 @@ function start_parameter_table( ) start_val = zeros(0) - for identifier_ram in ram_matrices.parameters + for param in ram_matrices.params found = false - for (i, identifier_table) in enumerate(parameter_table.identifier) - if identifier_ram == identifier_table + for (i, param_table) in enumerate(parameter_table.params) + if param == param_table push!(start_val, parameter_table.start[i]) found = true break diff --git a/src/additional_functions/start_val/start_simple.jl b/src/additional_functions/start_val/start_simple.jl index 4c4645256..2c4f661c1 100644 --- a/src/additional_functions/start_val/start_simple.jl +++ b/src/additional_functions/start_val/start_simple.jl @@ -10,7 +10,7 @@ start_covariances_obs_lat = 0.0, start_means = 0.0, kwargs...) - + Return a vector of simple starting values. """ function start_simple end @@ -62,13 +62,13 @@ function start_simple( start_means = 0.0, kwargs..., ) - A_ind, S_ind, F_ind, M_ind, parameters = ram_matrices.A_ind, + A_ind, S_ind, F_ind, M_ind, params = ram_matrices.A_ind, ram_matrices.S_ind, ram_matrices.F_ind, ram_matrices.M_ind, - ram_matrices.parameters + ram_matrices.params - n_par = length(parameters) + n_par = length(params) start_val = zeros(n_par) n_var, n_nod = ram_matrices.size_F diff --git a/src/frontend/fit/fitmeasures/n_par.jl b/src/frontend/fit/fitmeasures/n_par.jl index 9cb2d3479..c8553572b 100644 --- a/src/frontend/fit/fitmeasures/n_par.jl +++ b/src/frontend/fit/fitmeasures/n_par.jl @@ -5,7 +5,7 @@ n_par(sem_fit::SemFit) n_par(model::AbstractSemSingle) n_par(model::SemEnsemble) - n_par(identifier::Dict) + n_par(param_indices::Dict) Return the number of parameters. """ @@ -15,6 +15,6 @@ n_par(fit::SemFit) = n_par(fit.model) n_par(model::AbstractSemSingle) = n_par(model.imply) -n_par(model::SemEnsemble) = n_par(model.identifier) +n_par(model::SemEnsemble) = n_par(model.param_indices) -n_par(identifier::Dict) = length(identifier) +n_par(param_indices::Dict) = length(param_indices) diff --git a/src/frontend/specification/EnsembleParameterTable.jl b/src/frontend/specification/EnsembleParameterTable.jl index 79283953f..29a6cf984 100644 --- a/src/frontend/specification/EnsembleParameterTable.jl +++ b/src/frontend/specification/EnsembleParameterTable.jl @@ -109,12 +109,12 @@ get_group(partable::EnsembleParameterTable, group) = get_group(partable.tables, # update generic --------------------------------------------------------------------------- function update_partable!( partable::EnsembleParameterTable, - model_identifier::AbstractDict, + param_indices::AbstractDict, vec, column, ) for k in keys(partable.tables) - update_partable!(partable.tables[k], model_identifier, vec, column) + update_partable!(partable.tables[k], param_indices, vec, column) end return partable end diff --git a/src/frontend/specification/ParameterTable.jl b/src/frontend/specification/ParameterTable.jl index 4fc9d1513..c0430625f 100644 --- a/src/frontend/specification/ParameterTable.jl +++ b/src/frontend/specification/ParameterTable.jl @@ -192,16 +192,16 @@ push!(partable::ParameterTable, d::Nothing) = nothing function update_partable!( partable::ParameterTable, - model_identifier::AbstractDict, - vec, + param_indices::AbstractDict, + values::AbstractVector, column, ) new_col = Vector{eltype(vec)}(undef, length(partable)) - for (i, identifier) in enumerate(partable.columns[:identifier]) - if !(identifier == :const) - new_col[i] = vec[model_identifier[identifier]] - elseif identifier == :const - new_col[i] = zero(eltype(vec)) + for (i, param) in enumerate(partable.columns[:identifier]) + if !(param == :const) + new_col[i] = values[param_indices[param]] + elseif param == :const + new_col[i] = zero(eltype(values)) end end push!(partable.columns, column => new_col) @@ -210,14 +210,14 @@ end """ update_partable!(partable::AbstractParameterTable, sem_fit::SemFit, vec, column) - + Write `vec` to `column` of `partable`. # Arguments - `vec::Vector`: has to be in the same order as the `model` parameters """ update_partable!(partable::AbstractParameterTable, sem_fit::SemFit, vec, column) = - update_partable!(partable, identifier(sem_fit), vec, column) + update_partable!(partable, param_indices(sem_fit), vec, column) # update estimates ------------------------------------------------------------------------- """ @@ -254,7 +254,7 @@ function update_start!( if !(start_val isa Vector) start_val = start_val(model; kwargs...) end - return update_partable!(partable, identifier(model), start_val, :start) + return update_partable!(partable, param_indices(model), start_val, :start) end # update partable standard errors ---------------------------------------------------------- diff --git a/src/frontend/specification/RAMMatrices.jl b/src/frontend/specification/RAMMatrices.jl index 2faf84a8f..eb92c889b 100644 --- a/src/frontend/specification/RAMMatrices.jl +++ b/src/frontend/specification/RAMMatrices.jl @@ -11,7 +11,7 @@ struct RAMMatrices <: SemSpecification S_ind::ArrayParamsMap F_ind::Vector{Int} M_ind::Union{ArrayParamsMap, Nothing} - parameters::Any + params::Any colnames::Any constants::Any size_F::Any @@ -21,10 +21,10 @@ end ### Constructor ############################################################################################ -function RAMMatrices(; A, S, F, M = nothing, parameters, colnames) - A_indices = array_parameters_map(parameters, A) - S_indices = array_parameters_map(parameters, S) - M_indices = !isnothing(M) ? array_parameters_map(parameters, M) : nothing +function RAMMatrices(; A, S, F, M = nothing, params, colnames) + A_indices = array_params_map(params, A) + S_indices = array_params_map(params, S) + M_indices = !isnothing(M) ? array_params_map(params, M) : nothing F_indices = findall([any(isone.(col)) for col in eachcol(F)]) constants = get_RAMConstants(A, S, M) return RAMMatrices( @@ -32,7 +32,7 @@ function RAMMatrices(; A, S, F, M = nothing, parameters, colnames) S_indices, F_indices, M_indices, - parameters, + params, colnames, constants, size(F), @@ -107,10 +107,10 @@ end function RAMMatrices(partable::ParameterTable; par_id = nothing) if isnothing(par_id) - parameters, n_par, par_positions = get_par_npar_identifier(partable) + params, n_par, par_positions = get_par_npar_indices(partable) else - parameters, n_par, par_positions = - par_id[:parameters], par_id[:n_par], par_id[:par_positions] + params, n_par, par_positions = + par_id[:params], par_id[:n_par], par_id[:par_positions] end n_observed = size(partable.variables[:observed_vars], 1) @@ -169,7 +169,7 @@ function RAMMatrices(partable::ParameterTable; par_id = nothing) constants = Vector{RAMConstant}() for i in 1:length(partable) - from, parameter_type, to, free, value_fixed, identifier = partable[i] + from, parameter_type, to, free, value_fixed, param = partable[i] row_ind = positions[to] if from != Symbol("1") @@ -191,7 +191,7 @@ function RAMMatrices(partable::ParameterTable; par_id = nothing) ) end else - par_ind = par_positions[identifier] + par_ind = par_positions[param] if (parameter_type == :→) && (from == Symbol("1")) push!(M_ind[par_ind], row_ind) elseif parameter_type == :→ @@ -210,7 +210,7 @@ function RAMMatrices(partable::ParameterTable; par_id = nothing) S_ind, F_ind, M_ind, - parameters, + params, colnames, constants, (n_observed, n_node), @@ -243,7 +243,7 @@ function ParameterTable(ram_matrices::RAMMatrices) end # parameters - for (i, par) in enumerate(ram_matrices.parameters) + for (i, par) in enumerate(ram_matrices.params) push_partable_rows!( partable, position_names, @@ -266,9 +266,9 @@ end function RAMMatrices(partable::EnsembleParameterTable) ram_matrices = Dict{Symbol, RAMMatrices}() - parameters, n_par, par_positions = get_par_npar_identifier(partable) + params, n_par, par_positions = get_par_npar_indices(partable) par_id = - Dict(:parameters => parameters, :n_par => n_par, :par_positions => par_positions) + Dict(:params => params, :n_par => n_par, :par_positions => par_positions) for key in keys(partable.tables) ram_mat = RAMMatrices(partable.tables[key]; par_id = par_id) @@ -291,27 +291,27 @@ end ### Additional Functions ############################################################################################ -function get_par_npar_identifier(partable::ParameterTable) - parameters = unique(partable.columns[:identifier]) - filter!(x -> x != :const, parameters) - n_par = length(parameters) - par_positions = Dict(parameters .=> 1:n_par) - return parameters, n_par, par_positions +function get_par_npar_indices(partable::ParameterTable) + params = unique(partable.columns[:identifier]) + filter!(x -> x != :const, params) + n_par = length(params) + par_positions = Dict(params .=> 1:n_par) + return params, n_par, par_positions end -function get_par_npar_identifier(partable::EnsembleParameterTable) - parameters = Vector{Symbol}() +function get_par_npar_indices(partable::EnsembleParameterTable) + params = Vector{Symbol}() for key in keys(partable.tables) - append!(parameters, partable.tables[key].columns[:identifier]) + append!(params, partable.tables[key].columns[:identifier]) end - parameters = unique(parameters) - filter!(x -> x != :const, parameters) + params = unique(params) + filter!(x -> x != :const, params) - n_par = length(parameters) + n_par = length(params) - par_positions = Dict(parameters .=> 1:n_par) + par_positions = Dict(params .=> 1:n_par) - return parameters, n_par, par_positions + return params, n_par, par_positions end function get_partable_row(c::RAMConstant, position_names) @@ -330,7 +330,7 @@ function get_partable_row(c::RAMConstant, position_names) value_fixed = c.value start = 0.0 estimate = 0.0 - identifier = :const + return Dict( :from => from, :parameter_type => parameter_type, @@ -339,7 +339,7 @@ function get_partable_row(c::RAMConstant, position_names) :value_fixed => value_fixed, :start => start, :estimate => estimate, - :identifier => identifier, + :identifier => :const, ) end @@ -355,7 +355,7 @@ end cartesian_is_known(index, known_indices::Nothing) = false -function get_partable_row(par, position_names, index, matrix, n_nod, known_indices) +function get_partable_row(param, position_names, index, matrix, n_nod, known_indices) # variable names if matrix == :M @@ -387,7 +387,6 @@ function get_partable_row(par, position_names, index, matrix, n_nod, known_indic value_fixed = 0.0 start = 0.0 estimate = 0.0 - identifier = par return Dict( :from => from, @@ -397,7 +396,7 @@ function get_partable_row(par, position_names, index, matrix, n_nod, known_indic :value_fixed => value_fixed, :start => start, :estimate => estimate, - :identifier => identifier, + :identifier => param, ) end @@ -433,7 +432,7 @@ function ==(mat1::RAMMatrices, mat2::RAMMatrices) (mat1.S_ind == mat2.S_ind) && (mat1.F_ind == mat2.F_ind) && (mat1.M_ind == mat2.M_ind) && - (mat1.parameters == mat2.parameters) && + (mat1.params == mat2.params) && (mat1.colnames == mat2.colnames) && (mat1.size_F == mat2.size_F) && (mat1.constants == mat2.constants) diff --git a/src/frontend/specification/StenoGraphs.jl b/src/frontend/specification/StenoGraphs.jl index 5c9ce7fdb..a581e9e5a 100644 --- a/src/frontend/specification/StenoGraphs.jl +++ b/src/frontend/specification/StenoGraphs.jl @@ -38,8 +38,8 @@ function ParameterTable(; graph, observed_vars, latent_vars, g = 1, parname = : value_fixed = zeros(n) start = zeros(n) estimate = zeros(n) - identifier = Vector{Symbol}(undef, n) - identifier .= Symbol("") + params = Vector{Symbol}(undef, n) + params .= Symbol("") # group = Vector{Symbol}(undef, n) # start_partable = zeros(Bool, n) @@ -80,7 +80,7 @@ function ParameterTable(; graph, observed_vars, latent_vars, g = 1, parname = : if modifier.value[g] == :NaN throw(DomainError(NaN, "NaN is not allowed as a parameter label.")) end - identifier[i] = modifier.value[g] + params[i] = modifier.value[g] end end end @@ -88,13 +88,13 @@ function ParameterTable(; graph, observed_vars, latent_vars, g = 1, parname = : # make identifiers for parameters that are not labeled current_id = 1 - for i in 1:length(identifier) - if (identifier[i] == Symbol("")) & free[i] - identifier[i] = Symbol(parname, :_, current_id) + for i in 1:length(params) + if (params[i] == Symbol("")) & free[i] + params[i] = Symbol(parname, :_, current_id) current_id += 1 - elseif (identifier[i] == Symbol("")) & !free[i] - identifier[i] = :const - elseif (identifier[i] != Symbol("")) & !free[i] + elseif (params[i] == Symbol("")) & !free[i] + params[i] = :const + elseif (params[i] != Symbol("")) & !free[i] @warn "You labeled a constant. Please check if the labels of your graph are correct." end end @@ -108,7 +108,7 @@ function ParameterTable(; graph, observed_vars, latent_vars, g = 1, parname = : :value_fixed => value_fixed, :start => start, :estimate => estimate, - :identifier => identifier, + :identifier => params, ), Dict( :latent_vars => latent_vars, diff --git a/src/frontend/specification/documentation.jl b/src/frontend/specification/documentation.jl index e3be49971..27bedfea1 100644 --- a/src/frontend/specification/documentation.jl +++ b/src/frontend/specification/documentation.jl @@ -14,7 +14,7 @@ Return a `ParameterTable` constructed from (1) a graph or (2) RAM matrices. - `observed_vars::Vector{Symbol}`: observed variable names - `latent_vars::Vector{Symbol}`: latent variable names - `ram_matrices::RAMMatrices`: a `RAMMatrices` object - + # Examples See the online documentation on [Model specification](@ref) and the [ParameterTable interface](@ref). @@ -54,11 +54,11 @@ function EnsembleParameterTable end (1) RAMMatrices(partable::ParameterTable) - (2) RAMMatrices(;A, S, F, M = nothing, parameters, colnames) + (2) RAMMatrices(;A, S, F, M = nothing, params, colnames) (3) RAMMatrices(partable::EnsembleParameterTable) - -Return `RAMMatrices` constructed from (1) a parameter table or (2) individual matrices. + +Return `RAMMatrices` constructed from (1) a parameter table or (2) individual matrices. (3) Return a dictionary of `RAMMatrices` from an `EnsembleParameterTable` (keys are the group names). @@ -68,7 +68,7 @@ Return `RAMMatrices` constructed from (1) a parameter table or (2) individual ma - `S`: matrix of undirected effects - `F`: filter matrix - `M`: vector of mean effects -- `parameters::Vector{Symbol}`: parameter labels +- `params::Vector{Symbol}`: parameter labels - `colnames::Vector{Symbol}`: variable names corresponding to the A, S and F matrix columns # Examples @@ -79,7 +79,7 @@ function RAMMatrices end """ fixed(args...) -Fix parameters to a certain value. +Fix parameters to a certain value. For ensemble models, multiple values (one for each submodel/group) are needed. # Examples @@ -94,7 +94,7 @@ function fixed end """ start(args...) -Define starting values for parameters. +Define starting values for parameters. For ensemble models, multiple values (one for each submodel/group) are needed. # Examples diff --git a/src/imply/RAM/generic.jl b/src/imply/RAM/generic.jl index c14121bb4..e934f8b84 100644 --- a/src/imply/RAM/generic.jl +++ b/src/imply/RAM/generic.jl @@ -34,7 +34,7 @@ and for models with a meanstructure, the model implied means are computed as ``` ## Interfaces -- `identifier(::RAM) `-> Dict containing the parameter labels and their position +- `params(::RAM) `-> Dict containing the parameter labels and their position - `n_par(::RAM)` -> Number of parameters - `Σ(::RAM)` -> model implied covariance matrix @@ -111,7 +111,7 @@ mutable struct RAM{ ∇S::S2 ∇M::S3 - identifier::D + param_indices::D end using StructuralEquationModels @@ -128,12 +128,11 @@ function RAM(; kwargs..., ) ram_matrices = RAMMatrices(specification) - identifier = StructuralEquationModels.identifier(ram_matrices) + param_indices = SEM.param_indices(ram_matrices) # get dimensions of the model - n_par = length(ram_matrices.parameters) + n_par = length(ram_matrices.params) n_var, n_nod = ram_matrices.size_F - parameters = ram_matrices.parameters F = zeros(ram_matrices.size_F) F[CartesianIndex.(1:n_var, ram_matrices.F_ind)] .= 1.0 @@ -198,7 +197,7 @@ function RAM(; ∇A, ∇S, ∇M, - identifier, + param_indices, ) end @@ -213,7 +212,7 @@ gradient!(imply::RAM, par, model::AbstractSemSingle) = gradient!(imply, par, model, imply.has_meanstructure) # objective and gradient -function objective!(imply::RAM, parameters, model, has_meanstructure::Val{T}) where {T} +function objective!(imply::RAM, params, model, has_meanstructure::Val{T}) where {T} fill_A_S_M!( imply.A, imply.S, @@ -221,7 +220,7 @@ function objective!(imply::RAM, parameters, model, has_meanstructure::Val{T}) wh imply.A_indices, imply.S_indices, imply.M_indices, - parameters, + params, ) @. imply.I_A = -imply.A @@ -239,7 +238,7 @@ end function gradient!( imply::RAM, - parameters, + params, model::AbstractSemSingle, has_meanstructure::Val{T}, ) where {T} @@ -250,7 +249,7 @@ function gradient!( imply.A_indices, imply.S_indices, imply.M_indices, - parameters, + params, ) @. imply.I_A = -imply.A @@ -281,7 +280,7 @@ objective_gradient_hessian!(imply::RAM, par, model::AbstractSemSingle, has_means ### Recommended methods ############################################################################################ -identifier(imply::RAM) = imply.identifier +param_indices(imply::RAM) = imply.param_indices n_par(imply::RAM) = imply.n_par function update_observed(imply::RAM, observed::SemObserved; kwargs...) diff --git a/src/imply/RAM/symbolic.jl b/src/imply/RAM/symbolic.jl index fae687e1d..ce381d2b4 100644 --- a/src/imply/RAM/symbolic.jl +++ b/src/imply/RAM/symbolic.jl @@ -29,7 +29,7 @@ Subtype of `SemImply` that implements the RAM notation with symbolic precomputat Subtype of `SemImply`. ## Interfaces -- `identifier(::RAMSymbolic) `-> Dict containing the parameter labels and their position +- `params(::RAMSymbolic) `-> Dict containing the parameter labels and their position - `n_par(::RAMSymbolic)` -> Number of parameters - `Σ(::RAMSymbolic)` -> model implied covariance matrix @@ -79,7 +79,7 @@ struct RAMSymbolic{F1, F2, F3, A1, A2, A3, S1, S2, S3, V, V2, F4, A4, F5, A5, D1 μ::A4 ∇μ_function::F5 ∇μ::A5 - identifier::D1 + param_indices::D1 has_meanstructure::B end @@ -98,9 +98,9 @@ function RAMSymbolic(; kwargs..., ) ram_matrices = RAMMatrices(specification) - identifier = StructuralEquationModels.identifier(ram_matrices) + param_indices = SEM.param_indices(ram_matrices) - n_par = length(ram_matrices.parameters) + n_par = length(ram_matrices.params) n_var, n_nod = ram_matrices.size_F par = (Symbolics.@variables θ[1:n_par])[1] @@ -201,7 +201,7 @@ function RAMSymbolic(; μ, ∇μ_function, ∇μ, - identifier, + param_indices, has_meanstructure, ) end @@ -240,7 +240,7 @@ objective_gradient_hessian!(imply::RAMSymbolic, par, model) = gradient!(imply, p ### Recommended methods ############################################################################################ -identifier(imply::RAMSymbolic) = imply.identifier +param_indices(imply::RAMSymbolic) = imply.param_indices n_par(imply::RAMSymbolic) = imply.n_par function update_observed(imply::RAMSymbolic, observed::SemObserved; kwargs...) diff --git a/src/imply/empty.jl b/src/imply/empty.jl index ba8580d16..1d0ea69ff 100644 --- a/src/imply/empty.jl +++ b/src/imply/empty.jl @@ -19,14 +19,14 @@ model per group and an additional model with `ImplyEmpty` and `SemRidge` for the # Extended help ## Interfaces -- `identifier(::RAMSymbolic) `-> Dict containing the parameter labels and their position +- `params(::RAMSymbolic) `-> Dict containing the parameter labels and their position - `n_par(::RAMSymbolic)` -> Number of parameters ## Implementation Subtype of `SemImply`. """ struct ImplyEmpty{V, V2} <: SemImply - identifier::V2 + param_indices::V2 n_par::V end @@ -37,9 +37,9 @@ end function ImplyEmpty(; specification, kwargs...) ram_matrices = RAMMatrices(specification) - n_par = length(ram_matrices.parameters) + n_par = length(ram_matrices.params) - return ImplyEmpty(identifier(ram_matrices), n_par) + return ImplyEmpty(param_indices(ram_matrices), n_par) end ############################################################################################ @@ -54,7 +54,7 @@ hessian!(imply::ImplyEmpty, par, model) = nothing ### Recommended methods ############################################################################################ -identifier(imply::ImplyEmpty) = imply.identifier +param_indices(imply::ImplyEmpty) = imply.param_indices n_par(imply::ImplyEmpty) = imply.n_par update_observed(imply::ImplyEmpty, observed::SemObserved; kwargs...) = imply diff --git a/src/loss/ML/FIML.jl b/src/loss/ML/FIML.jl index 1cc7c123c..18cc88289 100644 --- a/src/loss/ML/FIML.jl +++ b/src/loss/ML/FIML.jl @@ -82,20 +82,20 @@ end ### methods ############################################################################################ -function objective!(semfiml::SemFIML, parameters, model) +function objective!(semfiml::SemFIML, params, model) if !check_fiml(semfiml, model) - return non_posdef_return(parameters) + return non_posdef_return(params) end prepare_SemFIML!(semfiml, model) - objective = F_FIML(rows(observed(model)), semfiml, model, parameters) + objective = F_FIML(rows(observed(model)), semfiml, model, params) return objective / n_obs(observed(model)) end -function gradient!(semfiml::SemFIML, parameters, model) +function gradient!(semfiml::SemFIML, params, model) if !check_fiml(semfiml, model) - return ones(eltype(parameters), size(parameters)) + return ones(eltype(params), size(params)) end prepare_SemFIML!(semfiml, model) @@ -104,15 +104,15 @@ function gradient!(semfiml::SemFIML, parameters, model) return gradient end -function objective_gradient!(semfiml::SemFIML, parameters, model) +function objective_gradient!(semfiml::SemFIML, params, model) if !check_fiml(semfiml, model) - return non_posdef_return(parameters), ones(eltype(parameters), size(parameters)) + return non_posdef_return(params), ones(eltype(params), size(params)) end prepare_SemFIML!(semfiml, model) objective = - F_FIML(rows(observed(model)), semfiml, model, parameters) / n_obs(observed(model)) + F_FIML(rows(observed(model)), semfiml, model, params) / n_obs(observed(model)) gradient = ∇F_FIML(rows(observed(model)), semfiml, model) / n_obs(observed(model)) return objective, gradient @@ -174,8 +174,8 @@ function ∇F_fiml_outer(JΣ, Jμ, imply, model, semfiml) return G end -function F_FIML(rows, semfiml, model, parameters) - F = zero(eltype(parameters)) +function F_FIML(rows, semfiml, model, params) + F = zero(eltype(params)) for i in 1:size(rows, 1) F += F_one_pattern( semfiml.meandiff[i], diff --git a/src/loss/regularization/ridge.jl b/src/loss/regularization/ridge.jl index ebf3e7bfe..0d9d10b4b 100644 --- a/src/loss/regularization/ridge.jl +++ b/src/loss/regularization/ridge.jl @@ -58,7 +58,7 @@ function SemRidge(; ), ) else - which_ridge = get_identifier_indices(which_ridge, imply) + which_ridge = params_to_indices(which_ridge, imply) end end which = [CartesianIndex(x) for x in which_ridge] diff --git a/src/objective_gradient_hessian.jl b/src/objective_gradient_hessian.jl index 53d68ec2c..61b78a54f 100644 --- a/src/objective_gradient_hessian.jl +++ b/src/objective_gradient_hessian.jl @@ -2,56 +2,56 @@ # methods for AbstractSem ############################################################################################ -function objective!(model::AbstractSemSingle, parameters) - objective!(imply(model), parameters, model) - return objective!(loss(model), parameters, model) +function objective!(model::AbstractSemSingle, params) + objective!(imply(model), params, model) + return objective!(loss(model), params, model) end -function gradient!(gradient, model::AbstractSemSingle, parameters) +function gradient!(gradient, model::AbstractSemSingle, params) fill!(gradient, zero(eltype(gradient))) - gradient!(imply(model), parameters, model) - gradient!(gradient, loss(model), parameters, model) + gradient!(imply(model), params, model) + gradient!(gradient, loss(model), params, model) end -function hessian!(hessian, model::AbstractSemSingle, parameters) +function hessian!(hessian, model::AbstractSemSingle, params) fill!(hessian, zero(eltype(hessian))) - hessian!(imply(model), parameters, model) - hessian!(hessian, loss(model), parameters, model) + hessian!(imply(model), params, model) + hessian!(hessian, loss(model), params, model) end -function objective_gradient!(gradient, model::AbstractSemSingle, parameters) +function objective_gradient!(gradient, model::AbstractSemSingle, params) fill!(gradient, zero(eltype(gradient))) - objective_gradient!(imply(model), parameters, model) - objective_gradient!(gradient, loss(model), parameters, model) + objective_gradient!(imply(model), params, model) + objective_gradient!(gradient, loss(model), params, model) end -function objective_hessian!(hessian, model::AbstractSemSingle, parameters) +function objective_hessian!(hessian, model::AbstractSemSingle, params) fill!(hessian, zero(eltype(hessian))) - objective_hessian!(imply(model), parameters, model) - objective_hessian!(hessian, loss(model), parameters, model) + objective_hessian!(imply(model), params, model) + objective_hessian!(hessian, loss(model), params, model) end -function gradient_hessian!(gradient, hessian, model::AbstractSemSingle, parameters) +function gradient_hessian!(gradient, hessian, model::AbstractSemSingle, params) fill!(gradient, zero(eltype(gradient))) fill!(hessian, zero(eltype(hessian))) - gradient_hessian!(imply(model), parameters, model) - gradient_hessian!(gradient, hessian, loss(model), parameters, model) + gradient_hessian!(imply(model), params, model) + gradient_hessian!(gradient, hessian, loss(model), params, model) end function objective_gradient_hessian!( gradient, hessian, model::AbstractSemSingle, - parameters, + params, ) fill!(gradient, zero(eltype(gradient))) fill!(hessian, zero(eltype(hessian))) - objective_gradient_hessian!(imply(model), parameters, model) - return objective_gradient_hessian!(gradient, hessian, loss(model), parameters, model) + objective_gradient_hessian!(imply(model), params, model) + return objective_gradient_hessian!(gradient, hessian, loss(model), params, model) end ############################################################################################ -# methods for SemFiniteDiff +# methods for SemFiniteDiff ############################################################################################ gradient!(gradient, model::SemFiniteDiff, par) = @@ -60,25 +60,25 @@ gradient!(gradient, model::SemFiniteDiff, par) = hessian!(hessian, model::SemFiniteDiff, par) = FiniteDiff.finite_difference_hessian!(hessian, x -> objective!(model, x), par) -function objective_gradient!(gradient, model::SemFiniteDiff, parameters) - gradient!(gradient, model, parameters) - return objective!(model, parameters) +function objective_gradient!(gradient, model::SemFiniteDiff, params) + gradient!(gradient, model, params) + return objective!(model, params) end # other methods -function gradient_hessian!(gradient, hessian, model::SemFiniteDiff, parameters) - gradient!(gradient, model, parameters) - hessian!(hessian, model, parameters) +function gradient_hessian!(gradient, hessian, model::SemFiniteDiff, params) + gradient!(gradient, model, params) + hessian!(hessian, model, params) end -function objective_hessian!(hessian, model::SemFiniteDiff, parameters) - hessian!(hessian, model, parameters) - return objective!(model, parameters) +function objective_hessian!(hessian, model::SemFiniteDiff, params) + hessian!(hessian, model, params) + return objective!(model, params) end -function objective_gradient_hessian!(gradient, hessian, model::SemFiniteDiff, parameters) - hessian!(hessian, model, parameters) - return objective_gradient!(gradient, model, parameters) +function objective_gradient_hessian!(gradient, hessian, model::SemFiniteDiff, params) + hessian!(hessian, model, params) + return objective_gradient!(gradient, model, params) end ############################################################################################ @@ -341,44 +341,44 @@ end # Documentation ############################################################################################ """ - objective!(model::AbstractSem, parameters) + objective!(model::AbstractSem, params) -Returns the objective value at `parameters`. +Returns the objective value at `params`. The model object can be modified. # Implementation To implement a new `SemImply` or `SemLossFunction` subtype, you need to add a method for - objective!(newtype::MyNewType, parameters, model::AbstractSemSingle) + objective!(newtype::MyNewType, params, model::AbstractSemSingle) To implement a new `AbstractSem` subtype, you need to add a method for - objective!(model::MyNewType, parameters) + objective!(model::MyNewType, params) """ function objective! end """ - gradient!(gradient, model::AbstractSem, parameters) + gradient!(gradient, model::AbstractSem, params) -Writes the gradient value at `parameters` to `gradient`. +Writes the gradient value at `params` to `gradient`. # Implementation To implement a new `SemImply` or `SemLossFunction` type, you can add a method for - gradient!(newtype::MyNewType, parameters, model::AbstractSemSingle) + gradient!(newtype::MyNewType, params, model::AbstractSemSingle) To implement a new `AbstractSem` subtype, you can add a method for - gradient!(gradient, model::MyNewType, parameters) + gradient!(gradient, model::MyNewType, params) """ function gradient! end """ - hessian!(hessian, model::AbstractSem, parameters) + hessian!(hessian, model::AbstractSem, params) -Writes the hessian value at `parameters` to `hessian`. +Writes the hessian value at `params` to `hessian`. # Implementation To implement a new `SemImply` or `SemLossFunction` type, you can add a method for - hessian!(newtype::MyNewType, parameters, model::AbstractSemSingle) + hessian!(newtype::MyNewType, params, model::AbstractSemSingle) To implement a new `AbstractSem` subtype, you can add a method for - hessian!(hessian, model::MyNewType, parameters) + hessian!(hessian, model::MyNewType, params) """ function hessian! end diff --git a/src/types.jl b/src/types.jl index 46b2781fb..f026b2cb0 100644 --- a/src/types.jl +++ b/src/types.jl @@ -153,14 +153,14 @@ Returns a SemEnsemble with fields - `sems::Tuple`: `AbstractSem`s. - `weights::Vector`: Weights for each model. - `optimizer::SemOptimizer`: Connects the model to the optimizer. See also [`SemOptimizer`](@ref). -- `identifier::Dict`: Stores parameter labels and their position. +- `param_indices::Dict`: Stores parameter labels and their position. """ struct SemEnsemble{N, T <: Tuple, V <: AbstractVector, D, I} <: AbstractSemCollection n::N sems::T weights::V optimizer::D - identifier::I + param_indices::I end function SemEnsemble(models...; optimizer = SemOptimizerOptim, weights = nothing, kwargs...) @@ -174,11 +174,11 @@ function SemEnsemble(models...; optimizer = SemOptimizerOptim, weights = nothing weights = [n_obs(model) / nobs_total for model in models] end - # check identifier equality - id = identifier(models[1]) + # check parameters equality + par_indices = param_indices(models[1]) for model in models - if id != identifier(model) - throw(ErrorException("The identifier of your models do not match. \n + if par_indices != param_indices(model) + throw(ErrorException("The parameters of your models do not match. \n Maybe you tried to specify models of an ensemble via ParameterTables. \n In that case, you may use RAMMatrices instead.")) end @@ -189,7 +189,7 @@ function SemEnsemble(models...; optimizer = SemOptimizerOptim, weights = nothing optimizer = optimizer(; kwargs...) end - return SemEnsemble(n, models, weights, optimizer, id) + return SemEnsemble(n, models, weights, optimizer, par_indices) end """ diff --git a/test/examples/helper.jl b/test/examples/helper.jl index 3bb4e217a..4ba264e37 100644 --- a/test/examples/helper.jl +++ b/test/examples/helper.jl @@ -1,43 +1,43 @@ -function test_gradient(model, parameters; rtol = 1e-10, atol = 0) +function test_gradient(model, params; rtol = 1e-10, atol = 0) true_grad = - FiniteDiff.finite_difference_gradient(Base.Fix1(objective!, model), parameters) - gradient = similar(parameters) + FiniteDiff.finite_difference_gradient(Base.Fix1(objective!, model), params) + gradient = similar(params) # F and G fill!(gradient, NaN) - gradient!(gradient, model, parameters) + gradient!(gradient, model, params) @test gradient ≈ true_grad rtol = rtol atol = atol # only G fill!(gradient, NaN) - objective_gradient!(gradient, model, parameters) + objective_gradient!(gradient, model, params) @test gradient ≈ true_grad rtol = rtol atol = atol end -function test_hessian(model, parameters; rtol = 1e-4, atol = 0) +function test_hessian(model, params; rtol = 1e-4, atol = 0) true_hessian = - FiniteDiff.finite_difference_hessian(Base.Fix1(objective!, model), parameters) - hessian = similar(parameters, size(true_hessian)) - gradient = similar(parameters) + FiniteDiff.finite_difference_hessian(Base.Fix1(objective!, model), params) + hessian = similar(params, size(true_hessian)) + gradient = similar(params) # H fill!(hessian, NaN) - hessian!(hessian, model, parameters) + hessian!(hessian, model, params) @test hessian ≈ true_hessian rtol = rtol atol = atol # F and H fill!(hessian, NaN) - objective_hessian!(hessian, model, parameters) + objective_hessian!(hessian, model, params) @test hessian ≈ true_hessian rtol = rtol atol = atol # G and H fill!(hessian, NaN) - gradient_hessian!(gradient, hessian, model, parameters) + gradient_hessian!(gradient, hessian, model, params) @test hessian ≈ true_hessian rtol = rtol atol = atol # F, G and H fill!(hessian, NaN) - objective_gradient_hessian!(gradient, hessian, model, parameters) + objective_gradient_hessian!(gradient, hessian, model, params) @test hessian ≈ true_hessian rtol = rtol atol = atol end diff --git a/test/examples/multigroup/build_models.jl b/test/examples/multigroup/build_models.jl index 9b97300df..8913860c8 100644 --- a/test/examples/multigroup/build_models.jl +++ b/test/examples/multigroup/build_models.jl @@ -122,7 +122,7 @@ struct UserSemML <: SemLossFunction end using LinearAlgebra: isposdef, logdet, tr, inv -function SEM.objective!(semml::UserSemML, parameters, model::AbstractSem) +function SEM.objective!(semml::UserSemML, params, model::AbstractSem) Σ = imply(model).Σ Σₒ = SEM.obs_cov(observed(model)) if !isposdef(Σ) diff --git a/test/examples/multigroup/multigroup.jl b/test/examples/multigroup/multigroup.jl index 759c24eda..818f9afdc 100644 --- a/test/examples/multigroup/multigroup.jl +++ b/test/examples/multigroup/multigroup.jl @@ -56,7 +56,7 @@ specification_g1 = RAMMatrices(; A = A, S = S1, F = F, - parameters = x, + params = x, colnames = [:x1, :x2, :x3, :x4, :x5, :x6, :x7, :x8, :x9, :visual, :textual, :speed], ) @@ -64,7 +64,7 @@ specification_g2 = RAMMatrices(; A = A, S = S2, F = F, - parameters = x, + params = x, colnames = [:x1, :x2, :x3, :x4, :x5, :x6, :x7, :x8, :x9, :visual, :textual, :speed], ) diff --git a/test/examples/political_democracy/political_democracy.jl b/test/examples/political_democracy/political_democracy.jl index 389800745..86b7e89bc 100644 --- a/test/examples/political_democracy/political_democracy.jl +++ b/test/examples/political_democracy/political_democracy.jl @@ -75,7 +75,7 @@ spec = RAMMatrices(; A = A, S = S, F = F, - parameters = x, + params = x, colnames = [ :x1, :x2, @@ -107,7 +107,7 @@ spec_mean = RAMMatrices(; S = S, F = F, M = M, - parameters = x, + params = x, colnames = [ :x1, :x2, diff --git a/test/examples/recover_parameters/recover_parameters_twofact.jl b/test/examples/recover_parameters/recover_parameters_twofact.jl index 68e44ce20..1bd7136bc 100644 --- a/test/examples/recover_parameters/recover_parameters_twofact.jl +++ b/test/examples/recover_parameters/recover_parameters_twofact.jl @@ -40,7 +40,7 @@ A = [ 0 0 0 0 0 0 0 0 ] -ram_matrices = RAMMatrices(; A = A, S = S, F = F, parameters = x, colnames = nothing) +ram_matrices = RAMMatrices(; A = A, S = S, F = F, params = x, colnames = nothing) true_val = [ repeat([1], 8) diff --git a/test/unit_tests/specification.jl b/test/unit_tests/specification.jl index 0bfc0de2d..42ad5e431 100644 --- a/test/unit_tests/specification.jl +++ b/test/unit_tests/specification.jl @@ -3,18 +3,18 @@ @test ram_matrices == RAMMatrices(partable) end -@test get_identifier_indices([:x2, :x10, :x28], model_ml) == [2, 10, 28] +@test params_to_indices([:x2, :x10, :x28], model_ml) == [2, 10, 28] -@testset "get_identifier_indices" begin +@testset "params_to_indices" begin pars = [:θ_1, :θ_7, :θ_21] - @test get_identifier_indices(pars, model_ml) == get_identifier_indices(pars, partable) - @test get_identifier_indices(pars, model_ml) == - get_identifier_indices(pars, RAMMatrices(partable)) + @test params_to_indices(pars, model_ml) == params_to_indices(pars, partable) + @test params_to_indices(pars, model_ml) == + params_to_indices(pars, RAMMatrices(partable)) end # from docstrings: -parameter_indices = get_identifier_indices([:λ₁, λ₂], my_fitted_sem) -values = solution(my_fitted_sem)[parameter_indices] +param_indices = params_to_indices([:λ₁, λ₂], my_fitted_sem) +values = solution(my_fitted_sem)[param_indices] graph = @StenoGraph begin # measurement model From 61282bc00e70da6cf2f928fec7210af7561fc44f Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sun, 28 Apr 2024 13:23:59 -0700 Subject: [PATCH 04/56] ParTable: columns[:identifier] => columns[:param] --- src/frontend/fit/summary.jl | 10 +++++----- src/frontend/specification/ParameterTable.jl | 8 ++++---- src/frontend/specification/RAMMatrices.jl | 8 ++++---- src/frontend/specification/StenoGraphs.jl | 2 +- 4 files changed, 14 insertions(+), 14 deletions(-) diff --git a/src/frontend/fit/summary.jl b/src/frontend/fit/summary.jl index 4cda902d7..1d75eb826 100644 --- a/src/frontend/fit/summary.jl +++ b/src/frontend/fit/summary.jl @@ -91,7 +91,7 @@ function sem_summary( printstyled("Loadings: \n"; color = color) print("\n") - sorted_columns = [:to, :estimate, :identifier, :value_fixed, :start] + sorted_columns = [:to, :estimate, :param, :value_fixed, :start] loading_columns = sort_partially(sorted_columns, columns) header_cols = copy(loading_columns) replace!(header_cols, :parameter_type => :type) @@ -140,7 +140,7 @@ function sem_summary( ) sorted_columns = - [:from, :parameter_type, :to, :estimate, :identifier, :value_fixed, :start] + [:from, :parameter_type, :to, :estimate, :param, :value_fixed, :start] regression_columns = sort_partially(sorted_columns, columns) regression_array = reduce( @@ -168,7 +168,7 @@ function sem_summary( ) sorted_columns = - [:from, :parameter_type, :to, :estimate, :identifier, :value_fixed, :start] + [:from, :parameter_type, :to, :estimate, :param, :value_fixed, :start] variance_columns = sort_partially(sorted_columns, columns) variance_array = reduce( @@ -196,7 +196,7 @@ function sem_summary( ) sorted_columns = - [:from, :parameter_type, :to, :estimate, :identifier, :value_fixed, :start] + [:from, :parameter_type, :to, :estimate, :param, :value_fixed, :start] variance_columns = sort_partially(sorted_columns, columns) variance_array = reduce( @@ -225,7 +225,7 @@ function sem_summary( printstyled("Means: \n"; color = color) sorted_columns = - [:from, :parameter_type, :to, :estimate, :identifier, :value_fixed, :start] + [:from, :parameter_type, :to, :estimate, :param, :value_fixed, :start] variance_columns = sort_partially(sorted_columns, columns) variance_array = reduce( diff --git a/src/frontend/specification/ParameterTable.jl b/src/frontend/specification/ParameterTable.jl index c0430625f..aff0380f9 100644 --- a/src/frontend/specification/ParameterTable.jl +++ b/src/frontend/specification/ParameterTable.jl @@ -21,7 +21,7 @@ function ParameterTable(::Nothing) :value_fixed => Vector{Float64}(), :start => Vector{Float64}(), :estimate => Vector{Float64}(), - :identifier => Vector{Symbol}(), + :param => Vector{Symbol}(), :start => Vector{Float64}(), ) @@ -66,7 +66,7 @@ function Base.show(io::IO, partable::ParameterTable) :start, :estimate, :se, - :identifier, + :param, ] existing_columns = [haskey(partable.columns, key) for key in relevant_columns] @@ -102,7 +102,7 @@ Base.getindex(partable::ParameterTable, i::Int) = ( partable.columns[:to][i], partable.columns[:free][i], partable.columns[:value_fixed][i], - partable.columns[:identifier][i], + partable.columns[:param][i], ) function Base.length(partable::ParameterTable) @@ -197,7 +197,7 @@ function update_partable!( column, ) new_col = Vector{eltype(vec)}(undef, length(partable)) - for (i, param) in enumerate(partable.columns[:identifier]) + for (i, param) in enumerate(partable.columns[:param]) if !(param == :const) new_col[i] = values[param_indices[param]] elseif param == :const diff --git a/src/frontend/specification/RAMMatrices.jl b/src/frontend/specification/RAMMatrices.jl index eb92c889b..f6e56b950 100644 --- a/src/frontend/specification/RAMMatrices.jl +++ b/src/frontend/specification/RAMMatrices.jl @@ -292,7 +292,7 @@ end ############################################################################################ function get_par_npar_indices(partable::ParameterTable) - params = unique(partable.columns[:identifier]) + params = unique(partable.columns[:param]) filter!(x -> x != :const, params) n_par = length(params) par_positions = Dict(params .=> 1:n_par) @@ -302,7 +302,7 @@ end function get_par_npar_indices(partable::EnsembleParameterTable) params = Vector{Symbol}() for key in keys(partable.tables) - append!(params, partable.tables[key].columns[:identifier]) + append!(params, partable.tables[key].columns[:param]) end params = unique(params) filter!(x -> x != :const, params) @@ -339,7 +339,7 @@ function get_partable_row(c::RAMConstant, position_names) :value_fixed => value_fixed, :start => start, :estimate => estimate, - :identifier => :const, + :param => :const, ) end @@ -396,7 +396,7 @@ function get_partable_row(param, position_names, index, matrix, n_nod, known_ind :value_fixed => value_fixed, :start => start, :estimate => estimate, - :identifier => param, + :param => param, ) end diff --git a/src/frontend/specification/StenoGraphs.jl b/src/frontend/specification/StenoGraphs.jl index a581e9e5a..bebbdeb2e 100644 --- a/src/frontend/specification/StenoGraphs.jl +++ b/src/frontend/specification/StenoGraphs.jl @@ -108,7 +108,7 @@ function ParameterTable(; graph, observed_vars, latent_vars, g = 1, parname = : :value_fixed => value_fixed, :start => start, :estimate => estimate, - :identifier => params, + :param => params, ), Dict( :latent_vars => latent_vars, From a0e76cc5ea6b0675568f4c1b8893561c0d72fbfc Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sat, 16 Mar 2024 18:41:13 -0700 Subject: [PATCH 05/56] getindex(EnsParTable, i) instead of get_group() --- src/frontend/specification/EnsembleParameterTable.jl | 4 +--- src/frontend/specification/RAMMatrices.jl | 4 ---- 2 files changed, 1 insertion(+), 7 deletions(-) diff --git a/src/frontend/specification/EnsembleParameterTable.jl b/src/frontend/specification/EnsembleParameterTable.jl index 29a6cf984..0865a4dcc 100644 --- a/src/frontend/specification/EnsembleParameterTable.jl +++ b/src/frontend/specification/EnsembleParameterTable.jl @@ -98,9 +98,7 @@ end push!(partable::EnsembleParameterTable, d::Nothing, group) = nothing -# get group -------------------------------------------------------------------------------- - -get_group(partable::EnsembleParameterTable, group) = get_group(partable.tables, group) +Base.getindex(partable::EnsembleParameterTable, group) = partable.tables[group] ############################################################################################ ### Update Partable from Fitted Model diff --git a/src/frontend/specification/RAMMatrices.jl b/src/frontend/specification/RAMMatrices.jl index f6e56b950..db52defcd 100644 --- a/src/frontend/specification/RAMMatrices.jl +++ b/src/frontend/specification/RAMMatrices.jl @@ -439,7 +439,3 @@ function ==(mat1::RAMMatrices, mat2::RAMMatrices) ) return res end - -function get_group(d::Dict, group) - return d[group] -end From e177e29c73194c423ac60ec3b912690fb44629d4 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sat, 16 Mar 2024 18:40:35 -0700 Subject: [PATCH 06/56] replace no-op ctors with convert(T, obj) convert() is a proper method to call to avoid unnecessary construction, ctor semantics requires that a new object is constructed --- src/frontend/specification/EnsembleParameterTable.jl | 6 ++---- src/frontend/specification/RAMMatrices.jl | 7 +++++-- src/imply/RAM/generic.jl | 2 +- src/imply/RAM/symbolic.jl | 2 +- 4 files changed, 9 insertions(+), 8 deletions(-) diff --git a/src/frontend/specification/EnsembleParameterTable.jl b/src/frontend/specification/EnsembleParameterTable.jl index 0865a4dcc..24a9a295a 100644 --- a/src/frontend/specification/EnsembleParameterTable.jl +++ b/src/frontend/specification/EnsembleParameterTable.jl @@ -20,10 +20,8 @@ end ### Convert to other types ############################################################################################ -import Base.Dict - -function Dict(partable::EnsembleParameterTable) - return partable.tables +function Base.convert(::Type{Dict}, partable::EnsembleParameterTable) + return convert(Dict, partable.tables) end #= function DataFrame( diff --git a/src/frontend/specification/RAMMatrices.jl b/src/frontend/specification/RAMMatrices.jl index db52defcd..6752d7c6b 100644 --- a/src/frontend/specification/RAMMatrices.jl +++ b/src/frontend/specification/RAMMatrices.jl @@ -39,8 +39,6 @@ function RAMMatrices(; A, S, F, M = nothing, params, colnames) ) end -RAMMatrices(a::RAMMatrices) = a - ############################################################################################ ### Constants ############################################################################################ @@ -217,6 +215,8 @@ function RAMMatrices(partable::ParameterTable; par_id = nothing) ) end +Base.convert(::Type{RAMMatrices}, partable::ParameterTable) = RAMMatrices(partable) + ############################################################################################ ### get parameter table from RAMMatrices ############################################################################################ @@ -259,6 +259,9 @@ function ParameterTable(ram_matrices::RAMMatrices) return partable end +Base.convert(::Type{<:ParameterTable}, ram_matrices::RAMMatrices) = + ParameterTable(ram_matrices) + ############################################################################################ ### get RAMMatrices from EnsembleParameterTable ############################################################################################ diff --git a/src/imply/RAM/generic.jl b/src/imply/RAM/generic.jl index e934f8b84..d43c8378e 100644 --- a/src/imply/RAM/generic.jl +++ b/src/imply/RAM/generic.jl @@ -127,7 +127,7 @@ function RAM(; meanstructure = false, kwargs..., ) - ram_matrices = RAMMatrices(specification) + ram_matrices = convert(RAMMatrices, specification) param_indices = SEM.param_indices(ram_matrices) # get dimensions of the model diff --git a/src/imply/RAM/symbolic.jl b/src/imply/RAM/symbolic.jl index ce381d2b4..0fe9c29bb 100644 --- a/src/imply/RAM/symbolic.jl +++ b/src/imply/RAM/symbolic.jl @@ -97,7 +97,7 @@ function RAMSymbolic(; approximate_hessian = false, kwargs..., ) - ram_matrices = RAMMatrices(specification) + ram_matrices = convert(RAMMatrices, specification) param_indices = SEM.param_indices(ram_matrices) n_par = length(ram_matrices.params) From 161dc584adbb5ea1fbc1a0de1583ebd92369058d Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Mon, 27 May 2024 14:29:44 -0700 Subject: [PATCH 07/56] ParamTable: convert vars from Dict to fields make the type immutable --- src/frontend/fit/summary.jl | 34 +++--- src/frontend/specification/ParameterTable.jl | 107 ++++++++----------- src/frontend/specification/RAMMatrices.jl | 83 ++++++-------- src/loss/ML/FIML.jl | 4 +- src/observed/get_colnames.jl | 15 +-- test/examples/helper.jl | 7 +- test/unit_tests/data_input_formats.jl | 7 +- 7 files changed, 106 insertions(+), 151 deletions(-) diff --git a/src/frontend/fit/summary.jl b/src/frontend/fit/summary.jl index 1d75eb826..4d6bc6181 100644 --- a/src/frontend/fit/summary.jl +++ b/src/frontend/fit/summary.jl @@ -60,19 +60,18 @@ function sem_summary( ) print("\n") printstyled("Latent variables: "; color = color) - for var in partable.variables[:latent_vars] + for var in partable.latent_vars print("$var ") end print("\n") printstyled("Observed variables: "; color = color) - for var in partable.variables[:observed_vars] + for var in partable.observed_vars print("$var ") end print("\n") - if haskey(partable.variables, :sorted_vars) && - (length(partable.variables[:sorted_vars]) > 0) + if length(partable.sorted_vars) > 0 printstyled("Sorted variables: "; color = color) - for var in partable.variables[:sorted_vars] + for var in partable.sorted_vars print("$var ") end print("\n") @@ -96,11 +95,11 @@ function sem_summary( header_cols = copy(loading_columns) replace!(header_cols, :parameter_type => :type) - for var in partable.variables[:latent_vars] + for var in partable.latent_vars indicator_indices = findall( (partable.columns[:from] .== var) .& (partable.columns[:parameter_type] .== :→) .& - (partable.columns[:to] .∈ [partable.variables[:observed_vars]]), + (partable.columns[:to] .∈ [partable.observed_vars]), ) loading_array = reduce( hcat, @@ -125,16 +124,16 @@ function sem_summary( regression_indices = findall( (partable.columns[:parameter_type] .== :→) .& ( ( - (partable.columns[:to] .∈ [partable.variables[:observed_vars]]) .& - (partable.columns[:from] .∈ [partable.variables[:observed_vars]]) + (partable.columns[:to] .∈ [partable.observed_vars]) .& + (partable.columns[:from] .∈ [partable.observed_vars]) ) .| ( - (partable.columns[:to] .∈ [partable.variables[:latent_vars]]) .& - (partable.columns[:from] .∈ [partable.variables[:observed_vars]]) + (partable.columns[:to] .∈ [partable.latent_vars]) .& + (partable.columns[:from] .∈ [partable.observed_vars]) ) .| ( - (partable.columns[:to] .∈ [partable.variables[:latent_vars]]) .& - (partable.columns[:from] .∈ [partable.variables[:latent_vars]]) + (partable.columns[:to] .∈ [partable.latent_vars]) .& + (partable.columns[:from] .∈ [partable.latent_vars]) ) ), ) @@ -266,19 +265,18 @@ function sem_summary( print("\n") let partable = partable.tables[[keys(partable.tables)...][1]] printstyled("Latent variables: "; color = color) - for var in partable.variables[:latent_vars] + for var in partable.latent_vars print("$var ") end print("\n") printstyled("Observed variables: "; color = color) - for var in partable.variables[:observed_vars] + for var in partable.observed_vars print("$var ") end print("\n") - if haskey(partable.variables, :sorted_vars) && - (length(partable.variables[:sorted_vars]) > 0) + if length(partable.sorted_vars) > 0 printstyled("Sorted variables: "; color = color) - for var in partable.variables[:sorted_vars] + for var in partable.sorted_vars print("$var ") end print("\n") diff --git a/src/frontend/specification/ParameterTable.jl b/src/frontend/specification/ParameterTable.jl index aff0380f9..49ea4664b 100644 --- a/src/frontend/specification/ParameterTable.jl +++ b/src/frontend/specification/ParameterTable.jl @@ -2,9 +2,11 @@ ### Types ############################################################################################ -mutable struct ParameterTable{C, V} <: AbstractParameterTable +struct ParameterTable{C} <: AbstractParameterTable columns::C - variables::V + observed_vars::Vector{Symbol} + latent_vars::Vector{Symbol} + sorted_vars::Vector{Symbol} end ############################################################################################ @@ -12,7 +14,10 @@ end ############################################################################################ # constuct an empty table -function ParameterTable(::Nothing) +function ParameterTable(; + observed_vars::Union{AbstractVector{Symbol}, Nothing} = nothing, + latent_vars::Union{AbstractVector{Symbol}, Nothing} = nothing, +) columns = Dict{Symbol, Any}( :from => Vector{Symbol}(), :parameter_type => Vector{Symbol}(), @@ -25,13 +30,10 @@ function ParameterTable(::Nothing) :start => Vector{Float64}(), ) - variables = Dict{Symbol, Any}( - :latent_vars => Vector{Symbol}(), - :observed_vars => Vector{Symbol}(), - :sorted_vars => Vector{Symbol}(), - ) - - return ParameterTable(columns, variables) + return ParameterTable(columns, + !isnothing(observed_vars) ? copy(observed_vars) : Vector{Symbol}(), + !isnothing(latent_vars) ? copy(latent_vars) : Vector{Symbol}(), + Vector{Symbol}()) end ############################################################################################ @@ -68,26 +70,21 @@ function Base.show(io::IO, partable::ParameterTable) :se, :param, ] - existing_columns = [haskey(partable.columns, key) for key in relevant_columns] + shown_columns = filter!( + col -> haskey(partable.columns, col) && length(partable.columns[col]) > 0, + relevant_columns, + ) - as_matrix = - hcat([partable.columns[key] for key in relevant_columns[existing_columns]]...) + as_matrix = mapreduce(col -> partable.columns[col], hcat, shown_columns) pretty_table( io, as_matrix, - header = ( - relevant_columns[existing_columns], - eltype.([partable.columns[key] for key in relevant_columns[existing_columns]]), - ), + header = (shown_columns, [eltype(partable.columns[col]) for col in shown_columns]), tf = PrettyTables.tf_compact, ) - if haskey(partable.variables, :latent_vars) - print(io, "Latent Variables: $(partable.variables[:latent_vars]) \n") - end - if haskey(partable.variables, :observed_vars) - print(io, "Observed Variables: $(partable.variables[:observed_vars]) \n") - end + print(io, "Latent Variables: $(partable.latent_vars) \n") + print(io, "Observed Variables: $(partable.observed_vars) \n") end ############################################################################################ @@ -96,7 +93,7 @@ end # Iteration -------------------------------------------------------------------------------- -Base.getindex(partable::ParameterTable, i::Int) = ( +Base.getindex(partable::ParameterTable, i::Integer) = ( partable.columns[:from][i], partable.columns[:parameter_type][i], partable.columns[:to][i], @@ -105,14 +102,7 @@ Base.getindex(partable::ParameterTable, i::Int) = ( partable.columns[:param][i], ) -function Base.length(partable::ParameterTable) - len = missing - for key in keys(partable.columns) - len = length(partable.columns[key]) - break - end - return len -end +Base.length(partable::ParameterTable) = length(first(values(partable.columns))) # Sorting ---------------------------------------------------------------------------------- @@ -122,51 +112,46 @@ end Base.showerror(io::IO, e::CyclicModelError) = print(io, e.msg) -import Base.sort!, Base.sort - -function sort!(partable::ParameterTable) - variables = [partable.variables[:latent_vars]; partable.variables[:observed_vars]] +function Base.sort!(partable::ParameterTable) + vars = [ + partable.latent_vars + partable.observed_vars + ] - is_regression = - (partable.columns[:parameter_type] .== :→) .& - (partable.columns[:from] .!= Symbol("1")) + is_regression = [ + (partype == :→) && (from != Symbol("1")) for + (partype, from) in zip(partable.columns[:parameter_type], partable.columns[:from]) + ] to = partable.columns[:to][is_regression] from = partable.columns[:from][is_regression] - sorted_variables = Vector{Symbol}() + sorted_vars = Vector{Symbol}() - sorted = false - while !sorted + while !isempty(vars) acyclic = false - for (i, variable) in enumerate(variables) - if !(variable ∈ to) - push!(sorted_variables, variable) - deleteat!(variables, i) - delete_edges = from .!= variable + for (i, var) in enumerate(vars) + if !(var ∈ to) + push!(sorted_vars, var) + deleteat!(vars, i) + delete_edges = from .!= var to = to[delete_edges] from = from[delete_edges] acyclic = true end end - if !acyclic + acyclic || throw(CyclicModelError("your model is cyclic and therefore can not be ordered")) - end - acyclic = false - - if length(variables) == 0 - sorted = true - end end - push!(partable.variables, :sorted_vars => sorted_variables) + copyto!(resize!(partable.sorted_vars, length(sorted_vars)), sorted_vars) return partable end -function sort(partable::ParameterTable) +function Base.sort(partable::ParameterTable) new_partable = deepcopy(partable) sort!(new_partable) return new_partable @@ -174,15 +159,13 @@ end # add a row -------------------------------------------------------------------------------- -import Base.push! - -function push!(partable::ParameterTable, d::AbstractDict) - for key in keys(d) - push!(partable.columns[key], d[key]) +function Base.push!(partable::ParameterTable, d::AbstractDict{Symbol}) + for (key, val) in pairs(d) + push!(partable.columns[key], val) end end -push!(partable::ParameterTable, d::Nothing) = nothing +Base.push!(partable::ParameterTable, d::Nothing) = nothing ############################################################################################ ### Update Partable from Fitted Model diff --git a/src/frontend/specification/RAMMatrices.jl b/src/frontend/specification/RAMMatrices.jl index 6752d7c6b..45bdfe57b 100644 --- a/src/frontend/specification/RAMMatrices.jl +++ b/src/frontend/specification/RAMMatrices.jl @@ -111,34 +111,24 @@ function RAMMatrices(partable::ParameterTable; par_id = nothing) par_id[:params], par_id[:n_par], par_id[:par_positions] end - n_observed = size(partable.variables[:observed_vars], 1) - n_latent = size(partable.variables[:latent_vars], 1) + n_observed = length(partable.observed_vars) + n_latent = length(partable.latent_vars) n_node = n_observed + n_latent # F indices - if length(partable.variables[:sorted_vars]) != 0 - F_ind = findall( - x -> x ∈ partable.variables[:observed_vars], - partable.variables[:sorted_vars], - ) - else - F_ind = 1:n_observed - end + F_ind = + length(partable.sorted_vars) != 0 ? + findall(∈(Set(partable.observed_vars)), partable.sorted_vars) : + 1:n_observed # indices of the colnames - if length(partable.variables[:sorted_vars]) != 0 - positions = - Dict(zip(partable.variables[:sorted_vars], collect(1:n_observed+n_latent))) - colnames = copy(partable.variables[:sorted_vars]) - else - positions = Dict( - zip( - [partable.variables[:observed_vars]; partable.variables[:latent_vars]], - collect(1:n_observed+n_latent), - ), - ) - colnames = [partable.variables[:observed_vars]; partable.variables[:latent_vars]] - end + colnames = + length(partable.sorted_vars) != 0 ? copy(partable.sorted_vars) : + [ + partable.observed_vars + partable.latent_vars + ] + col_indices = Dict(col => i for (i, col) in enumerate(colnames)) # fill Matrices # known_labels = Dict{Symbol, Int64}() @@ -154,51 +144,48 @@ function RAMMatrices(partable::ParameterTable; par_id = nothing) end # is there a meanstructure? - if any(partable.columns[:from] .== Symbol("1")) - M_ind = Vector{Vector{Int64}}(undef, n_par) - for i in 1:length(M_ind) - M_ind[i] = Vector{Int64}() - end - else - M_ind = nothing - end + M_ind = + any(==(Symbol("1")), partable.columns[:from]) ? [Vector{Int64}() for _ in 1:n_par] : + nothing - # handel constants + # handle constants constants = Vector{RAMConstant}() for i in 1:length(partable) from, parameter_type, to, free, value_fixed, param = partable[i] - row_ind = positions[to] - if from != Symbol("1") - col_ind = positions[from] - end + row_ind = col_indices[to] + col_ind = from != Symbol("1") ? col_indices[from] : nothing if !free - if (parameter_type == :→) & (from == Symbol("1")) + if (parameter_type == :→) && (from == Symbol("1")) push!(constants, RAMConstant(:M, row_ind, value_fixed)) elseif (parameter_type == :→) push!( constants, RAMConstant(:A, CartesianIndex(row_ind, col_ind), value_fixed), ) - else + elseif (parameter_type == :↔) push!( constants, RAMConstant(:S, CartesianIndex(row_ind, col_ind), value_fixed), ) + else + error("Unsupported parameter type: $(parameter_type)") end else par_ind = par_positions[param] if (parameter_type == :→) && (from == Symbol("1")) push!(M_ind[par_ind], row_ind) elseif parameter_type == :→ - push!(A_ind[par_ind], (row_ind + (col_ind - 1) * n_node)) - else + push!(A_ind[par_ind], row_ind + (col_ind - 1) * n_node) + elseif parameter_type == :↔ push!(S_ind[par_ind], row_ind + (col_ind - 1) * n_node) if row_ind != col_ind push!(S_ind[par_ind], col_ind + (row_ind - 1) * n_node) end + else + error("Unsupported parameter type: $(parameter_type)") end end end @@ -222,21 +209,15 @@ Base.convert(::Type{RAMMatrices}, partable::ParameterTable) = RAMMatrices(partab ############################################################################################ function ParameterTable(ram_matrices::RAMMatrices) - partable = ParameterTable(nothing) - colnames = ram_matrices.colnames - position_names = Dict{Int64, Symbol}(1:length(colnames) .=> colnames) - - # observed and latent variables - names_obs = colnames[ram_matrices.F_ind] - names_lat = colnames[findall(x -> !(x ∈ ram_matrices.F_ind), 1:length(colnames))] - partable.variables = Dict( - :sorted_vars => Vector{Symbol}(), - :observed_vars => names_obs, - :latent_vars => names_lat, + partable = ParameterTable( + observed_vars = colnames[ram_matrices.F_ind], + latent_vars = colnames[setdiff(eachindex(colnames), ram_matrices.F_ind)], ) + position_names = Dict{Int64, Symbol}(1:length(colnames) .=> colnames) + # constants for c in ram_matrices.constants push!(partable, get_partable_row(c, position_names)) diff --git a/src/loss/ML/FIML.jl b/src/loss/ML/FIML.jl index 18cc88289..9f5103f1f 100644 --- a/src/loss/ML/FIML.jl +++ b/src/loss/ML/FIML.jl @@ -252,5 +252,5 @@ end get_n_nodes(specification::RAMMatrices) = specification.size_F[2] get_n_nodes(specification::ParameterTable) = - length(specification.variables[:observed_vars]) + - length(specification.variables[:latent_vars]) + length(specification.observed_vars) + + length(specification.latent_vars) diff --git a/src/observed/get_colnames.jl b/src/observed/get_colnames.jl index d620de659..b8d89c3d0 100644 --- a/src/observed/get_colnames.jl +++ b/src/observed/get_colnames.jl @@ -1,15 +1,8 @@ -# specification colnames +# specification colnames (only observed) function get_colnames(specification::ParameterTable) - if !haskey(specification.variables, :sorted_vars) || - (length(specification.variables[:sorted_vars]) == 0) - colnames = specification.variables[:observed_vars] - else - is_obs = [ - var ∈ specification.variables[:observed_vars] for - var in specification.variables[:sorted_vars] - ] - colnames = specification.variables[:sorted_vars][is_obs] - end + colnames = + isempty(specification.sorted_vars) ? specification.observed_vars : + filter(in(Set(specification.observed_vars)), specification.sorted_vars) return colnames end diff --git a/test/examples/helper.jl b/test/examples/helper.jl index 4ba264e37..e93a1437d 100644 --- a/test/examples/helper.jl +++ b/test/examples/helper.jl @@ -118,8 +118,7 @@ function compare_estimates( if type == :↔ type = "~~" elseif type == :→ - if (from ∈ partable.variables[:latent_vars]) & - (to ∈ partable.variables[:observed_vars]) + if (from ∈ partable.latent_vars) && (to ∈ partable.observed_vars) type = "=~" else type = "~" @@ -251,8 +250,8 @@ function compare_estimates( if type == :↔ type = "~~" elseif type == :→ - if (from ∈ partable.variables[:latent_vars]) & - (to ∈ partable.variables[:observed_vars]) + if (from ∈ partable.latent_vars) && + (to ∈ partable.observed_vars) type = "=~" else type = "~" diff --git a/test/unit_tests/data_input_formats.jl b/test/unit_tests/data_input_formats.jl index 485cf82d2..7a048b280 100644 --- a/test/unit_tests/data_input_formats.jl +++ b/test/unit_tests/data_input_formats.jl @@ -2,9 +2,10 @@ using StructuralEquationModels, Test, Statistics using StructuralEquationModels: obs_cov, obs_mean, get_data ### model specification -------------------------------------------------------------------- -spec = ParameterTable(nothing) -spec.variables[:observed_vars] = [:x1, :x2, :x3, :y1, :y2, :y3, :y4, :y5, :y6, :y7, :y8] -spec.variables[:latent_vars] = [:ind60, :dem60, :dem65] +spec = ParameterTable( + observed_vars = [:x1, :x2, :x3, :y1, :y2, :y3, :y4, :y5, :y6, :y7, :y8], + latent_vars = [:ind60, :dem60, :dem65], +) ### data ----------------------------------------------------------------------------------- From afaff2c9512722342330e6931ed636539e054a14 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sun, 26 May 2024 21:23:55 -0700 Subject: [PATCH 08/56] ParamTable: update StenGraph-based ctor * use graph as a main parameter * simplify rows processing * don't reallocate table.columns Co-authored-by: Maximilian-Stefan-Ernst <34346372+Maximilian-Stefan-Ernst@users.noreply.github.com> --- src/frontend/specification/StenoGraphs.jl | 124 +++++++++--------- test/examples/multigroup/multigroup.jl | 8 +- .../political_democracy.jl | 7 +- 3 files changed, 66 insertions(+), 73 deletions(-) diff --git a/src/frontend/specification/StenoGraphs.jl b/src/frontend/specification/StenoGraphs.jl index bebbdeb2e..1d3332cb9 100644 --- a/src/frontend/specification/StenoGraphs.jl +++ b/src/frontend/specification/StenoGraphs.jl @@ -4,6 +4,9 @@ ### Define Modifiers ############################################################################################ +#FIXME: remove when StenoGraphs.jl will provide AbstractStenoGraph +const AbstractStenoGraph = AbstractArray{T, 1} where {T <: StenoGraphs.AbstractEdge} + # fixed parameter values struct Fixed{N} <: EdgeModifier value::N @@ -28,59 +31,59 @@ label(args...) = Label(args) ### constructor for parameter table from graph ############################################################################################ -function ParameterTable(; graph, observed_vars, latent_vars, g = 1, parname = :θ) +function ParameterTable( + graph::AbstractStenoGraph; + observed_vars, + latent_vars, + group::Integer = 1, + param_prefix = :θ, +) graph = unique(graph) n = length(graph) - from = Vector{Symbol}(undef, n) - parameter_type = Vector{Symbol}(undef, n) - to = Vector{Symbol}(undef, n) - free = ones(Bool, n) - value_fixed = zeros(n) - start = zeros(n) - estimate = zeros(n) - params = Vector{Symbol}(undef, n) - params .= Symbol("") - # group = Vector{Symbol}(undef, n) - # start_partable = zeros(Bool, n) - sorted_vars = Vector{Symbol}() + partable = ParameterTable(latent_vars = latent_vars, observed_vars = observed_vars) + from = resize!(partable.columns[:from], n) + parameter_type = resize!(partable.columns[:parameter_type], n) + to = resize!(partable.columns[:to], n) + free = fill!(resize!(partable.columns[:free], n), true) + value_fixed = fill!(resize!(partable.columns[:value_fixed], n), NaN) + start = fill!(resize!(partable.columns[:start], n), NaN) + param_refs = fill!(resize!(partable.columns[:param], n), Symbol("")) + # group = Vector{Symbol}(undef, n) for (i, element) in enumerate(graph) - if element isa DirectedEdge - from[i] = element.src.node - to[i] = element.dst.node + edge = element isa ModifiedEdge ? element.edge : element + from[i] = edge.src.node + to[i] = edge.dst.node + if edge isa DirectedEdge parameter_type[i] = :→ - elseif element isa UndirectedEdge - from[i] = element.src.node - to[i] = element.dst.node + elseif edge isa UndirectedEdge parameter_type[i] = :↔ - elseif element isa ModifiedEdge - if element.edge isa DirectedEdge - from[i] = element.edge.src.node - to[i] = element.edge.dst.node - parameter_type[i] = :→ - elseif element.edge isa UndirectedEdge - from[i] = element.edge.src.node - to[i] = element.edge.dst.node - parameter_type[i] = :↔ - end + else + throw( + ArgumentError( + "The graph contains an unsupported edge of type $(typeof(edge)).", + ), + ) + end + if element isa ModifiedEdge for modifier in values(element.modifiers) + modval = modifier.value[group] if modifier isa Fixed - if modifier.value[g] == :NaN + if modval == :NaN free[i] = true value_fixed[i] = 0.0 else free[i] = false - value_fixed[i] = modifier.value[g] + value_fixed[i] = modval end elseif modifier isa Start - start_partable[i] = modifier.value[g] == :NaN - start[i] = modifier.value[g] + start[i] = modval elseif modifier isa Label - if modifier.value[g] == :NaN + if modval == :NaN throw(DomainError(NaN, "NaN is not allowed as a parameter label.")) end - params[i] = modifier.value[g] + param_refs[i] = modval end end end @@ -88,41 +91,32 @@ function ParameterTable(; graph, observed_vars, latent_vars, g = 1, parname = : # make identifiers for parameters that are not labeled current_id = 1 - for i in 1:length(params) - if (params[i] == Symbol("")) & free[i] - params[i] = Symbol(parname, :_, current_id) - current_id += 1 - elseif (params[i] == Symbol("")) & !free[i] - params[i] = :const - elseif (params[i] != Symbol("")) & !free[i] - @warn "You labeled a constant. Please check if the labels of your graph are correct." + for i in eachindex(param_refs) + if param_refs[i] == Symbol("") + if free[i] + param_refs[i] = Symbol(param_prefix, :_, current_id) + current_id += 1 + else + param_refs[i] = :const + end + elseif !free[i] + @warn "You labeled a constant ($(param_refs[i])=$(value_fixed[i])). Please check if the labels of your graph are correct." end end - return StructuralEquationModels.ParameterTable( - Dict( - :from => from, - :parameter_type => parameter_type, - :to => to, - :free => free, - :value_fixed => value_fixed, - :start => start, - :estimate => estimate, - :param => params, - ), - Dict( - :latent_vars => latent_vars, - :observed_vars => observed_vars, - :sorted_vars => sorted_vars, - ), - ) + return partable end ############################################################################################ ### constructor for EnsembleParameterTable from graph ############################################################################################ -function EnsembleParameterTable(; graph, observed_vars, latent_vars, groups) +function EnsembleParameterTable( + graph::AbstractStenoGraph; + observed_vars, + latent_vars, + groups +) graph = unique(graph) partable = EnsembleParameterTable(nothing) @@ -130,12 +124,12 @@ function EnsembleParameterTable(; graph, observed_vars, latent_vars, groups) for (i, group) in enumerate(groups) push!( partable.tables, - Symbol(group) => ParameterTable(; - graph = graph, + Symbol(group) => ParameterTable( + graph; observed_vars = observed_vars, latent_vars = latent_vars, - g = i, - parname = Symbol(:g, i), + group = i, + param_prefix = Symbol(:g, i), ), ) end diff --git a/test/examples/multigroup/multigroup.jl b/test/examples/multigroup/multigroup.jl index 818f9afdc..0a648f2dc 100644 --- a/test/examples/multigroup/multigroup.jl +++ b/test/examples/multigroup/multigroup.jl @@ -111,8 +111,8 @@ graph = @StenoGraph begin _(latent_vars) ⇔ _(latent_vars) end -partable = EnsembleParameterTable(; - graph = graph, +partable = EnsembleParameterTable( + graph; observed_vars = observed_vars, latent_vars = latent_vars, groups = [:Pasteur, :Grant_White], @@ -140,8 +140,8 @@ graph = @StenoGraph begin Symbol("1") → _(observed_vars) end -partable_miss = EnsembleParameterTable(; - graph = graph, +partable_miss = EnsembleParameterTable( + graph; observed_vars = observed_vars, latent_vars = latent_vars, groups = [:Pasteur, :Grant_White], diff --git a/test/examples/political_democracy/political_democracy.jl b/test/examples/political_democracy/political_democracy.jl index 86b7e89bc..9085531b0 100644 --- a/test/examples/political_democracy/political_democracy.jl +++ b/test/examples/political_democracy/political_democracy.jl @@ -136,6 +136,7 @@ semoptimizer = SemOptimizerOptim @testset "RAMMatrices | constructor | Optim" begin include("constructor.jl") end + semoptimizer = SemOptimizerNLopt @testset "RAMMatrices | constructor | NLopt" begin include("constructor.jl") @@ -212,8 +213,7 @@ graph = @StenoGraph begin y8 ↔ y4 + y6 end -spec = - ParameterTable(latent_vars = latent_vars, observed_vars = observed_vars, graph = graph) +spec = ParameterTable(graph, latent_vars = latent_vars, observed_vars = observed_vars) sort!(spec) @@ -244,8 +244,7 @@ graph = @StenoGraph begin Symbol("1") → fixed(0) * ind60 end -spec_mean = - ParameterTable(latent_vars = latent_vars, observed_vars = observed_vars, graph = graph) +spec_mean = ParameterTable(graph, latent_vars = latent_vars, observed_vars = observed_vars) sort!(spec_mean) From 6308cdda46914428bcf2e67454bb560d4d79bec4 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sun, 26 May 2024 21:28:45 -0700 Subject: [PATCH 09/56] rename Base.sort() to sort_vars() because the ParTable contains rows and columns, it is not clear, what sort() actually sorts. Co-authored-by: Maximilian-Stefan-Ernst <34346372+Maximilian-Stefan-Ernst@users.noreply.github.com> --- src/StructuralEquationModels.jl | 2 ++ .../specification/EnsembleParameterTable.jl | 18 ++++-------- src/frontend/specification/ParameterTable.jl | 29 +++++++++++++++---- test/examples/multigroup/build_models.jl | 2 +- .../political_democracy.jl | 4 +-- test/unit_tests/sorting.jl | 4 +-- 6 files changed, 36 insertions(+), 23 deletions(-) diff --git a/src/StructuralEquationModels.jl b/src/StructuralEquationModels.jl index 3319f049b..5f8bd070f 100644 --- a/src/StructuralEquationModels.jl +++ b/src/StructuralEquationModels.jl @@ -154,6 +154,8 @@ export AbstractSem, Label, label, params_to_indices, + sort_vars!, + sort_vars, RAMMatrices, param_indices, fit_measures, diff --git a/src/frontend/specification/EnsembleParameterTable.jl b/src/frontend/specification/EnsembleParameterTable.jl index 24a9a295a..4653022cc 100644 --- a/src/frontend/specification/EnsembleParameterTable.jl +++ b/src/frontend/specification/EnsembleParameterTable.jl @@ -67,23 +67,17 @@ end ### Additional Methods ############################################################################################ -# Sorting ---------------------------------------------------------------------------------- +# Variables Sorting ------------------------------------------------------------------------ -# Sorting ---------------------------------------------------------------------------------- - -function sort!(ensemble_partable::EnsembleParameterTable) - for partable in values(ensemble_partable.tables) - sort!(partable) +function sort_vars!(partables::EnsembleParameterTable) + for partable in values(partables.tables) + sort_vars!(partable) end - return ensemble_partable + return partables end -function sort(partable::EnsembleParameterTable) - new_partable = deepcopy(partable) - sort!(new_partable) - return new_partable -end +sort_vars(partables::EnsembleParameterTable) = sort_vars!(deepcopy(partables)) # add a row -------------------------------------------------------------------------------- diff --git a/src/frontend/specification/ParameterTable.jl b/src/frontend/specification/ParameterTable.jl index 49ea4664b..d82b22f1d 100644 --- a/src/frontend/specification/ParameterTable.jl +++ b/src/frontend/specification/ParameterTable.jl @@ -112,7 +112,18 @@ end Base.showerror(io::IO, e::CyclicModelError) = print(io, e.msg) -function Base.sort!(partable::ParameterTable) +""" + sort_vars!(partable::ParameterTable) + sort_vars!(partables::EnsembleParameterTable) + +Sort variables in `partable` so that all independent variables are +before the dependent variables and store it in `partable.sorted_vars`. + +If the relations between the variables are acyclic, sorting will +make the resulting `A` matrix in the *RAM* model lower triangular +and allow faster calculations. +""" +function sort_vars!(partable::ParameterTable) vars = [ partable.latent_vars partable.observed_vars @@ -151,11 +162,17 @@ function Base.sort!(partable::ParameterTable) return partable end -function Base.sort(partable::ParameterTable) - new_partable = deepcopy(partable) - sort!(new_partable) - return new_partable -end +""" + sort_vars(partable::ParameterTable) + sort_vars(partables::EnsembleParameterTable) + +Sort variables in `partable` so that all independent variables are +before the dependent variables, and return a copy of `partable` +where the sorted variables are in `partable.sorted_vars`. + +See [sort_vars!](@ref) for in-place version. +""" +sort_vars(partable::ParameterTable) = sort_vars!(deepcopy(partable)) # add a row -------------------------------------------------------------------------------- diff --git a/test/examples/multigroup/build_models.jl b/test/examples/multigroup/build_models.jl index 8913860c8..e26facc5e 100644 --- a/test/examples/multigroup/build_models.jl +++ b/test/examples/multigroup/build_models.jl @@ -49,7 +49,7 @@ end # ML estimation - sorted ############################################################################################ -partable_s = sort(partable) +partable_s = sort_vars(partable) specification_s = RAMMatrices(partable_s) diff --git a/test/examples/political_democracy/political_democracy.jl b/test/examples/political_democracy/political_democracy.jl index 9085531b0..d7fbb8f2c 100644 --- a/test/examples/political_democracy/political_democracy.jl +++ b/test/examples/political_democracy/political_democracy.jl @@ -215,7 +215,7 @@ end spec = ParameterTable(graph, latent_vars = latent_vars, observed_vars = observed_vars) -sort!(spec) +sort_vars!(spec) partable = spec @@ -246,7 +246,7 @@ end spec_mean = ParameterTable(graph, latent_vars = latent_vars, observed_vars = observed_vars) -sort!(spec_mean) +sort_vars!(spec_mean) partable_mean = spec_mean diff --git a/test/unit_tests/sorting.jl b/test/unit_tests/sorting.jl index e573c6d22..5ca890c51 100644 --- a/test/unit_tests/sorting.jl +++ b/test/unit_tests/sorting.jl @@ -1,8 +1,8 @@ ############################################################################ -### test sorting +### test variables sorting ############################################################################ -sort!(partable) +sort_vars!(partable) model_ml_sorted = Sem(specification = partable, data = dat) From 92730f8c616c165b16a52261b8a3e017729128a9 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sun, 10 Mar 2024 21:51:05 -0700 Subject: [PATCH 10/56] don't import == --- src/frontend/specification/RAMMatrices.jl | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/src/frontend/specification/RAMMatrices.jl b/src/frontend/specification/RAMMatrices.jl index 45bdfe57b..f3c8f11fc 100644 --- a/src/frontend/specification/RAMMatrices.jl +++ b/src/frontend/specification/RAMMatrices.jl @@ -49,9 +49,7 @@ struct RAMConstant value::Any end -import Base.== - -function ==(c1::RAMConstant, c2::RAMConstant) +function Base.:(==)(c1::RAMConstant, c2::RAMConstant) res = ((c1.matrix == c2.matrix) && (c1.index == c2.index) && (c1.value == c2.value)) return res end @@ -410,7 +408,7 @@ function push_partable_rows!(partable, position_names, par, i, A_ind, S_ind, M_i return nothing end -function ==(mat1::RAMMatrices, mat2::RAMMatrices) +function Base.:(==)(mat1::RAMMatrices, mat2::RAMMatrices) res = ( (mat1.A_ind == mat2.A_ind) && (mat1.S_ind == mat2.S_ind) && From 83e782823685b3defa4b3b01c167f8e565e47a67 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sun, 10 Mar 2024 21:51:51 -0700 Subject: [PATCH 11/56] don't import push!() --- src/frontend/specification/EnsembleParameterTable.jl | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/src/frontend/specification/EnsembleParameterTable.jl b/src/frontend/specification/EnsembleParameterTable.jl index 4653022cc..1ce7a0d59 100644 --- a/src/frontend/specification/EnsembleParameterTable.jl +++ b/src/frontend/specification/EnsembleParameterTable.jl @@ -82,13 +82,11 @@ sort_vars(partables::EnsembleParameterTable) = sort_vars!(deepcopy(partables)) # add a row -------------------------------------------------------------------------------- # do we really need this? -import Base.push! - -function push!(partable::EnsembleParameterTable, d::AbstractDict, group) +function Base.push!(partable::EnsembleParameterTable, d::AbstractDict, group) push!(partable.tables[group], d) end -push!(partable::EnsembleParameterTable, d::Nothing, group) = nothing +Base.push!(partable::EnsembleParameterTable, d::Nothing, group) = nothing Base.getindex(partable::EnsembleParameterTable, group) = partable.tables[group] From a2eed98e66235429d49462dc5bea27f46ea41c9a Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Wed, 1 May 2024 22:42:28 -0700 Subject: [PATCH 12/56] don't import DataFrame --- src/StructuralEquationModels.jl | 1 - src/frontend/specification/ParameterTable.jl | 8 +++++--- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/src/StructuralEquationModels.jl b/src/StructuralEquationModels.jl index 5f8bd070f..be9eccc07 100644 --- a/src/StructuralEquationModels.jl +++ b/src/StructuralEquationModels.jl @@ -16,7 +16,6 @@ using LinearAlgebra, DelimitedFiles, DataFrames -import DataFrames: DataFrame export StenoGraphs, @StenoGraph, meld const SEM = StructuralEquationModels diff --git a/src/frontend/specification/ParameterTable.jl b/src/frontend/specification/ParameterTable.jl index d82b22f1d..c4ebefef3 100644 --- a/src/frontend/specification/ParameterTable.jl +++ b/src/frontend/specification/ParameterTable.jl @@ -46,12 +46,14 @@ function Dict(partable::ParameterTable) return partable.columns end -function DataFrame(partable::ParameterTable; columns = nothing) +function DataFrames.DataFrame( + partable::ParameterTable; + columns::Union{AbstractVector{Symbol}, Nothing} = nothing, +) if isnothing(columns) columns = keys(partable.columns) end - out = DataFrame([key => partable.columns[key] for key in columns]) - return DataFrame(out) + return DataFrame([col => partable.columns[col] for col in columns]) end ############################################################################################ From 458f1cfd397d60f32822532655c859c27383da7a Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sat, 23 Mar 2024 14:50:28 -0700 Subject: [PATCH 13/56] remove no-op push!() --- src/frontend/specification/EnsembleParameterTable.jl | 2 -- src/frontend/specification/ParameterTable.jl | 2 -- 2 files changed, 4 deletions(-) diff --git a/src/frontend/specification/EnsembleParameterTable.jl b/src/frontend/specification/EnsembleParameterTable.jl index 1ce7a0d59..672b9c25b 100644 --- a/src/frontend/specification/EnsembleParameterTable.jl +++ b/src/frontend/specification/EnsembleParameterTable.jl @@ -86,8 +86,6 @@ function Base.push!(partable::EnsembleParameterTable, d::AbstractDict, group) push!(partable.tables[group], d) end -Base.push!(partable::EnsembleParameterTable, d::Nothing, group) = nothing - Base.getindex(partable::EnsembleParameterTable, group) = partable.tables[group] ############################################################################################ diff --git a/src/frontend/specification/ParameterTable.jl b/src/frontend/specification/ParameterTable.jl index c4ebefef3..41f02da9e 100644 --- a/src/frontend/specification/ParameterTable.jl +++ b/src/frontend/specification/ParameterTable.jl @@ -184,8 +184,6 @@ function Base.push!(partable::ParameterTable, d::AbstractDict{Symbol}) end end -Base.push!(partable::ParameterTable, d::Nothing) = nothing - ############################################################################################ ### Update Partable from Fitted Model ############################################################################################ From ff2140a9442bc31d84851f86d93512b8b5f0b8e7 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Mon, 27 May 2024 14:30:12 -0700 Subject: [PATCH 14/56] ParTable ctor: simplify rows code * use named tuples * reduce code duplication * use colnames vector instead of position_names Dict --- src/frontend/specification/ParameterTable.jl | 2 +- src/frontend/specification/RAMMatrices.jl | 168 +++++++++---------- 2 files changed, 77 insertions(+), 93 deletions(-) diff --git a/src/frontend/specification/ParameterTable.jl b/src/frontend/specification/ParameterTable.jl index 41f02da9e..f73f80b77 100644 --- a/src/frontend/specification/ParameterTable.jl +++ b/src/frontend/specification/ParameterTable.jl @@ -178,7 +178,7 @@ sort_vars(partable::ParameterTable) = sort_vars!(deepcopy(partable)) # add a row -------------------------------------------------------------------------------- -function Base.push!(partable::ParameterTable, d::AbstractDict{Symbol}) +function Base.push!(partable::ParameterTable, d::Union{AbstractDict{Symbol}, NamedTuple}) for (key, val) in pairs(d) push!(partable.columns[key], val) end diff --git a/src/frontend/specification/RAMMatrices.jl b/src/frontend/specification/RAMMatrices.jl index f3c8f11fc..ba5af9243 100644 --- a/src/frontend/specification/RAMMatrices.jl +++ b/src/frontend/specification/RAMMatrices.jl @@ -214,18 +214,16 @@ function ParameterTable(ram_matrices::RAMMatrices) latent_vars = colnames[setdiff(eachindex(colnames), ram_matrices.F_ind)], ) - position_names = Dict{Int64, Symbol}(1:length(colnames) .=> colnames) - # constants for c in ram_matrices.constants - push!(partable, get_partable_row(c, position_names)) + push!(partable, partable_row(c, colnames)) end # parameters for (i, par) in enumerate(ram_matrices.params) - push_partable_rows!( + append_partable_rows!( partable, - position_names, + colnames, par, i, ram_matrices.A_ind, @@ -296,112 +294,98 @@ function get_par_npar_indices(partable::EnsembleParameterTable) return params, n_par, par_positions end -function get_partable_row(c::RAMConstant, position_names) - # variable names - from = position_names[c.index[2]] - to = position_names[c.index[1]] - # parameter type - if c.matrix == :A - parameter_type = :→ - elseif c.matrix == :S - parameter_type = :↔ - elseif c.matrix == :M - parameter_type = :→ - end - free = false - value_fixed = c.value - start = 0.0 - estimate = 0.0 - - return Dict( - :from => from, - :parameter_type => parameter_type, - :to => to, - :free => free, - :value_fixed => value_fixed, - :start => start, - :estimate => estimate, - :param => :const, - ) -end - -function cartesian_is_known(index, known_indices) - known = false - for k_in in known_indices - if (index == k_in) | ((index[1] == k_in[2]) & (index[2] == k_in[1])) - known = true - end +function matrix_to_parameter_type(matrix::Symbol) + if matrix == :A + return :→ + elseif matrix == :S + return :↔ + elseif matrix == :M + return :→ + else + throw( + ArgumentError( + "Unsupported matrix $matrix, supported matrices are :A, :S and :M", + ), + ) end - return known end -cartesian_is_known(index, known_indices::Nothing) = false - -function get_partable_row(param, position_names, index, matrix, n_nod, known_indices) +partable_row(c::RAMConstant, varnames::AbstractVector{Symbol}) = ( + from = varnames[c.index[2]], + parameter_type = matrix_to_parameter_type(c.matrix), + to = varnames[c.index[1]], + free = false, + value_fixed = c.value, + start = 0.0, + estimate = 0.0, + param = :const, +) + +function partable_row( + par::Symbol, + varnames::AbstractVector{Symbol}, + index::Integer, + matrix::Symbol, + n_nod::Integer, +) # variable names if matrix == :M from = Symbol("1") - to = position_names[index] + to = varnames[index] else - index = linear2cartesian(index, (n_nod, n_nod)) - - if (matrix == :S) & (cartesian_is_known(index, known_indices)) - return nothing - elseif matrix == :S - push!(known_indices, index) - end + cart_index = linear2cartesian(index, (n_nod, n_nod)) - from = position_names[index[2]] - to = position_names[index[1]] - end - - # parameter type - if matrix == :A - parameter_type = :→ - elseif matrix == :S - parameter_type = :↔ - elseif matrix == :M - parameter_type = :→ + from = varnames[cart_index[2]] + to = varnames[cart_index[1]] end - free = true - value_fixed = 0.0 - start = 0.0 - estimate = 0.0 - - return Dict( - :from => from, - :parameter_type => parameter_type, - :to => to, - :free => free, - :value_fixed => value_fixed, - :start => start, - :estimate => estimate, - :param => param, + return ( + from = from, + parameter_type = matrix_to_parameter_type(matrix), + to = to, + free = true, + value_fixed = 0.0, + start = 0.0, + estimate = 0.0, + param = par, ) end -function push_partable_rows!(partable, position_names, par, i, A_ind, S_ind, M_ind, n_nod) - A_ind = A_ind[i] - S_ind = S_ind[i] - isnothing(M_ind) || (M_ind = M_ind[i]) - - for ind in A_ind - push!(partable, get_partable_row(par, position_names, ind, :A, n_nod, nothing)) +function append_partable_rows!( + partable::ParameterTable, + varnames::AbstractVector{Symbol}, + par::Symbol, + par_index::Integer, + A_ind, + S_ind, + M_ind, + n_nod::Integer, +) + for ind in A_ind[par_index] + push!(partable, partable_row(par, varnames, ind, :A, n_nod)) end - known_indices = Vector{CartesianIndex}() - for ind in S_ind - push!( - partable, - get_partable_row(par, position_names, ind, :S, n_nod, known_indices), - ) + visited_S_indices = Set{Int}() + for ind in S_ind[par_index] + if ind ∉ visited_S_indices + push!(partable, partable_row(par, varnames, ind, :S, n_nod)) + # mark index and its symmetric as visited + push!(visited_S_indices, ind) + cart_index = linear2cartesian(ind, (n_nod, n_nod)) + push!( + visited_S_indices, + cartesian2linear( + CartesianIndex(cart_index[2], cart_index[1]), + (n_nod, n_nod), + ), + ) + end end if !isnothing(M_ind) - for ind in M_ind - push!(partable, get_partable_row(par, position_names, ind, :M, n_nod, nothing)) + for ind in M_ind[par_index] + push!(partable, partable_row(par, varnames, ind, :M, n_nod)) end end From edc8443e39338c32af665a5130ec713c27346558 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Mon, 27 May 2024 00:52:46 -0700 Subject: [PATCH 15/56] ParTable: full support for Iterator iface --- src/frontend/specification/ParameterTable.jl | 15 ++++++++++++++- 1 file changed, 14 insertions(+), 1 deletion(-) diff --git a/src/frontend/specification/ParameterTable.jl b/src/frontend/specification/ParameterTable.jl index f73f80b77..13d6448df 100644 --- a/src/frontend/specification/ParameterTable.jl +++ b/src/frontend/specification/ParameterTable.jl @@ -94,6 +94,14 @@ end ############################################################################################ # Iteration -------------------------------------------------------------------------------- +ParameterTableRow = @NamedTuple begin + from::Symbol + parameter_type::Symbol + to::Symbol + free::Bool + value_fixed::Any + param::Symbol +end Base.getindex(partable::ParameterTable, i::Integer) = ( partable.columns[:from][i], @@ -104,7 +112,12 @@ Base.getindex(partable::ParameterTable, i::Integer) = ( partable.columns[:param][i], ) -Base.length(partable::ParameterTable) = length(first(values(partable.columns))) +Base.length(partable::ParameterTable) = length(partable.columns[:param]) +Base.eachindex(partable::ParameterTable) = Base.OneTo(length(partable)) + +Base.eltype(::Type{<:ParameterTable}) = ParameterTableRow +Base.iterate(partable::ParameterTable, i::Integer = 1) = + i > length(partable) ? nothing : (partable[i], i + 1) # Sorting ---------------------------------------------------------------------------------- From 018b077a73806c39670b1d0b44468e8a3aa8356f Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sat, 9 Mar 2024 15:23:26 -0800 Subject: [PATCH 16/56] RAMConstant: simplify * declare RAMConstant field types * refactor constants collection to avoid code duplication --- src/frontend/specification/RAMMatrices.jl | 110 ++++++++++------------ 1 file changed, 50 insertions(+), 60 deletions(-) diff --git a/src/frontend/specification/RAMMatrices.jl b/src/frontend/specification/RAMMatrices.jl index ba5af9243..8c6314316 100644 --- a/src/frontend/specification/RAMMatrices.jl +++ b/src/frontend/specification/RAMMatrices.jl @@ -1,3 +1,48 @@ +############################################################################################ +### Constants +############################################################################################ + +struct RAMConstant + matrix::Symbol + index::Union{Int, CartesianIndex{2}} + value::Any +end + +function Base.:(==)(c1::RAMConstant, c2::RAMConstant) + res = ((c1.matrix == c2.matrix) && (c1.index == c2.index) && (c1.value == c2.value)) + return res +end + +function append_RAMConstants!( + constants::AbstractVector{RAMConstant}, + mtx_name::Symbol, + mtx::AbstractArray, +) + for (index, val) in pairs(mtx) + if isa(val, Number) && !iszero(val) + push!(constants, RAMConstant(mtx_name, index, val)) + end + end + return constants +end + +function set_RAMConstant!(A, S, M, rc::RAMConstant) + if rc.matrix == :A + A[rc.index] = rc.value + elseif rc.matrix == :S + S[rc.index] = rc.value + S[rc.index[2], rc.index[1]] = rc.value # symmetric + elseif rc.matrix == :M + M[rc.index] = rc.value + end +end + +function set_RAMConstants!(A, S, M, rc_vec::Vector{RAMConstant}) + for rc in rc_vec + set_RAMConstant!(A, S, M, rc) + end +end + ############################################################################################ ### Type ############################################################################################ @@ -13,7 +58,7 @@ struct RAMMatrices <: SemSpecification M_ind::Union{ArrayParamsMap, Nothing} params::Any colnames::Any - constants::Any + constants::Vector{RAMConstant} size_F::Any end @@ -26,7 +71,10 @@ function RAMMatrices(; A, S, F, M = nothing, params, colnames) S_indices = array_params_map(params, S) M_indices = !isnothing(M) ? array_params_map(params, M) : nothing F_indices = findall([any(isone.(col)) for col in eachcol(F)]) - constants = get_RAMConstants(A, S, M) + constants = Vector{RAMConstant}() + append_RAMConstants!(constants, :A, A) + append_RAMConstants!(constants, :S, S) + isnothing(M) || append_RAMConstants!(constants, :M, M) return RAMMatrices( A_indices, S_indices, @@ -39,64 +87,6 @@ function RAMMatrices(; A, S, F, M = nothing, params, colnames) ) end -############################################################################################ -### Constants -############################################################################################ - -struct RAMConstant - matrix::Any - index::Any - value::Any -end - -function Base.:(==)(c1::RAMConstant, c2::RAMConstant) - res = ((c1.matrix == c2.matrix) && (c1.index == c2.index) && (c1.value == c2.value)) - return res -end - -function get_RAMConstants(A, S, M) - constants = Vector{RAMConstant}() - - for index in CartesianIndices(A) - if (A[index] isa Number) && !iszero(A[index]) - push!(constants, RAMConstant(:A, index, A[index])) - end - end - - for index in CartesianIndices(S) - if (S[index] isa Number) && !iszero(S[index]) - push!(constants, RAMConstant(:S, index, S[index])) - end - end - - if !isnothing(M) - for index in CartesianIndices(M) - if (M[index] isa Number) && !iszero(M[index]) - push!(constants, RAMConstant(:M, index, M[index])) - end - end - end - - return constants -end - -function set_RAMConstant!(A, S, M, rc::RAMConstant) - if rc.matrix == :A - A[rc.index] = rc.value - elseif rc.matrix == :S - S[rc.index] = rc.value - S[rc.index[2], rc.index[1]] = rc.value - elseif rc.matrix == :M - M[rc.index] = rc.value - end -end - -function set_RAMConstants!(A, S, M, rc_vec::Vector{RAMConstant}) - for rc in rc_vec - set_RAMConstant!(A, S, M, rc) - end -end - ############################################################################################ ### get RAMMatrices from parameter table ############################################################################################ From 70e6199b8c03949c0d429d7127df9261eb0abb2b Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sat, 9 Mar 2024 15:24:39 -0800 Subject: [PATCH 17/56] RAMMatrices: optimize F_indices init --- src/frontend/specification/RAMMatrices.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/frontend/specification/RAMMatrices.jl b/src/frontend/specification/RAMMatrices.jl index 8c6314316..718e319a1 100644 --- a/src/frontend/specification/RAMMatrices.jl +++ b/src/frontend/specification/RAMMatrices.jl @@ -70,7 +70,7 @@ function RAMMatrices(; A, S, F, M = nothing, params, colnames) A_indices = array_params_map(params, A) S_indices = array_params_map(params, S) M_indices = !isnothing(M) ? array_params_map(params, M) : nothing - F_indices = findall([any(isone.(col)) for col in eachcol(F)]) + F_indices = [i for (i, col) in zip(axes(F, 2), eachcol(F)) if any(isone, col)] constants = Vector{RAMConstant}() append_RAMConstants!(constants, :A, A) append_RAMConstants!(constants, :S, S) From 9bf094e1780b6f63a12f41cd33f248fd1d7a7573 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sat, 9 Mar 2024 15:26:01 -0800 Subject: [PATCH 18/56] RAMMatrices: declare types for all fields --- src/frontend/specification/RAMMatrices.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/frontend/specification/RAMMatrices.jl b/src/frontend/specification/RAMMatrices.jl index 718e319a1..6f674b6bd 100644 --- a/src/frontend/specification/RAMMatrices.jl +++ b/src/frontend/specification/RAMMatrices.jl @@ -56,10 +56,10 @@ struct RAMMatrices <: SemSpecification S_ind::ArrayParamsMap F_ind::Vector{Int} M_ind::Union{ArrayParamsMap, Nothing} - params::Any - colnames::Any + params::Vector{Symbol} + colnames::Union{Vector{Symbol}, Nothing} constants::Vector{RAMConstant} - size_F::Any + size_F::Tuple{Int, Int} end ############################################################################################ From 9a5842ee287064763050363dfce88efe88d49419 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sun, 10 Mar 2024 12:14:33 -0700 Subject: [PATCH 19/56] RAMMatrices: option to keep zero constants --- src/frontend/specification/RAMMatrices.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/frontend/specification/RAMMatrices.jl b/src/frontend/specification/RAMMatrices.jl index 6f674b6bd..e605d432a 100644 --- a/src/frontend/specification/RAMMatrices.jl +++ b/src/frontend/specification/RAMMatrices.jl @@ -16,10 +16,11 @@ end function append_RAMConstants!( constants::AbstractVector{RAMConstant}, mtx_name::Symbol, - mtx::AbstractArray, + mtx::AbstractArray; + skip_zeros::Bool = true, ) for (index, val) in pairs(mtx) - if isa(val, Number) && !iszero(val) + if isa(val, Number) && !(skip_zeros && iszero(val)) push!(constants, RAMConstant(mtx_name, index, val)) end end From 296d827d9d5adf76b20a001f6df0092946d4f6ee Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Wed, 3 Apr 2024 22:49:34 -0700 Subject: [PATCH 20/56] nonunique() helper function --- src/additional_functions/helper.jl | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/src/additional_functions/helper.jl b/src/additional_functions/helper.jl index b96813dc3..3e614e57b 100644 --- a/src/additional_functions/helper.jl +++ b/src/additional_functions/helper.jl @@ -142,3 +142,18 @@ function elimination_matrix(n::Integer) end return L end + +# returns the vector of non-unique values in the order of appearance +# each non-unique values is reported once +function nonunique(values::AbstractVector) + value_counts = Dict{eltype(values), Int}() + res = similar(values, 0) + for v in values + n = get!(value_counts, v, 0) + if n == 1 # second encounter + push!(res, v) + end + value_counts[v] = n + 1 + end + return res +end From dec1f4d533df23821bf6f35bfe69db5bcf74acbd Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Fri, 3 May 2024 07:37:14 -0700 Subject: [PATCH 21/56] add check_vars() and check_params() --- src/StructuralEquationModels.jl | 1 + src/frontend/specification/checks.jl | 42 ++++++++++++++++++++++++++++ 2 files changed, 43 insertions(+) create mode 100644 src/frontend/specification/checks.jl diff --git a/src/StructuralEquationModels.jl b/src/StructuralEquationModels.jl index be9eccc07..8fe2ff90e 100644 --- a/src/StructuralEquationModels.jl +++ b/src/StructuralEquationModels.jl @@ -30,6 +30,7 @@ include("additional_functions/commutation_matrix.jl") # fitted objects include("frontend/fit/SemFit.jl") # specification of models +include("frontend/specification/checks.jl") include("frontend/specification/ParameterTable.jl") include("frontend/specification/EnsembleParameterTable.jl") include("frontend/specification/RAMMatrices.jl") diff --git a/src/frontend/specification/checks.jl b/src/frontend/specification/checks.jl new file mode 100644 index 000000000..5326e535f --- /dev/null +++ b/src/frontend/specification/checks.jl @@ -0,0 +1,42 @@ +# check if params vector correctly matches the parameter references (from the ParameterTable) +function check_params( + params::AbstractVector{Symbol}, + param_refs::Union{AbstractVector{Symbol}, Nothing}, +) + dup_params = nonunique(params) + isempty(dup_params) || + throw(ArgumentError("Duplicate parameters detected: $(join(dup_params, ", "))")) + any(==(:const), params) && + throw(ArgumentError("Parameters constain reserved :const name")) + + if !isnothing(param_refs) + # check if all references parameters are present + all_refs = Set(id for id in param_refs if id != :const) + undecl_params = setdiff(all_refs, params) + if !isempty(undecl_params) + throw( + ArgumentError( + "The following $(length(undecl_params)) parameters present in the table, but are not declared: " * + join(sort!(collect(undecl_params))), + ), + ) + end + end + + return nothing +end + +function check_vars(vars::AbstractVector{Symbol}, nvars::Union{Integer, Nothing}) + isnothing(nvars) || + length(vars) == nvars || + throw( + DimensionMismatch( + "variables length ($(length(vars))) does not match the number of columns in A matrix ($nvars)", + ), + ) + dup_vars = nonunique(vars) + isempty(dup_vars) || + throw(ArgumentError("Duplicate variables detected: $(join(dup_vars, ", "))")) + + return nothing +end From a7a17dfed691b98d3c03f90088e716bcf66070cc Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Fri, 3 May 2024 07:43:25 -0700 Subject: [PATCH 22/56] RAMMatrices ctor: dims and vars checks --- src/frontend/specification/RAMMatrices.jl | 36 ++++++++++++++++++++++- 1 file changed, 35 insertions(+), 1 deletion(-) diff --git a/src/frontend/specification/RAMMatrices.jl b/src/frontend/specification/RAMMatrices.jl index e605d432a..551cf4c60 100644 --- a/src/frontend/specification/RAMMatrices.jl +++ b/src/frontend/specification/RAMMatrices.jl @@ -67,7 +67,41 @@ end ### Constructor ############################################################################################ -function RAMMatrices(; A, S, F, M = nothing, params, colnames) +function RAMMatrices(; + A::AbstractMatrix, + S::AbstractMatrix, + F::AbstractMatrix, + M::Union{AbstractVector, Nothing} = nothing, + params::AbstractVector{Symbol}, + colnames::Union{AbstractVector{Symbol}, Nothing} = nothing, +) + ncols = size(A, 2) + isnothing(colnames) || check_vars(colnames, ncols) + + size(A, 1) == size(A, 2) || throw(DimensionMismatch("A must be a square matrix")) + size(S, 1) == size(S, 2) || throw(DimensionMismatch("S must be a square matrix")) + size(A, 2) == ncols || throw( + DimensionMismatch( + "A should have as many rows and columns as colnames length ($ncols), $(size(A)) found", + ), + ) + size(S, 2) == ncols || throw( + DimensionMismatch( + "S should have as many rows and columns as colnames length ($ncols), $(size(S)) found", + ), + ) + size(F, 2) == ncols || throw( + DimensionMismatch( + "F should have as many columns as colnames length ($ncols), $(size(F, 2)) found", + ), + ) + if !isnothing(M) + length(M) == ncols || throw( + DimensionMismatch( + "M should have as many elements as colnames length ($ncols), $(length(M)) found", + ), + ) + end A_indices = array_params_map(params, A) S_indices = array_params_map(params, S) M_indices = !isnothing(M) ? array_params_map(params, M) : nothing From 25cd574d6495cd0df5246b02c9e72d63766c1889 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Fri, 3 May 2024 07:43:44 -0700 Subject: [PATCH 23/56] RAMMatrices: cleanup params index * simplify parameters() function to return just a vector of params * RAMMatrices ctor: use check_params() --- .../specification/EnsembleParameterTable.jl | 9 +++ src/frontend/specification/ParameterTable.jl | 8 +- src/frontend/specification/RAMMatrices.jl | 76 ++++--------------- 3 files changed, 31 insertions(+), 62 deletions(-) diff --git a/src/frontend/specification/EnsembleParameterTable.jl b/src/frontend/specification/EnsembleParameterTable.jl index 672b9c25b..6d8961523 100644 --- a/src/frontend/specification/EnsembleParameterTable.jl +++ b/src/frontend/specification/EnsembleParameterTable.jl @@ -67,6 +67,15 @@ end ### Additional Methods ############################################################################################ +# get the vector of all parameters in the table +# the position of the parameter is based on its first appearance in the table (and the ensemble) +function params(partable::EnsembleParameterTable) + params = mapreduce(vcat, values(partable.tables)) do tbl + tbl.columns[:param] + end + return filter!(!=(:const), unique!(params)) # exclude constants +end + # Variables Sorting ------------------------------------------------------------------------ function sort_vars!(partables::EnsembleParameterTable) diff --git a/src/frontend/specification/ParameterTable.jl b/src/frontend/specification/ParameterTable.jl index 13d6448df..5efa2905f 100644 --- a/src/frontend/specification/ParameterTable.jl +++ b/src/frontend/specification/ParameterTable.jl @@ -119,8 +119,14 @@ Base.eltype(::Type{<:ParameterTable}) = ParameterTableRow Base.iterate(partable::ParameterTable, i::Integer = 1) = i > length(partable) ? nothing : (partable[i], i + 1) -# Sorting ---------------------------------------------------------------------------------- +# get the vector of all parameters in the table +# the position of the parameter is based on its first appearance in the table (and the ensemble) +params(partable::ParameterTable) = + filter!(!=(:const), unique(partable.columns[:param])) + + +# Sorting ---------------------------------------------------------------------------------- struct CyclicModelError <: Exception msg::AbstractString end diff --git a/src/frontend/specification/RAMMatrices.jl b/src/frontend/specification/RAMMatrices.jl index 551cf4c60..023f46342 100644 --- a/src/frontend/specification/RAMMatrices.jl +++ b/src/frontend/specification/RAMMatrices.jl @@ -102,6 +102,9 @@ function RAMMatrices(; ), ) end + + check_params(params, nothing) + A_indices = array_params_map(params, A) S_indices = array_params_map(params, S) M_indices = !isnothing(M) ? array_params_map(params, M) : nothing @@ -126,13 +129,13 @@ end ### get RAMMatrices from parameter table ############################################################################################ -function RAMMatrices(partable::ParameterTable; par_id = nothing) - if isnothing(par_id) - params, n_par, par_positions = get_par_npar_indices(partable) - else - params, n_par, par_positions = - par_id[:params], par_id[:n_par], par_id[:par_positions] - end +function RAMMatrices( + partable::ParameterTable; + params::Union{AbstractVector{Symbol}, Nothing} = nothing, +) + params = copy(isnothing(params) ? SEM.params(partable) : params) + check_params(params, partable.columns[:param]) + params_index = Dict(param => i for (i, param) in enumerate(params)) n_observed = length(partable.observed_vars) n_latent = length(partable.latent_vars) @@ -156,20 +159,13 @@ function RAMMatrices(partable::ParameterTable; par_id = nothing) # fill Matrices # known_labels = Dict{Symbol, Int64}() - A_ind = Vector{Vector{Int64}}(undef, n_par) - for i in 1:length(A_ind) - A_ind[i] = Vector{Int64}() - end - S_ind = Vector{Vector{Int64}}(undef, n_par) - S_ind .= [Vector{Int64}()] - for i in 1:length(S_ind) - S_ind[i] = Vector{Int64}() - end + A_ind = [Vector{Int64}() for _ in 1:length(params)] + S_ind = [Vector{Int64}() for _ in 1:length(params)] # is there a meanstructure? M_ind = - any(==(Symbol("1")), partable.columns[:from]) ? [Vector{Int64}() for _ in 1:n_par] : - nothing + any(==(Symbol("1")), partable.columns[:from]) ? + [Vector{Int64}() for _ in 1:length(params)] : nothing # handle constants constants = Vector{RAMConstant}() @@ -197,7 +193,7 @@ function RAMMatrices(partable::ParameterTable; par_id = nothing) error("Unsupported parameter type: $(parameter_type)") end else - par_ind = par_positions[param] + par_ind = params_index[param] if (parameter_type == :→) && (from == Symbol("1")) push!(M_ind[par_ind], row_ind) elseif parameter_type == :→ @@ -264,25 +260,6 @@ end Base.convert(::Type{<:ParameterTable}, ram_matrices::RAMMatrices) = ParameterTable(ram_matrices) -############################################################################################ -### get RAMMatrices from EnsembleParameterTable -############################################################################################ - -function RAMMatrices(partable::EnsembleParameterTable) - ram_matrices = Dict{Symbol, RAMMatrices}() - - params, n_par, par_positions = get_par_npar_indices(partable) - par_id = - Dict(:params => params, :n_par => n_par, :par_positions => par_positions) - - for key in keys(partable.tables) - ram_mat = RAMMatrices(partable.tables[key]; par_id = par_id) - push!(ram_matrices, key => ram_mat) - end - - return ram_matrices -end - ############################################################################################ ### Pretty Printing ############################################################################################ @@ -296,29 +273,6 @@ end ### Additional Functions ############################################################################################ -function get_par_npar_indices(partable::ParameterTable) - params = unique(partable.columns[:param]) - filter!(x -> x != :const, params) - n_par = length(params) - par_positions = Dict(params .=> 1:n_par) - return params, n_par, par_positions -end - -function get_par_npar_indices(partable::EnsembleParameterTable) - params = Vector{Symbol}() - for key in keys(partable.tables) - append!(params, partable.tables[key].columns[:param]) - end - params = unique(params) - filter!(x -> x != :const, params) - - n_par = length(params) - - par_positions = Dict(params .=> 1:n_par) - - return params, n_par, par_positions -end - function matrix_to_parameter_type(matrix::Symbol) if matrix == :A return :→ From 98b74d9cb7f1d38ac3a54921e4917338005ab9e3 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sun, 17 Mar 2024 17:19:57 -0700 Subject: [PATCH 24/56] include RAMMatrices before EnsParTable --- src/StructuralEquationModels.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/StructuralEquationModels.jl b/src/StructuralEquationModels.jl index 8fe2ff90e..a8598ebc7 100644 --- a/src/StructuralEquationModels.jl +++ b/src/StructuralEquationModels.jl @@ -32,8 +32,8 @@ include("frontend/fit/SemFit.jl") # specification of models include("frontend/specification/checks.jl") include("frontend/specification/ParameterTable.jl") -include("frontend/specification/EnsembleParameterTable.jl") include("frontend/specification/RAMMatrices.jl") +include("frontend/specification/EnsembleParameterTable.jl") include("frontend/specification/StenoGraphs.jl") include("frontend/fit/summary.jl") # pretty printing From 5cdbe7c8fd0d1ca9d711713a022af040832c50ad Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sun, 17 Mar 2024 00:35:14 -0700 Subject: [PATCH 25/56] fix EnsParTable to Dict{RAMMatrices} convert * this method is not RAMMatrices ctor, it is Dict{K, RAMMatrices} convert * use comprehension to construct dict --- .../specification/EnsembleParameterTable.jl | 13 +++++++++++++ src/frontend/specification/ParameterTable.jl | 4 +--- test/examples/multigroup/build_models.jl | 2 +- test/examples/multigroup/multigroup.jl | 4 ++-- 4 files changed, 17 insertions(+), 6 deletions(-) diff --git a/src/frontend/specification/EnsembleParameterTable.jl b/src/frontend/specification/EnsembleParameterTable.jl index 6d8961523..8192ab6f8 100644 --- a/src/frontend/specification/EnsembleParameterTable.jl +++ b/src/frontend/specification/EnsembleParameterTable.jl @@ -24,6 +24,19 @@ function Base.convert(::Type{Dict}, partable::EnsembleParameterTable) return convert(Dict, partable.tables) end +function Base.convert( + ::Type{Dict{K, RAMMatrices}}, + partables::EnsembleParameterTable; + params::Union{AbstractVector{Symbol}, Nothing} = nothing, +) where {K} + isnothing(params) || (params = SEM.params(partables)) + + return Dict{K, RAMMatrices}( + K(key) => RAMMatrices(partable; params = params) for + (key, partable) in pairs(partables.tables) + ) +end + #= function DataFrame( partable::ParameterTable; columns = nothing) diff --git a/src/frontend/specification/ParameterTable.jl b/src/frontend/specification/ParameterTable.jl index 5efa2905f..5e4fe157f 100644 --- a/src/frontend/specification/ParameterTable.jl +++ b/src/frontend/specification/ParameterTable.jl @@ -40,9 +40,7 @@ end ### Convert to other types ############################################################################################ -import Base.Dict - -function Dict(partable::ParameterTable) +function Base.convert(::Type{Dict}, partable::ParameterTable) return partable.columns end diff --git a/test/examples/multigroup/build_models.jl b/test/examples/multigroup/build_models.jl index e26facc5e..23d429796 100644 --- a/test/examples/multigroup/build_models.jl +++ b/test/examples/multigroup/build_models.jl @@ -51,7 +51,7 @@ end partable_s = sort_vars(partable) -specification_s = RAMMatrices(partable_s) +specification_s = convert(Dict{Symbol, RAMMatrices}, partable_s) specification_g1_s = specification_s[:Pasteur] specification_g2_s = specification_s[:Grant_White] diff --git a/test/examples/multigroup/multigroup.jl b/test/examples/multigroup/multigroup.jl index 0a648f2dc..552a65cfb 100644 --- a/test/examples/multigroup/multigroup.jl +++ b/test/examples/multigroup/multigroup.jl @@ -118,7 +118,7 @@ partable = EnsembleParameterTable( groups = [:Pasteur, :Grant_White], ) -specification = RAMMatrices(partable) +specification = convert(Dict{Symbol, RAMMatrices}, partable) specification_g1 = specification[:Pasteur] specification_g2 = specification[:Grant_White] @@ -147,7 +147,7 @@ partable_miss = EnsembleParameterTable( groups = [:Pasteur, :Grant_White], ) -specification_miss = RAMMatrices(partable_miss) +specification_miss = convert(Dict{Symbol, RAMMatrices}, partable_miss) specification_miss_g1 = specification_miss[:Pasteur] specification_miss_g2 = specification_miss[:Grant_White] From e71573b8ac67a3a7abbf37510201dcd049670a45 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Wed, 20 Mar 2024 16:40:51 -0700 Subject: [PATCH 26/56] DataFrame(EnsParTable) --- .../specification/EnsembleParameterTable.jl | 17 ++++++++++------- 1 file changed, 10 insertions(+), 7 deletions(-) diff --git a/src/frontend/specification/EnsembleParameterTable.jl b/src/frontend/specification/EnsembleParameterTable.jl index 8192ab6f8..4430651e9 100644 --- a/src/frontend/specification/EnsembleParameterTable.jl +++ b/src/frontend/specification/EnsembleParameterTable.jl @@ -37,13 +37,16 @@ function Base.convert( ) end -#= function DataFrame( - partable::ParameterTable; - columns = nothing) - if isnothing(columns) columns = keys(partable.columns) end - out = DataFrame([key => partable.columns[key] for key in columns]) - return DataFrame(out) -end =# +function DataFrames.DataFrame( + partables::EnsembleParameterTable; + columns::Union{AbstractVector{Symbol}, Nothing} = nothing, +) + mapreduce(vcat, pairs(partables.tables)) do (key, partable) + df = DataFrame(partable; columns = columns) + df[!, :group] .= key + return df + end +end ############################################################################################ ### get parameter table from RAMMatrices From d5ac0d19c49aa7aa60510121d5a581b28f41f569 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Wed, 1 May 2024 16:45:38 -0700 Subject: [PATCH 27/56] params() API method * remove n_par.jl * remove identifier.jl --- src/StructuralEquationModels.jl | 6 +- src/additional_functions/identifier.jl | 59 ------------------- src/frontend/fit/SemFit.jl | 3 + src/frontend/fit/fitmeasures/n_par.jl | 20 ------- .../specification/EnsembleParameterTable.jl | 6 +- src/frontend/specification/ParameterTable.jl | 11 ++-- src/frontend/specification/RAMMatrices.jl | 2 + src/imply/RAM/generic.jl | 11 ---- src/imply/RAM/symbolic.jl | 10 +--- src/imply/empty.jl | 14 +---- src/loss/regularization/ridge.jl | 3 +- src/types.jl | 36 +++++++++-- test/unit_tests/specification.jl | 15 ++--- 13 files changed, 56 insertions(+), 140 deletions(-) delete mode 100644 src/additional_functions/identifier.jl delete mode 100644 src/frontend/fit/fitmeasures/n_par.jl diff --git a/src/StructuralEquationModels.jl b/src/StructuralEquationModels.jl index a8598ebc7..109a913e7 100644 --- a/src/StructuralEquationModels.jl +++ b/src/StructuralEquationModels.jl @@ -74,15 +74,12 @@ include("additional_functions/start_val/start_partable.jl") include("additional_functions/start_val/start_simple.jl") include("additional_functions/artifacts.jl") include("additional_functions/simulation.jl") -# identifier -include("additional_functions/identifier.jl") # fit measures include("frontend/fit/fitmeasures/AIC.jl") include("frontend/fit/fitmeasures/BIC.jl") include("frontend/fit/fitmeasures/chi2.jl") include("frontend/fit/fitmeasures/df.jl") include("frontend/fit/fitmeasures/minus2ll.jl") -include("frontend/fit/fitmeasures/n_par.jl") include("frontend/fit/fitmeasures/n_obs.jl") include("frontend/fit/fitmeasures/p.jl") include("frontend/fit/fitmeasures/RMSEA.jl") @@ -153,11 +150,10 @@ export AbstractSem, start, Label, label, - params_to_indices, sort_vars!, sort_vars, RAMMatrices, - param_indices, + params, fit_measures, AIC, BIC, diff --git a/src/additional_functions/identifier.jl b/src/additional_functions/identifier.jl deleted file mode 100644 index 1b10357d6..000000000 --- a/src/additional_functions/identifier.jl +++ /dev/null @@ -1,59 +0,0 @@ -############################################################################################ -# get a map from parameters to their indices -############################################################################################ - -param_indices(sem_fit::SemFit) = param_indices(sem_fit.model) -param_indices(model::AbstractSemSingle) = param_indices(model.imply) -param_indices(model::SemEnsemble) = model.param_indices - -############################################################################################ -# construct a map from parameters to indices -############################################################################################ - -param_indices(ram_matrices::RAMMatrices) = - Dict(par => i for (i, par) in enumerate(ram_matrices.params)) -function param_indices(partable::ParameterTable) - _, _, param_indices = get_par_npar_indices(partable) - return param_indices -end - -############################################################################################ -# get indices of a Vector of parameter labels -############################################################################################ - -params_to_indices(params, param_indices::Dict{Symbol, Int}) = - [param_indices[par] for par in params] - -params_to_indices( - params, - obj::Union{SemFit, AbstractSemSingle, SemEnsemble, SemImply}, -) = params_to_indices(params, params(obj)) - -function params_to_indices(params, obj::Union{ParameterTable, RAMMatrices}) - @warn "You are trying to find parameter indices from a ParameterTable or RAMMatrices object. \n - If your model contains user-defined types, this may lead to wrong results. \n - To be on the safe side, try to reference parameters by labels or query the indices from - the constructed model (`params_to_indices(params, model)`)." maxlog = 1 - return params_to_indices(params, params(obj)) -end - -############################################################################################ -# documentation -############################################################################################ -""" - params_to_indices(params, model) - -Returns the indices of `params`. - -# Arguments -- `params::Vector{Symbol}`: parameter labels -- `model`: either a SEM or a fitted SEM - -# Examples -```julia -parameter_indices = params_to_indices([:λ₁, λ₂], my_fitted_sem) - -values = solution(my_fitted_sem)[parameter_indices] -``` -""" -function params_to_indices end diff --git a/src/frontend/fit/SemFit.jl b/src/frontend/fit/SemFit.jl index 97cd9c5a6..19a6d4441 100644 --- a/src/frontend/fit/SemFit.jl +++ b/src/frontend/fit/SemFit.jl @@ -46,6 +46,9 @@ end # additional methods ############################################################################################ +params(fit::SemFit) = params(fit.model) +n_par(fit::SemFit) = n_par(fit.model) + # access fields minimum(sem_fit::SemFit) = sem_fit.minimum solution(sem_fit::SemFit) = sem_fit.solution diff --git a/src/frontend/fit/fitmeasures/n_par.jl b/src/frontend/fit/fitmeasures/n_par.jl deleted file mode 100644 index c8553572b..000000000 --- a/src/frontend/fit/fitmeasures/n_par.jl +++ /dev/null @@ -1,20 +0,0 @@ -############################################################################################ -### get number of parameters -############################################################################################ -""" - n_par(sem_fit::SemFit) - n_par(model::AbstractSemSingle) - n_par(model::SemEnsemble) - n_par(param_indices::Dict) - -Return the number of parameters. -""" -function n_par end - -n_par(fit::SemFit) = n_par(fit.model) - -n_par(model::AbstractSemSingle) = n_par(model.imply) - -n_par(model::SemEnsemble) = n_par(model.param_indices) - -n_par(param_indices::Dict) = length(param_indices) diff --git a/src/frontend/specification/EnsembleParameterTable.jl b/src/frontend/specification/EnsembleParameterTable.jl index 4430651e9..37d7ef15c 100644 --- a/src/frontend/specification/EnsembleParameterTable.jl +++ b/src/frontend/specification/EnsembleParameterTable.jl @@ -120,12 +120,12 @@ Base.getindex(partable::EnsembleParameterTable, group) = partable.tables[group] # update generic --------------------------------------------------------------------------- function update_partable!( partable::EnsembleParameterTable, - param_indices::AbstractDict, - vec, + params::AbstractVector{Symbol}, + values::AbstractVector, column, ) for k in keys(partable.tables) - update_partable!(partable.tables[k], param_indices, vec, column) + update_partable!(partable.tables[k], params, values, column) end return partable end diff --git a/src/frontend/specification/ParameterTable.jl b/src/frontend/specification/ParameterTable.jl index 5e4fe157f..60cc93ec9 100644 --- a/src/frontend/specification/ParameterTable.jl +++ b/src/frontend/specification/ParameterTable.jl @@ -209,14 +209,15 @@ end function update_partable!( partable::ParameterTable, - param_indices::AbstractDict, + params::AbstractVector{Symbol}, values::AbstractVector, column, ) new_col = Vector{eltype(vec)}(undef, length(partable)) + params_index = Dict(param => i for (i, param) in enumerate(params)) for (i, param) in enumerate(partable.columns[:param]) if !(param == :const) - new_col[i] = values[param_indices[param]] + new_col[i] = values[params_index[param]] elseif param == :const new_col[i] = zero(eltype(values)) end @@ -233,8 +234,8 @@ Write `vec` to `column` of `partable`. # Arguments - `vec::Vector`: has to be in the same order as the `model` parameters """ -update_partable!(partable::AbstractParameterTable, sem_fit::SemFit, vec, column) = - update_partable!(partable, param_indices(sem_fit), vec, column) +update_partable!(partable::AbstractParameterTable, sem_fit::SemFit, values, column) = + update_partable!(partable, params(sem_fit), values, column) # update estimates ------------------------------------------------------------------------- """ @@ -271,7 +272,7 @@ function update_start!( if !(start_val isa Vector) start_val = start_val(model; kwargs...) end - return update_partable!(partable, param_indices(model), start_val, :start) + return update_partable!(partable, params(model), start_val, :start) end # update partable standard errors ---------------------------------------------------------- diff --git a/src/frontend/specification/RAMMatrices.jl b/src/frontend/specification/RAMMatrices.jl index 023f46342..5d3abe295 100644 --- a/src/frontend/specification/RAMMatrices.jl +++ b/src/frontend/specification/RAMMatrices.jl @@ -63,6 +63,8 @@ struct RAMMatrices <: SemSpecification size_F::Tuple{Int, Int} end +n_par(ram::RAMMatrices) = length(ram.A_ind) + ############################################################################################ ### Constructor ############################################################################################ diff --git a/src/imply/RAM/generic.jl b/src/imply/RAM/generic.jl index d43c8378e..0988d8e99 100644 --- a/src/imply/RAM/generic.jl +++ b/src/imply/RAM/generic.jl @@ -72,7 +72,6 @@ mutable struct RAM{ A4, A5, A6, - V, V2, I1, I2, @@ -85,7 +84,6 @@ mutable struct RAM{ S2, S3, B, - D, } <: SemImply Σ::A1 A::A2 @@ -94,7 +92,6 @@ mutable struct RAM{ μ::A5 M::A6 - n_par::V ram_matrices::V2 has_meanstructure::B @@ -110,8 +107,6 @@ mutable struct RAM{ ∇A::S1 ∇S::S2 ∇M::S3 - - param_indices::D end using StructuralEquationModels @@ -128,7 +123,6 @@ function RAM(; kwargs..., ) ram_matrices = convert(RAMMatrices, specification) - param_indices = SEM.param_indices(ram_matrices) # get dimensions of the model n_par = length(ram_matrices.params) @@ -184,7 +178,6 @@ function RAM(; F, μ, M_pre, - n_par, ram_matrices, has_meanstructure, A_indices, @@ -197,7 +190,6 @@ function RAM(; ∇A, ∇S, ∇M, - param_indices, ) end @@ -280,9 +272,6 @@ objective_gradient_hessian!(imply::RAM, par, model::AbstractSemSingle, has_means ### Recommended methods ############################################################################################ -param_indices(imply::RAM) = imply.param_indices -n_par(imply::RAM) = imply.n_par - function update_observed(imply::RAM, observed::SemObserved; kwargs...) if n_man(observed) == size(imply.Σ, 1) return imply diff --git a/src/imply/RAM/symbolic.jl b/src/imply/RAM/symbolic.jl index 0fe9c29bb..db5997d2b 100644 --- a/src/imply/RAM/symbolic.jl +++ b/src/imply/RAM/symbolic.jl @@ -62,7 +62,7 @@ and for models with a meanstructure, the model implied means are computed as \mu = F(I-A)^{-1}M ``` """ -struct RAMSymbolic{F1, F2, F3, A1, A2, A3, S1, S2, S3, V, V2, F4, A4, F5, A5, D1, B} <: +struct RAMSymbolic{F1, F2, F3, A1, A2, A3, S1, S2, S3, V2, F4, A4, F5, A5, B} <: SemImplySymbolic Σ_function::F1 ∇Σ_function::F2 @@ -73,13 +73,11 @@ struct RAMSymbolic{F1, F2, F3, A1, A2, A3, S1, S2, S3, V, V2, F4, A4, F5, A5, D1 Σ_symbolic::S1 ∇Σ_symbolic::S2 ∇²Σ_symbolic::S3 - n_par::V ram_matrices::V2 μ_function::F4 μ::A4 ∇μ_function::F5 ∇μ::A5 - param_indices::D1 has_meanstructure::B end @@ -98,7 +96,6 @@ function RAMSymbolic(; kwargs..., ) ram_matrices = convert(RAMMatrices, specification) - param_indices = SEM.param_indices(ram_matrices) n_par = length(ram_matrices.params) n_var, n_nod = ram_matrices.size_F @@ -195,13 +192,11 @@ function RAMSymbolic(; Σ_symbolic, ∇Σ_symbolic, ∇²Σ_symbolic, - n_par, ram_matrices, μ_function, μ, ∇μ_function, ∇μ, - param_indices, has_meanstructure, ) end @@ -240,9 +235,6 @@ objective_gradient_hessian!(imply::RAMSymbolic, par, model) = gradient!(imply, p ### Recommended methods ############################################################################################ -param_indices(imply::RAMSymbolic) = imply.param_indices -n_par(imply::RAMSymbolic) = imply.n_par - function update_observed(imply::RAMSymbolic, observed::SemObserved; kwargs...) if n_man(observed) == size(imply.Σ, 1) return imply diff --git a/src/imply/empty.jl b/src/imply/empty.jl index 1d0ea69ff..56297ea06 100644 --- a/src/imply/empty.jl +++ b/src/imply/empty.jl @@ -25,9 +25,8 @@ model per group and an additional model with `ImplyEmpty` and `SemRidge` for the ## Implementation Subtype of `SemImply`. """ -struct ImplyEmpty{V, V2} <: SemImply - param_indices::V2 - n_par::V +struct ImplyEmpty{V2} <: SemImply + ram_matrices::V2 end ############################################################################################ @@ -35,11 +34,7 @@ end ############################################################################################ function ImplyEmpty(; specification, kwargs...) - ram_matrices = RAMMatrices(specification) - - n_par = length(ram_matrices.params) - - return ImplyEmpty(param_indices(ram_matrices), n_par) + return ImplyEmpty(convert(RAMMatrices, specification)) end ############################################################################################ @@ -54,7 +49,4 @@ hessian!(imply::ImplyEmpty, par, model) = nothing ### Recommended methods ############################################################################################ -param_indices(imply::ImplyEmpty) = imply.param_indices -n_par(imply::ImplyEmpty) = imply.n_par - update_observed(imply::ImplyEmpty, observed::SemObserved; kwargs...) = imply diff --git a/src/loss/regularization/ridge.jl b/src/loss/regularization/ridge.jl index 0d9d10b4b..2d098d550 100644 --- a/src/loss/regularization/ridge.jl +++ b/src/loss/regularization/ridge.jl @@ -58,7 +58,8 @@ function SemRidge(; ), ) else - which_ridge = params_to_indices(which_ridge, imply) + par2ind = Dict(par => ind for (ind, par) in enumerate(params(imply))) + which_ridge = getindex.(Ref(par2ind), which_ridge) end end which = [CartesianIndex(x) for x in which_ridge] diff --git a/src/types.jl b/src/types.jl index f026b2cb0..70e0b4fbf 100644 --- a/src/types.jl +++ b/src/types.jl @@ -13,6 +13,23 @@ abstract type AbstractSemCollection <: AbstractSem end "Supertype for all loss functions of SEMs. If you want to implement a custom loss function, it should be a subtype of `SemLossFunction`." abstract type SemLossFunction end +""" + params(semobj) + +Return the vector of SEM model parameters. +""" +params(model::AbstractSem) = model.params + +""" + n_par(semobj) + +Return the number of SEM model parameters. +""" +n_par(model::AbstractSem) = length(params(model)) + +params(model::AbstractSemSingle) = params(model.imply) +n_par(model::AbstractSemSingle) = n_par(model.imply) + """ SemLoss(args...; loss_weights = nothing, ...) @@ -75,6 +92,9 @@ If you would like to implement a different notation, e.g. LISREL, you should imp """ abstract type SemImply end +params(imply::SemImply) = params(imply.ram_matrices) +n_par(imply::SemImply) = n_par(imply.ram_matrices) + "Subtype of SemImply for all objects that can serve as the imply field of a SEM and use some form of symbolic precomputation." abstract type SemImplySymbolic <: SemImply end @@ -153,14 +173,14 @@ Returns a SemEnsemble with fields - `sems::Tuple`: `AbstractSem`s. - `weights::Vector`: Weights for each model. - `optimizer::SemOptimizer`: Connects the model to the optimizer. See also [`SemOptimizer`](@ref). -- `param_indices::Dict`: Stores parameter labels and their position. +- `params::Vector`: Stores parameter labels and their position. """ struct SemEnsemble{N, T <: Tuple, V <: AbstractVector, D, I} <: AbstractSemCollection n::N sems::T weights::V optimizer::D - param_indices::I + params::I end function SemEnsemble(models...; optimizer = SemOptimizerOptim, weights = nothing, kwargs...) @@ -175,9 +195,9 @@ function SemEnsemble(models...; optimizer = SemOptimizerOptim, weights = nothing end # check parameters equality - par_indices = param_indices(models[1]) + params = SEM.params(models[1]) for model in models - if par_indices != param_indices(model) + if params != SEM.params(model) throw(ErrorException("The parameters of your models do not match. \n Maybe you tried to specify models of an ensemble via ParameterTables. \n In that case, you may use RAMMatrices instead.")) @@ -189,9 +209,12 @@ function SemEnsemble(models...; optimizer = SemOptimizerOptim, weights = nothing optimizer = optimizer(; kwargs...) end - return SemEnsemble(n, models, weights, optimizer, par_indices) + return SemEnsemble(n, models, weights, optimizer, params) end +params(ensemble::SemEnsemble) = ensemble.params +n_par(ensemble::SemEnsemble) = length(ensemble.params) + """ n_models(ensemble::SemEnsemble) -> Integer @@ -253,4 +276,7 @@ Base type for all SEM specifications. """ abstract type SemSpecification end +params(spec::SemSpecification) = spec.params +n_par(spec::SemSpecification) = length(params(spec)) + abstract type AbstractParameterTable <: SemSpecification end diff --git a/test/unit_tests/specification.jl b/test/unit_tests/specification.jl index 42ad5e431..c081dc0f9 100644 --- a/test/unit_tests/specification.jl +++ b/test/unit_tests/specification.jl @@ -3,19 +3,12 @@ @test ram_matrices == RAMMatrices(partable) end -@test params_to_indices([:x2, :x10, :x28], model_ml) == [2, 10, 28] - -@testset "params_to_indices" begin - pars = [:θ_1, :θ_7, :θ_21] - @test params_to_indices(pars, model_ml) == params_to_indices(pars, partable) - @test params_to_indices(pars, model_ml) == - params_to_indices(pars, RAMMatrices(partable)) +@testset "params()" begin + @test params(model_ml)[2, 10, 28] == [:x2, :x10, :x28] + @test params(model_ml) == params(partable) + @test params(model_ml) == params(RAMMatrices(partable)) end -# from docstrings: -param_indices = params_to_indices([:λ₁, λ₂], my_fitted_sem) -values = solution(my_fitted_sem)[param_indices] - graph = @StenoGraph begin # measurement model visual → fixed(1.0, 1.0) * x1 + fixed(0.5, 0.5) * x2 + fixed(0.6, 0.8) * x3 From 957aa8c5b1732cb31e648c94f4a3c38062744279 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Thu, 9 May 2024 09:41:38 -0700 Subject: [PATCH 28/56] EnsParTable ctor: enforce same params in tables * fix EnsParTable container to Dict{Symbol, ParTable} * don't use keywords for main params as it complicates dispatch Co-authored-by: Maximilian-Stefan-Ernst <34346372+Maximilian-Stefan-Ernst@users.noreply.github.com> --- .../specification/EnsembleParameterTable.jl | 62 ++++++++++--------- src/frontend/specification/StenoGraphs.jl | 28 ++++----- test/examples/multigroup/multigroup.jl | 4 +- 3 files changed, 46 insertions(+), 48 deletions(-) diff --git a/src/frontend/specification/EnsembleParameterTable.jl b/src/frontend/specification/EnsembleParameterTable.jl index 37d7ef15c..0161cc81c 100644 --- a/src/frontend/specification/EnsembleParameterTable.jl +++ b/src/frontend/specification/EnsembleParameterTable.jl @@ -2,8 +2,9 @@ ### Types ############################################################################################ -mutable struct EnsembleParameterTable{C} <: AbstractParameterTable - tables::C +struct EnsembleParameterTable <: AbstractParameterTable + tables::Dict{Symbol, ParameterTable} + params::Vector{Symbol} end ############################################################################################ @@ -11,9 +12,37 @@ end ############################################################################################ # constuct an empty table -function EnsembleParameterTable(::Nothing) - tables = Dict{Symbol, ParameterTable}() - return EnsembleParameterTable(tables) +EnsembleParameterTable(::Nothing; params::Union{Nothing, Vector{Symbol}} = nothing) = + EnsembleParameterTable( + Dict{Symbol, ParameterTable}(), + isnothing(params) ? Symbol[] : copy(params), + ) + +# dictionary of SEM specifications +function EnsembleParameterTable( + spec_ensemble::AbstractDict{K, V}; + params::Union{Nothing, Vector{Symbol}} = nothing, +) where {K, V <: SemSpecification} + partables = Dict{Symbol, ParameterTable}( + Symbol(group) => convert(ParameterTable, spec; params = params) for + (group, spec) in pairs(spec_ensemble) + ) + + if isnothing(params) + # collect all SEM parameters in ensemble if not specified + # and apply the set to all partables + params = + unique(mapreduce(SEM.params, vcat, values(partables), init = Vector{Symbol}())) + for partable in values(partables) + if partable.params != params + copyto!(resize!(partable.params, length(params)), params) + #throw(ArgumentError("The parameter sets of the SEM specifications in the ensemble do not match.")) + end + end + else + params = copy(params) + end + return EnsembleParameterTable(partables, params) end ############################################################################################ @@ -48,20 +77,6 @@ function DataFrames.DataFrame( end end -############################################################################################ -### get parameter table from RAMMatrices -############################################################################################ - -function EnsembleParameterTable(args...; groups) - partable = EnsembleParameterTable(nothing) - - for (group, ram_matrices) in zip(groups, args) - push!(partable.tables, group => ParameterTable(ram_matrices)) - end - - return partable -end - ############################################################################################ ### Pretty Printing ############################################################################################ @@ -83,15 +98,6 @@ end ### Additional Methods ############################################################################################ -# get the vector of all parameters in the table -# the position of the parameter is based on its first appearance in the table (and the ensemble) -function params(partable::EnsembleParameterTable) - params = mapreduce(vcat, values(partable.tables)) do tbl - tbl.columns[:param] - end - return filter!(!=(:const), unique!(params)) # exclude constants -end - # Variables Sorting ------------------------------------------------------------------------ function sort_vars!(partables::EnsembleParameterTable) diff --git a/src/frontend/specification/StenoGraphs.jl b/src/frontend/specification/StenoGraphs.jl index 1d3332cb9..424878f2c 100644 --- a/src/frontend/specification/StenoGraphs.jl +++ b/src/frontend/specification/StenoGraphs.jl @@ -115,24 +115,18 @@ function EnsembleParameterTable( graph::AbstractStenoGraph; observed_vars, latent_vars, - groups + groups, ) graph = unique(graph) - partable = EnsembleParameterTable(nothing) - - for (i, group) in enumerate(groups) - push!( - partable.tables, - Symbol(group) => ParameterTable( - graph; - observed_vars = observed_vars, - latent_vars = latent_vars, - group = i, - param_prefix = Symbol(:g, i), - ), - ) - end - - return partable + partables = Dict( + group => ParameterTable( + graph; + observed_vars = observed_vars, + latent_vars = latent_vars, + group = i, + param_prefix = Symbol(:g, group), + ) for (i, group) in enumerate(groups) + ) + return EnsembleParameterTable(partables) end diff --git a/test/examples/multigroup/multigroup.jl b/test/examples/multigroup/multigroup.jl index 552a65cfb..e428eba1d 100644 --- a/test/examples/multigroup/multigroup.jl +++ b/test/examples/multigroup/multigroup.jl @@ -69,9 +69,7 @@ specification_g2 = RAMMatrices(; ) partable = EnsembleParameterTable( - specification_g1, - specification_g2; - groups = [:Pasteur, :Grant_White], + Dict(:Pasteur => specification_g1, :Grant_White => specification_g2), ) specification_miss_g1 = nothing From 7def9bfb431543edc8636f0cd62209b0199fd4b1 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Wed, 1 May 2024 15:19:09 -0700 Subject: [PATCH 29/56] formatting fixes --- src/frontend/fit/summary.jl | 9 +++----- src/frontend/specification/ParameterTable.jl | 18 +++++----------- src/frontend/specification/RAMMatrices.jl | 3 +-- src/imply/RAM/generic.jl | 22 ++------------------ src/loss/ML/FIML.jl | 3 +-- src/objective_gradient_hessian.jl | 7 +------ test/examples/helper.jl | 6 ++---- 7 files changed, 15 insertions(+), 53 deletions(-) diff --git a/src/frontend/fit/summary.jl b/src/frontend/fit/summary.jl index 4d6bc6181..85c09590b 100644 --- a/src/frontend/fit/summary.jl +++ b/src/frontend/fit/summary.jl @@ -138,8 +138,7 @@ function sem_summary( ), ) - sorted_columns = - [:from, :parameter_type, :to, :estimate, :param, :value_fixed, :start] + sorted_columns = [:from, :parameter_type, :to, :estimate, :param, :value_fixed, :start] regression_columns = sort_partially(sorted_columns, columns) regression_array = reduce( @@ -166,8 +165,7 @@ function sem_summary( (partable.columns[:to] .== partable.columns[:from]), ) - sorted_columns = - [:from, :parameter_type, :to, :estimate, :param, :value_fixed, :start] + sorted_columns = [:from, :parameter_type, :to, :estimate, :param, :value_fixed, :start] variance_columns = sort_partially(sorted_columns, columns) variance_array = reduce( @@ -194,8 +192,7 @@ function sem_summary( (partable.columns[:to] .!= partable.columns[:from]), ) - sorted_columns = - [:from, :parameter_type, :to, :estimate, :param, :value_fixed, :start] + sorted_columns = [:from, :parameter_type, :to, :estimate, :param, :value_fixed, :start] variance_columns = sort_partially(sorted_columns, columns) variance_array = reduce( diff --git a/src/frontend/specification/ParameterTable.jl b/src/frontend/specification/ParameterTable.jl index 60cc93ec9..3b6ff7ed6 100644 --- a/src/frontend/specification/ParameterTable.jl +++ b/src/frontend/specification/ParameterTable.jl @@ -59,17 +59,8 @@ end ############################################################################################ function Base.show(io::IO, partable::ParameterTable) - relevant_columns = [ - :from, - :parameter_type, - :to, - :free, - :value_fixed, - :start, - :estimate, - :se, - :param, - ] + relevant_columns = + [:from, :parameter_type, :to, :free, :value_fixed, :start, :estimate, :se, :param] shown_columns = filter!( col -> haskey(partable.columns, col) && length(partable.columns[col]) > 0, relevant_columns, @@ -125,6 +116,7 @@ params(partable::ParameterTable) = # Sorting ---------------------------------------------------------------------------------- + struct CyclicModelError <: Exception msg::AbstractString end @@ -149,8 +141,8 @@ function sort_vars!(partable::ParameterTable) ] is_regression = [ - (partype == :→) && (from != Symbol("1")) for - (partype, from) in zip(partable.columns[:parameter_type], partable.columns[:from]) + (partype == :→) && (from != Symbol("1")) for (partype, from) in + zip(partable.columns[:parameter_type], partable.columns[:from]) ] to = partable.columns[:to][is_regression] diff --git a/src/frontend/specification/RAMMatrices.jl b/src/frontend/specification/RAMMatrices.jl index 5d3abe295..138bc1c9b 100644 --- a/src/frontend/specification/RAMMatrices.jl +++ b/src/frontend/specification/RAMMatrices.jl @@ -146,8 +146,7 @@ function RAMMatrices( # F indices F_ind = length(partable.sorted_vars) != 0 ? - findall(∈(Set(partable.observed_vars)), partable.sorted_vars) : - 1:n_observed + findall(∈(Set(partable.observed_vars)), partable.sorted_vars) : 1:n_observed # indices of the colnames colnames = diff --git a/src/imply/RAM/generic.jl b/src/imply/RAM/generic.jl index 0988d8e99..07aec6648 100644 --- a/src/imply/RAM/generic.jl +++ b/src/imply/RAM/generic.jl @@ -65,26 +65,8 @@ Additional interfaces Only available in gradient! calls: - `I_A⁻¹(::RAM)` -> ``(I-A)^{-1}`` """ -mutable struct RAM{ - A1, - A2, - A3, - A4, - A5, - A6, - V2, - I1, - I2, - I3, - M1, - M2, - M3, - M4, - S1, - S2, - S3, - B, -} <: SemImply +mutable struct RAM{A1, A2, A3, A4, A5, A6, V2, I1, I2, I3, M1, M2, M3, M4, S1, S2, S3, B} <: + SemImply Σ::A1 A::A2 S::A3 diff --git a/src/loss/ML/FIML.jl b/src/loss/ML/FIML.jl index 9f5103f1f..d4870ac1b 100644 --- a/src/loss/ML/FIML.jl +++ b/src/loss/ML/FIML.jl @@ -252,5 +252,4 @@ end get_n_nodes(specification::RAMMatrices) = specification.size_F[2] get_n_nodes(specification::ParameterTable) = - length(specification.observed_vars) + - length(specification.latent_vars) + length(specification.observed_vars) + length(specification.latent_vars) diff --git a/src/objective_gradient_hessian.jl b/src/objective_gradient_hessian.jl index 61b78a54f..2debbcd40 100644 --- a/src/objective_gradient_hessian.jl +++ b/src/objective_gradient_hessian.jl @@ -38,12 +38,7 @@ function gradient_hessian!(gradient, hessian, model::AbstractSemSingle, params) gradient_hessian!(gradient, hessian, loss(model), params, model) end -function objective_gradient_hessian!( - gradient, - hessian, - model::AbstractSemSingle, - params, -) +function objective_gradient_hessian!(gradient, hessian, model::AbstractSemSingle, params) fill!(gradient, zero(eltype(gradient))) fill!(hessian, zero(eltype(hessian))) objective_gradient_hessian!(imply(model), params, model) diff --git a/test/examples/helper.jl b/test/examples/helper.jl index e93a1437d..0f10ce838 100644 --- a/test/examples/helper.jl +++ b/test/examples/helper.jl @@ -1,6 +1,5 @@ function test_gradient(model, params; rtol = 1e-10, atol = 0) - true_grad = - FiniteDiff.finite_difference_gradient(Base.Fix1(objective!, model), params) + true_grad = FiniteDiff.finite_difference_gradient(Base.Fix1(objective!, model), params) gradient = similar(params) # F and G @@ -250,8 +249,7 @@ function compare_estimates( if type == :↔ type = "~~" elseif type == :→ - if (from ∈ partable.latent_vars) && - (to ∈ partable.observed_vars) + if (from ∈ partable.latent_vars) && (to ∈ partable.observed_vars) type = "=~" else type = "~" From becc6b52a4ed34262bd9a9b3971f35173a0da963 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Wed, 1 May 2024 15:15:09 -0700 Subject: [PATCH 30/56] ParTable ctor: allow providing columns data --- src/frontend/specification/ParameterTable.jl | 30 +++++++++++--------- src/frontend/specification/StenoGraphs.jl | 18 ++++++------ 2 files changed, 25 insertions(+), 23 deletions(-) diff --git a/src/frontend/specification/ParameterTable.jl b/src/frontend/specification/ParameterTable.jl index 3b6ff7ed6..38f84917b 100644 --- a/src/frontend/specification/ParameterTable.jl +++ b/src/frontend/specification/ParameterTable.jl @@ -13,23 +13,25 @@ end ### Constructors ############################################################################################ -# constuct an empty table -function ParameterTable(; +# construct a dictionary with the default partable columns +# optionally pre-allocate data for nrows +empty_partable_columns(nrows::Integer = 0) = Dict{Symbol, Vector}( + :from => fill(Symbol(), nrows), + :parameter_type => fill(Symbol(), nrows), + :to => fill(Symbol(), nrows), + :free => fill(true, nrows), + :value_fixed => fill(NaN, nrows), + :start => fill(NaN, nrows), + :estimate => fill(NaN, nrows), + :param => fill(Symbol(), nrows), +) + +# construct using the provided columns data or create and empty table +function ParameterTable( + columns::Dict{Symbol, Vector} = empty_partable_columns(); observed_vars::Union{AbstractVector{Symbol}, Nothing} = nothing, latent_vars::Union{AbstractVector{Symbol}, Nothing} = nothing, ) - columns = Dict{Symbol, Any}( - :from => Vector{Symbol}(), - :parameter_type => Vector{Symbol}(), - :to => Vector{Symbol}(), - :free => Vector{Bool}(), - :value_fixed => Vector{Float64}(), - :start => Vector{Float64}(), - :estimate => Vector{Float64}(), - :param => Vector{Symbol}(), - :start => Vector{Float64}(), - ) - return ParameterTable(columns, !isnothing(observed_vars) ? copy(observed_vars) : Vector{Symbol}(), !isnothing(latent_vars) ? copy(latent_vars) : Vector{Symbol}(), diff --git a/src/frontend/specification/StenoGraphs.jl b/src/frontend/specification/StenoGraphs.jl index 424878f2c..42edd6a13 100644 --- a/src/frontend/specification/StenoGraphs.jl +++ b/src/frontend/specification/StenoGraphs.jl @@ -41,14 +41,14 @@ function ParameterTable( graph = unique(graph) n = length(graph) - partable = ParameterTable(latent_vars = latent_vars, observed_vars = observed_vars) - from = resize!(partable.columns[:from], n) - parameter_type = resize!(partable.columns[:parameter_type], n) - to = resize!(partable.columns[:to], n) - free = fill!(resize!(partable.columns[:free], n), true) - value_fixed = fill!(resize!(partable.columns[:value_fixed], n), NaN) - start = fill!(resize!(partable.columns[:start], n), NaN) - param_refs = fill!(resize!(partable.columns[:param], n), Symbol("")) + columns = empty_partable_columns(n) + from = columns[:from] + parameter_type = columns[:parameter_type] + to = columns[:to] + free = columns[:free] + value_fixed = columns[:value_fixed] + start = columns[:start] + param_refs = columns[:param] # group = Vector{Symbol}(undef, n) for (i, element) in enumerate(graph) @@ -104,7 +104,7 @@ function ParameterTable( end end - return partable + return ParameterTable(columns; latent_vars, observed_vars) end ############################################################################################ From c29b0c753a2e3a7cbce4fd3d2f19e917ee52d9da Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Wed, 1 May 2024 15:28:45 -0700 Subject: [PATCH 31/56] update_partable!() cleanup + docstring --- src/frontend/specification/ParameterTable.jl | 39 ++++++++++---------- 1 file changed, 19 insertions(+), 20 deletions(-) diff --git a/src/frontend/specification/ParameterTable.jl b/src/frontend/specification/ParameterTable.jl index 38f84917b..52344d72e 100644 --- a/src/frontend/specification/ParameterTable.jl +++ b/src/frontend/specification/ParameterTable.jl @@ -201,36 +201,35 @@ end # update generic --------------------------------------------------------------------------- +""" + update_partable!(partable::AbstractParameterTable, params::Vector{Symbol}, values, column) + +Write parameter `values` into `column` of `partable`. + +The `params` and `values` vectors define the pairs of value +parameters, which are being matched to the `:param` column +of the `partable`. +""" function update_partable!( partable::ParameterTable, params::AbstractVector{Symbol}, values::AbstractVector, - column, + column::Symbol, ) - new_col = Vector{eltype(vec)}(undef, length(partable)) - params_index = Dict(param => i for (i, param) in enumerate(params)) + length(params) == length(values) || throw( + ArgumentError( + "The length of `params` ($(length(params))) and their `values` ($(length(values))) must be the same", + ), + ) + coldata = get!(() -> Vector{eltype(values)}(), partable.columns, column) + resize!(coldata, length(partable)) + params_index = Dict(zip(params, eachindex(params))) for (i, param) in enumerate(partable.columns[:param]) - if !(param == :const) - new_col[i] = values[params_index[param]] - elseif param == :const - new_col[i] = zero(eltype(values)) - end + coldata[i] = param != :const ? values[params_index[param]] : zero(eltype(values)) end - push!(partable.columns, column => new_col) return partable end -""" - update_partable!(partable::AbstractParameterTable, sem_fit::SemFit, vec, column) - -Write `vec` to `column` of `partable`. - -# Arguments -- `vec::Vector`: has to be in the same order as the `model` parameters -""" -update_partable!(partable::AbstractParameterTable, sem_fit::SemFit, values, column) = - update_partable!(partable, params(sem_fit), values, column) - # update estimates ------------------------------------------------------------------------- """ update_estimate!( From b9d8bc5d5957121ef9c6cdd4fa17ab4c99bb3288 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Wed, 1 May 2024 15:29:48 -0700 Subject: [PATCH 32/56] update_partable!(): SemFit methods use basic one --- src/frontend/specification/ParameterTable.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/frontend/specification/ParameterTable.jl b/src/frontend/specification/ParameterTable.jl index 52344d72e..4749ecc8b 100644 --- a/src/frontend/specification/ParameterTable.jl +++ b/src/frontend/specification/ParameterTable.jl @@ -239,7 +239,7 @@ end Write parameter estimates from `sem_fit` to the `:estimate` column of `partable` """ update_estimate!(partable::AbstractParameterTable, sem_fit::SemFit) = - update_partable!(partable, sem_fit, sem_fit.solution, :estimate) + update_partable!(partable, params(sem_fit), sem_fit.solution, :estimate) # update starting values ------------------------------------------------------------------- """ @@ -254,7 +254,7 @@ Write starting values from `sem_fit` or `start_val` to the `:estimate` column of - `kwargs...`: are passed to `start_val` """ update_start!(partable::AbstractParameterTable, sem_fit::SemFit) = - update_partable!(partable, sem_fit, sem_fit.start_val, :start) + update_partable!(partable, params(sem_fit), sem_fit.start_val, :start) function update_start!( partable::AbstractParameterTable, @@ -289,5 +289,5 @@ function update_se_hessian!( hessian = :finitediff, ) se = se_hessian(sem_fit; hessian = hessian) - return update_partable!(partable, sem_fit, se, :se) + return update_partable!(partable, params(sem_fit), se, :se) end From 6a484f1338ba010dfac6db641036237505581aee Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Wed, 1 May 2024 15:34:16 -0700 Subject: [PATCH 33/56] ParTable: add explicit params field --- .../specification/EnsembleParameterTable.jl | 24 ++++------- src/frontend/specification/ParameterTable.jl | 43 +++++++++++++++---- src/frontend/specification/RAMMatrices.jl | 40 ++++++++++++++--- 3 files changed, 77 insertions(+), 30 deletions(-) diff --git a/src/frontend/specification/EnsembleParameterTable.jl b/src/frontend/specification/EnsembleParameterTable.jl index 0161cc81c..10b59fa15 100644 --- a/src/frontend/specification/EnsembleParameterTable.jl +++ b/src/frontend/specification/EnsembleParameterTable.jl @@ -23,25 +23,19 @@ function EnsembleParameterTable( spec_ensemble::AbstractDict{K, V}; params::Union{Nothing, Vector{Symbol}} = nothing, ) where {K, V <: SemSpecification} - partables = Dict{Symbol, ParameterTable}( - Symbol(group) => convert(ParameterTable, spec; params = params) for - (group, spec) in pairs(spec_ensemble) - ) - - if isnothing(params) + params = if isnothing(params) # collect all SEM parameters in ensemble if not specified # and apply the set to all partables - params = - unique(mapreduce(SEM.params, vcat, values(partables), init = Vector{Symbol}())) - for partable in values(partables) - if partable.params != params - copyto!(resize!(partable.params, length(params)), params) - #throw(ArgumentError("The parameter sets of the SEM specifications in the ensemble do not match.")) - end - end + unique(mapreduce(SEM.params, vcat, values(spec_ensemble), init = Vector{Symbol}())) else - params = copy(params) + copy(params) end + + # convert each model specification to ParameterTable + partables = Dict{Symbol, ParameterTable}( + Symbol(group) => convert(ParameterTable, spec; params) for + (group, spec) in pairs(spec_ensemble) + ) return EnsembleParameterTable(partables, params) end diff --git a/src/frontend/specification/ParameterTable.jl b/src/frontend/specification/ParameterTable.jl index 4749ecc8b..86811b0f4 100644 --- a/src/frontend/specification/ParameterTable.jl +++ b/src/frontend/specification/ParameterTable.jl @@ -7,6 +7,7 @@ struct ParameterTable{C} <: AbstractParameterTable observed_vars::Vector{Symbol} latent_vars::Vector{Symbol} sorted_vars::Vector{Symbol} + params::Vector{Symbol} end ############################################################################################ @@ -31,11 +32,32 @@ function ParameterTable( columns::Dict{Symbol, Vector} = empty_partable_columns(); observed_vars::Union{AbstractVector{Symbol}, Nothing} = nothing, latent_vars::Union{AbstractVector{Symbol}, Nothing} = nothing, + params::Union{AbstractVector{Symbol}, Nothing} = nothing, ) - return ParameterTable(columns, + params = isnothing(params) ? unique!(filter(!=(:const), columns[:param])) : copy(params) + check_params(params, columns[:param]) + return ParameterTable( + columns, !isnothing(observed_vars) ? copy(observed_vars) : Vector{Symbol}(), !isnothing(latent_vars) ? copy(latent_vars) : Vector{Symbol}(), - Vector{Symbol}()) + Vector{Symbol}(), + params, + ) +end + +# new parameter table with different parameters order +function ParameterTable( + partable::ParameterTable; + params::Union{AbstractVector{Symbol}, Nothing} = nothing, +) + isnothing(params) || check_params(params, partable.columns[:param]) + + return ParameterTable( + Dict(col => copy(values) for (col, values) in pairs(partable.columns)), + observed_vars = copy(partable.observed_vars), + latent_vars = copy(partable.latent_vars), + params = params, + ) end ############################################################################################ @@ -46,6 +68,15 @@ function Base.convert(::Type{Dict}, partable::ParameterTable) return partable.columns end +function Base.convert( + ::Type{ParameterTable}, + partable::ParameterTable; + params::Union{AbstractVector{Symbol}, Nothing} = nothing, +) + return isnothing(params) || partable.params == params ? partable : + ParameterTable(partable; params) +end + function DataFrames.DataFrame( partable::ParameterTable; columns::Union{AbstractVector{Symbol}, Nothing} = nothing, @@ -110,12 +141,8 @@ Base.eltype(::Type{<:ParameterTable}) = ParameterTableRow Base.iterate(partable::ParameterTable, i::Integer = 1) = i > length(partable) ? nothing : (partable[i], i + 1) - -# get the vector of all parameters in the table -# the position of the parameter is based on its first appearance in the table (and the ensemble) -params(partable::ParameterTable) = - filter!(!=(:const), unique(partable.columns[:param])) - +params(partable::ParameterTable) = partable.params +n_par(partable::ParameterTable) = length(params(partable)) # Sorting ---------------------------------------------------------------------------------- diff --git a/src/frontend/specification/RAMMatrices.jl b/src/frontend/specification/RAMMatrices.jl index 138bc1c9b..0f79d592b 100644 --- a/src/frontend/specification/RAMMatrices.jl +++ b/src/frontend/specification/RAMMatrices.jl @@ -222,18 +222,40 @@ function RAMMatrices( ) end -Base.convert(::Type{RAMMatrices}, partable::ParameterTable) = RAMMatrices(partable) +Base.convert( + ::Type{RAMMatrices}, + partable::ParameterTable; + params::Union{AbstractVector{Symbol}, Nothing} = nothing, +) = RAMMatrices(partable; params) ############################################################################################ ### get parameter table from RAMMatrices ############################################################################################ -function ParameterTable(ram_matrices::RAMMatrices) - colnames = ram_matrices.colnames +function ParameterTable( + ram_matrices::RAMMatrices; + params::Union{AbstractVector{Symbol}, Nothing} = nothing, + observed_var_prefix::Symbol = :obs, + latent_var_prefix::Symbol = :var, +) + # defer parameter checks until we know which ones are used + if !isnothing(ram_matrices.colnames) + colnames = ram_matrices.colnames + observed_vars = colnames[ram_matrices.F_ind] + latent_vars = colnames[setdiff(eachindex(colnames), ram_matrices.F_ind)] + else + observed_vars = + [Symbol("$(observed_var_prefix)_$i") for i in 1:nobserved_vars(ram_matrices)] + latent_vars = + [Symbol("$(latent_var_prefix)_$i") for i in 1:nlatent_vars(ram_matrices)] + colnames = vcat(observed_vars, latent_vars) + end + # construct an empty table partable = ParameterTable( - observed_vars = colnames[ram_matrices.F_ind], - latent_vars = colnames[setdiff(eachindex(colnames), ram_matrices.F_ind)], + observed_vars = observed_vars, + latent_vars = latent_vars, + params = isnothing(params) ? SEM.params(ram_matrices) : params, ) # constants @@ -254,12 +276,16 @@ function ParameterTable(ram_matrices::RAMMatrices) ram_matrices.size_F[2], ) end + check_params(SEM.params(partable), partable.columns[:param]) return partable end -Base.convert(::Type{<:ParameterTable}, ram_matrices::RAMMatrices) = - ParameterTable(ram_matrices) +Base.convert( + ::Type{<:ParameterTable}, + ram::RAMMatrices; + params::Union{AbstractVector{Symbol}, Nothing} = nothing, +) = ParameterTable(ram; params) ############################################################################################ ### Pretty Printing From 3804cc3814f6e1dc70bad65254bc46a2fdde32a0 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sun, 17 Mar 2024 00:02:00 -0700 Subject: [PATCH 34/56] n_par() -> nparams() for clarity and aligning to Julia naming conventions --- docs/src/developer/imply.md | 4 ++-- src/StructuralEquationModels.jl | 2 +- src/additional_functions/simulation.jl | 2 +- src/frontend/fit/SemFit.jl | 2 +- src/frontend/fit/fitmeasures/AIC.jl | 2 +- src/frontend/fit/fitmeasures/BIC.jl | 2 +- src/frontend/fit/fitmeasures/df.jl | 2 +- src/frontend/fit/fitmeasures/fit_measures.jl | 2 +- src/frontend/fit/summary.jl | 2 +- src/frontend/specification/ParameterTable.jl | 2 +- src/frontend/specification/RAMMatrices.jl | 2 +- src/frontend/specification/Sem.jl | 2 +- src/imply/RAM/generic.jl | 6 +++--- src/imply/RAM/symbolic.jl | 7 +++---- src/imply/empty.jl | 4 ++-- src/loss/regularization/ridge.jl | 12 ++++++------ src/types.jl | 13 ++++++------- test/examples/helper.jl | 4 ++-- test/examples/political_democracy/by_parts.jl | 2 +- .../recover_parameters_twofact.jl | 2 +- 20 files changed, 37 insertions(+), 39 deletions(-) diff --git a/docs/src/developer/imply.md b/docs/src/developer/imply.md index 44e0f6ff4..cb30e40fe 100644 --- a/docs/src/developer/imply.md +++ b/docs/src/developer/imply.md @@ -30,10 +30,10 @@ To make stored computations available to loss functions, simply write a function Additionally, you can specify methods for `gradient` and `hessian` as well as the combinations described in [Custom loss functions](@ref). -The last thing nedded to make it work is a method for `n_par` that takes your imply type and returns the number of parameters of the model: +The last thing nedded to make it work is a method for `nparams` that takes your imply type and returns the number of parameters of the model: ```julia -n_par(imply::MyImply) = ... +nparams(imply::MyImply) = ... ``` Just as described in [Custom loss functions](@ref), you may define a constructor. Typically, this will depend on the `specification = ...` argument that can be a `ParameterTable` or a `RAMMatrices` object. diff --git a/src/StructuralEquationModels.jl b/src/StructuralEquationModels.jl index 109a913e7..7812fa819 100644 --- a/src/StructuralEquationModels.jl +++ b/src/StructuralEquationModels.jl @@ -154,6 +154,7 @@ export AbstractSem, sort_vars, RAMMatrices, params, + nparams, fit_measures, AIC, BIC, @@ -161,7 +162,6 @@ export AbstractSem, df, fit_measures, minus2ll, - n_par, n_obs, p_value, RMSEA, diff --git a/src/additional_functions/simulation.jl b/src/additional_functions/simulation.jl index 58e8432e1..0dda725c6 100644 --- a/src/additional_functions/simulation.jl +++ b/src/additional_functions/simulation.jl @@ -73,7 +73,7 @@ function swap_observed( # update imply imply = update_observed(imply, new_observed; kwargs...) kwargs[:imply] = imply - kwargs[:n_par] = n_par(imply) + kwargs[:nparams] = nparams(imply) # update loss loss = update_observed(loss, new_observed; kwargs...) diff --git a/src/frontend/fit/SemFit.jl b/src/frontend/fit/SemFit.jl index 19a6d4441..ace9ed320 100644 --- a/src/frontend/fit/SemFit.jl +++ b/src/frontend/fit/SemFit.jl @@ -47,7 +47,7 @@ end ############################################################################################ params(fit::SemFit) = params(fit.model) -n_par(fit::SemFit) = n_par(fit.model) +nparams(fit::SemFit) = nparams(fit.model) # access fields minimum(sem_fit::SemFit) = sem_fit.minimum diff --git a/src/frontend/fit/fitmeasures/AIC.jl b/src/frontend/fit/fitmeasures/AIC.jl index 519f7beb7..f26f1f4dc 100644 --- a/src/frontend/fit/fitmeasures/AIC.jl +++ b/src/frontend/fit/fitmeasures/AIC.jl @@ -3,4 +3,4 @@ Return the akaike information criterion. """ -AIC(sem_fit::SemFit) = minus2ll(sem_fit) + 2n_par(sem_fit) +AIC(sem_fit::SemFit) = minus2ll(sem_fit) + 2nparams(sem_fit) diff --git a/src/frontend/fit/fitmeasures/BIC.jl b/src/frontend/fit/fitmeasures/BIC.jl index 56200f32b..47bd12f1b 100644 --- a/src/frontend/fit/fitmeasures/BIC.jl +++ b/src/frontend/fit/fitmeasures/BIC.jl @@ -3,4 +3,4 @@ Return the bayesian information criterion. """ -BIC(sem_fit::SemFit) = minus2ll(sem_fit) + log(n_obs(sem_fit)) * n_par(sem_fit) +BIC(sem_fit::SemFit) = minus2ll(sem_fit) + log(n_obs(sem_fit)) * nparams(sem_fit) diff --git a/src/frontend/fit/fitmeasures/df.jl b/src/frontend/fit/fitmeasures/df.jl index f546bb000..d4a4376dd 100644 --- a/src/frontend/fit/fitmeasures/df.jl +++ b/src/frontend/fit/fitmeasures/df.jl @@ -8,7 +8,7 @@ function df end df(sem_fit::SemFit) = df(sem_fit.model) -df(model::AbstractSem) = n_dp(model) - n_par(model) +df(model::AbstractSem) = n_dp(model) - nparams(model) function n_dp(model::AbstractSemSingle) nman = n_man(model) diff --git a/src/frontend/fit/fitmeasures/fit_measures.jl b/src/frontend/fit/fitmeasures/fit_measures.jl index e3f85a0f2..40e3caae0 100644 --- a/src/frontend/fit/fitmeasures/fit_measures.jl +++ b/src/frontend/fit/fitmeasures/fit_measures.jl @@ -1,5 +1,5 @@ fit_measures(sem_fit) = - fit_measures(sem_fit, n_par, df, AIC, BIC, RMSEA, χ², p_value, minus2ll) + fit_measures(sem_fit, nparams, df, AIC, BIC, RMSEA, χ², p_value, minus2ll) function fit_measures(sem_fit, args...) measures = Dict{Symbol, Union{Float64, Missing}}() diff --git a/src/frontend/fit/summary.jl b/src/frontend/fit/summary.jl index 85c09590b..a31d4796f 100644 --- a/src/frontend/fit/summary.jl +++ b/src/frontend/fit/summary.jl @@ -16,7 +16,7 @@ function sem_summary( println("Convergence: $(convergence(sem_fit))") println("No. iterations/evaluations: $(n_iterations(sem_fit))") print("\n") - println("Number of parameters: $(n_par(sem_fit))") + println("Number of parameters: $(nparams(sem_fit))") println("Number of observations: $(n_obs(sem_fit))") print("\n") printstyled( diff --git a/src/frontend/specification/ParameterTable.jl b/src/frontend/specification/ParameterTable.jl index 86811b0f4..fd34570fb 100644 --- a/src/frontend/specification/ParameterTable.jl +++ b/src/frontend/specification/ParameterTable.jl @@ -142,7 +142,7 @@ Base.iterate(partable::ParameterTable, i::Integer = 1) = i > length(partable) ? nothing : (partable[i], i + 1) params(partable::ParameterTable) = partable.params -n_par(partable::ParameterTable) = length(params(partable)) +nparams(partable::ParameterTable) = length(params(partable)) # Sorting ---------------------------------------------------------------------------------- diff --git a/src/frontend/specification/RAMMatrices.jl b/src/frontend/specification/RAMMatrices.jl index 0f79d592b..80ac34cf9 100644 --- a/src/frontend/specification/RAMMatrices.jl +++ b/src/frontend/specification/RAMMatrices.jl @@ -63,7 +63,7 @@ struct RAMMatrices <: SemSpecification size_F::Tuple{Int, Int} end -n_par(ram::RAMMatrices) = length(ram.A_ind) +nparams(ram::RAMMatrices) = length(ram.A_ind) ############################################################################################ ### Constructor diff --git a/src/frontend/specification/Sem.jl b/src/frontend/specification/Sem.jl index 208ef3000..73d4e81da 100644 --- a/src/frontend/specification/Sem.jl +++ b/src/frontend/specification/Sem.jl @@ -73,7 +73,7 @@ function get_fields!(kwargs, observed, imply, loss, optimizer) end kwargs[:imply] = imply - kwargs[:n_par] = n_par(imply) + kwargs[:nparams] = nparams(imply) # loss loss = get_SemLoss(loss; kwargs...) diff --git a/src/imply/RAM/generic.jl b/src/imply/RAM/generic.jl index 07aec6648..9eb694d51 100644 --- a/src/imply/RAM/generic.jl +++ b/src/imply/RAM/generic.jl @@ -34,8 +34,8 @@ and for models with a meanstructure, the model implied means are computed as ``` ## Interfaces -- `params(::RAM) `-> Dict containing the parameter labels and their position -- `n_par(::RAM)` -> Number of parameters +- `params(::RAM) `-> vector of parameter labels +- `nparams(::RAM)` -> number of parameters - `Σ(::RAM)` -> model implied covariance matrix - `μ(::RAM)` -> model implied mean vector @@ -107,7 +107,7 @@ function RAM(; ram_matrices = convert(RAMMatrices, specification) # get dimensions of the model - n_par = length(ram_matrices.params) + n_par = nparams(ram_matrices) n_var, n_nod = ram_matrices.size_F F = zeros(ram_matrices.size_F) F[CartesianIndex.(1:n_var, ram_matrices.F_ind)] .= 1.0 diff --git a/src/imply/RAM/symbolic.jl b/src/imply/RAM/symbolic.jl index db5997d2b..6eb372d4d 100644 --- a/src/imply/RAM/symbolic.jl +++ b/src/imply/RAM/symbolic.jl @@ -29,8 +29,8 @@ Subtype of `SemImply` that implements the RAM notation with symbolic precomputat Subtype of `SemImply`. ## Interfaces -- `params(::RAMSymbolic) `-> Dict containing the parameter labels and their position -- `n_par(::RAMSymbolic)` -> Number of parameters +- `params(::RAMSymbolic) `-> vector of parameter ids +- `nparams(::RAMSymbolic)` -> number of parameters - `Σ(::RAMSymbolic)` -> model implied covariance matrix - `μ(::RAMSymbolic)` -> model implied mean vector @@ -97,7 +97,7 @@ function RAMSymbolic(; ) ram_matrices = convert(RAMMatrices, specification) - n_par = length(ram_matrices.params) + n_par = nparams(ram_matrices) n_var, n_nod = ram_matrices.size_F par = (Symbolics.@variables θ[1:n_par])[1] @@ -141,7 +141,6 @@ function RAMSymbolic(; if hessian & !approximate_hessian n_sig = length(Σ_symbolic) - n_par = size(par, 1) ∇²Σ_symbolic_vec = [Symbolics.sparsehessian(σᵢ, [par...]) for σᵢ in vec(Σ_symbolic)] @variables J[1:n_sig] diff --git a/src/imply/empty.jl b/src/imply/empty.jl index 56297ea06..f1af2ec42 100644 --- a/src/imply/empty.jl +++ b/src/imply/empty.jl @@ -19,8 +19,8 @@ model per group and an additional model with `ImplyEmpty` and `SemRidge` for the # Extended help ## Interfaces -- `params(::RAMSymbolic) `-> Dict containing the parameter labels and their position -- `n_par(::RAMSymbolic)` -> Number of parameters +- `params(::RAMSymbolic) `-> Vector of parameter labels +- `nparams(::RAMSymbolic)` -> Number of parameters ## Implementation Subtype of `SemImply`. diff --git a/src/loss/regularization/ridge.jl b/src/loss/regularization/ridge.jl index 2d098d550..a61dd2af0 100644 --- a/src/loss/regularization/ridge.jl +++ b/src/loss/regularization/ridge.jl @@ -8,18 +8,18 @@ Ridge regularization. # Constructor - SemRidge(;α_ridge, which_ridge, n_par, parameter_type = Float64, imply = nothing, kwargs...) + SemRidge(;α_ridge, which_ridge, nparams, parameter_type = Float64, imply = nothing, kwargs...) # Arguments - `α_ridge`: hyperparameter for penalty term - `which_ridge::Vector`: Vector of parameter labels (Symbols) or indices that indicate which parameters should be regularized. -- `n_par::Int`: number of parameters of the model +- `nparams::Int`: number of parameters of the model - `imply::SemImply`: imply part of the model - `parameter_type`: type of the parameters # Examples ```julia -my_ridge = SemRidge(;α_ridge = 0.02, which_ridge = [:λ₁, :λ₂, :ω₂₃], n_par = 30, imply = my_imply) +my_ridge = SemRidge(;α_ridge = 0.02, which_ridge = [:λ₁, :λ₂, :ω₂₃], nparams = 30, imply = my_imply) ``` # Interfaces @@ -45,7 +45,7 @@ end function SemRidge(; α_ridge, which_ridge, - n_par, + nparams, parameter_type = Float64, imply = nothing, kwargs..., @@ -68,8 +68,8 @@ function SemRidge(; α_ridge, which, which_H, - zeros(parameter_type, n_par), - zeros(parameter_type, n_par, n_par), + zeros(parameter_type, nparams), + zeros(parameter_type, nparams, nparams), ) end diff --git a/src/types.jl b/src/types.jl index 70e0b4fbf..9db48b6db 100644 --- a/src/types.jl +++ b/src/types.jl @@ -21,14 +21,14 @@ Return the vector of SEM model parameters. params(model::AbstractSem) = model.params """ - n_par(semobj) + nparams(semobj) Return the number of SEM model parameters. """ -n_par(model::AbstractSem) = length(params(model)) +nparams(model::AbstractSem) = length(params(model)) params(model::AbstractSemSingle) = params(model.imply) -n_par(model::AbstractSemSingle) = n_par(model.imply) +nparams(model::AbstractSemSingle) = nparams(model.imply) """ SemLoss(args...; loss_weights = nothing, ...) @@ -93,7 +93,7 @@ If you would like to implement a different notation, e.g. LISREL, you should imp abstract type SemImply end params(imply::SemImply) = params(imply.ram_matrices) -n_par(imply::SemImply) = n_par(imply.ram_matrices) +nparams(imply::SemImply) = nparams(imply.ram_matrices) "Subtype of SemImply for all objects that can serve as the imply field of a SEM and use some form of symbolic precomputation." abstract type SemImplySymbolic <: SemImply end @@ -185,7 +185,6 @@ end function SemEnsemble(models...; optimizer = SemOptimizerOptim, weights = nothing, kwargs...) n = length(models) - npar = n_par(models[1]) # default weights @@ -213,7 +212,7 @@ function SemEnsemble(models...; optimizer = SemOptimizerOptim, weights = nothing end params(ensemble::SemEnsemble) = ensemble.params -n_par(ensemble::SemEnsemble) = length(ensemble.params) +nparams(ensemble::SemEnsemble) = length(ensemble.params) """ n_models(ensemble::SemEnsemble) -> Integer @@ -277,6 +276,6 @@ Base type for all SEM specifications. abstract type SemSpecification end params(spec::SemSpecification) = spec.params -n_par(spec::SemSpecification) = length(params(spec)) +nparams(spec::SemSpecification) = length(params(spec)) abstract type AbstractParameterTable <: SemSpecification end diff --git a/test/examples/helper.jl b/test/examples/helper.jl index 0f10ce838..0230dd497 100644 --- a/test/examples/helper.jl +++ b/test/examples/helper.jl @@ -46,7 +46,7 @@ fitmeasure_names_ml = Dict( :df => "df", :χ² => "chisq", :p_value => "pvalue", - :n_par => "npar", + :nparams => "npar", :RMSEA => "rmsea", ) @@ -54,7 +54,7 @@ fitmeasure_names_ls = Dict( :df => "df", :χ² => "chisq", :p_value => "pvalue", - :n_par => "npar", + :nparams => "npar", :RMSEA => "rmsea", ) diff --git a/test/examples/political_democracy/by_parts.jl b/test/examples/political_democracy/by_parts.jl index c071d9e00..c06c9929d 100644 --- a/test/examples/political_democracy/by_parts.jl +++ b/test/examples/political_democracy/by_parts.jl @@ -15,7 +15,7 @@ ml = SemML(observed = observed) wls = SemWLS(observed = observed) -ridge = SemRidge(α_ridge = 0.001, which_ridge = 16:20, n_par = 31) +ridge = SemRidge(α_ridge = 0.001, which_ridge = 16:20, nparams = 31) constant = SemConstant(constant_loss = 3.465) diff --git a/test/examples/recover_parameters/recover_parameters_twofact.jl b/test/examples/recover_parameters/recover_parameters_twofact.jl index 1bd7136bc..5aa79842c 100644 --- a/test/examples/recover_parameters/recover_parameters_twofact.jl +++ b/test/examples/recover_parameters/recover_parameters_twofact.jl @@ -63,7 +63,7 @@ Random.seed!(1234) x = transpose(rand(true_dist, 100000)) semobserved = SemObservedData(data = x, specification = nothing) -loss_ml = SemLoss(SemML(; observed = semobserved, n_par = length(start))) +loss_ml = SemLoss(SemML(; observed = semobserved, nparams = length(start))) optimizer = SemOptimizerOptim( BFGS(; linesearch = BackTracking(order = 3), alphaguess = InitialHagerZhang()),# m = 100), From ca9f860bb680e190ea992634ac1a7c31ad7037ed Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sun, 26 May 2024 21:24:19 -0700 Subject: [PATCH 35/56] param_values(ParTable) Co-authored-by: Maximilian-Stefan-Ernst <34346372+Maximilian-Stefan-Ernst@users.noreply.github.com> --- src/frontend/specification/ParameterTable.jl | 56 ++++++++++++++++++++ 1 file changed, 56 insertions(+) diff --git a/src/frontend/specification/ParameterTable.jl b/src/frontend/specification/ParameterTable.jl index fd34570fb..0e3125665 100644 --- a/src/frontend/specification/ParameterTable.jl +++ b/src/frontend/specification/ParameterTable.jl @@ -318,3 +318,59 @@ function update_se_hessian!( se = se_hessian(sem_fit; hessian = hessian) return update_partable!(partable, params(sem_fit), se, :se) end + +""" + param_values!(out::AbstractVector, partable::ParameterTable, + col::Symbol = :estimate) + +Extract parameter values from the `col` column of `partable` +into the `out` vector. + +The `out` vector should be of `nparams(partable)` length. +The *i*-th element of the `out` vector will contain the +value of the *i*-th parameter from `params(partable)`. + +Note that the function combines the duplicate occurences of the +same parameter in `partable` and will raise an error if the +values do not match. +""" +function param_values!( + out::AbstractVector, + partable::ParameterTable, + col::Symbol = :estimate, +) + (length(out) == nparams(partable)) || throw( + DimensionMismatch( + "The length of parameter values vector ($(length(out))) does not match the number of parameters ($(nparams(partable)))", + ), + ) + param_index = Dict(param => i for (i, param) in enumerate(params(partable))) + param_values_col = partable.columns[col] + for (i, param) in enumerate(partable.columns[:param]) + (param == :const) && continue + param_ind = get(param_index, param, nothing) + @assert !isnothing(param_ind) "Parameter table contains unregistered parameter :$param at row #$i" + val = param_values_col[i] + if !isnan(out[param_ind]) + @assert isequal(out[param_ind], val) "Parameter :$param value at row #$i ($val) differs from the earlier encountered value ($(out[param_ind]))" + else + out[param_ind] = val + end + end + return out +end + +""" + param_values(out::AbstractVector, col::Symbol = :estimate) + +Extract parameter values from the `col` column of `partable`. + +Returns the values vector. The *i*-th element corresponds to +the value of *i*-th parameter from `params(partable)`. + +Note that the function combines the duplicate occurences of the +same parameter in `partable` and will raise an error if the +values do not match. +""" +param_values(partable::ParameterTable, col::Symbol = :estimate) = + param_values!(fill(NaN, nparams(partable)), partable, col) From b117b34cb03b42fc848dd9886ec7b64ffbbd497f Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Wed, 1 May 2024 16:46:34 -0700 Subject: [PATCH 36/56] lavaan_param_values(lav_fit, partable) --- src/frontend/specification/ParameterTable.jl | 142 +++++++++++++++++++ 1 file changed, 142 insertions(+) diff --git a/src/frontend/specification/ParameterTable.jl b/src/frontend/specification/ParameterTable.jl index 0e3125665..43d0e1b11 100644 --- a/src/frontend/specification/ParameterTable.jl +++ b/src/frontend/specification/ParameterTable.jl @@ -374,3 +374,145 @@ values do not match. """ param_values(partable::ParameterTable, col::Symbol = :estimate) = param_values!(fill(NaN, nparams(partable)), partable, col) + +""" + lavaan_param_values!(out::AbstractVector, partable_lav, + partable::ParameterTable, + lav_col::Symbol = :est, lav_group = nothing) + +Extract parameter values from the `partable_lav` lavaan model that +match the parameters of `partable` into the `out` vector. + +The method sets the *i*-th element of the `out` vector to +the value of *i*-th parameter from `params(partable)`. + +Note that the lavaan and `partable` models are matched by the +the names of variables in the tables (`from` and `to` columns) +as well as the type of their relationship (`relation` column), +and not by the names of the model parameters. +""" +function lavaan_param_values!( + out::AbstractVector, + partable_lav, + partable::ParameterTable, + lav_col::Symbol = :est, + lav_group = nothing, +) + + # find indices of all df row where f is true + findallrows(f::Function, df) = findall(f(r) for r in eachrow(df)) + + (length(out) == nparams(partable)) || throw( + DimensionMismatch( + "The length of parameter values vector ($(length(out))) does not match the number of parameters ($(nparams(partable)))", + ), + ) + partable_mask = findall(partable.columns[:free]) + param_index = Dict(param => i for (i, param) in enumerate(params(partable))) + + lav_values = partable_lav[:, lav_col] + for (from, to, type, id) in zip( + [ + view(partable.columns[k], partable_mask) for + k in [:from, :to, :parameter_type, :param] + ]..., + ) + lav_ind = nothing + + if from == Symbol("1") + lav_ind = findallrows( + r -> + r[:lhs] == String(to) && + r[:op] == "~1" && + (isnothing(lav_group) || r[:group] == lav_group), + partable_lav, + ) + else + if type == :↔ + lav_type = "~~" + elseif type == :→ + if (from ∈ partable.latent_vars) && (to ∈ partable.observed_vars) + lav_type = "=~" + else + lav_type = "~" + from, to = to, from + end + end + + if lav_type == "~~" + lav_ind = findallrows( + r -> + ( + (r[:lhs] == String(from) && r[:rhs] == String(to)) || + (r[:lhs] == String(to) && r[:rhs] == String(from)) + ) && + r[:op] == lav_type && + (isnothing(lav_group) || r[:group] == lav_group), + partable_lav, + ) + else + lav_ind = findallrows( + r -> + r[:lhs] == String(from) && + r[:rhs] == String(to) && + r[:op] == lav_type && + (isnothing(lav_group) || r[:group] == lav_group), + partable_lav, + ) + end + end + + if length(lav_ind) == 0 + throw( + ErrorException( + "Parameter $id ($from $type $to) could not be found in the lavaan solution", + ), + ) + elseif length(lav_ind) > 1 + throw( + ErrorException( + "At least one parameter was found twice in the lavaan solution", + ), + ) + end + + param_ind = param_index[id] + param_val = lav_values[lav_ind[1]] + if isnan(out[param_ind]) + out[param_ind] = param_val + else + @assert out[param_ind] ≈ param_val atol = 1E-10 "Parameter :$id value at row #$lav_ind ($param_val) differs from the earlier encountered value ($(out[param_ind]))" + end + end + + return out +end + +""" + lavaan_param_values(partable_lav, partable::ParameterTable, + lav_col::Symbol = :est, lav_group = nothing) + +Extract parameter values from the `partable_lav` lavaan model that +match the parameters of `partable`. + +The `out` vector should be of `nparams(partable)` length. +The *i*-th element of the `out` vector will contain the +value of the *i*-th parameter from `params(partable)`. + +Note that the lavaan and `partable` models are matched by the +the names of variables in the tables (`from` and `to` columns), +and the type of their relationship (`relation` column), +but not by the ids of the model parameters. +""" +lavaan_param_values( + partable_lav, + partable::ParameterTable, + lav_col::Symbol = :est, + lav_group = nothing, +) = lavaan_param_values!( + fill(NaN, nparams(partable)), + partable_lav, + partable, + lav_col, + lav_group, +) From 398870b23d9617752885b00eb26c87313e225ab0 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Wed, 1 May 2024 16:46:54 -0700 Subject: [PATCH 37/56] compare_estimates() -> test_estimates() * do tests inside * use param_values()/lavaan_param_values() --- test/examples/helper.jl | 283 +++--------------- test/examples/multigroup/build_models.jl | 18 +- test/examples/political_democracy/by_parts.jl | 37 +-- .../political_democracy/constructor.jl | 37 +-- test/unit_tests/sorting.jl | 2 +- 5 files changed, 79 insertions(+), 298 deletions(-) diff --git a/test/examples/helper.jl b/test/examples/helper.jl index 0230dd497..d4c140d67 100644 --- a/test/examples/helper.jl +++ b/test/examples/helper.jl @@ -1,3 +1,5 @@ +using LinearAlgebra: norm + function test_gradient(model, params; rtol = 1e-10, atol = 0) true_grad = FiniteDiff.finite_difference_gradient(Base.Fix1(objective!, model), params) gradient = similar(params) @@ -71,262 +73,63 @@ function test_fitmeasures( end end -function compare_estimates( +function test_estimates( partable::ParameterTable, partable_lav; rtol = 1e-10, atol = 0, col = :estimate, lav_col = :est, + lav_group = nothing, + skip::Bool = false, ) - correct = [] - - for i in findall(partable.columns[:free]) - from = partable.columns[:from][i] - to = partable.columns[:to][i] - type = partable.columns[:parameter_type][i] - estimate = partable.columns[col][i] - - if from == Symbol("1") - lav_ind = - findall((partable_lav.lhs .== String(to)) .& (partable_lav.op .== "~1")) - - if length(lav_ind) == 0 - throw( - ErrorException( - "Parameter from: $from, to: $to, type: $type, could not be found in the lavaan solution", - ), - ) - elseif length(lav_ind) > 1 - throw( - ErrorException( - "At least one parameter was found twice in the lavaan solution", - ), - ) - else - is_correct = isapprox( - estimate, - partable_lav[:, lav_col][lav_ind[1]]; - rtol = rtol, - atol = atol, - ) - push!(correct, is_correct) - end - - else - if type == :↔ - type = "~~" - elseif type == :→ - if (from ∈ partable.latent_vars) && (to ∈ partable.observed_vars) - type = "=~" - else - type = "~" - from, to = to, from - end - end - - if type == "~~" - lav_ind = findall( - ( - ( - (partable_lav.lhs .== String(from)) .& - (partable_lav.rhs .== String(to)) - ) .| ( - (partable_lav.lhs .== String(to)) .& - (partable_lav.rhs .== String(from)) - ) - ) .& (partable_lav.op .== type), - ) - - if length(lav_ind) == 0 - throw( - ErrorException( - "Parameter from: $from, to: $to, type: $type, could not be found in the lavaan solution", - ), - ) - elseif length(lav_ind) > 1 - throw( - ErrorException( - "At least one parameter was found twice in the lavaan solution", - ), - ) - else - is_correct = isapprox( - estimate, - partable_lav[:, lav_col][lav_ind[1]]; - rtol = rtol, - atol = atol, - ) - push!(correct, is_correct) - end - - else - lav_ind = findall( - (partable_lav.lhs .== String(from)) .& - (partable_lav.rhs .== String(to)) .& - (partable_lav.op .== type), - ) - - if length(lav_ind) == 0 - throw( - ErrorException( - "Parameter from: $from, to: $to, type: $type, could not be found in the lavaan solution", - ), - ) - elseif length(lav_ind) > 1 - throw( - ErrorException( - "At least one parameter was found twice in the lavaan solution", - ), - ) - else - is_correct = isapprox( - estimate, - partable_lav[:, lav_col][lav_ind[1]]; - rtol = rtol, - atol = atol, - ) - push!(correct, is_correct) - end - end - end + actual = StructuralEquationModels.param_values(partable, col) + expected = StructuralEquationModels.lavaan_param_values( + partable_lav, + partable, + lav_col, + lav_group, + ) + @test !any(isnan, actual) + @test !any(isnan, expected) + + if skip # workaround skip=false not supported in earlier versions + @test actual ≈ expected rtol = rtol atol = atol norm = Base.Fix2(norm, Inf) skip = + skip + else + @test actual ≈ expected rtol = rtol atol = atol norm = Base.Fix2(norm, Inf) end - - return all(correct) end -function compare_estimates( +function test_estimates( ens_partable::EnsembleParameterTable, partable_lav; rtol = 1e-10, atol = 0, col = :estimate, lav_col = :est, - lav_groups, + lav_groups::AbstractDict, + skip::Bool = false, ) - correct = [] - - for key in keys(ens_partable.tables) - group = lav_groups[key] - partable = ens_partable.tables[key] - - for i in findall(partable.columns[:free]) - from = partable.columns[:from][i] - to = partable.columns[:to][i] - type = partable.columns[:parameter_type][i] - estimate = partable.columns[col][i] - - if from == Symbol("1") - lav_ind = findall( - (partable_lav.lhs .== String(to)) .& - (partable_lav.op .== "~1") .& - (partable_lav.group .== group), - ) - - if length(lav_ind) == 0 - throw( - ErrorException( - "Mean parameter of variable $to could not be found in the lavaan solution", - ), - ) - elseif length(lav_ind) > 1 - throw( - ErrorException( - "At least one parameter was found twice in the lavaan solution", - ), - ) - else - is_correct = isapprox( - estimate, - partable_lav[:, lav_col][lav_ind[1]]; - rtol = rtol, - atol = atol, - ) - push!(correct, is_correct) - end - - else - if type == :↔ - type = "~~" - elseif type == :→ - if (from ∈ partable.latent_vars) && (to ∈ partable.observed_vars) - type = "=~" - else - type = "~" - from, to = to, from - end - end - - if type == "~~" - lav_ind = findall( - ( - ( - (partable_lav.lhs .== String(from)) .& - (partable_lav.rhs .== String(to)) - ) .| ( - (partable_lav.lhs .== String(to)) .& - (partable_lav.rhs .== String(from)) - ) - ) .& - (partable_lav.op .== type) .& - (partable_lav.group .== group), - ) - - if length(lav_ind) == 0 - throw( - ErrorException( - "Parameter from: $from, to: $to, type: $type, could not be found in the lavaan solution", - ), - ) - elseif length(lav_ind) > 1 - throw( - ErrorException( - "At least one parameter was found twice in the lavaan solution", - ), - ) - else - is_correct = isapprox( - estimate, - partable_lav[:, lav_col][lav_ind[1]]; - rtol = rtol, - atol = atol, - ) - push!(correct, is_correct) - end - - else - lav_ind = findall( - (partable_lav.lhs .== String(from)) .& - (partable_lav.rhs .== String(to)) .& - (partable_lav.op .== type) .& - (partable_lav.group .== group), - ) - - if length(lav_ind) == 0 - throw( - ErrorException( - "Parameter $from $type $to could not be found in the lavaan solution", - ), - ) - elseif length(lav_ind) > 1 - throw( - ErrorException( - "At least one parameter was found twice in the lavaan solution", - ), - ) - else - is_correct = isapprox( - estimate, - partable_lav[:, lav_col][lav_ind[1]]; - rtol = rtol, - atol = atol, - ) - push!(correct, is_correct) - end - end - end - end + actual = fill(NaN, nparams(ens_partable)) + expected = fill(NaN, nparams(ens_partable)) + for (key, partable) in pairs(ens_partable.tables) + StructuralEquationModels.param_values!(actual, partable, col) + StructuralEquationModels.lavaan_param_values!( + expected, + partable_lav, + partable, + lav_col, + lav_groups[key], + ) + end + @test !any(isnan, actual) + @test !any(isnan, expected) + + if skip # workaround skip=false not supported in earlier versions + @test actual ≈ expected rtol = rtol atol = atol norm = Base.Fix2(norm, Inf) skip = + skip + else + @test actual ≈ expected rtol = rtol atol = atol norm = Base.Fix2(norm, Inf) end - - return all(correct) end diff --git a/test/examples/multigroup/build_models.jl b/test/examples/multigroup/build_models.jl index 23d429796..70d2bb914 100644 --- a/test/examples/multigroup/build_models.jl +++ b/test/examples/multigroup/build_models.jl @@ -17,7 +17,7 @@ end @testset "ml_solution_multigroup" begin solution = sem_fit(model_ml_multigroup) update_estimate!(partable, solution) - @test compare_estimates( + test_estimates( partable, solution_lav[:parameter_estimates_ml]; atol = 1e-4, @@ -35,7 +35,7 @@ end ) update_se_hessian!(partable, solution_ml) - @test compare_estimates( + test_estimates( partable, solution_lav[:parameter_estimates_ml]; atol = 1e-3, @@ -78,7 +78,7 @@ grad_fd = FiniteDiff.finite_difference_gradient( @testset "ml_solution_multigroup | sorted" begin solution = sem_fit(model_ml_multigroup) update_estimate!(partable_s, solution) - @test compare_estimates( + test_estimates( partable_s, solution_lav[:parameter_estimates_ml]; atol = 1e-4, @@ -96,7 +96,7 @@ end ) update_se_hessian!(partable_s, solution_ml) - @test compare_estimates( + test_estimates( partable_s, solution_lav[:parameter_estimates_ml]; atol = 1e-3, @@ -152,7 +152,7 @@ end @testset "solution_user_defined_loss" begin solution = sem_fit(model_ml_multigroup) update_estimate!(partable, solution) - @test compare_estimates( + test_estimates( partable, solution_lav[:parameter_estimates_ml]; atol = 1e-4, @@ -179,7 +179,7 @@ end @testset "ls_solution_multigroup" begin solution = sem_fit(model_ls_multigroup) update_estimate!(partable, solution) - @test compare_estimates( + test_estimates( partable, solution_lav[:parameter_estimates_ls]; atol = 1e-4, @@ -198,7 +198,7 @@ end ) update_se_hessian!(partable, solution_ls) - @test compare_estimates( + test_estimates( partable, solution_lav[:parameter_estimates_ls]; atol = 1e-2, @@ -266,7 +266,7 @@ if !isnothing(specification_miss_g1) @testset "fiml_solution_multigroup" begin solution = sem_fit(model_ml_multigroup) update_estimate!(partable_miss, solution) - @test compare_estimates( + test_estimates( partable_miss, solution_lav[:parameter_estimates_fiml]; atol = 1e-4, @@ -284,7 +284,7 @@ if !isnothing(specification_miss_g1) ) update_se_hessian!(partable_miss, solution) - @test compare_estimates( + test_estimates( partable_miss, solution_lav[:parameter_estimates_fiml]; atol = 1e-3, diff --git a/test/examples/political_democracy/by_parts.jl b/test/examples/political_democracy/by_parts.jl index c06c9929d..11953ccb6 100644 --- a/test/examples/political_democracy/by_parts.jl +++ b/test/examples/political_democracy/by_parts.jl @@ -73,7 +73,7 @@ for (model, name, solution_name) in zip(models, model_names, solution_names) @testset "$(name)_solution" begin solution = sem_fit(model) update_estimate!(partable, solution) - @test compare_estimates(partable, solution_lav[solution_name]; atol = 1e-2) + test_estimates(partable, solution_lav[solution_name]; atol = 1e-2) end catch end @@ -114,7 +114,7 @@ end test_fitmeasures(fit_measures(solution_ml), solution_lav[:fitmeasures_ml]; atol = 1e-3) update_se_hessian!(partable, solution_ml) - @test compare_estimates( + test_estimates( partable, solution_lav[:parameter_estimates_ml]; atol = 1e-3, @@ -135,7 +135,7 @@ end @test (fm[:AIC] === missing) & (fm[:BIC] === missing) & (fm[:minus2ll] === missing) update_se_hessian!(partable, solution_ls) - @test compare_estimates( + test_estimates( partable, solution_lav[:parameter_estimates_ls]; atol = 1e-2, @@ -178,21 +178,18 @@ if semoptimizer == SemOptimizerOptim @testset "ml_solution_hessian" begin solution = sem_fit(model_ml) update_estimate!(partable, solution) - @test compare_estimates( - partable, - solution_lav[:parameter_estimates_ml]; - atol = 1e-3, - ) + test_estimates(partable, solution_lav[:parameter_estimates_ml]; atol = 1e-3) end @testset "ls_solution_hessian" begin solution = sem_fit(model_ls) update_estimate!(partable, solution) - @test compare_estimates( + test_estimates( partable, solution_lav[:parameter_estimates_ls]; atol = 1e-3, - ) skip = true + skip = true, + ) end end @@ -260,7 +257,7 @@ for (model, name, solution_name) in zip(models, model_names, solution_names) @testset "$(name)_solution_mean" begin solution = sem_fit(model) update_estimate!(partable_mean, solution) - @test compare_estimates(partable_mean, solution_lav[solution_name]; atol = 1e-2) + test_estimates(partable_mean, solution_lav[solution_name]; atol = 1e-2) end catch end @@ -279,7 +276,7 @@ end ) update_se_hessian!(partable_mean, solution_ml) - @test compare_estimates( + test_estimates( partable_mean, solution_lav[:parameter_estimates_ml_mean]; atol = 0.002, @@ -300,7 +297,7 @@ end @test (fm[:AIC] === missing) & (fm[:BIC] === missing) & (fm[:minus2ll] === missing) update_se_hessian!(partable_mean, solution_ls) - @test compare_estimates( + test_estimates( partable_mean, solution_lav[:parameter_estimates_ls_mean]; atol = 1e-2, @@ -342,21 +339,13 @@ end @testset "fiml_solution" begin solution = sem_fit(model_ml) update_estimate!(partable_mean, solution) - @test compare_estimates( - partable_mean, - solution_lav[:parameter_estimates_fiml]; - atol = 1e-2, - ) + test_estimates(partable_mean, solution_lav[:parameter_estimates_fiml]; atol = 1e-2) end @testset "fiml_solution_symbolic" begin solution = sem_fit(model_ml_sym) update_estimate!(partable_mean, solution) - @test compare_estimates( - partable_mean, - solution_lav[:parameter_estimates_fiml]; - atol = 1e-2, - ) + test_estimates(partable_mean, solution_lav[:parameter_estimates_fiml]; atol = 1e-2) end ############################################################################################ @@ -372,7 +361,7 @@ end ) update_se_hessian!(partable_mean, solution_ml) - @test compare_estimates( + test_estimates( partable_mean, solution_lav[:parameter_estimates_fiml]; atol = 1e-3, diff --git a/test/examples/political_democracy/constructor.jl b/test/examples/political_democracy/constructor.jl index 4ca1994bd..5f1c838e8 100644 --- a/test/examples/political_democracy/constructor.jl +++ b/test/examples/political_democracy/constructor.jl @@ -87,7 +87,7 @@ for (model, name, solution_name) in zip(models, model_names, solution_names) @testset "$(name)_solution" begin solution = sem_fit(model) update_estimate!(partable, solution) - @test compare_estimates(partable, solution_lav[solution_name]; atol = 1e-2) + test_estimates(partable, solution_lav[solution_name]; atol = 1e-2) end catch end @@ -131,7 +131,7 @@ end test_fitmeasures(fit_measures(solution_ml), solution_lav[:fitmeasures_ml]; atol = 1e-3) update_se_hessian!(partable, solution_ml) - @test compare_estimates( + test_estimates( partable, solution_lav[:parameter_estimates_ml]; atol = 1e-3, @@ -152,7 +152,7 @@ end @test ismissing(fm[:AIC]) && ismissing(fm[:BIC]) && ismissing(fm[:minus2ll]) update_se_hessian!(partable, solution_ls) - @test compare_estimates( + test_estimates( partable, solution_lav[:parameter_estimates_ls]; atol = 1e-2, @@ -199,22 +199,19 @@ if semoptimizer == SemOptimizerOptim @testset "ml_solution_hessian" begin solution = sem_fit(model_ml) update_estimate!(partable, solution) - @test compare_estimates( - partable, - solution_lav[:parameter_estimates_ml]; - atol = 1e-3, - ) + test_estimates(partable, solution_lav[:parameter_estimates_ml]; atol = 1e-3) end @testset "ls_solution_hessian" begin solution = sem_fit(model_ls) update_estimate!(partable, solution) - @test compare_estimates( + test_estimates( partable, solution_lav[:parameter_estimates_ls]; atol = 0.002, rtol = 0.0, - ) skip = true + skip = true, + ) end end @@ -286,7 +283,7 @@ for (model, name, solution_name) in zip(models, model_names, solution_names) @testset "$(name)_solution_mean" begin solution = sem_fit(model) update_estimate!(partable_mean, solution) - @test compare_estimates(partable_mean, solution_lav[solution_name]; atol = 1e-2) + test_estimates(partable_mean, solution_lav[solution_name]; atol = 1e-2) end catch end @@ -305,7 +302,7 @@ end ) update_se_hessian!(partable_mean, solution_ml) - @test compare_estimates( + test_estimates( partable_mean, solution_lav[:parameter_estimates_ml_mean]; atol = 0.002, @@ -326,7 +323,7 @@ end @test ismissing(fm[:AIC]) && ismissing(fm[:BIC]) && ismissing(fm[:minus2ll]) update_se_hessian!(partable_mean, solution_ls) - @test compare_estimates( + test_estimates( partable_mean, solution_lav[:parameter_estimates_ls_mean]; atol = 1e-2, @@ -379,21 +376,13 @@ end @testset "fiml_solution" begin solution = sem_fit(model_ml) update_estimate!(partable_mean, solution) - @test compare_estimates( - partable_mean, - solution_lav[:parameter_estimates_fiml]; - atol = 1e-2, - ) + test_estimates(partable_mean, solution_lav[:parameter_estimates_fiml]; atol = 1e-2) end @testset "fiml_solution_symbolic" begin solution = sem_fit(model_ml_sym) update_estimate!(partable_mean, solution) - @test compare_estimates( - partable_mean, - solution_lav[:parameter_estimates_fiml]; - atol = 1e-2, - ) + test_estimates(partable_mean, solution_lav[:parameter_estimates_fiml]; atol = 1e-2) end ############################################################################################ @@ -409,7 +398,7 @@ end ) update_se_hessian!(partable_mean, solution_ml) - @test compare_estimates( + test_estimates( partable_mean, solution_lav[:parameter_estimates_fiml]; atol = 0.002, diff --git a/test/unit_tests/sorting.jl b/test/unit_tests/sorting.jl index 5ca890c51..f5bc38ae0 100644 --- a/test/unit_tests/sorting.jl +++ b/test/unit_tests/sorting.jl @@ -13,5 +13,5 @@ end @testset "ml_solution_sorted" begin solution_ml_sorted = sem_fit(model_ml_sorted) update_estimate!(partable, solution_ml_sorted) - @test SEM.compare_estimates(par_ml, partable, 0.01) + @test test_estimates(par_ml, partable, 0.01) end From a2518b718a1dbf9bf0592efc2113ff1e4c8195d8 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sun, 26 May 2024 21:25:03 -0700 Subject: [PATCH 38/56] update_partable!(): dict-based generic version Co-authored-by: Maximilian-Stefan-Ernst <34346372+Maximilian-Stefan-Ernst@users.noreply.github.com> --- .../specification/EnsembleParameterTable.jl | 23 ++++-- src/frontend/specification/ParameterTable.jl | 81 ++++++++++++++----- 2 files changed, 76 insertions(+), 28 deletions(-) diff --git a/src/frontend/specification/EnsembleParameterTable.jl b/src/frontend/specification/EnsembleParameterTable.jl index 10b59fa15..b0b50448b 100644 --- a/src/frontend/specification/EnsembleParameterTable.jl +++ b/src/frontend/specification/EnsembleParameterTable.jl @@ -117,15 +117,24 @@ Base.getindex(partable::EnsembleParameterTable, group) = partable.tables[group] ### Update Partable from Fitted Model ############################################################################################ -# update generic --------------------------------------------------------------------------- function update_partable!( - partable::EnsembleParameterTable, + partables::EnsembleParameterTable, + column::Symbol, + param_values::AbstractDict{Symbol}, + default::Any = nothing, +) + for partable in values(partables.tables) + update_partable!(partable, column, param_values, default) + end + return partables +end + +function update_partable!( + partables::EnsembleParameterTable, + column::Symbol, params::AbstractVector{Symbol}, values::AbstractVector, - column, + default::Any = nothing, ) - for k in keys(partable.tables) - update_partable!(partable.tables[k], params, values, column) - end - return partable + return update_partable!(partables, column, Dict(zip(params, values)), default) end diff --git a/src/frontend/specification/ParameterTable.jl b/src/frontend/specification/ParameterTable.jl index 43d0e1b11..6dcbbbf84 100644 --- a/src/frontend/specification/ParameterTable.jl +++ b/src/frontend/specification/ParameterTable.jl @@ -227,6 +227,32 @@ end ############################################################################################ # update generic --------------------------------------------------------------------------- +function update_partable!( + partable::ParameterTable, + column::Symbol, + param_values::AbstractDict{Symbol, T}, + default::Any = nothing, +) where {T} + coldata = get!(() -> Vector{T}(undef, length(partable)), partable.columns, column) + + isvec_def = (default isa AbstractVector) && (length(default) == length(partable)) + + for (i, par) in enumerate(partable.columns[:param]) + if par == :const + coldata[i] = !isnothing(default) ? (isvec_def ? default[i] : default) : zero(T) + elseif haskey(param_values, par) + coldata[i] = param_values[par] + else + if isnothing(default) + throw(KeyError(par)) + else + coldata[i] = isvec_def ? default[i] : default + end + end + end + + return partable +end """ update_partable!(partable::AbstractParameterTable, params::Vector{Symbol}, values, column) @@ -239,49 +265,62 @@ of the `partable`. """ function update_partable!( partable::ParameterTable, + column::Symbol, params::AbstractVector{Symbol}, values::AbstractVector, - column::Symbol, + default::Any = nothing, ) length(params) == length(values) || throw( ArgumentError( "The length of `params` ($(length(params))) and their `values` ($(length(values))) must be the same", ), ) - coldata = get!(() -> Vector{eltype(values)}(), partable.columns, column) - resize!(coldata, length(partable)) - params_index = Dict(zip(params, eachindex(params))) - for (i, param) in enumerate(partable.columns[:param]) - coldata[i] = param != :const ? values[params_index[param]] : zero(eltype(values)) + param_values = Dict(zip(params, values)) + if length(param_values) != length(params) + throw(ArgumentError("Duplicate parameter names in `params`")) end - return partable + update_partable!(partable, column, param_values, default) end # update estimates ------------------------------------------------------------------------- """ update_estimate!( partable::AbstractParameterTable, - sem_fit::SemFit) + fit::SemFit) -Write parameter estimates from `sem_fit` to the `:estimate` column of `partable` +Write parameter estimates from `fit` to the `:estimate` column of `partable` """ -update_estimate!(partable::AbstractParameterTable, sem_fit::SemFit) = - update_partable!(partable, params(sem_fit), sem_fit.solution, :estimate) +update_estimate!(partable::ParameterTable, fit::SemFit) = update_partable!( + partable, + :estimate, + params(fit), + fit.solution, + partable.columns[:value_fixed], +) + +# fallback method for ensemble +update_estimate!(partable::AbstractParameterTable, fit::SemFit) = + update_partable!(partable, :estimate, params(fit), fit.solution) # update starting values ------------------------------------------------------------------- """ - update_start!(partable::AbstractParameterTable, sem_fit::SemFit) + update_start!(partable::AbstractParameterTable, 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`. +Write starting values from `fit` or `start_val` to the `:estimate` column of `partable`. # Arguments - `start_val`: either a vector of starting values or a function to compute starting values from `model` - `kwargs...`: are passed to `start_val` """ -update_start!(partable::AbstractParameterTable, sem_fit::SemFit) = - update_partable!(partable, params(sem_fit), sem_fit.start_val, :start) +update_start!(partable::AbstractParameterTable, fit::SemFit) = update_partable!( + partable, + :start, + params(fit), + fit.start_val, + partable.columns[:value_fixed], +) function update_start!( partable::AbstractParameterTable, @@ -292,17 +331,17 @@ function update_start!( if !(start_val isa Vector) start_val = start_val(model; kwargs...) end - return update_partable!(partable, params(model), start_val, :start) + return update_partable!(partable, :start, params(model), start_val) end # update partable standard errors ---------------------------------------------------------- """ update_se_hessian!( partable::AbstractParameterTable, - sem_fit::SemFit; + fit::SemFit; hessian = :finitediff) -Write hessian standard errors computed for `sem_fit` to the `:se` column of `partable` +Write hessian standard errors computed for `fit` to the `:se` column of `partable` # Arguments - `hessian::Symbol`: how to compute the hessian, see [se_hessian](@ref) for more information. @@ -312,11 +351,11 @@ Write hessian standard errors computed for `sem_fit` to the `:se` column of `par """ function update_se_hessian!( partable::AbstractParameterTable, - sem_fit::SemFit; + fit::SemFit; hessian = :finitediff, ) - se = se_hessian(sem_fit; hessian = hessian) - return update_partable!(partable, params(sem_fit), se, :se) + se = se_hessian(fit; hessian = hessian) + return update_partable!(partable, :se, params(fit), se) end """ From 84945fa2dfdab324c8c0d65b8e806d45dabe91b3 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Wed, 1 May 2024 21:19:28 -0700 Subject: [PATCH 39/56] ParTable: getindex() returns NamedTuple so the downstream code doesn't rely on the order of tuple elements --- src/frontend/specification/ParameterTable.jl | 12 +++---- src/frontend/specification/RAMMatrices.jl | 36 +++++++++----------- 2 files changed, 23 insertions(+), 25 deletions(-) diff --git a/src/frontend/specification/ParameterTable.jl b/src/frontend/specification/ParameterTable.jl index 6dcbbbf84..d1590269c 100644 --- a/src/frontend/specification/ParameterTable.jl +++ b/src/frontend/specification/ParameterTable.jl @@ -126,12 +126,12 @@ ParameterTableRow = @NamedTuple begin end Base.getindex(partable::ParameterTable, i::Integer) = ( - partable.columns[:from][i], - partable.columns[:parameter_type][i], - partable.columns[:to][i], - partable.columns[:free][i], - partable.columns[:value_fixed][i], - partable.columns[:param][i], + from = partable.columns[:from][i], + parameter_type = partable.columns[:parameter_type][i], + to = partable.columns[:to][i], + free = partable.columns[:free][i], + value_fixed = partable.columns[:value_fixed][i], + param = partable.columns[:param][i], ) Base.length(partable::ParameterTable) = length(partable.columns[:param]) diff --git a/src/frontend/specification/RAMMatrices.jl b/src/frontend/specification/RAMMatrices.jl index 80ac34cf9..a900f53a9 100644 --- a/src/frontend/specification/RAMMatrices.jl +++ b/src/frontend/specification/RAMMatrices.jl @@ -171,41 +171,39 @@ function RAMMatrices( # handle constants constants = Vector{RAMConstant}() - for i in 1:length(partable) - from, parameter_type, to, free, value_fixed, param = partable[i] - - row_ind = col_indices[to] - col_ind = from != Symbol("1") ? col_indices[from] : nothing - - if !free - if (parameter_type == :→) && (from == Symbol("1")) - push!(constants, RAMConstant(:M, row_ind, value_fixed)) - elseif (parameter_type == :→) + for r in partable + row_ind = col_indices[r.to] + col_ind = r.from != Symbol("1") ? col_indices[r.from] : nothing + + if !r.free + if (r.parameter_type == :→) && (r.from == Symbol("1")) + push!(constants, RAMConstant(:M, row_ind, r.value_fixed)) + elseif r.parameter_type == :→ push!( constants, - RAMConstant(:A, CartesianIndex(row_ind, col_ind), value_fixed), + RAMConstant(:A, CartesianIndex(row_ind, col_ind), r.value_fixed), ) - elseif (parameter_type == :↔) + elseif r.parameter_type == :↔ push!( constants, - RAMConstant(:S, CartesianIndex(row_ind, col_ind), value_fixed), + RAMConstant(:S, CartesianIndex(row_ind, col_ind), r.value_fixed), ) else - error("Unsupported parameter type: $(parameter_type)") + error("Unsupported parameter type: $(r.parameter_type)") end else - par_ind = params_index[param] - if (parameter_type == :→) && (from == Symbol("1")) + par_ind = params_index[r.param] + if (r.parameter_type == :→) && (r.from == Symbol("1")) push!(M_ind[par_ind], row_ind) - elseif parameter_type == :→ + elseif r.parameter_type == :→ push!(A_ind[par_ind], row_ind + (col_ind - 1) * n_node) - elseif parameter_type == :↔ + elseif r.parameter_type == :↔ push!(S_ind[par_ind], row_ind + (col_ind - 1) * n_node) if row_ind != col_ind push!(S_ind[par_ind], col_ind + (row_ind - 1) * n_node) end else - error("Unsupported parameter type: $(parameter_type)") + error("Unsupported parameter type: $(r.parameter_type)") end end end From 60e864bb5258822d70f65fd2bca11f8c43c5c3e7 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Wed, 1 May 2024 21:42:33 -0700 Subject: [PATCH 40/56] ParTable: graph-based ctor supports params= kw --- src/frontend/specification/StenoGraphs.jl | 20 ++++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/src/frontend/specification/StenoGraphs.jl b/src/frontend/specification/StenoGraphs.jl index 42edd6a13..69b91eabf 100644 --- a/src/frontend/specification/StenoGraphs.jl +++ b/src/frontend/specification/StenoGraphs.jl @@ -33,8 +33,9 @@ label(args...) = Label(args) function ParameterTable( graph::AbstractStenoGraph; - observed_vars, - latent_vars, + observed_vars::AbstractVector{Symbol}, + latent_vars::AbstractVector{Symbol}, + params::Union{AbstractVector{Symbol}, Nothing} = nothing, group::Integer = 1, param_prefix = :θ, ) @@ -104,7 +105,7 @@ function ParameterTable( end end - return ParameterTable(columns; latent_vars, observed_vars) + return ParameterTable(columns; latent_vars, observed_vars, params) end ############################################################################################ @@ -113,8 +114,9 @@ end function EnsembleParameterTable( graph::AbstractStenoGraph; - observed_vars, - latent_vars, + observed_vars::AbstractVector{Symbol}, + latent_vars::AbstractVector{Symbol}, + params::Union{AbstractVector{Symbol}, Nothing} = nothing, groups, ) graph = unique(graph) @@ -122,11 +124,13 @@ function EnsembleParameterTable( partables = Dict( group => ParameterTable( graph; - observed_vars = observed_vars, - latent_vars = latent_vars, + observed_vars, + latent_vars, + params, group = i, param_prefix = Symbol(:g, group), ) for (i, group) in enumerate(groups) ) - return EnsembleParameterTable(partables) + + return EnsembleParameterTable(partables; params) end From 88ea6b6d4a31027beec92189d28de4cfb89a79c5 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Wed, 1 May 2024 22:36:44 -0700 Subject: [PATCH 41/56] rename parameter_type to relation for clarity --- src/frontend/fit/summary.jl | 30 ++++++-------------- src/frontend/specification/ParameterTable.jl | 14 ++++----- src/frontend/specification/RAMMatrices.jl | 23 ++++++++------- src/frontend/specification/StenoGraphs.jl | 6 ++-- 4 files changed, 31 insertions(+), 42 deletions(-) diff --git a/src/frontend/fit/summary.jl b/src/frontend/fit/summary.jl index a31d4796f..f7ecdb331 100644 --- a/src/frontend/fit/summary.jl +++ b/src/frontend/fit/summary.jl @@ -93,19 +93,13 @@ function sem_summary( sorted_columns = [:to, :estimate, :param, :value_fixed, :start] loading_columns = sort_partially(sorted_columns, columns) header_cols = copy(loading_columns) - replace!(header_cols, :parameter_type => :type) for var in partable.latent_vars indicator_indices = findall( (partable.columns[:from] .== var) .& - (partable.columns[:parameter_type] .== :→) .& + (partable.columns[:relation] .== :→) .& (partable.columns[:to] .∈ [partable.observed_vars]), ) - loading_array = reduce( - hcat, - check_round(partable.columns[c][indicator_indices]; digits = digits) for - c in loading_columns - ) printstyled(var; color = secondary_color) print("\n") @@ -122,7 +116,7 @@ function sem_summary( printstyled("Directed Effects: \n"; color = color) regression_indices = findall( - (partable.columns[:parameter_type] .== :→) .& ( + (partable.columns[:relation] .== :→) .& ( ( (partable.columns[:to] .∈ [partable.observed_vars]) .& (partable.columns[:from] .∈ [partable.observed_vars]) @@ -138,7 +132,7 @@ function sem_summary( ), ) - sorted_columns = [:from, :parameter_type, :to, :estimate, :param, :value_fixed, :start] + sorted_columns = [:from, :relation, :to, :estimate, :param, :value_fixed, :start] regression_columns = sort_partially(sorted_columns, columns) regression_array = reduce( @@ -147,7 +141,6 @@ function sem_summary( c in regression_columns ) regression_columns[2] = Symbol("") - replace!(regression_columns, :parameter_type => :type) print("\n") pretty_table( @@ -161,11 +154,11 @@ function sem_summary( printstyled("Variances: \n"; color = color) variance_indices = findall( - (partable.columns[:parameter_type] .== :↔) .& + (partable.columns[:relation] .== :↔) .& (partable.columns[:to] .== partable.columns[:from]), ) - sorted_columns = [:from, :parameter_type, :to, :estimate, :param, :value_fixed, :start] + sorted_columns = [:from, :relation, :to, :estimate, :param, :value_fixed, :start] variance_columns = sort_partially(sorted_columns, columns) variance_array = reduce( @@ -174,7 +167,6 @@ function sem_summary( c in variance_columns ) variance_columns[2] = Symbol("") - replace!(variance_columns, :parameter_type => :type) print("\n") pretty_table( @@ -188,11 +180,11 @@ function sem_summary( printstyled("Covariances: \n"; color = color) variance_indices = findall( - (partable.columns[:parameter_type] .== :↔) .& + (partable.columns[:relation] .== :↔) .& (partable.columns[:to] .!= partable.columns[:from]), ) - sorted_columns = [:from, :parameter_type, :to, :estimate, :param, :value_fixed, :start] + sorted_columns = [:from, :relation, :to, :estimate, :param, :value_fixed, :start] variance_columns = sort_partially(sorted_columns, columns) variance_array = reduce( @@ -201,7 +193,6 @@ function sem_summary( c in variance_columns ) variance_columns[2] = Symbol("") - replace!(variance_columns, :parameter_type => :type) print("\n") pretty_table( @@ -213,15 +204,13 @@ function sem_summary( print("\n") mean_indices = findall( - (partable.columns[:parameter_type] .== :→) .& - (partable.columns[:from] .== Symbol("1")), + (partable.columns[:relation] .== :→) .& (partable.columns[:from] .== Symbol("1")), ) if length(mean_indices) > 0 printstyled("Means: \n"; color = color) - sorted_columns = - [:from, :parameter_type, :to, :estimate, :param, :value_fixed, :start] + sorted_columns = [:from, :relation, :to, :estimate, :param, :value_fixed, :start] variance_columns = sort_partially(sorted_columns, columns) variance_array = reduce( @@ -230,7 +219,6 @@ function sem_summary( c in variance_columns ) variance_columns[2] = Symbol("") - replace!(variance_columns, :parameter_type => :type) print("\n") pretty_table( diff --git a/src/frontend/specification/ParameterTable.jl b/src/frontend/specification/ParameterTable.jl index d1590269c..eeb749401 100644 --- a/src/frontend/specification/ParameterTable.jl +++ b/src/frontend/specification/ParameterTable.jl @@ -18,7 +18,7 @@ end # optionally pre-allocate data for nrows empty_partable_columns(nrows::Integer = 0) = Dict{Symbol, Vector}( :from => fill(Symbol(), nrows), - :parameter_type => fill(Symbol(), nrows), + :relation => fill(Symbol(), nrows), :to => fill(Symbol(), nrows), :free => fill(true, nrows), :value_fixed => fill(NaN, nrows), @@ -93,7 +93,7 @@ end function Base.show(io::IO, partable::ParameterTable) relevant_columns = - [:from, :parameter_type, :to, :free, :value_fixed, :start, :estimate, :se, :param] + [:from, :relation, :to, :free, :value_fixed, :start, :estimate, :se, :param] shown_columns = filter!( col -> haskey(partable.columns, col) && length(partable.columns[col]) > 0, relevant_columns, @@ -118,7 +118,7 @@ end # Iteration -------------------------------------------------------------------------------- ParameterTableRow = @NamedTuple begin from::Symbol - parameter_type::Symbol + relation::Symbol to::Symbol free::Bool value_fixed::Any @@ -127,7 +127,7 @@ end Base.getindex(partable::ParameterTable, i::Integer) = ( from = partable.columns[:from][i], - parameter_type = partable.columns[:parameter_type][i], + relation = partable.columns[:relation][i], to = partable.columns[:to][i], free = partable.columns[:free][i], value_fixed = partable.columns[:value_fixed][i], @@ -170,8 +170,8 @@ function sort_vars!(partable::ParameterTable) ] is_regression = [ - (partype == :→) && (from != Symbol("1")) for (partype, from) in - zip(partable.columns[:parameter_type], partable.columns[:from]) + (rel == :→) && (from != Symbol("1")) for + (rel, from) in zip(partable.columns[:relation], partable.columns[:from]) ] to = partable.columns[:to][is_regression] @@ -453,7 +453,7 @@ function lavaan_param_values!( for (from, to, type, id) in zip( [ view(partable.columns[k], partable_mask) for - k in [:from, :to, :parameter_type, :param] + k in [:from, :to, :relation, :param] ]..., ) lav_ind = nothing diff --git a/src/frontend/specification/RAMMatrices.jl b/src/frontend/specification/RAMMatrices.jl index a900f53a9..b52c15a51 100644 --- a/src/frontend/specification/RAMMatrices.jl +++ b/src/frontend/specification/RAMMatrices.jl @@ -176,34 +176,34 @@ function RAMMatrices( col_ind = r.from != Symbol("1") ? col_indices[r.from] : nothing if !r.free - if (r.parameter_type == :→) && (r.from == Symbol("1")) + if (r.relation == :→) && (r.from == Symbol("1")) push!(constants, RAMConstant(:M, row_ind, r.value_fixed)) - elseif r.parameter_type == :→ + elseif r.relation == :→ push!( constants, RAMConstant(:A, CartesianIndex(row_ind, col_ind), r.value_fixed), ) - elseif r.parameter_type == :↔ + elseif r.relation == :↔ push!( constants, RAMConstant(:S, CartesianIndex(row_ind, col_ind), r.value_fixed), ) else - error("Unsupported parameter type: $(r.parameter_type)") + error("Unsupported parameter type: $(r.relation)") end else par_ind = params_index[r.param] - if (r.parameter_type == :→) && (r.from == Symbol("1")) + if (r.relation == :→) && (r.from == Symbol("1")) push!(M_ind[par_ind], row_ind) - elseif r.parameter_type == :→ + elseif r.relation == :→ push!(A_ind[par_ind], row_ind + (col_ind - 1) * n_node) - elseif r.parameter_type == :↔ + elseif r.relation == :↔ push!(S_ind[par_ind], row_ind + (col_ind - 1) * n_node) if row_ind != col_ind push!(S_ind[par_ind], col_ind + (row_ind - 1) * n_node) end else - error("Unsupported parameter type: $(r.parameter_type)") + error("Unsupported parameter type: $(r.relation)") end end end @@ -298,7 +298,8 @@ end ### Additional Functions ############################################################################################ -function matrix_to_parameter_type(matrix::Symbol) +# return the `from □ to` variables relation symbol (□) given the name of the source RAM matrix +function matrix_to_relation(matrix::Symbol) if matrix == :A return :→ elseif matrix == :S @@ -316,7 +317,7 @@ end partable_row(c::RAMConstant, varnames::AbstractVector{Symbol}) = ( from = varnames[c.index[2]], - parameter_type = matrix_to_parameter_type(c.matrix), + relation = matrix_to_relation(c.matrix), to = varnames[c.index[1]], free = false, value_fixed = c.value, @@ -346,7 +347,7 @@ function partable_row( return ( from = from, - parameter_type = matrix_to_parameter_type(matrix), + relation = matrix_to_relation(matrix), to = to, free = true, value_fixed = 0.0, diff --git a/src/frontend/specification/StenoGraphs.jl b/src/frontend/specification/StenoGraphs.jl index 69b91eabf..67bb7973c 100644 --- a/src/frontend/specification/StenoGraphs.jl +++ b/src/frontend/specification/StenoGraphs.jl @@ -44,7 +44,7 @@ function ParameterTable( columns = empty_partable_columns(n) from = columns[:from] - parameter_type = columns[:parameter_type] + relation = columns[:relation] to = columns[:to] free = columns[:free] value_fixed = columns[:value_fixed] @@ -57,9 +57,9 @@ function ParameterTable( from[i] = edge.src.node to[i] = edge.dst.node if edge isa DirectedEdge - parameter_type[i] = :→ + relation[i] = :→ elseif edge isa UndirectedEdge - parameter_type[i] = :↔ + relation[i] = :↔ else throw( ArgumentError( From d57c08be5f7291341bb834aee227ffbf42f66348 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Wed, 1 May 2024 22:40:59 -0700 Subject: [PATCH 42/56] sem_summary(): cleanup filters --- src/frontend/fit/summary.jl | 80 +++++++++++++++---------------------- 1 file changed, 32 insertions(+), 48 deletions(-) diff --git a/src/frontend/fit/summary.jl b/src/frontend/fit/summary.jl index f7ecdb331..214eb58f1 100644 --- a/src/frontend/fit/summary.jl +++ b/src/frontend/fit/summary.jl @@ -96,9 +96,9 @@ function sem_summary( for var in partable.latent_vars indicator_indices = findall( - (partable.columns[:from] .== var) .& - (partable.columns[:relation] .== :→) .& - (partable.columns[:to] .∈ [partable.observed_vars]), + r -> + (r.from == var) && (r.relation == :→) && (r.to ∈ partable.observed_vars), + partable, ) printstyled(var; color = secondary_color) @@ -116,20 +116,13 @@ function sem_summary( printstyled("Directed Effects: \n"; color = color) regression_indices = findall( - (partable.columns[:relation] .== :→) .& ( - ( - (partable.columns[:to] .∈ [partable.observed_vars]) .& - (partable.columns[:from] .∈ [partable.observed_vars]) - ) .| - ( - (partable.columns[:to] .∈ [partable.latent_vars]) .& - (partable.columns[:from] .∈ [partable.observed_vars]) - ) .| - ( - (partable.columns[:to] .∈ [partable.latent_vars]) .& - (partable.columns[:from] .∈ [partable.latent_vars]) - ) - ), + r -> + (r.relation == :→) && ( + ((r.to ∈ partable.observed_vars) && (r.from ∈ partable.observed_vars)) || + ((r.to ∈ partable.latent_vars) && (r.from ∈ partable.observed_vars)) || + ((r.to ∈ partable.latent_vars) && (r.from ∈ partable.latent_vars)) + ), + partable, ) sorted_columns = [:from, :relation, :to, :estimate, :param, :value_fixed, :start] @@ -153,25 +146,22 @@ function sem_summary( printstyled("Variances: \n"; color = color) - variance_indices = findall( - (partable.columns[:relation] .== :↔) .& - (partable.columns[:to] .== partable.columns[:from]), - ) + var_indices = findall(r -> r.relation == :↔ && r.to == r.from, partable) sorted_columns = [:from, :relation, :to, :estimate, :param, :value_fixed, :start] - variance_columns = sort_partially(sorted_columns, columns) + var_columns = sort_partially(sorted_columns, columns) - variance_array = reduce( + var_array = reduce( hcat, - check_round(partable.columns[c][variance_indices]; digits = digits) for + check_round(partable.columns[c][var_indices]; digits = digits) for c in variance_columns ) - variance_columns[2] = Symbol("") + var_columns[2] = Symbol("") print("\n") pretty_table( - variance_array; - header = variance_columns, + var_array; + header = var_columns, tf = PrettyTables.tf_borderless, alignment = :l, ) @@ -179,51 +169,45 @@ function sem_summary( printstyled("Covariances: \n"; color = color) - variance_indices = findall( - (partable.columns[:relation] .== :↔) .& - (partable.columns[:to] .!= partable.columns[:from]), - ) + covar_indices = findall(r -> r.relation == :↔ && r.to != r.from, partable) - sorted_columns = [:from, :relation, :to, :estimate, :param, :value_fixed, :start] - variance_columns = sort_partially(sorted_columns, columns) + covar_columns = sort_partially(sorted_columns, columns) - variance_array = reduce( + covar_array = reduce( hcat, - check_round(partable.columns[c][variance_indices]; digits = digits) for - c in variance_columns + check_round(partable.columns[c][covar_indices]; digits = digits) for + c in covar_columns ) - variance_columns[2] = Symbol("") + covar_columns[2] = Symbol("") print("\n") pretty_table( - variance_array; - header = variance_columns, + covar_array; + header = covar_columns, tf = PrettyTables.tf_borderless, alignment = :l, ) print("\n") - mean_indices = findall( - (partable.columns[:relation] .== :→) .& (partable.columns[:from] .== Symbol("1")), - ) + mean_indices = findall(r -> (r.relation == :→) && (r.from == Symbol("1")), partable) if length(mean_indices) > 0 printstyled("Means: \n"; color = color) sorted_columns = [:from, :relation, :to, :estimate, :param, :value_fixed, :start] - variance_columns = sort_partially(sorted_columns, columns) + mean_columns = sort_partially(sorted_columns, columns) - variance_array = reduce( + mean_array = reduce( hcat, check_round(partable.columns[c][mean_indices]; digits = digits) for - c in variance_columns + c in mean_columns ) - variance_columns[2] = Symbol("") + mean_columns[2] = Symbol("") print("\n") pretty_table( - variance_array; - header = variance_columns, + mean_array; + header = mean_columns, tf = PrettyTables.tf_borderless, alignment = :l, ) From 14970a2086c9c86a3c6e1c488eaf5c996b6b32a7 Mon Sep 17 00:00:00 2001 From: Maximilian Ernst Date: Sun, 19 May 2024 21:16:40 +0200 Subject: [PATCH 43/56] fix sem_summary method for partable --- src/frontend/fit/summary.jl | 20 +++++++++++++++++++- 1 file changed, 19 insertions(+), 1 deletion(-) diff --git a/src/frontend/fit/summary.jl b/src/frontend/fit/summary.jl index 214eb58f1..621791211 100644 --- a/src/frontend/fit/summary.jl +++ b/src/frontend/fit/summary.jl @@ -100,6 +100,11 @@ function sem_summary( (r.from == var) && (r.relation == :→) && (r.to ∈ partable.observed_vars), partable, ) + loading_array = reduce( + hcat, + check_round(partable.columns[c][indicator_indices]; digits = digits) for + c in loading_columns + ) printstyled(var; color = secondary_color) print("\n") @@ -109,6 +114,7 @@ function sem_summary( header = header_cols, tf = PrettyTables.tf_borderless, alignment = :l, + formatters = (v, i, j) -> isa(v, Number) && isnan(v) ? "" : v, ) print("\n") end @@ -141,6 +147,7 @@ function sem_summary( header = regression_columns, tf = PrettyTables.tf_borderless, alignment = :l, + formatters = (v, i, j) -> isa(v, Number) && isnan(v) ? "" : v, ) print("\n") @@ -154,7 +161,7 @@ function sem_summary( var_array = reduce( hcat, check_round(partable.columns[c][var_indices]; digits = digits) for - c in variance_columns + c in var_columns ) var_columns[2] = Symbol("") @@ -164,6 +171,7 @@ function sem_summary( header = var_columns, tf = PrettyTables.tf_borderless, alignment = :l, + formatters = (v, i, j) -> isa(v, Number) && isnan(v) ? "" : v, ) print("\n") @@ -186,6 +194,7 @@ function sem_summary( header = covar_columns, tf = PrettyTables.tf_borderless, alignment = :l, + formatters = (v, i, j) -> isa(v, Number) && isnan(v) ? "" : v, ) print("\n") @@ -210,6 +219,7 @@ function sem_summary( header = mean_columns, tf = PrettyTables.tf_borderless, alignment = :l, + formatters = (v, i, j) -> isa(v, Number) && isnan(v) ? "" : v, ) print("\n") end @@ -290,6 +300,14 @@ function sort_partially(sorted, to_sort) return out end +function Base.findall(fun::Function, partable::ParameterTable) + rows = Int[] + for (i, r) in enumerate(partable) + fun(r) ? push!(rows, i) : nothing + end + return rows +end + """ (1) sem_summary(sem_fit::SemFit; show_fitmeasures = false) From 4ba5fa8081b4d4aa6c71c9bfd63cf48b171b832b Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Mon, 27 May 2024 02:26:15 -0700 Subject: [PATCH 44/56] show(ParTable): suppress NaNs --- src/frontend/specification/ParameterTable.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/frontend/specification/ParameterTable.jl b/src/frontend/specification/ParameterTable.jl index eeb749401..6455f3332 100644 --- a/src/frontend/specification/ParameterTable.jl +++ b/src/frontend/specification/ParameterTable.jl @@ -105,6 +105,8 @@ function Base.show(io::IO, partable::ParameterTable) as_matrix, header = (shown_columns, [eltype(partable.columns[col]) for col in shown_columns]), tf = PrettyTables.tf_compact, + # TODO switch to `missing` as non-specified values and suppress printing of `missing` instead + formatters = (v, i, j) -> isa(v, Number) && isnan(v) ? "" : v, ) print(io, "Latent Variables: $(partable.latent_vars) \n") From 10ada011d571e7cdb779361bc9add4048c3be374 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sun, 19 May 2024 19:39:01 +0200 Subject: [PATCH 45/56] sort_vars!(ParTable): cleanup --- src/frontend/specification/ParameterTable.jl | 28 +++++++++++++------- 1 file changed, 18 insertions(+), 10 deletions(-) diff --git a/src/frontend/specification/ParameterTable.jl b/src/frontend/specification/ParameterTable.jl index 6455f3332..7f4eab486 100644 --- a/src/frontend/specification/ParameterTable.jl +++ b/src/frontend/specification/ParameterTable.jl @@ -171,13 +171,15 @@ function sort_vars!(partable::ParameterTable) partable.observed_vars ] - is_regression = [ - (rel == :→) && (from != Symbol("1")) for - (rel, from) in zip(partable.columns[:relation], partable.columns[:from]) + # regression edges (excluding intercept) + edges = [ + (from, to) for (rel, from, to) in zip( + partable.columns[:relation], + partable.columns[:from], + partable.columns[:to], + ) if (rel == :→) && (from != Symbol("1")) ] - - to = partable.columns[:to][is_regression] - from = partable.columns[:from][is_regression] + sort!(edges, by = last) # sort edges by target sorted_vars = Vector{Symbol}() @@ -185,21 +187,27 @@ function sort_vars!(partable::ParameterTable) acyclic = false for (i, var) in enumerate(vars) - if !(var ∈ to) + # check if var has any incoming edge + eix = searchsortedfirst(edges, (var, var), by = last) + if !(eix <= length(edges) && last(edges[eix]) == var) + # var is source, no edges to it push!(sorted_vars, var) deleteat!(vars, i) - delete_edges = from .!= var - to = to[delete_edges] - from = from[delete_edges] + # remove var outgoing edges + filter!(e -> e[1] != var, edges) acyclic = true + break end end + # if acyclic is false, all vars have incoming edge acyclic || throw(CyclicModelError("your model is cyclic and therefore can not be ordered")) end copyto!(resize!(partable.sorted_vars, length(sorted_vars)), sorted_vars) + @assert length(partable.sorted_vars) == + length(partable.observed_vars) + length(partable.latent_vars) return partable end From 42da8a8bd35b1c1f79c15c5d7cd60547e29dbfdc Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sat, 11 May 2024 13:21:57 -0700 Subject: [PATCH 46/56] Project.toml: disable SymbolicUtils 1.6 causes problems with sparsehessian(). It is a temporary fix until the compatibility issues are resolved in Symbolics.jl --- Project.toml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/Project.toml b/Project.toml index 2fa168a9e..3a44943a6 100644 --- a/Project.toml +++ b/Project.toml @@ -22,6 +22,7 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" StenoGraphs = "78862bba-adae-4a83-bb4d-33c106177f81" Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" +SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b" [compat] julia = "1.9, 1.10" @@ -36,6 +37,7 @@ Optim = "1" PrettyTables = "2" StatsBase = "0.33, 0.34" Symbolics = "4, 5" +SymbolicUtils = "1.4 - 1.5" [extras] Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" From 7892b6d6b19dd99c8033ebec706de418373fef45 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sat, 11 May 2024 13:22:15 -0700 Subject: [PATCH 47/56] Project.toml: support StenoGraphs 0.3 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 3a44943a6..9edeb4536 100644 --- a/Project.toml +++ b/Project.toml @@ -26,7 +26,7 @@ SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b" [compat] julia = "1.9, 1.10" -StenoGraphs = "0.2" +StenoGraphs = "0.2, 0.3" DataFrames = "1" Distributions = "0.25" FiniteDiff = "2" From 63ae853d34a9ded04be2f866b86f78f69f8457fb Mon Sep 17 00:00:00 2001 From: Maximilian Ernst Date: Sun, 19 May 2024 22:15:15 +0200 Subject: [PATCH 48/56] RAM ctor: better error for missing meanstruct --- src/imply/RAM/generic.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/imply/RAM/generic.jl b/src/imply/RAM/generic.jl index 9eb694d51..bc81a71a0 100644 --- a/src/imply/RAM/generic.jl +++ b/src/imply/RAM/generic.jl @@ -143,6 +143,7 @@ function RAM(; # μ if meanstructure has_meanstructure = Val(true) + !isnothing(M_indices) || throw(ArgumentError("You set `meanstructure = true`, but your model specification contains no mean parameters.")) ∇M = gradient ? matrix_gradient(M_indices, n_nod) : nothing μ = zeros(n_var) else From 7a39a0dc249843abc6603fcca743133153e46172 Mon Sep 17 00:00:00 2001 From: Maximilian Ernst Date: Sun, 19 May 2024 23:43:54 +0200 Subject: [PATCH 49/56] add function param_indices --- src/StructuralEquationModels.jl | 1 + src/types.jl | 17 +++++++++++++++++ 2 files changed, 18 insertions(+) diff --git a/src/StructuralEquationModels.jl b/src/StructuralEquationModels.jl index 7812fa819..19c7653bc 100644 --- a/src/StructuralEquationModels.jl +++ b/src/StructuralEquationModels.jl @@ -155,6 +155,7 @@ export AbstractSem, RAMMatrices, params, nparams, + param_indices, fit_measures, AIC, BIC, diff --git a/src/types.jl b/src/types.jl index 9db48b6db..0be745e80 100644 --- a/src/types.jl +++ b/src/types.jl @@ -29,6 +29,23 @@ nparams(model::AbstractSem) = length(params(model)) params(model::AbstractSemSingle) = params(model.imply) nparams(model::AbstractSemSingle) = nparams(model.imply) +""" + param_indices(semobj) + param_indices(param_names, semobj) + +Returns either a dict of parameter names and their indices in `semobj`. +If `param_names` are provided, returns a vector their indices in `semobj` instead. + +# Examples +```julia +parind = param_indices(my_fitted_sem) +parind[:param_name] + +parind = param_indices([:param_name_1, param_name_2], my_fitted_sem) +``` +""" +param_indices(semobj) = Dict(params(semobj) .=> 1:nparams(semobj)) +param_indices(param_names, semobj) = getindex.([Dict(params(semobj) .=> 1:nparams(semobj))], param_names) """ SemLoss(args...; loss_weights = nothing, ...) From 79af8105eeae66519a060124207e9e9249cdc9b7 Mon Sep 17 00:00:00 2001 From: Maximilian Ernst Date: Sun, 19 May 2024 23:45:48 +0200 Subject: [PATCH 50/56] start fixing docs --- docs/src/developer/loss.md | 9 ++--- docs/src/developer/sem.md | 3 +- docs/src/performance/simulation.md | 5 +-- docs/src/performance/sorting.md | 2 +- docs/src/tutorials/collection/multigroup.md | 6 ++-- docs/src/tutorials/constraints/constraints.md | 35 +++++++++++-------- .../tutorials/construction/build_by_parts.md | 4 +-- .../construction/outer_constructor.md | 2 +- docs/src/tutorials/first_model.md | 4 +-- docs/src/tutorials/inspection/inspection.md | 10 +++--- docs/src/tutorials/meanstructure.md | 8 ++--- .../regularization/regularization.md | 7 ++-- .../specification/graph_interface.md | 8 ++--- .../tutorials/specification/ram_matrices.md | 4 +-- .../tutorials/specification/specification.md | 4 +-- 15 files changed, 61 insertions(+), 50 deletions(-) diff --git a/docs/src/developer/loss.md b/docs/src/developer/loss.md index 8bd654bf1..e1137dbf1 100644 --- a/docs/src/developer/loss.md +++ b/docs/src/developer/loss.md @@ -60,11 +60,12 @@ graph = @StenoGraph begin end partable = ParameterTable( + graph, latent_vars = latent_vars, - observed_vars = observed_vars, - graph = graph) + observed_vars = observed_vars +) -parameter_indices = get_identifier_indices([:a, :b, :c], partable) +parameter_indices = param_indices([:a, :b, :c], partable) myridge = Ridge(0.01, parameter_indices) model = SemFiniteDiff( @@ -269,4 +270,4 @@ model_ml = SemFiniteDiff( model_fit = sem_fit(model_ml) ``` -If you want to differentiate your own loss functions via automatic differentiation, check out the [AutoDiffSEM](https://github.com/StructuralEquationModels/AutoDiffSEM) package (spoiler allert: it's really easy). +If you want to differentiate your own loss functions via automatic differentiation, check out the [AutoDiffSEM](https://github.com/StructuralEquationModels/AutoDiffSEM) package. diff --git a/docs/src/developer/sem.md b/docs/src/developer/sem.md index c6b9f0523..528da88b8 100644 --- a/docs/src/developer/sem.md +++ b/docs/src/developer/sem.md @@ -11,7 +11,8 @@ struct SemFiniteDiff{ observed::O imply::I loss::L - optimizer::Dend + optimizer::D +end ``` Additionally, we need to define a method to compute at least the objective value, and if you want to use gradient based optimizers (which you most probably will), we need also to define a method to compute the gradient. For example, the respective fallback methods for all `AbstractSemSingle` models are defined as diff --git a/docs/src/performance/simulation.md b/docs/src/performance/simulation.md index 4b00df6a4..b8a5081fe 100644 --- a/docs/src/performance/simulation.md +++ b/docs/src/performance/simulation.md @@ -43,9 +43,10 @@ graph = @StenoGraph begin end partable = ParameterTable( + graph, latent_vars = latent_vars, - observed_vars = observed_vars, - graph = graph) + observed_vars = observed_vars +) ``` ```@example swap_observed diff --git a/docs/src/performance/sorting.md b/docs/src/performance/sorting.md index 802720099..78fd09411 100644 --- a/docs/src/performance/sorting.md +++ b/docs/src/performance/sorting.md @@ -13,7 +13,7 @@ To automatically reorder your variables in a way that makes this optimization po We use it as ```julia -sort!(parameter_table) +sort_vars!(parameter_table) model = Sem( specification = parameter_table, diff --git a/docs/src/tutorials/collection/multigroup.md b/docs/src/tutorials/collection/multigroup.md index 4e6105128..399d89760 100644 --- a/docs/src/tutorials/collection/multigroup.md +++ b/docs/src/tutorials/collection/multigroup.md @@ -61,8 +61,8 @@ You can then use the resulting graph to specify an `EnsembleParameterTable` ```@example mg; ansicolor = true groups = [:Pasteur, :Grant_White] -partable = EnsembleParameterTable(; - graph = graph, +partable = EnsembleParameterTable( + graph, observed_vars = observed_vars, latent_vars = latent_vars, groups = groups) @@ -71,7 +71,7 @@ partable = EnsembleParameterTable(; The parameter table can be used to create a `Dict` of RAMMatrices with keys equal to the group names and parameter tables as values: ```@example mg; ansicolor = true -specification = RAMMatrices(partable) +specification = convert(Dict{Symbol, RAMMatrices}, partable) ``` That is, you can asses the group-specific `RAMMatrices` as `specification[:group_name]`. diff --git a/docs/src/tutorials/constraints/constraints.md b/docs/src/tutorials/constraints/constraints.md index 7e2ec53e1..a67ad7372 100644 --- a/docs/src/tutorials/constraints/constraints.md +++ b/docs/src/tutorials/constraints/constraints.md @@ -16,13 +16,13 @@ graph = @StenoGraph begin # loadings ind60 → fixed(1)*x1 + x2 + x3 - dem60 → fixed(1)*y1 + y2 + y3 + y4 + dem60 → fixed(1)*y1 + label(:λ₂)*y2 + label(:λ₃)*y3 + y4 dem65 → fixed(1)*y5 + y6 + y7 + y8 # latent regressions ind60 → dem60 dem60 → dem65 - ind60 → dem65 + ind60 → label(:λₗ)*dem65 # variances _(observed_vars) ↔ _(observed_vars) @@ -31,15 +31,15 @@ graph = @StenoGraph begin # covariances y1 ↔ y5 y2 ↔ y4 + y6 - y3 ↔ y7 - y8 ↔ y4 + y6 + y3 ↔ label(:y3y7)*y7 + y8 ↔ label(:y8y4)*y4 + y6 end partable = ParameterTable( + graph, latent_vars = latent_vars, - observed_vars = observed_vars, - graph = graph) + observed_vars = observed_vars) data = example_data("political_democracy") @@ -64,17 +64,19 @@ Let's introduce some constraints: (Of course those constaints only serve an illustratory purpose.) -We first need to get the indices of the respective parameters that are invoved in the constraints. We can look up their labels in the output above, and retrieve their indices as +We first need to get the indices of the respective parameters that are invoved in the constraints. +We can look up their labels in the output above, and retrieve their indices as ```@example constraints -parameter_indices = get_identifier_indices([:θ_29, :θ_30, :θ_3, :θ_4, :θ_11], model) +parind = param_indices(model) +parind[:y3y7] # 29 ``` -The bound constraint is easy to specify: Just give a vector of upper or lower bounds that contains the bound for each parameter. In our example, only parameter number 11 has an upper bound, and the number of total parameters is `n_par(model) = 31`, so we define +The bound constraint is easy to specify: Just give a vector of upper or lower bounds that contains the bound for each parameter. In our example, only the parameter labeled `:λₗ` has an upper bound, and the number of total parameters is `n_par(model) = 31`, so we define ```@example constraints upper_bounds = fill(Inf, 31) -upper_bounds[11] = 0.5 +upper_bounds[parind[:λₗ]] = 0.5 ``` The equailty and inequality constraints have to be reformulated to be of the form `x = 0` or `x ≤ 0`: @@ -84,6 +86,8 @@ The equailty and inequality constraints have to be reformulated to be of the for Now they can be defined as functions of the parameter vector: ```@example constraints +parind[:y3y7] # 29 +parind[:y8y4] # 30 # θ[29] + θ[30] - 1 = 0.0 function eq_constraint(θ, gradient) if length(gradient) > 0 @@ -94,6 +98,8 @@ function eq_constraint(θ, gradient) return θ[29] + θ[30] - 1 end +parind[:λ₂] # 3 +parind[:λ₃] # 4 # θ[3] - θ[4] - 0.1 ≤ 0 function ineq_constraint(θ, gradient) if length(gradient) > 0 @@ -109,7 +115,7 @@ If the algorithm needs gradients at an iteration, it will pass the vector `gradi With `if length(gradient) > 0` we check if the algorithm needs gradients, and if it does, we fill the `gradient` vector with the gradients of the constraint w.r.t. the parameters. -In NLopt, vector-valued constraints are also possible, but we refer to the documentation fot that. +In NLopt, vector-valued constraints are also possible, but we refer to the documentation for that. ### Fit the model @@ -153,10 +159,11 @@ As you can see, the optimizer converged (`:XTOL_REACHED`) and investigating the ```@example constraints update_partable!( - partable, - model_fit_constrained, + partable, + :estimate_constr, + params(model_fit_constrained), solution(model_fit_constrained), - :estimate_constr) + ) sem_summary(partable) ``` diff --git a/docs/src/tutorials/construction/build_by_parts.md b/docs/src/tutorials/construction/build_by_parts.md index 5a56f1ccf..779949d98 100644 --- a/docs/src/tutorials/construction/build_by_parts.md +++ b/docs/src/tutorials/construction/build_by_parts.md @@ -39,9 +39,9 @@ graph = @StenoGraph begin end partable = ParameterTable( + graph, latent_vars = latent_vars, - observed_vars = observed_vars, - graph = graph) + observed_vars = observed_vars) ``` Now, we construct the different parts: diff --git a/docs/src/tutorials/construction/outer_constructor.md b/docs/src/tutorials/construction/outer_constructor.md index 21f6bfd3f..f072b80bc 100644 --- a/docs/src/tutorials/construction/outer_constructor.md +++ b/docs/src/tutorials/construction/outer_constructor.md @@ -74,7 +74,7 @@ model = Sem( specification = partable, data = data, imply = RAMSymbolic, - loss = SemWLS + loss = SemWLS, wls_weight_matrix = W ) diff --git a/docs/src/tutorials/first_model.md b/docs/src/tutorials/first_model.md index b19a22200..7568a5917 100644 --- a/docs/src/tutorials/first_model.md +++ b/docs/src/tutorials/first_model.md @@ -83,9 +83,9 @@ We then use this graph to define a `ParameterTable` object ```@example high_level; ansicolor = true partable = ParameterTable( + graph, latent_vars = latent_vars, - observed_vars = observed_vars, - graph = graph) + observed_vars = observed_vars) ``` load the example data diff --git a/docs/src/tutorials/inspection/inspection.md b/docs/src/tutorials/inspection/inspection.md index 5bc7946ba..b2eefadb2 100644 --- a/docs/src/tutorials/inspection/inspection.md +++ b/docs/src/tutorials/inspection/inspection.md @@ -31,9 +31,9 @@ graph = @StenoGraph begin end partable = ParameterTable( + graph, latent_vars = latent_vars, - observed_vars = observed_vars, - graph = graph) + observed_vars = observed_vars) data = example_data("political_democracy") @@ -87,8 +87,8 @@ We can also update the `ParameterTable` object with other information via [`upda se_bs = se_bootstrap(model_fit; n_boot = 20) se_he = se_hessian(model_fit) -update_partable!(partable, model_fit, se_he, :se_hessian) -update_partable!(partable, model_fit, se_bs, :se_bootstrap) +update_partable!(partable, :se_hessian, params(model_fit), se_he) +update_partable!(partable, :se_bootstrap, params(model_fit), se_bs) sem_summary(partable) ``` @@ -130,7 +130,7 @@ df minus2ll n_man n_obs -n_par +nparams p_value RMSEA ``` diff --git a/docs/src/tutorials/meanstructure.md b/docs/src/tutorials/meanstructure.md index 9f2c167df..c6ad692b6 100644 --- a/docs/src/tutorials/meanstructure.md +++ b/docs/src/tutorials/meanstructure.md @@ -39,9 +39,9 @@ graph = @StenoGraph begin end partable = ParameterTable( + graph, latent_vars = latent_vars, - observed_vars = observed_vars, - graph = graph) + observed_vars = observed_vars) ``` ```julia @@ -77,9 +77,9 @@ graph = @StenoGraph begin end partable = ParameterTable( + graph, latent_vars = latent_vars, - observed_vars = observed_vars, - graph = graph) + observed_vars = observed_vars) ``` that is, all observed variable means are estimated freely. diff --git a/docs/src/tutorials/regularization/regularization.md b/docs/src/tutorials/regularization/regularization.md index b7d9affab..ce554a91d 100644 --- a/docs/src/tutorials/regularization/regularization.md +++ b/docs/src/tutorials/regularization/regularization.md @@ -86,9 +86,10 @@ graph = @StenoGraph begin end partable = ParameterTable( + graph, latent_vars = latent_vars, - observed_vars = observed_vars, - graph = graph) + observed_vars = observed_vars +) data = example_data("political_democracy") @@ -101,7 +102,7 @@ model = Sem( We labeled the covariances between the items because we want to regularize those: ```@example reg -ind = get_identifier_indices([:cov_15, :cov_24, :cov_26, :cov_37, :cov_48, :cov_68], model) +ind = param_indices([:cov_15, :cov_24, :cov_26, :cov_37, :cov_48, :cov_68], model) ``` In the following, we fit the same model with lasso regularization of those covariances. diff --git a/docs/src/tutorials/specification/graph_interface.md b/docs/src/tutorials/specification/graph_interface.md index 7a03083c7..609c844c3 100644 --- a/docs/src/tutorials/specification/graph_interface.md +++ b/docs/src/tutorials/specification/graph_interface.md @@ -16,9 +16,9 @@ observed_vars = ... latent_vars = ... partable = ParameterTable( + graph, latent_vars = latent_vars, - observed_vars = observed_vars, - graph = graph) + observed_vars = observed_vars) model = Sem( specification = partable, @@ -65,9 +65,9 @@ As you saw above and in the [A first model](@ref) example, the graph object need ```julia partable = ParameterTable( + graph, latent_vars = latent_vars, - observed_vars = observed_vars, - graph = graph) + observed_vars = observed_vars) ``` The `ParameterTable` constructor also needs you to specify a vector of observed and latent variables, in the example above this would correspond to diff --git a/docs/src/tutorials/specification/ram_matrices.md b/docs/src/tutorials/specification/ram_matrices.md index 8eea6967c..5f0757238 100644 --- a/docs/src/tutorials/specification/ram_matrices.md +++ b/docs/src/tutorials/specification/ram_matrices.md @@ -59,7 +59,7 @@ spec = RAMMatrices(; A = A, S = S, F = F, - parameters = θ, + params = θ, colnames = [:x1, :x2, :x3, :y1, :y2, :y3, :y4, :y5, :y6, :y7, :y8, :ind60, :dem60, :dem65] ) @@ -90,7 +90,7 @@ spec = RAMMatrices(; A = A, S = S, F = F, - parameters = θ, + params = θ, colnames = [:x1, :x2, :x3, :y1, :y2, :y3, :y4, :y5, :y6, :y7, :y8, :ind60, :dem60, :dem65] ) ``` diff --git a/docs/src/tutorials/specification/specification.md b/docs/src/tutorials/specification/specification.md index 88f19ce3d..c426443f4 100644 --- a/docs/src/tutorials/specification/specification.md +++ b/docs/src/tutorials/specification/specification.md @@ -18,9 +18,9 @@ graph = @StenoGraph begin end partable = ParameterTable( + graph, latent_vars = latent_vars, - observed_vars = observed_vars, - graph = graph) + observed_vars = observed_vars) model = Sem( specification = partable, From 7ba872abe81795fb0f092e5a3ed0260e384c52ee Mon Sep 17 00:00:00 2001 From: Maximilian Ernst Date: Mon, 20 May 2024 00:04:28 +0200 Subject: [PATCH 51/56] fix regularization docs --- docs/src/tutorials/regularization/regularization.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/src/tutorials/regularization/regularization.md b/docs/src/tutorials/regularization/regularization.md index ce554a91d..4aaff1d0a 100644 --- a/docs/src/tutorials/regularization/regularization.md +++ b/docs/src/tutorials/regularization/regularization.md @@ -146,7 +146,7 @@ fit = sem_fit(model) update_estimate!(partable, fit) -update_partable!(partable, fit_lasso, solution(fit_lasso), :estimate_lasso) +update_partable!(partable, :estimate_lasso, params(fit_lasso), solution(fit_lasso)) sem_summary(partable) ``` @@ -180,7 +180,7 @@ fit_mixed = sem_fit(model_mixed) Let's again compare the different results: ```@example reg -update_partable!(partable, fit_mixed, solution(fit_mixed), :estimate_mixed) +update_partable!(partable, :estimate_mixed, params(fit_mixed), solution(fit_mixed)) sem_summary(partable) ``` \ No newline at end of file From c3ee769fa9f12808f27806d90f3030b4e7760aab Mon Sep 17 00:00:00 2001 From: Maximilian Ernst Date: Mon, 20 May 2024 00:07:06 +0200 Subject: [PATCH 52/56] introduce formatting error --- src/additional_functions/helper.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/additional_functions/helper.jl b/src/additional_functions/helper.jl index 3e614e57b..138ae431e 100644 --- a/src/additional_functions/helper.jl +++ b/src/additional_functions/helper.jl @@ -41,7 +41,7 @@ function get_observed(rowind, data, semobserved; args = (), kwargs = NamedTuple( return observed_vec end -skipmissing_mean(mat::AbstractMatrix) = +skipmissing_mean(mat::AbstractMatrix) = [mean(skipmissing(coldata)) for coldata in eachcol(mat)] function F_one_person(imp_mean, meandiff, inverse, data, logdet) From 5bb89a0400440697aa94ace4e54e91c7f8287b6a Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sun, 26 May 2024 21:25:20 -0700 Subject: [PATCH 53/56] update_start(): fix docstring typo Co-authored-by: Maximilian-Stefan-Ernst <34346372+Maximilian-Stefan-Ernst@users.noreply.github.com> --- src/frontend/specification/ParameterTable.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/frontend/specification/ParameterTable.jl b/src/frontend/specification/ParameterTable.jl index 7f4eab486..4f1bf747d 100644 --- a/src/frontend/specification/ParameterTable.jl +++ b/src/frontend/specification/ParameterTable.jl @@ -317,7 +317,7 @@ update_estimate!(partable::AbstractParameterTable, fit::SemFit) = update_start!(partable::AbstractParameterTable, fit::SemFit) update_start!(partable::AbstractParameterTable, model::AbstractSem, start_val; kwargs...) -Write starting values from `fit` or `start_val` to the `:estimate` column of `partable`. +Write starting values from `fit` or `start_val` to the `:start` column of `partable`. # Arguments - `start_val`: either a vector of starting values or a function to compute starting values From 43dceff22f33d835296c7780e212d0bd7f2b2e26 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sun, 26 May 2024 21:28:14 -0700 Subject: [PATCH 54/56] push!(::ParTable, Tuple): check keys compat Co-authored-by: Maximilian-Stefan-Ernst <34346372+Maximilian-Stefan-Ernst@users.noreply.github.com> --- src/frontend/specification/ParameterTable.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/frontend/specification/ParameterTable.jl b/src/frontend/specification/ParameterTable.jl index 4f1bf747d..03b111e87 100644 --- a/src/frontend/specification/ParameterTable.jl +++ b/src/frontend/specification/ParameterTable.jl @@ -227,6 +227,8 @@ sort_vars(partable::ParameterTable) = sort_vars!(deepcopy(partable)) # add a row -------------------------------------------------------------------------------- function Base.push!(partable::ParameterTable, d::Union{AbstractDict{Symbol}, NamedTuple}) + issetequal(keys(partable.columns), keys(d)) || + throw(ArgumentError("The new row needs to have the same keys as the columns of the parameter table.")) for (key, val) in pairs(d) push!(partable.columns[key], val) end From 165474ca7a1c7f6a4616d03d9e61b4b825b0e551 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Mon, 27 May 2024 14:38:02 -0700 Subject: [PATCH 55/56] SemObsCov ctor: restrict n_obs to integer don't allow missing n_obs --- src/observed/covariance.jl | 4 ++-- test/unit_tests/data_input_formats.jl | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/observed/covariance.jl b/src/observed/covariance.jl index 8d73b1a99..9be35e510 100644 --- a/src/observed/covariance.jl +++ b/src/observed/covariance.jl @@ -47,13 +47,13 @@ struct SemObservedCovariance{B, C} <: SemObserved end function SemObservedCovariance(; - specification::Union{SemSpecification, Nothing}, + specification::Union{SemSpecification, Nothing} = nothing, obs_cov, obs_colnames = nothing, spec_colnames = nothing, obs_mean = nothing, meanstructure = false, - n_obs = nothing, + n_obs::Integer, kwargs..., ) if !meanstructure & !isnothing(obs_mean) diff --git a/test/unit_tests/data_input_formats.jl b/test/unit_tests/data_input_formats.jl index 7a048b280..44656c331 100644 --- a/test/unit_tests/data_input_formats.jl +++ b/test/unit_tests/data_input_formats.jl @@ -209,7 +209,7 @@ end ) end -@test_throws UndefKeywordError(:specification) SemObservedCovariance(obs_cov = dat_cov) +@test_throws UndefKeywordError(:n_obs) SemObservedCovariance(obs_cov = dat_cov) @test_throws ArgumentError("no `obs_colnames` were specified") begin SemObservedCovariance(specification = spec, obs_cov = dat_cov, n_obs = 75) From b2012b09f5d4b2bac0b69381eb1d8878e25924fa Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Mon, 27 May 2024 14:04:48 -0700 Subject: [PATCH 56/56] fixup param_indices() --- src/types.jl | 10 +++------- 1 file changed, 3 insertions(+), 7 deletions(-) diff --git a/src/types.jl b/src/types.jl index 0be745e80..d194709b8 100644 --- a/src/types.jl +++ b/src/types.jl @@ -29,23 +29,19 @@ nparams(model::AbstractSem) = length(params(model)) params(model::AbstractSemSingle) = params(model.imply) nparams(model::AbstractSemSingle) = nparams(model.imply) + """ param_indices(semobj) - param_indices(param_names, semobj) -Returns either a dict of parameter names and their indices in `semobj`. -If `param_names` are provided, returns a vector their indices in `semobj` instead. +Returns a dict of parameter names and their indices in `semobj`. # Examples ```julia parind = param_indices(my_fitted_sem) parind[:param_name] - -parind = param_indices([:param_name_1, param_name_2], my_fitted_sem) ``` """ -param_indices(semobj) = Dict(params(semobj) .=> 1:nparams(semobj)) -param_indices(param_names, semobj) = getindex.([Dict(params(semobj) .=> 1:nparams(semobj))], param_names) +param_indices(semobj) = Dict(par => i for (i, par) in enumerate(params(semobj))) """ SemLoss(args...; loss_weights = nothing, ...)