diff --git a/Project.toml b/Project.toml index 2fa168a9e..9edeb4536 100644 --- a/Project.toml +++ b/Project.toml @@ -22,10 +22,11 @@ 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" -StenoGraphs = "0.2" +StenoGraphs = "0.2, 0.3" DataFrames = "1" Distributions = "0.25" FiniteDiff = "2" @@ -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" 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/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..4aaff1d0a 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. @@ -145,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) ``` @@ -179,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 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, diff --git a/src/StructuralEquationModels.jl b/src/StructuralEquationModels.jl index 113022960..19c7653bc 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 @@ -31,9 +30,10 @@ 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") +include("frontend/specification/EnsembleParameterTable.jl") include("frontend/specification/StenoGraphs.jl") include("frontend/fit/summary.jl") # pretty printing @@ -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,9 +150,12 @@ export AbstractSem, start, Label, label, - get_identifier_indices, + sort_vars!, + sort_vars, RAMMatrices, - identifier, + params, + nparams, + param_indices, fit_measures, AIC, BIC, @@ -163,7 +163,6 @@ export AbstractSem, df, fit_measures, minus2ll, - n_par, n_obs, p_value, RMSEA, diff --git a/src/additional_functions/helper.jl b/src/additional_functions/helper.jl index b96813dc3..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) @@ -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 diff --git a/src/additional_functions/identifier.jl b/src/additional_functions/identifier.jl deleted file mode 100644 index fefcc1be5..000000000 --- a/src/additional_functions/identifier.jl +++ /dev/null @@ -1,59 +0,0 @@ -############################################################################################ -# get parameter identifier -############################################################################################ - -identifier(sem_fit::SemFit) = identifier(sem_fit.model) -identifier(model::AbstractSemSingle) = identifier(model.imply) -identifier(model::SemEnsemble) = model.identifier - -############################################################################################ -# construct identifier -############################################################################################ - -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 -end - -############################################################################################ -# get indices of a Vector of parameter labels -############################################################################################ - -get_identifier_indices(parameters, identifier::Dict{Symbol, Int}) = - [identifier[par] for par in parameters] - -get_identifier_indices( - parameters, - obj::Union{SemFit, AbstractSemSingle, SemEnsemble, SemImply}, -) = get_identifier_indices(parameters, identifier(obj)) - -function get_identifier_indices(parameters, 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)) -end - -############################################################################################ -# documentation -############################################################################################ -""" - get_identifier_indices(parameters, model) - -Returns the indices of `parameters`. - -# Arguments -- `parameters::Vector{Symbol}`: parameter labels -- `model`: either a SEM or a fitted SEM - -# Examples -```julia -parameter_indices = get_identifier_indices([:λ₁, λ₂], my_fitted_sem) - -values = solution(my_fitted_sem)[parameter_indices] -``` -""" -function get_identifier_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/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/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/SemFit.jl b/src/frontend/fit/SemFit.jl index 97cd9c5a6..ace9ed320 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) +nparams(fit::SemFit) = nparams(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/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/fitmeasures/n_par.jl b/src/frontend/fit/fitmeasures/n_par.jl deleted file mode 100644 index 9cb2d3479..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(identifier::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.identifier) - -n_par(identifier::Dict) = length(identifier) diff --git a/src/frontend/fit/summary.jl b/src/frontend/fit/summary.jl index 4cda902d7..621791211 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( @@ -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") @@ -91,16 +90,15 @@ 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) - 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]]), + r -> + (r.from == var) && (r.relation == :→) && (r.to ∈ partable.observed_vars), + partable, ) loading_array = reduce( hcat, @@ -116,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 @@ -123,24 +122,16 @@ function sem_summary( printstyled("Directed Effects: \n"; color = color) 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.variables[:latent_vars]]) .& - (partable.columns[:from] .∈ [partable.variables[:observed_vars]]) - ) .| - ( - (partable.columns[:to] .∈ [partable.variables[:latent_vars]]) .& - (partable.columns[:from] .∈ [partable.variables[: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, :parameter_type, :to, :estimate, :identifier, :value_fixed, :start] + sorted_columns = [:from, :relation, :to, :estimate, :param, :value_fixed, :start] regression_columns = sort_partially(sorted_columns, columns) regression_array = reduce( @@ -149,7 +140,6 @@ function sem_summary( c in regression_columns ) regression_columns[2] = Symbol("") - replace!(regression_columns, :parameter_type => :type) print("\n") pretty_table( @@ -157,91 +147,79 @@ function sem_summary( header = regression_columns, tf = PrettyTables.tf_borderless, alignment = :l, + formatters = (v, i, j) -> isa(v, Number) && isnan(v) ? "" : v, ) print("\n") printstyled("Variances: \n"; color = color) - variance_indices = findall( - (partable.columns[:parameter_type] .== :↔) .& - (partable.columns[:to] .== partable.columns[:from]), - ) + var_indices = findall(r -> r.relation == :↔ && r.to == r.from, partable) - sorted_columns = - [:from, :parameter_type, :to, :estimate, :identifier, :value_fixed, :start] - variance_columns = sort_partially(sorted_columns, columns) + sorted_columns = [:from, :relation, :to, :estimate, :param, :value_fixed, :start] + var_columns = sort_partially(sorted_columns, columns) - variance_array = reduce( + var_array = reduce( hcat, - check_round(partable.columns[c][variance_indices]; digits = digits) for - c in variance_columns + check_round(partable.columns[c][var_indices]; digits = digits) for + c in var_columns ) - variance_columns[2] = Symbol("") - replace!(variance_columns, :parameter_type => :type) + var_columns[2] = Symbol("") print("\n") pretty_table( - variance_array; - header = variance_columns, + var_array; + header = var_columns, tf = PrettyTables.tf_borderless, alignment = :l, + formatters = (v, i, j) -> isa(v, Number) && isnan(v) ? "" : v, ) print("\n") printstyled("Covariances: \n"; color = color) - variance_indices = findall( - (partable.columns[:parameter_type] .== :↔) .& - (partable.columns[:to] .!= partable.columns[:from]), - ) + covar_indices = findall(r -> r.relation == :↔ && r.to != r.from, partable) - sorted_columns = - [:from, :parameter_type, :to, :estimate, :identifier, :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("") - replace!(variance_columns, :parameter_type => :type) + covar_columns[2] = Symbol("") print("\n") pretty_table( - variance_array; - header = variance_columns, + covar_array; + header = covar_columns, tf = PrettyTables.tf_borderless, alignment = :l, + formatters = (v, i, j) -> isa(v, Number) && isnan(v) ? "" : v, ) print("\n") - mean_indices = findall( - (partable.columns[:parameter_type] .== :→) .& - (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, :parameter_type, :to, :estimate, :identifier, :value_fixed, :start] - variance_columns = sort_partially(sorted_columns, columns) + sorted_columns = [:from, :relation, :to, :estimate, :param, :value_fixed, :start] + 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("") - replace!(variance_columns, :parameter_type => :type) + mean_columns[2] = Symbol("") print("\n") pretty_table( - variance_array; - header = variance_columns, + mean_array; + header = mean_columns, tf = PrettyTables.tf_borderless, alignment = :l, + formatters = (v, i, j) -> isa(v, Number) && isnan(v) ? "" : v, ) print("\n") end @@ -266,19 +244,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") @@ -323,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) diff --git a/src/frontend/specification/EnsembleParameterTable.jl b/src/frontend/specification/EnsembleParameterTable.jl index 79283953f..b0b50448b 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,41 +12,63 @@ 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} + params = if isnothing(params) + # collect all SEM parameters in ensemble if not specified + # and apply the set to all partables + unique(mapreduce(SEM.params, vcat, values(spec_ensemble), init = Vector{Symbol}())) + else + 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 ############################################################################################ ### 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( - 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 =# - -############################################################################################ -### get parameter table from RAMMatrices -############################################################################################ - -function EnsembleParameterTable(args...; groups) - partable = EnsembleParameterTable(nothing) +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 - for (group, ram_matrices) in zip(groups, args) - push!(partable.tables, group => ParameterTable(ram_matrices)) +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 - - return partable end ############################################################################################ @@ -69,52 +92,49 @@ 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 -------------------------------------------------------------------------------- # 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 - -# 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 ############################################################################################ -# update generic --------------------------------------------------------------------------- function update_partable!( - partable::EnsembleParameterTable, - model_identifier::AbstractDict, - vec, - column, + partables::EnsembleParameterTable, + column::Symbol, + param_values::AbstractDict{Symbol}, + default::Any = nothing, ) - for k in keys(partable.tables) - update_partable!(partable.tables[k], model_identifier, vec, column) + for partable in values(partables.tables) + update_partable!(partable, column, param_values, default) end - return partable + return partables +end + +function update_partable!( + partables::EnsembleParameterTable, + column::Symbol, + params::AbstractVector{Symbol}, + values::AbstractVector, + default::Any = nothing, +) + 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 1910d666e..03b111e87 100644 --- a/src/frontend/specification/ParameterTable.jl +++ b/src/frontend/specification/ParameterTable.jl @@ -1,57 +1,90 @@ -abstract type AbstractParameterTable end - ############################################################################################ ### 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} + params::Vector{Symbol} end ############################################################################################ ### Constructors ############################################################################################ -# constuct an empty table -function ParameterTable(::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}(), - :identifier => Vector{Symbol}(), - :start => Vector{Float64}(), - ) +# 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), + :relation => 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), +) - variables = Dict{Symbol, Any}( - :latent_vars => Vector{Symbol}(), - :observed_vars => Vector{Symbol}(), - :sorted_vars => Vector{Symbol}(), +# 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, + params::Union{AbstractVector{Symbol}, Nothing} = nothing, +) + 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}(), + 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(columns, variables) + 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 ############################################################################################ ### Convert to other types ############################################################################################ -import Base.Dict - -function Dict(partable::ParameterTable) +function Base.convert(::Type{Dict}, partable::ParameterTable) return partable.columns end -function DataFrame(partable::ParameterTable; columns = nothing) +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, +) 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 ############################################################################################ @@ -59,37 +92,25 @@ end ############################################################################################ function Base.show(io::IO, partable::ParameterTable) - relevant_columns = [ - :from, - :parameter_type, - :to, - :free, - :value_fixed, - :start, - :estimate, - :se, - :identifier, - ] - existing_columns = [haskey(partable.columns, key) for key in relevant_columns] + relevant_columns = + [: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, + ) - 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, + # TODO switch to `missing` as non-specified values and suppress printing of `missing` instead + formatters = (v, i, j) -> isa(v, Number) && isnan(v) ? "" : v, ) - 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 ############################################################################################ @@ -97,24 +118,33 @@ end ############################################################################################ # Iteration -------------------------------------------------------------------------------- +ParameterTableRow = @NamedTuple begin + from::Symbol + relation::Symbol + to::Symbol + free::Bool + value_fixed::Any + param::Symbol +end -Base.getindex(partable::ParameterTable, i::Int) = ( - partable.columns[:from][i], - partable.columns[:parameter_type][i], - partable.columns[:to][i], - partable.columns[:free][i], - partable.columns[:value_fixed][i], - partable.columns[:identifier][i], +Base.getindex(partable::ParameterTable, i::Integer) = ( + from = partable.columns[:from][i], + relation = partable.columns[:relation][i], + to = partable.columns[:to][i], + free = partable.columns[:free][i], + value_fixed = partable.columns[:value_fixed][i], + param = 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(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) + +params(partable::ParameterTable) = partable.params +nparams(partable::ParameterTable) = length(params(partable)) # Sorting ---------------------------------------------------------------------------------- @@ -124,128 +154,185 @@ end Base.showerror(io::IO, e::CyclicModelError) = print(io, e.msg) -import Base.sort!, Base.sort +""" + sort_vars!(partable::ParameterTable) + sort_vars!(partables::EnsembleParameterTable) -function sort!(partable::ParameterTable) - variables = [partable.variables[:latent_vars]; partable.variables[:observed_vars]] +Sort variables in `partable` so that all independent variables are +before the dependent variables and store it in `partable.sorted_vars`. - is_regression = - (partable.columns[:parameter_type] .== :→) .& - (partable.columns[:from] .!= Symbol("1")) +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 + ] - to = partable.columns[:to][is_regression] - from = partable.columns[:from][is_regression] + # 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")) + ] + sort!(edges, by = last) # sort edges by target - 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 - to = to[delete_edges] - from = from[delete_edges] + for (i, var) in enumerate(vars) + # 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) + # remove var outgoing edges + filter!(e -> e[1] != var, edges) acyclic = true + break end end - if !acyclic + # if acyclic is false, all vars have incoming edge + 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) + @assert length(partable.sorted_vars) == + length(partable.observed_vars) + length(partable.latent_vars) return partable end -function sort(partable::ParameterTable) - new_partable = deepcopy(partable) - sort!(new_partable) - return new_partable -end +""" + sort_vars(partable::ParameterTable) + sort_vars(partables::EnsembleParameterTable) -# add a row -------------------------------------------------------------------------------- +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`. -import Base.push! +See [sort_vars!](@ref) for in-place version. +""" +sort_vars(partable::ParameterTable) = sort_vars!(deepcopy(partable)) + +# add a row -------------------------------------------------------------------------------- -function push!(partable::ParameterTable, d::AbstractDict) - for key in keys(d) - push!(partable.columns[key], d[key]) +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 end -push!(partable::ParameterTable, d::Nothing) = nothing - ############################################################################################ ### Update Partable from Fitted Model ############################################################################################ # update generic --------------------------------------------------------------------------- - function update_partable!( partable::ParameterTable, - model_identifier::AbstractDict, - vec, - 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)) + 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 - push!(partable.columns, column => new_col) + return partable end """ - update_partable!(partable::AbstractParameterTable, sem_fit::SemFit, vec, column) - -Write `vec` to `column` of `partable`. + update_partable!(partable::AbstractParameterTable, params::Vector{Symbol}, values, column) -# Arguments -- `vec::Vector`: has to be in the same order as the `model` parameters +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`. """ -update_partable!(partable::AbstractParameterTable, sem_fit::SemFit, vec, column) = - update_partable!(partable, identifier(sem_fit), vec, column) +function update_partable!( + partable::ParameterTable, + column::Symbol, + params::AbstractVector{Symbol}, + values::AbstractVector, + default::Any = nothing, +) + length(params) == length(values) || throw( + ArgumentError( + "The length of `params` ($(length(params))) and their `values` ($(length(values))) must be the same", + ), + ) + param_values = Dict(zip(params, values)) + if length(param_values) != length(params) + throw(ArgumentError("Duplicate parameter names in `params`")) + end + 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, 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 `:start` 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, 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, @@ -256,17 +343,17 @@ 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, :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. @@ -276,9 +363,207 @@ 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, sem_fit, se, :se) + se = se_hessian(fit; hessian = hessian) + return update_partable!(partable, :se, params(fit), 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) + +""" + 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, :relation, :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, +) diff --git a/src/frontend/specification/RAMMatrices.jl b/src/frontend/specification/RAMMatrices.jl index e0fcc575c..b52c15a51 100644 --- a/src/frontend/specification/RAMMatrices.jl +++ b/src/frontend/specification/RAMMatrices.jl @@ -1,86 +1,29 @@ -############################################################################################ -### Type -############################################################################################ - -# map from parameter index to linear indices of matrix/vector positions where it occurs -AbstractArrayParamsMap = AbstractVector{<:AbstractVector{<:Integer}} -ArrayParamsMap = Vector{Vector{Int}} - -struct RAMMatrices - A_ind::ArrayParamsMap - S_ind::ArrayParamsMap - F_ind::Vector{Int} - M_ind::Union{ArrayParamsMap, Nothing} - parameters::Any - colnames::Any - constants::Any - size_F::Any -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 - F_indices = findall([any(isone.(col)) for col in eachcol(F)]) - constants = get_RAMConstants(A, S, M) - return RAMMatrices( - A_indices, - S_indices, - F_indices, - M_indices, - parameters, - colnames, - constants, - size(F), - ) -end - -RAMMatrices(a::RAMMatrices) = a - ############################################################################################ ### Constants ############################################################################################ struct RAMConstant - matrix::Any - index::Any + matrix::Symbol + index::Union{Int, CartesianIndex{2}} 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 -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])) +function append_RAMConstants!( + constants::AbstractVector{RAMConstant}, + mtx_name::Symbol, + mtx::AbstractArray; + skip_zeros::Bool = true, +) + for (index, val) in pairs(mtx) + if isa(val, Number) && !(skip_zeros && iszero(val)) + push!(constants, RAMConstant(mtx_name, index, val)) 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 @@ -89,7 +32,7 @@ function set_RAMConstant!(A, S, M, rc::RAMConstant) A[rc.index] = rc.value elseif rc.matrix == :S S[rc.index] = rc.value - S[rc.index[2], rc.index[1]] = rc.value + S[rc.index[2], rc.index[1]] = rc.value # symmetric elseif rc.matrix == :M M[rc.index] = rc.value end @@ -102,105 +45,165 @@ function set_RAMConstants!(A, S, M, rc_vec::Vector{RAMConstant}) end ############################################################################################ -### get RAMMatrices from parameter table +### Type ############################################################################################ -function RAMMatrices(partable::ParameterTable; par_id = nothing) - if isnothing(par_id) - parameters, n_par, par_positions = get_par_npar_identifier(partable) - else - parameters, n_par, par_positions = - par_id[:parameters], par_id[:n_par], par_id[:par_positions] +# map from parameter index to linear indices of matrix/vector positions where it occurs +AbstractArrayParamsMap = AbstractVector{<:AbstractVector{<:Integer}} +ArrayParamsMap = Vector{Vector{Int}} + +struct RAMMatrices <: SemSpecification + A_ind::ArrayParamsMap + S_ind::ArrayParamsMap + F_ind::Vector{Int} + M_ind::Union{ArrayParamsMap, Nothing} + params::Vector{Symbol} + colnames::Union{Vector{Symbol}, Nothing} + constants::Vector{RAMConstant} + size_F::Tuple{Int, Int} +end + +nparams(ram::RAMMatrices) = length(ram.A_ind) + +############################################################################################ +### Constructor +############################################################################################ + +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 - n_observed = size(partable.variables[:observed_vars], 1) - n_latent = size(partable.variables[:latent_vars], 1) + 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 + 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) + isnothing(M) || append_RAMConstants!(constants, :M, M) + return RAMMatrices( + A_indices, + S_indices, + F_indices, + M_indices, + params, + colnames, + constants, + size(F), + ) +end + +############################################################################################ +### get RAMMatrices from parameter table +############################################################################################ + +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) 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}() - 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? - 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:length(params)] : nothing - # handel constants + # handle constants constants = Vector{RAMConstant}() - for i in 1:length(partable) - from, parameter_type, to, free, value_fixed, identifier = partable[i] + for r in partable + row_ind = col_indices[r.to] + col_ind = r.from != Symbol("1") ? col_indices[r.from] : nothing - row_ind = positions[to] - if from != Symbol("1") - col_ind = positions[from] - end - - if !free - if (parameter_type == :→) & (from == Symbol("1")) - push!(constants, RAMConstant(:M, row_ind, value_fixed)) - elseif (parameter_type == :→) + if !r.free + if (r.relation == :→) && (r.from == Symbol("1")) + push!(constants, RAMConstant(:M, row_ind, r.value_fixed)) + elseif r.relation == :→ push!( constants, - RAMConstant(:A, CartesianIndex(row_ind, col_ind), value_fixed), + RAMConstant(:A, CartesianIndex(row_ind, col_ind), r.value_fixed), ) - else + elseif r.relation == :↔ 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: $(r.relation)") end else - par_ind = par_positions[identifier] - if (parameter_type == :→) && (from == Symbol("1")) + par_ind = params_index[r.param] + if (r.relation == :→) && (r.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 + elseif r.relation == :→ + push!(A_ind[par_ind], row_ind + (col_ind - 1) * n_node) + 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.relation)") end end end @@ -210,43 +213,59 @@ function RAMMatrices(partable::ParameterTable; par_id = nothing) S_ind, F_ind, M_ind, - parameters, + params, colnames, constants, (n_observed, n_node), ) end +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) - 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))] +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 - partable.variables = Dict( - :sorted_vars => Vector{Symbol}(), - :observed_vars => names_obs, - :latent_vars => names_lat, + # construct an empty table + partable = ParameterTable( + observed_vars = observed_vars, + latent_vars = latent_vars, + params = isnothing(params) ? SEM.params(ram_matrices) : params, ) # 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.parameters) - push_partable_rows!( + for (i, par) in enumerate(ram_matrices.params) + append_partable_rows!( partable, - position_names, + colnames, par, i, ram_matrices.A_ind, @@ -255,28 +274,16 @@ function ParameterTable(ram_matrices::RAMMatrices) ram_matrices.size_F[2], ) end + check_params(SEM.params(partable), partable.columns[:param]) return partable end -############################################################################################ -### get RAMMatrices from EnsembleParameterTable -############################################################################################ - -function RAMMatrices(partable::EnsembleParameterTable) - ram_matrices = Dict{Symbol, RAMMatrices}() - - parameters, n_par, par_positions = get_par_npar_identifier(partable) - par_id = - Dict(:parameters => parameters, :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 +Base.convert( + ::Type{<:ParameterTable}, + ram::RAMMatrices; + params::Union{AbstractVector{Symbol}, Nothing} = nothing, +) = ParameterTable(ram; params) ############################################################################################ ### Pretty Printing @@ -291,156 +298,115 @@ 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 -end - -function get_par_npar_identifier(partable::EnsembleParameterTable) - parameters = Vector{Symbol}() - for key in keys(partable.tables) - append!(parameters, partable.tables[key].columns[:identifier]) - end - parameters = unique(parameters) - filter!(x -> x != :const, parameters) - - n_par = length(parameters) - - par_positions = Dict(parameters .=> 1:n_par) - - return parameters, 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 - identifier = :const - return Dict( - :from => from, - :parameter_type => parameter_type, - :to => to, - :free => free, - :value_fixed => value_fixed, - :start => start, - :estimate => estimate, - :identifier => identifier, - ) -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 +# 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 + 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(par, position_names, index, matrix, n_nod, known_indices) +partable_row(c::RAMConstant, varnames::AbstractVector{Symbol}) = ( + from = varnames[c.index[2]], + relation = matrix_to_relation(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 - - from = position_names[index[2]] - to = position_names[index[1]] - end + cart_index = linear2cartesian(index, (n_nod, n_nod)) - # 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 - identifier = par - - return Dict( - :from => from, - :parameter_type => parameter_type, - :to => to, - :free => free, - :value_fixed => value_fixed, - :start => start, - :estimate => estimate, - :identifier => identifier, + return ( + from = from, + relation = matrix_to_relation(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 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) && (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) ) return res end - -function get_group(d::Dict, group) - return d[group] -end 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/frontend/specification/StenoGraphs.jl b/src/frontend/specification/StenoGraphs.jl index 5c9ce7fdb..67bb7973c 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,60 @@ 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::AbstractVector{Symbol}, + latent_vars::AbstractVector{Symbol}, + params::Union{AbstractVector{Symbol}, Nothing} = nothing, + 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) - identifier = Vector{Symbol}(undef, n) - identifier .= Symbol("") - # group = Vector{Symbol}(undef, n) - # start_partable = zeros(Bool, n) - sorted_vars = Vector{Symbol}() + columns = empty_partable_columns(n) + from = columns[:from] + relation = columns[:relation] + 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) - if element isa DirectedEdge - from[i] = element.src.node - to[i] = element.dst.node - parameter_type[i] = :→ - elseif element isa UndirectedEdge - from[i] = element.src.node - to[i] = element.dst.node - 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 + edge = element isa ModifiedEdge ? element.edge : element + from[i] = edge.src.node + to[i] = edge.dst.node + if edge isa DirectedEdge + relation[i] = :→ + elseif edge isa UndirectedEdge + relation[i] = :↔ + 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 - identifier[i] = modifier.value[g] + param_refs[i] = modval end end end @@ -88,57 +92,45 @@ 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) - current_id += 1 - elseif (identifier[i] == Symbol("")) & !free[i] - identifier[i] = :const - elseif (identifier[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, - :identifier => identifier, - ), - Dict( - :latent_vars => latent_vars, - :observed_vars => observed_vars, - :sorted_vars => sorted_vars, - ), - ) + return ParameterTable(columns; latent_vars, observed_vars, params) end ############################################################################################ ### constructor for EnsembleParameterTable from graph ############################################################################################ -function EnsembleParameterTable(; graph, observed_vars, latent_vars, groups) +function EnsembleParameterTable( + graph::AbstractStenoGraph; + observed_vars::AbstractVector{Symbol}, + latent_vars::AbstractVector{Symbol}, + params::Union{AbstractVector{Symbol}, Nothing} = nothing, + groups, +) graph = unique(graph) - partable = EnsembleParameterTable(nothing) - - for (i, group) in enumerate(groups) - push!( - partable.tables, - Symbol(group) => ParameterTable(; - graph = graph, - observed_vars = observed_vars, - latent_vars = latent_vars, - g = i, - parname = Symbol(:g, i), - ), - ) - end + partables = Dict( + group => ParameterTable( + graph; + observed_vars, + latent_vars, + params, + group = i, + param_prefix = Symbol(:g, group), + ) for (i, group) in enumerate(groups) + ) - return partable + return EnsembleParameterTable(partables; params) end 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 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 00c0d0ef9..bc81a71a0 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 -- `identifier(::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 @@ -65,28 +65,8 @@ Additional interfaces Only available in gradient! calls: - `I_A⁻¹(::RAM)` -> ``(I-A)^{-1}`` """ -mutable struct RAM{ - A1, - A2, - A3, - A4, - A5, - A6, - V, - V2, - I1, - I2, - I3, - M1, - M2, - M3, - M4, - S1, - S2, - S3, - B, - D, -} <: 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 @@ -94,7 +74,6 @@ mutable struct RAM{ μ::A5 M::A6 - n_par::V ram_matrices::V2 has_meanstructure::B @@ -110,8 +89,6 @@ mutable struct RAM{ ∇A::S1 ∇S::S2 ∇M::S3 - - identifier::D end using StructuralEquationModels @@ -121,19 +98,17 @@ using StructuralEquationModels ############################################################################################ function RAM(; - specification, + specification::SemSpecification, #vech = false, gradient = true, meanstructure = false, kwargs..., ) - ram_matrices = RAMMatrices(specification) - identifier = StructuralEquationModels.identifier(ram_matrices) + ram_matrices = convert(RAMMatrices, specification) # get dimensions of the model - n_par = length(ram_matrices.parameters) + n_par = nparams(ram_matrices) 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 @@ -168,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 @@ -185,7 +161,6 @@ function RAM(; F, μ, M_pre, - n_par, ram_matrices, has_meanstructure, A_indices, @@ -198,7 +173,6 @@ function RAM(; ∇A, ∇S, ∇M, - identifier, ) end @@ -213,7 +187,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 +195,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 +213,7 @@ end function gradient!( imply::RAM, - parameters, + params, model::AbstractSemSingle, has_meanstructure::Val{T}, ) where {T} @@ -250,7 +224,7 @@ function gradient!( imply.A_indices, imply.S_indices, imply.M_indices, - parameters, + params, ) @. imply.I_A = -imply.A @@ -281,9 +255,6 @@ objective_gradient_hessian!(imply::RAM, par, model::AbstractSemSingle, has_means ### Recommended methods ############################################################################################ -identifier(imply::RAM) = imply.identifier -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 5c5e52112..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 -- `identifier(::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 @@ -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 - identifier::D1 has_meanstructure::B end @@ -88,7 +86,7 @@ end ############################################################################################ function RAMSymbolic(; - specification, + specification::SemSpecification, loss_types = nothing, vech = false, gradient = true, @@ -97,10 +95,9 @@ function RAMSymbolic(; approximate_hessian = false, kwargs..., ) - ram_matrices = RAMMatrices(specification) - identifier = StructuralEquationModels.identifier(ram_matrices) + ram_matrices = convert(RAMMatrices, specification) - n_par = length(ram_matrices.parameters) + n_par = nparams(ram_matrices) n_var, n_nod = ram_matrices.size_F par = (Symbolics.@variables θ[1:n_par])[1] @@ -144,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] @@ -195,13 +191,11 @@ function RAMSymbolic(; Σ_symbolic, ∇Σ_symbolic, ∇²Σ_symbolic, - n_par, ram_matrices, μ_function, μ, ∇μ_function, ∇μ, - identifier, has_meanstructure, ) end @@ -240,9 +234,6 @@ objective_gradient_hessian!(imply::RAMSymbolic, par, model) = gradient!(imply, p ### Recommended methods ############################################################################################ -identifier(imply::RAMSymbolic) = imply.identifier -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 ba8580d16..f1af2ec42 100644 --- a/src/imply/empty.jl +++ b/src/imply/empty.jl @@ -19,15 +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 -- `n_par(::RAMSymbolic)` -> Number of parameters +- `params(::RAMSymbolic) `-> Vector of parameter labels +- `nparams(::RAMSymbolic)` -> Number of parameters ## Implementation Subtype of `SemImply`. """ -struct ImplyEmpty{V, V2} <: SemImply - identifier::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.parameters) - - return ImplyEmpty(identifier(ram_matrices), n_par) + return ImplyEmpty(convert(RAMMatrices, specification)) end ############################################################################################ @@ -54,7 +49,4 @@ hessian!(imply::ImplyEmpty, par, model) = nothing ### Recommended methods ############################################################################################ -identifier(imply::ImplyEmpty) = imply.identifier -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..d4870ac1b 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], @@ -252,5 +252,4 @@ 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/loss/regularization/ridge.jl b/src/loss/regularization/ridge.jl index ebf3e7bfe..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..., @@ -58,7 +58,8 @@ function SemRidge(; ), ) else - which_ridge = get_identifier_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] @@ -67,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/objective_gradient_hessian.jl b/src/objective_gradient_hessian.jl index 53d68ec2c..2debbcd40 100644 --- a/src/objective_gradient_hessian.jl +++ b/src/objective_gradient_hessian.jl @@ -2,56 +2,51 @@ # 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, -) +function objective_gradient_hessian!(gradient, hessian, model::AbstractSemSingle, 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 +55,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 +336,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/observed/covariance.jl b/src/observed/covariance.jl index 1b5de9fc2..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, + 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/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/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/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, diff --git a/src/types.jl b/src/types.jl index 803bc733a..d194709b8 100644 --- a/src/types.jl +++ b/src/types.jl @@ -13,6 +13,36 @@ 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 + +""" + nparams(semobj) + +Return the number of SEM model parameters. +""" +nparams(model::AbstractSem) = length(params(model)) + +params(model::AbstractSemSingle) = params(model.imply) +nparams(model::AbstractSemSingle) = nparams(model.imply) + +""" + param_indices(semobj) + +Returns a dict of parameter names and their indices in `semobj`. + +# Examples +```julia +parind = param_indices(my_fitted_sem) +parind[:param_name] +``` +""" +param_indices(semobj) = Dict(par => i for (i, par) in enumerate(params(semobj))) + """ SemLoss(args...; loss_weights = nothing, ...) @@ -75,6 +105,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) +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 @@ -153,19 +186,18 @@ 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. +- `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 - identifier::I + params::I end function SemEnsemble(models...; optimizer = SemOptimizerOptim, weights = nothing, kwargs...) n = length(models) - npar = n_par(models[1]) # default weights @@ -174,11 +206,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 + params = SEM.params(models[1]) for model in models - if id != identifier(model) - throw(ErrorException("The identifier of your models do not match. \n + 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.")) end @@ -189,9 +221,12 @@ function SemEnsemble(models...; optimizer = SemOptimizerOptim, weights = nothing optimizer = optimizer(; kwargs...) end - return SemEnsemble(n, models, weights, optimizer, id) + return SemEnsemble(n, models, weights, optimizer, params) end +params(ensemble::SemEnsemble) = ensemble.params +nparams(ensemble::SemEnsemble) = length(ensemble.params) + """ n_models(ensemble::SemEnsemble) -> Integer @@ -247,3 +282,13 @@ 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 + +params(spec::SemSpecification) = spec.params +nparams(spec::SemSpecification) = length(params(spec)) + +abstract type AbstractParameterTable <: SemSpecification end diff --git a/test/examples/helper.jl b/test/examples/helper.jl index 3bb4e217a..d4c140d67 100644 --- a/test/examples/helper.jl +++ b/test/examples/helper.jl @@ -1,43 +1,44 @@ -function test_gradient(model, parameters; rtol = 1e-10, atol = 0) - true_grad = - FiniteDiff.finite_difference_gradient(Base.Fix1(objective!, model), parameters) - gradient = similar(parameters) +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) # 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 @@ -47,7 +48,7 @@ fitmeasure_names_ml = Dict( :df => "df", :χ² => "chisq", :p_value => "pvalue", - :n_par => "npar", + :nparams => "npar", :RMSEA => "rmsea", ) @@ -55,7 +56,7 @@ fitmeasure_names_ls = Dict( :df => "df", :χ² => "chisq", :p_value => "pvalue", - :n_par => "npar", + :nparams => "npar", :RMSEA => "rmsea", ) @@ -72,264 +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.variables[:latent_vars]) & - (to ∈ partable.variables[: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.variables[:latent_vars]) & - (to ∈ partable.variables[: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 9b97300df..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, @@ -49,9 +49,9 @@ end # ML estimation - sorted ############################################################################################ -partable_s = sort(partable) +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] @@ -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, @@ -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(Σ) @@ -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/multigroup/multigroup.jl b/test/examples/multigroup/multigroup.jl index 759c24eda..e428eba1d 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,14 +64,12 @@ 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], ) partable = EnsembleParameterTable( - specification_g1, - specification_g2; - groups = [:Pasteur, :Grant_White], + Dict(:Pasteur => specification_g1, :Grant_White => specification_g2), ) specification_miss_g1 = nothing @@ -111,14 +109,14 @@ 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], ) -specification = RAMMatrices(partable) +specification = convert(Dict{Symbol, RAMMatrices}, partable) specification_g1 = specification[:Pasteur] specification_g2 = specification[:Grant_White] @@ -140,14 +138,14 @@ 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], ) -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] diff --git a/test/examples/political_democracy/by_parts.jl b/test/examples/political_democracy/by_parts.jl index c071d9e00..11953ccb6 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) @@ -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/examples/political_democracy/political_democracy.jl b/test/examples/political_democracy/political_democracy.jl index 389800745..d7fbb8f2c 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, @@ -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,10 +213,9 @@ 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) +sort_vars!(spec) partable = spec @@ -244,10 +244,9 @@ 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) +sort_vars!(spec_mean) partable_mean = spec_mean diff --git a/test/examples/recover_parameters/recover_parameters_twofact.jl b/test/examples/recover_parameters/recover_parameters_twofact.jl index 68e44ce20..5aa79842c 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) @@ -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), diff --git a/test/unit_tests/data_input_formats.jl b/test/unit_tests/data_input_formats.jl index 485cf82d2..44656c331 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 ----------------------------------------------------------------------------------- @@ -208,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) diff --git a/test/unit_tests/sorting.jl b/test/unit_tests/sorting.jl index e573c6d22..f5bc38ae0 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) @@ -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 diff --git a/test/unit_tests/specification.jl b/test/unit_tests/specification.jl index 0bfc0de2d..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 get_identifier_indices([:x2, :x10, :x28], model_ml) == [2, 10, 28] - -@testset "get_identifier_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)) +@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: -parameter_indices = get_identifier_indices([:λ₁, λ₂], my_fitted_sem) -values = solution(my_fitted_sem)[parameter_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