From 3bda87d61858cc890ac19ebfc8f37b83417142d7 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Fri, 13 Sep 2024 17:55:58 -0700 Subject: [PATCH 01/12] Project.toml: update Symbolics deps --- Project.toml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index 9edeb4536..b038c3364 100644 --- a/Project.toml +++ b/Project.toml @@ -36,8 +36,8 @@ NLopt = "0.6, 1" Optim = "1" PrettyTables = "2" StatsBase = "0.33, 0.34" -Symbolics = "4, 5" -SymbolicUtils = "1.4 - 1.5" +Symbolics = "4, 5, 6" +SymbolicUtils = "1.4 - 1.5, 1.7, 2, 3" [extras] Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" From 70a4e1f58f56651450fe57d3ed3ba8a75410e76d Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Fri, 15 Mar 2024 08:36:18 -0700 Subject: [PATCH 02/12] tests/examples: import -> using no declarations, so import is not required --- test/examples/multigroup/build_models.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/test/examples/multigroup/build_models.jl b/test/examples/multigroup/build_models.jl index 4b5afd58e..3f29a6898 100644 --- a/test/examples/multigroup/build_models.jl +++ b/test/examples/multigroup/build_models.jl @@ -1,3 +1,5 @@ +const SEM = StructuralEquationModels + ############################################################################################ # ML estimation ############################################################################################ From 56a1b0426461a7ec4767f4a73ddfd17f32b575e4 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Fri, 22 Nov 2024 11:28:41 -0800 Subject: [PATCH 03/12] add ParamsArray replaces RAMMatrices indices and constants vectors with dedicated class that incapsulate this logic, resulting in overall cleaner interface A_ind, S_ind, M_ind become ParamsArray F_ind becomes SparseMatrixCSC parameters.jl is not longer required and is removed --- src/StructuralEquationModels.jl | 2 +- src/additional_functions/parameters.jl | 137 ------- src/additional_functions/params_array.jl | 204 ++++++++++ .../start_val/start_fabin3.jl | 152 ++++--- .../start_val/start_simple.jl | 20 +- src/frontend/specification/RAMMatrices.jl | 381 ++++++++---------- src/imply/RAM/generic.jl | 61 +-- src/imply/RAM/symbolic.jl | 21 +- 8 files changed, 486 insertions(+), 492 deletions(-) delete mode 100644 src/additional_functions/parameters.jl create mode 100644 src/additional_functions/params_array.jl diff --git a/src/StructuralEquationModels.jl b/src/StructuralEquationModels.jl index 944542379..6172af1ea 100644 --- a/src/StructuralEquationModels.jl +++ b/src/StructuralEquationModels.jl @@ -26,6 +26,7 @@ include("objective_gradient_hessian.jl") # helper objects and functions include("additional_functions/commutation_matrix.jl") +include("additional_functions/params_array.jl") # fitted objects include("frontend/fit/SemFit.jl") @@ -69,7 +70,6 @@ include("optimizer/optim.jl") include("optimizer/NLopt.jl") # helper functions include("additional_functions/helper.jl") -include("additional_functions/parameters.jl") include("additional_functions/start_val/start_val.jl") include("additional_functions/start_val/start_fabin3.jl") include("additional_functions/start_val/start_partable.jl") diff --git a/src/additional_functions/parameters.jl b/src/additional_functions/parameters.jl deleted file mode 100644 index d6e8eb535..000000000 --- a/src/additional_functions/parameters.jl +++ /dev/null @@ -1,137 +0,0 @@ -# fill A, S, and M matrices with the parameter values according to the parameters map -function fill_A_S_M!( - A::AbstractMatrix, - S::AbstractMatrix, - M::Union{AbstractVector, Nothing}, - A_indices::AbstractArrayParamsMap, - S_indices::AbstractArrayParamsMap, - M_indices::Union{AbstractArrayParamsMap, Nothing}, - params::AbstractVector, -) - @inbounds for (iA, iS, par) in zip(A_indices, S_indices, params) - for index_A in iA - A[index_A] = par - end - - for index_S in iS - S[index_S] = par - end - end - - if !isnothing(M) - @inbounds for (iM, par) in zip(M_indices, params) - for index_M in iM - M[index_M] = par - end - end - end -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_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(params)] - for (i, val) in enumerate(M) - par_ind = get(params_index, val, nothing) - if !isnothing(par_ind) - push!(res[par_ind], i) - end - end - return res -end - -function eachindex_lower(M; linear_indices = false, kwargs...) - indices = CartesianIndices(M) - indices = filter(x -> (x[1] >= x[2]), indices) - - if linear_indices - indices = cartesian2linear(indices, M) - end - - return indices -end - -function cartesian2linear(ind_cart, dims) - ind_lin = LinearIndices(dims)[ind_cart] - return ind_lin -end - -function linear2cartesian(ind_lin, dims) - ind_cart = CartesianIndices(dims)[ind_lin] - return ind_cart -end - -function set_constants!(M, M_pre) - for index in eachindex(M) - δ = tryparse(Float64, string(M[index])) - - if !iszero(M[index]) & (δ !== nothing) - M_pre[index] = δ - end - end -end - -function check_constants(M) - for index in eachindex(M) - δ = tryparse(Float64, string(M[index])) - - if !iszero(M[index]) & (δ !== nothing) - return true - end - end - - return false -end - -# construct length(M)×length(parameters) sparse matrix of 1s at the positions, -# where the corresponding parameter occurs in the M matrix -function matrix_gradient(M_indices::ArrayParamsMap, M_length::Integer) - rowval = reduce(vcat, M_indices) - colptr = - pushfirst!(accumulate((ptr, M_ind) -> ptr + length(M_ind), M_indices, init = 1), 1) - return SparseMatrixCSC( - M_length, - length(M_indices), - colptr, - rowval, - ones(length(rowval)), - ) -end - -# fill M with parameters -function fill_matrix!( - M::AbstractMatrix, - M_indices::AbstractArrayParamsMap, - params::AbstractVector, -) - for (iM, par) in zip(M_indices, params) - for index_M in iM - M[index_M] = par - end - end - return M -end - -# range of parameters that are referenced in the matrix -function param_range(mtx_indices::AbstractArrayParamsMap) - first_i = findfirst(!isempty, mtx_indices) - last_i = findlast(!isempty, mtx_indices) - - if !isnothing(first_i) && !isnothing(last_i) - for i in first_i:last_i - if isempty(mtx_indices[i]) - # TODO show which parameter is missing in which matrix - throw( - ErrorException( - "Your parameter vector is not partitioned into directed and undirected effects", - ), - ) - end - end - end - - return first_i:last_i -end diff --git a/src/additional_functions/params_array.jl b/src/additional_functions/params_array.jl new file mode 100644 index 000000000..f20a6518b --- /dev/null +++ b/src/additional_functions/params_array.jl @@ -0,0 +1,204 @@ +""" +Array with partially parameterized elements. +""" +struct ParamsArray{T, N} <: AbstractArray{T, N} + linear_indices::Vector{Int} + param_ptr::Vector{Int} + constants::Vector{Pair{Int, T}} + size::NTuple{N, Int} +end + +ParamsVector{T} = ParamsArray{T, 1} +ParamsMatrix{T} = ParamsArray{T, 2} + +function ParamsArray{T, N}( + params_map::AbstractVector{<:AbstractVector{Int}}, + constants::Vector{Pair{Int, T}}, + size::NTuple{N, Int}, +) where {T, N} + params_ptr = + pushfirst!(accumulate((ptr, inds) -> ptr + length(inds), params_map, init = 1), 1) + return ParamsArray{T, N}( + reduce(vcat, params_map, init = Vector{Int}()), + params_ptr, + constants, + size, + ) +end + +function ParamsArray{T, N}( + arr::AbstractArray{<:Any, N}, + params::AbstractVector{Symbol}; + skip_zeros::Bool = true, +) where {T, N} + params_index = Dict(param => i for (i, param) in enumerate(params)) + constants = Vector{Pair{Int, T}}() + params_map = [Vector{Int}() for _ in eachindex(params)] + arr_ixs = CartesianIndices(arr) + for (i, val) in pairs(vec(arr)) + ismissing(val) && continue + if isa(val, Number) + (skip_zeros && iszero(val)) || push!(constants, i => val) + else + par_ind = get(params_index, val, nothing) + if !isnothing(par_ind) + push!(params_map[par_ind], i) + else + throw(KeyError("Unrecognized parameter $val at position $(arr_ixs[i])")) + end + end + end + return ParamsArray{T, N}(params_map, constants, size(arr)) +end + +ParamsArray{T}( + arr::AbstractArray{<:Any, N}, + params::AbstractVector{Symbol}; + kwargs..., +) where {T, N} = ParamsArray{T, N}(arr, params; kwargs...) + +nparams(arr::ParamsArray) = length(arr.param_ptr) - 1 + +Base.size(arr::ParamsArray) = arr.size +Base.size(arr::ParamsArray, i::Integer) = arr.size[i] + +Base.:(==)(a::ParamsArray, b::ParamsArray) = return eltype(a) == eltype(b) && + size(a) == size(b) && + a.constants == b.constants && + a.param_ptr == b.param_ptr && + a.linear_indices == b.linear_indices + +# the range of arr.param_ptr indices that correspond to i-th parameter +param_occurences_range(arr::ParamsArray, i::Integer) = + arr.param_ptr[i]:(arr.param_ptr[i+1]-1) + +""" + param_occurences(arr::ParamsArray, i::Integer) + +Get the linear indices of the elements in `arr` that correspond to the +`i`-th parameter. +""" +param_occurences(arr::ParamsArray, i::Integer) = + view(arr.linear_indices, arr.param_ptr[i]:(arr.param_ptr[i+1]-1)) + +""" + materialize!(dest::AbstractArray{<:Any, N}, src::ParamsArray{<:Any, N}, + param_values::AbstractVector; + set_constants::Bool = true, + set_zeros::Bool = false) + +Materialize the parameterized array `src` into `dest` by substituting the parameter +references with the parameter values from `param_values`. +""" +function materialize!( + dest::AbstractArray{<:Any, N}, + src::ParamsArray{<:Any, N}, + param_values::AbstractVector; + set_constants::Bool = true, + set_zeros::Bool = false, +) where {N} + size(dest) == size(src) || throw( + DimensionMismatch( + "Parameters ($(size(params_arr))) and destination ($(size(dest))) array sizes don't match", + ), + ) + nparams(src) == length(param_values) || throw( + DimensionMismatch( + "Number of values ($(length(param_values))) does not match the number of parameters ($(nparams(src)))", + ), + ) + Z = eltype(dest) <: Number ? eltype(dest) : eltype(src) + set_zeros && fill!(dest, zero(Z)) + if set_constants + @inbounds for (i, val) in src.constants + dest[i] = val + end + end + @inbounds for (i, val) in enumerate(param_values) + for j in param_occurences_range(src, i) + dest[src.linear_indices[j]] = val + end + end + return dest +end + +""" + materialize([T], src::ParamsArray{<:Any, N}, + param_values::AbstractVector{T}) where T + +Materialize the parameterized array `src` into a new array of type `T` +by substituting the parameter references with the parameter values from `param_values`. +""" +materialize(::Type{T}, arr::ParamsArray, param_values::AbstractVector) where {T} = + materialize!(similar(arr, T), arr, param_values, set_constants = true, set_zeros = true) + +materialize(arr::ParamsArray, param_values::AbstractVector{T}) where {T} = + materialize(Union{T, eltype(arr)}, arr, param_values) + +function sparse_materialize( + ::Type{T}, + arr::ParamsMatrix, + param_values::AbstractVector, +) where {T} + nparams(arr) == length(param_values) || throw( + DimensionMismatch( + "Number of values ($(length(param)values))) does not match the number of parameter ($(nparams(arr)))", + ), + ) + # constant values in sparse matrix + cvals = [T(v) for (_, v) in arr.constants] + # parameter values in sparse matrix + parvals = Vector{T}(undef, length(arr.linear_indices)) + @inbounds for (i, val) in enumerate(param_values) + for j in param_occurences_range(arr, i) + parvals[j] = val + end + end + nzixs = [first.(arr.constants); arr.linear_indices] + ixorder = sortperm(nzixs) + nzixs = nzixs[ixorder] + nzvals = [cvals; parvals][ixorder] + arr_ixs = CartesianIndices(size(arr)) + return sparse( + [arr_ixs[i][1] for i in nzixs], + [arr_ixs[i][2] for i in nzixs], + nzvals, + size(arr)..., + ) +end + +sparse_materialize(arr::ParamsArray, params::AbstractVector{T}) where {T} = + sparse_materialize(Union{T, eltype(arr)}, arr, params) + +# construct length(M)×length(params) sparse matrix of 1s at the positions, +# where the corresponding parameter occurs in the arr +sparse_gradient(::Type{T}, arr::ParamsArray) where {T} = SparseMatrixCSC( + length(arr), + nparams(arr), + arr.param_ptr, + arr.linear_indices, + ones(T, length(arr.linear_indices)), +) + +sparse_gradient(arr::ParamsArray{T}) where {T} = sparse_gradient(T, arr) + +# range of parameters that are referenced in the matrix +function params_range(arr::ParamsArray; allow_gaps::Bool = false) + first_i = findfirst(i -> arr.param_ptr[i+1] > arr.param_ptr[i], 1:nparams(arr)-1) + last_i = findlast(i -> arr.param_ptr[i+1] > arr.param_ptr[i], 1:nparams(arr)-1) + + if !allow_gaps && !isnothing(first_i) && !isnothing(last_i) + for i in first_i:last_i + if isempty(param_occurences_range(arr, i)) + # TODO show which parameter is missing in which matrix + throw( + ErrorException( + "Parameter vector is not partitioned into directed and undirected effects", + ), + ) + end + end + end + + return first_i:last_i +end diff --git a/src/additional_functions/start_val/start_fabin3.jl b/src/additional_functions/start_val/start_fabin3.jl index 081af3ba1..9d692437e 100644 --- a/src/additional_functions/start_val/start_fabin3.jl +++ b/src/additional_functions/start_val/start_fabin3.jl @@ -31,23 +31,20 @@ function start_fabin3(observed::SemObservedMissing, imply, optimizer, args...; k end function start_fabin3(ram_matrices::RAMMatrices, Σ, μ) - A_ind, S_ind, F_ind, M_ind, n_par = ram_matrices.A_ind, - ram_matrices.S_ind, - ram_matrices.F_ind, - ram_matrices.M_ind, + A, S, F, M, n_par = ram_matrices.A, + ram_matrices.S, + ram_matrices.F, + ram_matrices.M, nparams(ram_matrices) start_val = zeros(n_par) - n_obs = nobserved_vars(ram_matrices) - n_var = nvars(ram_matrices) - n_latent = nlatent_vars(ram_matrices) - - C_indices = CartesianIndices((n_var, n_var)) + F_var2obs = Dict( + i => F.rowval[F.colptr[i]] for i in axes(F, 2) if isobserved_var(ram_matrices, i) + ) + @assert length(F_var2obs) == size(F, 1) # check in which matrix each parameter appears - indices = Vector{CartesianIndex{2}}(undef, n_par) - #= in_S = length.(S_ind) .!= 0 in_A = length.(A_ind) .!= 0 A_ind_c = [linear2cartesian(ind, (n_var, n_var)) for ind in A_ind] @@ -65,26 +62,53 @@ function start_fabin3(ram_matrices::RAMMatrices, Σ, μ) end =# # set undirected parameters in S - for (i, S_ind) in enumerate(S_ind) - for c_ind in C_indices[S_ind] - (c_ind[1] == c_ind[2]) || continue # covariances stay 0 - pos = searchsortedfirst(F_ind, c_ind[1]) - start_val[i] = - (pos <= length(F_ind)) && (F_ind[pos] == c_ind[1]) ? Σ[pos, pos] / 2 : 0.05 - break # i-th parameter initialized + S_indices = CartesianIndices(S) + for j in 1:nparams(S) + for lin_ind in param_occurences(S, j) + to, from = Tuple(S_indices[lin_ind]) + if (to == from) # covariances start with 0 + # half of observed variance for observed, 0.05 for latent + obs = get(F_var2obs, to, nothing) + start_val[j] = !isnothing(obs) ? Σ[obs, obs] / 2 : 0.05 + break # j-th parameter initialized + end end end # set loadings - constants = ram_matrices.constants - A_ind_c = [linear2cartesian(ind, (n_var, n_var)) for ind in A_ind] + A_indices = CartesianIndices(A) # ind_Λ = findall([is_in_Λ(ind_vec, F_ind) for ind_vec in A_ind_c]) - function calculate_lambda( - ref::Integer, - indicator::Integer, - indicators::AbstractVector{<:Integer}, - ) + # collect latent variable indicators in A + # maps latent parameter to the vector of dependent vars + # the 2nd index in the pair specified the parameter index, + # 0 if no parameter (constant), -1 if constant=1 + var2indicators = Dict{Int, Vector{Pair{Int, Int}}}() + for j in 1:nparams(A) + for lin_ind in param_occurences(A, j) + to, from = Tuple(A_indices[lin_ind]) + haskey(F_var2obs, from) && continue # skip observed + obs = get(F_var2obs, to, nothing) + if !isnothing(obs) + indicators = get!(() -> Vector{Pair{Int, Int}}(), var2indicators, from) + push!(indicators, obs => j) + end + end + end + + for (lin_ind, val) in A.constants + iszero(val) && continue # only non-zero loadings + to, from = Tuple(A_indices[lin_ind]) + haskey(F_var2obs, from) && continue # skip observed + obs = get(F_var2obs, to, nothing) + if !isnothing(obs) + indicators = get!(() -> Vector{Pair{Int, Int}}(), var2indicators, from) + push!(indicators, obs => ifelse(isone(val), -1, 0)) # no parameter associated, -1 = reference, 0 = indicator + end + end + + # calculate starting values for parameters of latent regression vars + function calculate_lambda(ref::Integer, indicator::Integer, indicators::AbstractVector) instruments = filter(i -> (i != ref) && (i != indicator), indicators) if length(instruments) == 1 s13 = Σ[ref, instruments[1]] @@ -99,61 +123,33 @@ function start_fabin3(ram_matrices::RAMMatrices, Σ, μ) end end - for i in setdiff(1:n_var, F_ind) - reference = Int64[] - indicators = Int64[] - indicator2parampos = Dict{Int, Int}() - - for (j, Aj_ind_c) in enumerate(A_ind_c) - for ind_c in Aj_ind_c - (ind_c[2] == i) || continue - ind_pos = searchsortedfirst(F_ind, ind_c[1]) - if (ind_pos <= length(F_ind)) && (F_ind[ind_pos] == ind_c[1]) - push!(indicators, ind_pos) - indicator2parampos[ind_pos] = j - end - end - end - - for ram_const in constants - if (ram_const.matrix == :A) && (ram_const.index[2] == i) - ind_pos = searchsortedfirst(F_ind, ram_const.index[1]) - if (ind_pos <= length(F_ind)) && (F_ind[ind_pos] == ram_const.index[1]) - if isone(ram_const.value) - push!(reference, ind_pos) - else - push!(indicators, ind_pos) - # no parameter associated - end - end - end - end - + for (i, indicators) in pairs(var2indicators) + reference = [obs for (obs, param) in indicators if param == -1] + indicator_obs = first.(indicators) # is there at least one reference indicator? if length(reference) > 0 - if (length(reference) > 1) && isempty(indicator2parampos) # don't warn if entire column is fixed + if (length(reference) > 1) && any(((obs, param),) -> param > 0, indicators) # don't warn if entire column is fixed @warn "You have more than 1 scaling indicator for $(ram_matrices.colnames[i])" end ref = reference[1] - for (j, indicator) in enumerate(indicators) - if (indicator != ref) && - (parampos = get(indicator2parampos, indicator, 0)) != 0 - start_val[parampos] = calculate_lambda(ref, indicator, indicators) + for (indicator, param) in indicators + if (indicator != ref) && (param > 0) + start_val[param] = calculate_lambda(ref, indicator, indicator_obs) end end # no reference indicator: - elseif length(indicators) > 0 - ref = indicators[1] - λ = Vector{Float64}(undef, length(indicators)) + else + ref = indicator_obs[1] + λ = Vector{Float64}(undef, length(indicator_obs)) λ[1] = 1.0 - for (j, indicator) in enumerate(indicators) + for (j, indicator) in enumerate(indicator_obs) if indicator != ref - λ[j] = calculate_lambda(ref, indicator, indicators) + λ[j] = calculate_lambda(ref, indicator, indicator_obs) end end - Σ_λ = Σ[indicators, indicators] + Σ_λ = Σ[indicator_obs, indicator_obs] l₂ = sum(abs2, λ) D = λ * λ' ./ l₂ θ = (I - D .^ 2) \ (diag(Σ_λ - D * Σ_λ * D)) @@ -164,24 +160,22 @@ function start_fabin3(ram_matrices::RAMMatrices, Σ, μ) λ .*= sign(Ψ) * sqrt(abs(Ψ)) - for (j, indicator) in enumerate(indicators) - if (parampos = get(indicator2parampos, indicator, 0)) != 0 - start_val[parampos] = λ[j] + for (j, (_, param)) in enumerate(indicators) + if param > 0 + start_val[param] = λ[j] end end - else - @warn "No scaling indicators for $(ram_matrices.colnames[i])" end end - # set means - if !isnothing(M_ind) - for (i, M_ind) in enumerate(M_ind) - if length(M_ind) != 0 - ind = M_ind[1] - pos = searchsortedfirst(F_ind, ind[1]) - if (pos <= length(F_ind)) && (F_ind[pos] == ind[1]) - start_val[i] = μ[pos] + if !isnothing(M) + # set starting values of the observed means + for j in 1:nparams(M) + M_ind = param_occurences(M, j) + if !isempty(M_ind) + obs = get(F_var2obs, M_ind[1], nothing) + if !isnothing(obs) + start_val[j] = μ[obs] end # latent means stay 0 end end diff --git a/src/additional_functions/start_val/start_simple.jl b/src/additional_functions/start_val/start_simple.jl index 3b29ec178..1f73a3583 100644 --- a/src/additional_functions/start_val/start_simple.jl +++ b/src/additional_functions/start_val/start_simple.jl @@ -62,10 +62,10 @@ function start_simple( start_means = 0.0, kwargs..., ) - A_ind, S_ind, F_ind, M_ind, n_par = ram_matrices.A_ind, - ram_matrices.S_ind, - ram_matrices.F_ind, - ram_matrices.M_ind, + A, S, F_ind, M, n_par = ram_matrices.A, + ram_matrices.S, + observed_var_indices(ram_matrices), + ram_matrices.M, nparams(ram_matrices) start_val = zeros(n_par) @@ -75,9 +75,11 @@ function start_simple( C_indices = CartesianIndices((n_var, n_var)) for i in 1:n_par - if length(S_ind[i]) != 0 + Si_ind = param_occurences(S, i) + Ai_ind = param_occurences(A, i) + if length(Si_ind) != 0 # use the first occurence of the parameter to determine starting value - c_ind = C_indices[S_ind[i][1]] + c_ind = C_indices[Si_ind[1]] if c_ind[1] == c_ind[2] if c_ind[1] ∈ F_ind start_val[i] = start_variances_observed @@ -95,14 +97,14 @@ function start_simple( start_val[i] = start_covariances_obs_lat end end - elseif length(A_ind[i]) != 0 - c_ind = C_indices[A_ind[i][1]] + elseif length(Ai_ind) != 0 + c_ind = C_indices[Ai_ind[1]] if (c_ind[1] ∈ F_ind) & !(c_ind[2] ∈ F_ind) start_val[i] = start_loadings else start_val[i] = start_regressions end - elseif !isnothing(M_ind) && (length(M_ind[i]) != 0) + elseif !isnothing(M) && (length(param_occurences(M, i)) != 0) start_val[i] = start_means end end diff --git a/src/frontend/specification/RAMMatrices.jl b/src/frontend/specification/RAMMatrices.jl index 6ba6be3d0..b140ae026 100644 --- a/src/frontend/specification/RAMMatrices.jl +++ b/src/frontend/specification/RAMMatrices.jl @@ -1,82 +1,56 @@ -############################################################################################ -### Constants -############################################################################################ - -struct RAMConstant - matrix::Symbol - index::Union{Int, CartesianIndex{2}} - value::Any -end - -function Base.:(==)(c1::RAMConstant, c2::RAMConstant) - res = ((c1.matrix == c2.matrix) && (c1.index == c2.index) && (c1.value == c2.value)) - return res -end - -function append_RAMConstants!( - constants::AbstractVector{RAMConstant}, - mtx_name::Symbol, - mtx::AbstractArray; - 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 - return constants -end - -function set_RAMConstant!(A, S, M, rc::RAMConstant) - if rc.matrix == :A - A[rc.index] = rc.value - elseif rc.matrix == :S - S[rc.index] = rc.value - S[rc.index[2], rc.index[1]] = rc.value # symmetric - elseif rc.matrix == :M - M[rc.index] = rc.value - end -end - -function set_RAMConstants!(A, S, M, rc_vec::Vector{RAMConstant}) - for rc in rc_vec - set_RAMConstant!(A, S, M, rc) - end -end ############################################################################################ ### Type ############################################################################################ -# 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} + A::ParamsMatrix{Float64} + S::ParamsMatrix{Float64} + F::SparseMatrixCSC{Float64} + M::Union{ParamsVector{Float64}, Nothing} params::Vector{Symbol} - colnames::Union{Vector{Symbol}, Nothing} - constants::Vector{RAMConstant} - size_F::Tuple{Int, Int} + colnames::Union{Vector{Symbol}, Nothing} # better call it "variables": it's a mixture of observed and latent (and it gets confusing with get_colnames()) end -nparams(ram::RAMMatrices) = length(ram.A_ind) - -nvars(ram::RAMMatrices) = ram.size_F[2] -nobserved_vars(ram::RAMMatrices) = ram.size_F[1] +nparams(ram::RAMMatrices) = nparams(ram.A) +nvars(ram::RAMMatrices) = size(ram.F, 2) +nobserved_vars(ram::RAMMatrices) = size(ram.F, 1) nlatent_vars(ram::RAMMatrices) = nvars(ram) - nobserved_vars(ram) vars(ram::RAMMatrices) = ram.colnames +isobserved_var(ram::RAMMatrices, i::Integer) = ram.F.colptr[i+1] > ram.F.colptr[i] +islatent_var(ram::RAMMatrices, i::Integer) = ram.F.colptr[i+1] == ram.F.colptr[i] + +# indices of observed variables in the order as they appear in ram.F rows +function observed_var_indices(ram::RAMMatrices) + obs_inds = Vector{Int}(undef, nobserved_vars(ram)) + @inbounds for i in 1:nvars(ram) + colptr = ram.F.colptr[i] + if ram.F.colptr[i+1] > colptr # is observed + obs_inds[ram.F.rowval[colptr]] = i + end + end + return obs_inds +end + +latent_var_indices(ram::RAMMatrices) = + [i for i in axes(ram.F, 2) if islatent_var(ram, i)] + +# observed variables in the order as they appear in ram.F rows function observed_vars(ram::RAMMatrices) if isnothing(ram.colnames) @warn "Your RAMMatrices do not contain column names. Please make sure the order of variables in your data is correct!" return nothing else - return view(ram.colnames, ram.F_ind) + obs_vars = Vector{Symbol}(undef, nobserved_vars(ram)) + @inbounds for (i, v) in enumerate(vars(ram)) + colptr = ram.F.colptr[i] + if ram.F.colptr[i+1] > colptr # is observed + obs_vars[ram.F.rowval[colptr]] = v + end + end + return obs_vars end end @@ -85,7 +59,7 @@ function latent_vars(ram::RAMMatrices) @warn "Your RAMMatrices do not contain column names. Please make sure the order of variables in your data is correct!" return nothing else - return view(ram.colnames, setdiff(eachindex(ram.colnames), ram.F_ind)) + return [col for (i, col) in enumerate(ram.colnames) if islatent_var(ram, i)] end end @@ -128,27 +102,16 @@ function RAMMatrices(; ), ) end - check_params(params, nothing) - A_indices = array_params_map(params, A) - S_indices = array_params_map(params, S) - M_indices = !isnothing(M) ? array_params_map(params, M) : nothing - 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), - ) + A = ParamsMatrix{Float64}(A, params) + S = ParamsMatrix{Float64}(S, params) + M = !isnothing(M) ? ParamsVector{Float64}(M, params) : nothing + spF = sparse(F) + if any(!isone, spF.nzval) + throw(ArgumentError("F should contain only 0s and 1s")) + end + return RAMMatrices(A, S, F, M, params, colnames) end ############################################################################################ @@ -165,83 +128,102 @@ function RAMMatrices( n_observed = length(partable.observed_vars) n_latent = length(partable.latent_vars) - n_node = n_observed + n_latent - - # F indices - F_ind = - length(partable.sorted_vars) != 0 ? - findall(∈(Set(partable.observed_vars)), partable.sorted_vars) : 1:n_observed - - # indices of the colnames - 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)) + n_vars = n_observed + n_latent + + if length(partable.sorted_vars) != 0 + @assert length(partable.sorted_vars) == nvars(partable) + vars_sorted = copy(partable.sorted_vars) + else + vars_sorted = [partable.observed_vars + partable.latent_vars] + end + + # indices of the vars (A/S/M rows or columns) + vars_index = Dict(col => i for (i, col) in enumerate(vars_sorted)) # fill Matrices # known_labels = Dict{Symbol, Int64}() - A_ind = [Vector{Int64}() for _ in 1:length(params)] - S_ind = [Vector{Int64}() for _ in 1:length(params)] - + T = nonmissingtype(eltype(partable.columns[:value_fixed])) + A_inds = [Vector{Int64}() for _ in 1:length(params)] + A_lin_ixs = LinearIndices((n_vars, n_vars)) + S_inds = [Vector{Int64}() for _ in 1:length(params)] + S_lin_ixs = LinearIndices((n_vars, n_vars)) + A_consts = Vector{Pair{Int, T}}() + S_consts = Vector{Pair{Int, T}}() # is there a meanstructure? - M_ind = + M_inds = any(==(Symbol("1")), partable.columns[:from]) ? [Vector{Int64}() for _ in 1:length(params)] : nothing - - # handle constants - constants = Vector{RAMConstant}() + M_consts = !isnothing(M_inds) ? Vector{Pair{Int, T}}() : nothing for r in partable - row_ind = col_indices[r.to] - col_ind = r.from != Symbol("1") ? col_indices[r.from] : nothing + row_ind = vars_index[r.to] + col_ind = r.from != Symbol("1") ? vars_index[r.from] : nothing if !r.free if (r.relation == :→) && (r.from == Symbol("1")) - push!(constants, RAMConstant(:M, row_ind, r.value_fixed)) + push!(M_consts, row_ind => r.value_fixed) elseif r.relation == :→ push!( - constants, - RAMConstant(:A, CartesianIndex(row_ind, col_ind), r.value_fixed), + A_consts, + A_lin_ixs[CartesianIndex(row_ind, col_ind)] => r.value_fixed, ) elseif r.relation == :↔ push!( - constants, - RAMConstant(:S, CartesianIndex(row_ind, col_ind), r.value_fixed), + S_consts, + S_lin_ixs[CartesianIndex(row_ind, col_ind)] => r.value_fixed, ) + if row_ind != col_ind # symmetric + push!( + S_consts, + S_lin_ixs[CartesianIndex(col_ind, row_ind)] => r.value_fixed, + ) + end else - error("Unsupported parameter type: $(r.relation)") + error("Unsupported relation: $(r.relation)") end else par_ind = params_index[r.param] if (r.relation == :→) && (r.from == Symbol("1")) - push!(M_ind[par_ind], row_ind) + push!(M_inds[par_ind], row_ind) elseif r.relation == :→ - push!(A_ind[par_ind], row_ind + (col_ind - 1) * n_node) + push!(A_inds[par_ind], A_lin_ixs[CartesianIndex(row_ind, col_ind)]) 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) + push!(S_inds[par_ind], S_lin_ixs[CartesianIndex(row_ind, col_ind)]) + if row_ind != col_ind # symmetric + push!(S_inds[par_ind], S_lin_ixs[CartesianIndex(col_ind, row_ind)]) end else - error("Unsupported parameter type: $(r.relation)") + error("Unsupported relation: $(r.relation)") end end end + # sort linear indices + for A_ind in A_inds + sort!(A_ind) + end + for S_ind in S_inds + unique!(sort!(S_ind)) # also symmetric duplicates + end + if !isnothing(M_inds) + for M_ind in M_inds + sort!(M_ind) + end + end + sort!(A_consts, by = first) + sort!(S_consts, by = first) + if !isnothing(M_consts) + sort!(M_consts, by = first) + end - return RAMMatrices( - A_ind, - S_ind, - F_ind, - M_ind, - params, - colnames, - constants, - (n_observed, n_node), - ) + return RAMMatrices(ParamsMatrix{T}(A_inds, A_consts, (n_vars, n_vars)), + ParamsMatrix{T}(S_inds, S_consts, (n_vars, n_vars)), + sparse(1:n_observed, + [vars_index[var] for var in partable.observed_vars], + ones(T, n_observed), n_observed, n_vars), + !isnothing(M_inds) ? ParamsVector{T}(M_inds, M_consts, (n_vars,)) : nothing, + params, vars_sorted) end Base.convert( @@ -255,21 +237,20 @@ Base.convert( ############################################################################################ function ParameterTable( - ram_matrices::RAMMatrices; + ram::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)] + + if !isnothing(ram.colnames) + latent_vars = SEM.latent_vars(ram) + observed_vars = SEM.observed_vars(ram) + colnames = ram.colnames 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)] + observed_vars = [Symbol("$(observed_var_prefix)_$i") for i in 1:nobserved_vars(ram)] + latent_vars = [Symbol("$(latent_var_prefix)_$i") for i in 1:nlatent_vars(ram)] colnames = vcat(observed_vars, latent_vars) end @@ -277,27 +258,16 @@ function ParameterTable( partable = ParameterTable( observed_vars = observed_vars, latent_vars = latent_vars, - params = isnothing(params) ? SEM.params(ram_matrices) : params, + params = isnothing(params) ? SEM.params(ram) : params, ) - # constants - for c in ram_matrices.constants - push!(partable, partable_row(c, colnames)) + # fill the table + append_rows!(partable, ram.S, :S, ram.params, colnames, skip_symmetric = true) + append_rows!(partable, ram.A, :A, ram.params, colnames) + if !isnothing(ram.M) + append_rows!(partable, ram.M, :M, ram.params, colnames) end - # parameters - for (i, par) in enumerate(ram_matrices.params) - append_partable_rows!( - partable, - colnames, - par, - i, - ram_matrices.A_ind, - ram_matrices.S_ind, - ram_matrices.M_ind, - ram_matrices.size_F[2], - ) - end check_params(SEM.params(partable), partable.columns[:param]) return partable @@ -339,23 +309,13 @@ function matrix_to_relation(matrix::Symbol) end end -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, -) - +# generates a ParTable row NamedTuple for a given element of RAM matrix function partable_row( - par::Symbol, - varnames::AbstractVector{Symbol}, - index::Integer, + val, + index, matrix::Symbol, - n_nod::Integer, + varnames::AbstractVector{Symbol}; + free::Bool = true, ) # variable names @@ -363,58 +323,65 @@ function partable_row( from = Symbol("1") to = varnames[index] else - cart_index = linear2cartesian(index, (n_nod, n_nod)) - - from = varnames[cart_index[2]] - to = varnames[cart_index[1]] + from = varnames[index[2]] + to = varnames[index[1]] end return ( from = from, relation = matrix_to_relation(matrix), to = to, - free = true, - value_fixed = 0.0, + free = free, + value_fixed = free ? 0.0 : val, start = 0.0, estimate = 0.0, - param = par, + param = free ? val : :const, ) end -function append_partable_rows!( +function append_rows!( partable::ParameterTable, - varnames::AbstractVector{Symbol}, - par::Symbol, - par_index::Integer, - A_ind, - S_ind, - M_ind, - n_nod::Integer, + arr::ParamsArray, + arr_name::Symbol, + params::AbstractVector, + varnames::AbstractVector{Symbol}; + skip_symmetric::Bool = false, ) - for ind in A_ind[par_index] - push!(partable, partable_row(par, varnames, ind, :A, n_nod)) - end + nparams(arr) == length(params) || throw( + ArgumentError( + "Length of parameters vector ($(length(params))) does not match the number of parameters in the matrix ($(nparams(arr)))", + ), + ) + arr_ixs = eachindex(arr) + + # add parameters + visited_indices = Set{eltype(arr_ixs)}() + for (i, par) in enumerate(params) + for j in param_occurences_range(arr, i) + arr_ix = arr_ixs[arr.linear_indices[j]] + skip_symmetric && (arr_ix ∈ visited_indices) && continue - 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), - ), + partable, + partable_row(par, arr_ix, arr_name, varnames, free = true), ) + if skip_symmetric + # mark index and its symmetric as visited + push!(visited_indices, arr_ix) + push!(visited_indices, CartesianIndex(arr_ix[2], arr_ix[1])) + end end end - if !isnothing(M_ind) - for ind in M_ind[par_index] - push!(partable, partable_row(par, varnames, ind, :M, n_nod)) + # add constants + for (i, val) in arr.constants + arr_ix = arr_ixs[i] + skip_symmetric && (arr_ix ∈ visited_indices) && continue + push!(partable, partable_row(val, arr_ix, arr_name, varnames, free = false)) + if skip_symmetric + # mark index and its symmetric as visited + push!(visited_indices, arr_ix) + push!(visited_indices, CartesianIndex(arr_ix[2], arr_ix[1])) end end @@ -423,14 +390,12 @@ end 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.A == mat2.A) && + (mat1.S == mat2.S) && + (mat1.F == mat2.F) && + (mat1.M == mat2.M) && (mat1.params == mat2.params) && - (mat1.colnames == mat2.colnames) && - (mat1.size_F == mat2.size_F) && - (mat1.constants == mat2.constants) + (mat1.colnames == mat2.colnames) ) return res end diff --git a/src/imply/RAM/generic.jl b/src/imply/RAM/generic.jl index e7e0b36f5..a16aac179 100644 --- a/src/imply/RAM/generic.jl +++ b/src/imply/RAM/generic.jl @@ -74,9 +74,6 @@ mutable struct RAM{ A5, A6, V2, - I1, - I2, - I3, M1, M2, M3, @@ -97,10 +94,6 @@ mutable struct RAM{ ram_matrices::V2 - A_indices::I1 - S_indices::I2 - M_indices::I3 - F⨉I_A⁻¹::M1 F⨉I_A⁻¹S::M2 I_A::M3 @@ -131,22 +124,14 @@ function RAM(; n_par = nparams(ram_matrices) n_obs = nobserved_vars(ram_matrices) n_var = nvars(ram_matrices) - F = zeros(ram_matrices.size_F) - F[CartesianIndex.(1:n_obs, ram_matrices.F_ind)] .= 1.0 - - # get indices - A_indices = copy(ram_matrices.A_ind) - S_indices = copy(ram_matrices.S_ind) - M_indices = !isnothing(ram_matrices.M_ind) ? copy(ram_matrices.M_ind) : nothing #preallocate arrays - A_pre = zeros(n_var, n_var) - S_pre = zeros(n_var, n_var) - M_pre = !isnothing(M_indices) ? zeros(n_var) : nothing - - set_RAMConstants!(A_pre, S_pre, M_pre, ram_matrices.constants) + nan_params = fill(NaN, n_par) + A_pre = materialize(ram_matrices.A, nan_params) + S_pre = materialize(ram_matrices.S, nan_params) + F = Matrix(ram_matrices.F) - A_pre = check_acyclic(A_pre, n_par, A_indices) + A_pre = check_acyclic(A_pre, ram_matrices.A) # pre-allocate some matrices Σ = zeros(n_obs, n_obs) @@ -155,8 +140,8 @@ function RAM(; I_A = similar(A_pre) if gradient_required - ∇A = matrix_gradient(A_indices, n_var^2) - ∇S = matrix_gradient(S_indices, n_var^2) + ∇A = sparse_gradient(ram_matrices.A) + ∇S = sparse_gradient(ram_matrices.S) else ∇A = nothing ∇S = nothing @@ -165,16 +150,16 @@ function RAM(; # μ if meanstructure MS = HasMeanStruct - !isnothing(M_indices) || throw( + !isnothing(ram_matrices.M) || throw( ArgumentError( "You set `meanstructure = true`, but your model specification contains no mean parameters.", ), ) - ∇M = gradient_required ? matrix_gradient(M_indices, n_var) : nothing + M_pre = materialize(ram_matrices.M, nan_params) + ∇M = gradient_required ? sparse_gradient(ram_matrices.M) : nothing μ = zeros(n_obs) else MS = NoMeanStruct - M_indices = nothing M_pre = nothing μ = nothing ∇M = nothing @@ -188,9 +173,6 @@ function RAM(; μ, M_pre, ram_matrices, - A_indices, - S_indices, - M_indices, F⨉I_A⁻¹, F⨉I_A⁻¹S, I_A, @@ -206,15 +188,11 @@ end ############################################################################################ function update!(targets::EvaluationTargets, imply::RAM, model::AbstractSemSingle, params) - fill_A_S_M!( - imply.A, - imply.S, - imply.M, - imply.A_indices, - imply.S_indices, - imply.M_indices, - params, - ) + materialize!(imply.A, imply.ram_matrices.A, params) + materialize!(imply.S, imply.ram_matrices.S, params) + if !isnothing(imply.M) + materialize!(imply.M, imply.ram_matrices.M, params) + end @. imply.I_A = -imply.A @view(imply.I_A[diagind(imply.I_A)]) .+= 1 @@ -251,12 +229,9 @@ end ### additional functions ############################################################################################ -function check_acyclic(A_pre, n_par, A_indices) - # fill copy of A-matrix with random parameters - A_rand = copy(A_pre) - randpar = rand(n_par) - - fill_matrix!(A_rand, A_indices, randpar) +function check_acyclic(A_pre::AbstractMatrix, A::ParamsMatrix) + # fill copy of A with random parameters + A_rand = materialize(A, rand(nparams(A))) # check if the model is acyclic acyclic = isone(det(I - A_rand)) diff --git a/src/imply/RAM/symbolic.jl b/src/imply/RAM/symbolic.jl index 9a96942ae..32ffcc068 100644 --- a/src/imply/RAM/symbolic.jl +++ b/src/imply/RAM/symbolic.jl @@ -102,24 +102,15 @@ function RAMSymbolic(; ram_matrices = convert(RAMMatrices, specification) n_par = nparams(ram_matrices) - n_obs = nobserved_vars(ram_matrices) - n_var = nvars(ram_matrices) - par = (Symbolics.@variables θ[1:n_par])[1] - A = zeros(Num, n_var, n_var) - S = zeros(Num, n_var, n_var) - !isnothing(ram_matrices.M_ind) ? M = zeros(Num, n_var) : M = nothing - F = zeros(ram_matrices.size_F) - F[CartesianIndex.(1:n_obs, ram_matrices.F_ind)] .= 1.0 - - set_RAMConstants!(A, S, M, ram_matrices.constants) - fill_A_S_M!(A, S, M, ram_matrices.A_ind, ram_matrices.S_ind, ram_matrices.M_ind, par) - - A, S, F = sparse(A), sparse(S), sparse(F) + A = sparse_materialize(Num, ram_matrices.A, par) + S = sparse_materialize(Num, ram_matrices.S, par) + M = !isnothing(ram_matrices.M) ? materialize(Num, ram_matrices.M, par) : nothing + F = ram_matrices.F - if !isnothing(loss_types) - any(loss_types .<: SemWLS) ? vech = true : nothing + if !isnothing(loss_types) && any(T -> T <: SemWLS, loss_types) + vech = true end I_A⁻¹ = neumann_series(A) From 81d0ab7f89bffb3e2792573a803dbc49abfc4629 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Fri, 22 Mar 2024 15:02:40 -0700 Subject: [PATCH 04/12] materialize!(Symm/LowTri/UpTri) --- src/additional_functions/params_array.jl | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/additional_functions/params_array.jl b/src/additional_functions/params_array.jl index f20a6518b..b79b6454f 100644 --- a/src/additional_functions/params_array.jl +++ b/src/additional_functions/params_array.jl @@ -135,6 +135,14 @@ materialize(::Type{T}, arr::ParamsArray, param_values::AbstractVector) where {T} materialize(arr::ParamsArray, param_values::AbstractVector{T}) where {T} = materialize(Union{T, eltype(arr)}, arr, param_values) +# the hack to update the structured matrix (should be fine since the structure is imposed by ParamsMatrix) +materialize!( + dest::Union{Symmetric, LowerTriangular, UpperTriangular}, + src::ParamsMatrix{<:Any}, + param_values::AbstractVector; + kwargs..., +) = materialize!(parent(dest), src, param_values; kwargs...) + function sparse_materialize( ::Type{T}, arr::ParamsMatrix, From fd13c740b6efd4c1d7406d4547ee93a0e4508acd Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Fri, 22 Mar 2024 16:13:34 -0700 Subject: [PATCH 05/12] ParamsArray: faster sparse materialize! --- src/additional_functions/params_array.jl | 91 ++++++++++++++++++----- src/frontend/specification/RAMMatrices.jl | 2 +- 2 files changed, 72 insertions(+), 21 deletions(-) diff --git a/src/additional_functions/params_array.jl b/src/additional_functions/params_array.jl index b79b6454f..13ae2eeaf 100644 --- a/src/additional_functions/params_array.jl +++ b/src/additional_functions/params_array.jl @@ -2,10 +2,14 @@ Array with partially parameterized elements. """ struct ParamsArray{T, N} <: AbstractArray{T, N} - linear_indices::Vector{Int} - param_ptr::Vector{Int} - constants::Vector{Pair{Int, T}} - size::NTuple{N, Int} + linear_indices::Vector{Int} # linear indices of the parameter refs in the destination array + nz_indices::Vector{Int} # indices of the parameters refs in nonzero elements vector + # (including the constants) ordered by the linear index + param_ptr::Vector{Int} # i-th element marks the start of the range in linear/nonzero + # indices arrays that corresponds to the i-th parameter + # (nparams + 1 elements) + constants::Vector{Tuple{Int, Int, T}} # linear index, index in nonzero vector, value + size::NTuple{N, Int} # size of the destination array end ParamsVector{T} = ParamsArray{T, 1} @@ -18,10 +22,16 @@ function ParamsArray{T, N}( ) where {T, N} params_ptr = pushfirst!(accumulate((ptr, inds) -> ptr + length(inds), params_map, init = 1), 1) + param_lin_inds = reduce(vcat, params_map, init = Vector{Int}()) + nz_lin_inds = unique!(sort!([param_lin_inds; first.(constants)])) + if length(nz_lin_inds) < length(param_lin_inds) + length(constants) + throw(ArgumentError("Duplicate linear indices in the parameterized array")) + end return ParamsArray{T, N}( - reduce(vcat, params_map, init = Vector{Int}()), + param_lin_inds, + searchsortedfirst.(Ref(nz_lin_inds), param_lin_inds), params_ptr, - constants, + [(c[1], searchsortedfirst(nz_lin_inds, c[1]), c[2]) for c in constants], size, ) end @@ -58,6 +68,7 @@ ParamsArray{T}( ) where {T, N} = ParamsArray{T, N}(arr, params; kwargs...) nparams(arr::ParamsArray) = length(arr.param_ptr) - 1 +SparseArrays.nnz(arr::ParamsArray) = length(arr.linear_indices) + length(arr.constants) Base.size(arr::ParamsArray) = arr.size Base.size(arr::ParamsArray, i::Integer) = arr.size[i] @@ -110,7 +121,7 @@ function materialize!( Z = eltype(dest) <: Number ? eltype(dest) : eltype(src) set_zeros && fill!(dest, zero(Z)) if set_constants - @inbounds for (i, val) in src.constants + @inbounds for (i, _, val) in src.constants dest[i] = val end end @@ -122,6 +133,43 @@ function materialize!( return dest end +function materialize!( + dest::SparseMatrixCSC, + src::ParamsMatrix, + param_values::AbstractVector; + set_constants::Bool = true, + set_zeros::Bool = false, +) + set_zeros && throw(ArgumentError("Cannot set zeros for sparse matrix")) + size(dest) == size(src) || throw( + DimensionMismatch( + "Parameters ($(size(params_arr))) and destination ($(size(dest))) array sizes don't match", + ), + ) + nparams(src) == length(param_values) || throw( + DimensionMismatch( + "Number of values ($(length(param_values))) does not match the number of parameters ($(nparams(src)))", + ), + ) + + nnz(dest) == nnz(src) || throw( + DimensionMismatch( + "Number of non-zero elements ($(nnz(dest))) does not match the number of parameter references and constants ($(nnz(src)))", + ), + ) + if set_constants + @inbounds for (_, j, val) in src.constants + dest.nzval[j] = val + end + end + @inbounds for (i, val) in enumerate(param_values) + for j in param_occurences_range(src, i) + dest.nzval[src.nz_indices[j]] = val + end + end + return dest +end + """ materialize([T], src::ParamsArray{<:Any, N}, param_values::AbstractVector{T}) where T @@ -150,27 +198,30 @@ function sparse_materialize( ) where {T} nparams(arr) == length(param_values) || throw( DimensionMismatch( - "Number of values ($(length(param)values))) does not match the number of parameter ($(nparams(arr)))", + "Number of values ($(length(param_values))) does not match the number of parameter ($(nparams(arr)))", ), ) - # constant values in sparse matrix - cvals = [T(v) for (_, v) in arr.constants] - # parameter values in sparse matrix - parvals = Vector{T}(undef, length(arr.linear_indices)) + + nz_vals = Vector{T}(undef, nnz(arr)) + nz_lininds = Vector{Int}(undef, nnz(arr)) + # fill constants + @inbounds for (lin_ind, nz_ind, val) in arr.constants + nz_vals[nz_ind] = val + nz_lininds[nz_ind] = lin_ind + end + # fill parameters @inbounds for (i, val) in enumerate(param_values) for j in param_occurences_range(arr, i) - parvals[j] = val + nz_ind = arr.nz_indices[j] + nz_vals[nz_ind] = val + nz_lininds[nz_ind] = arr.linear_indices[j] end end - nzixs = [first.(arr.constants); arr.linear_indices] - ixorder = sortperm(nzixs) - nzixs = nzixs[ixorder] - nzvals = [cvals; parvals][ixorder] arr_ixs = CartesianIndices(size(arr)) return sparse( - [arr_ixs[i][1] for i in nzixs], - [arr_ixs[i][2] for i in nzixs], - nzvals, + [arr_ixs[i][1] for i in nz_lininds], + [arr_ixs[i][2] for i in nz_lininds], + nz_vals, size(arr)..., ) end diff --git a/src/frontend/specification/RAMMatrices.jl b/src/frontend/specification/RAMMatrices.jl index b140ae026..ee487c25e 100644 --- a/src/frontend/specification/RAMMatrices.jl +++ b/src/frontend/specification/RAMMatrices.jl @@ -374,7 +374,7 @@ function append_rows!( end # add constants - for (i, val) in arr.constants + for (i, _, val) in arr.constants arr_ix = arr_ixs[i] skip_symmetric && (arr_ix ∈ visited_indices) && continue push!(partable, partable_row(val, arr_ix, arr_name, varnames, free = false)) From 19497d5bd7fbc8450ea861a80bb0d9c8321c8d2a Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Tue, 2 Jul 2024 17:23:27 -0700 Subject: [PATCH 06/12] ParamsArray: use Iterators.flatten() (faster) --- src/additional_functions/params_array.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/additional_functions/params_array.jl b/src/additional_functions/params_array.jl index 13ae2eeaf..bbee1dcf9 100644 --- a/src/additional_functions/params_array.jl +++ b/src/additional_functions/params_array.jl @@ -22,7 +22,7 @@ function ParamsArray{T, N}( ) where {T, N} params_ptr = pushfirst!(accumulate((ptr, inds) -> ptr + length(inds), params_map, init = 1), 1) - param_lin_inds = reduce(vcat, params_map, init = Vector{Int}()) + param_lin_inds = collect(Iterators.flatten(params_map)) nz_lin_inds = unique!(sort!([param_lin_inds; first.(constants)])) if length(nz_lin_inds) < length(param_lin_inds) + length(constants) throw(ArgumentError("Duplicate linear indices in the parameterized array")) From 139338d6884e45d898b6624cf43a0bbbe84d6eb7 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sun, 11 Aug 2024 13:07:50 -0700 Subject: [PATCH 07/12] Base.hash(::ParamsArray) --- src/additional_functions/params_array.jl | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/additional_functions/params_array.jl b/src/additional_functions/params_array.jl index bbee1dcf9..3a58171aa 100644 --- a/src/additional_functions/params_array.jl +++ b/src/additional_functions/params_array.jl @@ -79,6 +79,14 @@ Base.:(==)(a::ParamsArray, b::ParamsArray) = return eltype(a) == eltype(b) && a.param_ptr == b.param_ptr && a.linear_indices == b.linear_indices +Base.hash(a::ParamsArray, h::UInt) = hash( + typeof(a), + hash( + eltype(a), + hash(size(a), hash(a.constants, hash(a.param_ptr, hash(a.linear_indices, h)))), + ), +) + # the range of arr.param_ptr indices that correspond to i-th parameter param_occurences_range(arr::ParamsArray, i::Integer) = arr.param_ptr[i]:(arr.param_ptr[i+1]-1) From 58507840210ed7d7cc31c4f9c7d9f7036d1ca8cc Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sun, 11 Aug 2024 13:07:38 -0700 Subject: [PATCH 08/12] colnames -> vars --- .../start_val/start_fabin3.jl | 2 +- src/frontend/specification/RAMMatrices.jl | 45 +++++++++---------- src/frontend/specification/documentation.jl | 4 +- test/examples/multigroup/multigroup.jl | 4 +- .../political_democracy.jl | 4 +- .../recover_parameters_twofact.jl | 2 +- 6 files changed, 30 insertions(+), 31 deletions(-) diff --git a/src/additional_functions/start_val/start_fabin3.jl b/src/additional_functions/start_val/start_fabin3.jl index 9d692437e..d86b992da 100644 --- a/src/additional_functions/start_val/start_fabin3.jl +++ b/src/additional_functions/start_val/start_fabin3.jl @@ -129,7 +129,7 @@ function start_fabin3(ram_matrices::RAMMatrices, Σ, μ) # is there at least one reference indicator? if length(reference) > 0 if (length(reference) > 1) && any(((obs, param),) -> param > 0, indicators) # don't warn if entire column is fixed - @warn "You have more than 1 scaling indicator for $(ram_matrices.colnames[i])" + @warn "You have more than 1 scaling indicator for $(ram_matrices.vars[i])" end ref = reference[1] diff --git a/src/frontend/specification/RAMMatrices.jl b/src/frontend/specification/RAMMatrices.jl index ee487c25e..451f5fd69 100644 --- a/src/frontend/specification/RAMMatrices.jl +++ b/src/frontend/specification/RAMMatrices.jl @@ -9,7 +9,7 @@ struct RAMMatrices <: SemSpecification F::SparseMatrixCSC{Float64} M::Union{ParamsVector{Float64}, Nothing} params::Vector{Symbol} - colnames::Union{Vector{Symbol}, Nothing} # better call it "variables": it's a mixture of observed and latent (and it gets confusing with get_colnames()) + vars::Union{Vector{Symbol}, Nothing} # better call it "variables": it's a mixture of observed and latent (and it gets confusing with get_vars()) end nparams(ram::RAMMatrices) = nparams(ram.A) @@ -17,7 +17,7 @@ nvars(ram::RAMMatrices) = size(ram.F, 2) nobserved_vars(ram::RAMMatrices) = size(ram.F, 1) nlatent_vars(ram::RAMMatrices) = nvars(ram) - nobserved_vars(ram) -vars(ram::RAMMatrices) = ram.colnames +vars(ram::RAMMatrices) = ram.vars isobserved_var(ram::RAMMatrices, i::Integer) = ram.F.colptr[i+1] > ram.F.colptr[i] islatent_var(ram::RAMMatrices, i::Integer) = ram.F.colptr[i+1] == ram.F.colptr[i] @@ -34,13 +34,12 @@ function observed_var_indices(ram::RAMMatrices) return obs_inds end -latent_var_indices(ram::RAMMatrices) = - [i for i in axes(ram.F, 2) if islatent_var(ram, i)] +latent_var_indices(ram::RAMMatrices) = [i for i in axes(ram.F, 2) if islatent_var(ram, i)] # observed variables in the order as they appear in ram.F rows function observed_vars(ram::RAMMatrices) - if isnothing(ram.colnames) - @warn "Your RAMMatrices do not contain column names. Please make sure the order of variables in your data is correct!" + if isnothing(ram.vars) + @warn "Your RAMMatrices do not contain variable names. Please make sure the order of variables in your data is correct!" return nothing else obs_vars = Vector{Symbol}(undef, nobserved_vars(ram)) @@ -55,11 +54,11 @@ function observed_vars(ram::RAMMatrices) end function latent_vars(ram::RAMMatrices) - if isnothing(ram.colnames) - @warn "Your RAMMatrices do not contain column names. Please make sure the order of variables in your data is correct!" + if isnothing(ram.vars) + @warn "Your RAMMatrices do not contain variable names. Please make sure the order of variables in your data is correct!" return nothing else - return [col for (i, col) in enumerate(ram.colnames) if islatent_var(ram, i)] + return [col for (i, col) in enumerate(ram.vars) if islatent_var(ram, i)] end end @@ -73,32 +72,32 @@ function RAMMatrices(; F::AbstractMatrix, M::Union{AbstractVector, Nothing} = nothing, params::AbstractVector{Symbol}, - colnames::Union{AbstractVector{Symbol}, Nothing} = nothing, + vars::Union{AbstractVector{Symbol}, Nothing} = nothing, ) ncols = size(A, 2) - isnothing(colnames) || check_vars(colnames, ncols) + isnothing(vars) || check_vars(vars, 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", + "A should have as many rows and columns as vars 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", + "S should have as many rows and columns as vars 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", + "F should have as many columns as vars 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", + "M should have as many elements as vars length ($ncols), $(length(M)) found", ), ) end @@ -111,7 +110,7 @@ function RAMMatrices(; if any(!isone, spF.nzval) throw(ArgumentError("F should contain only 0s and 1s")) end - return RAMMatrices(A, S, F, M, params, colnames) + return RAMMatrices(A, S, F, M, params, vars) end ############################################################################################ @@ -244,14 +243,14 @@ function ParameterTable( ) # defer parameter checks until we know which ones are used - if !isnothing(ram.colnames) + if !isnothing(ram.vars) latent_vars = SEM.latent_vars(ram) observed_vars = SEM.observed_vars(ram) - colnames = ram.colnames + vars = ram.vars else observed_vars = [Symbol("$(observed_var_prefix)_$i") for i in 1:nobserved_vars(ram)] latent_vars = [Symbol("$(latent_var_prefix)_$i") for i in 1:nlatent_vars(ram)] - colnames = vcat(observed_vars, latent_vars) + vars = vcat(observed_vars, latent_vars) end # construct an empty table @@ -262,10 +261,10 @@ function ParameterTable( ) # fill the table - append_rows!(partable, ram.S, :S, ram.params, colnames, skip_symmetric = true) - append_rows!(partable, ram.A, :A, ram.params, colnames) + append_rows!(partable, ram.S, :S, ram.params, vars, skip_symmetric = true) + append_rows!(partable, ram.A, :A, ram.params, vars) if !isnothing(ram.M) - append_rows!(partable, ram.M, :M, ram.params, colnames) + append_rows!(partable, ram.M, :M, ram.params, vars) end check_params(SEM.params(partable), partable.columns[:param]) @@ -395,7 +394,7 @@ function Base.:(==)(mat1::RAMMatrices, mat2::RAMMatrices) (mat1.F == mat2.F) && (mat1.M == mat2.M) && (mat1.params == mat2.params) && - (mat1.colnames == mat2.colnames) + (mat1.vars == mat2.vars) ) return res end diff --git a/src/frontend/specification/documentation.jl b/src/frontend/specification/documentation.jl index e869dd43f..46135ead0 100644 --- a/src/frontend/specification/documentation.jl +++ b/src/frontend/specification/documentation.jl @@ -95,7 +95,7 @@ function EnsembleParameterTable end (1) RAMMatrices(partable::ParameterTable) - (2) RAMMatrices(;A, S, F, M = nothing, params, colnames) + (2) RAMMatrices(;A, S, F, M = nothing, params, vars) (3) RAMMatrices(partable::EnsembleParameterTable) @@ -110,7 +110,7 @@ Return `RAMMatrices` constructed from (1) a parameter table or (2) individual ma - `F`: filter matrix - `M`: vector of mean effects - `params::Vector{Symbol}`: parameter labels -- `colnames::Vector{Symbol}`: variable names corresponding to the A, S and F matrix columns +- `vars::Vector{Symbol}`: variable names corresponding to the A, S and F matrix columns # Examples See the online documentation on [Model specification](@ref) and the [RAMMatrices interface](@ref). diff --git a/test/examples/multigroup/multigroup.jl b/test/examples/multigroup/multigroup.jl index a2f277d91..caaa5c3f7 100644 --- a/test/examples/multigroup/multigroup.jl +++ b/test/examples/multigroup/multigroup.jl @@ -60,7 +60,7 @@ specification_g1 = RAMMatrices(; S = S1, F = F, params = x, - colnames = [:x1, :x2, :x3, :x4, :x5, :x6, :x7, :x8, :x9, :visual, :textual, :speed], + vars = [:x1, :x2, :x3, :x4, :x5, :x6, :x7, :x8, :x9, :visual, :textual, :speed], ) specification_g2 = RAMMatrices(; @@ -68,7 +68,7 @@ specification_g2 = RAMMatrices(; S = S2, F = F, params = x, - colnames = [:x1, :x2, :x3, :x4, :x5, :x6, :x7, :x8, :x9, :visual, :textual, :speed], + vars = [:x1, :x2, :x3, :x4, :x5, :x6, :x7, :x8, :x9, :visual, :textual, :speed], ) partable = EnsembleParameterTable( diff --git a/test/examples/political_democracy/political_democracy.jl b/test/examples/political_democracy/political_democracy.jl index d7fbb8f2c..2f570302a 100644 --- a/test/examples/political_democracy/political_democracy.jl +++ b/test/examples/political_democracy/political_democracy.jl @@ -76,7 +76,7 @@ spec = RAMMatrices(; S = S, F = F, params = x, - colnames = [ + vars = [ :x1, :x2, :x3, @@ -108,7 +108,7 @@ spec_mean = RAMMatrices(; F = F, M = M, params = x, - colnames = [ + vars = [ :x1, :x2, :x3, diff --git a/test/examples/recover_parameters/recover_parameters_twofact.jl b/test/examples/recover_parameters/recover_parameters_twofact.jl index f00187fac..89c1225e2 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, params = x, colnames = nothing) +ram_matrices = RAMMatrices(; A = A, S = S, F = F, params = x, vars = nothing) true_val = [ repeat([1], 8) From 0f747b7ec789d5aceb4eed3b2b4612de2f7c60f8 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Wed, 3 Apr 2024 22:50:18 -0700 Subject: [PATCH 09/12] update_partable!(): better params unique check --- src/frontend/specification/ParameterTable.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/frontend/specification/ParameterTable.jl b/src/frontend/specification/ParameterTable.jl index df2cc165b..05350fb12 100644 --- a/src/frontend/specification/ParameterTable.jl +++ b/src/frontend/specification/ParameterTable.jl @@ -309,10 +309,10 @@ function update_partable!( "The length of `params` ($(length(params))) and their `values` ($(length(values))) must be the same", ), ) + dup_params = nonunique(params) + isempty(dup_params) || + throw(ArgumentError("Duplicate parameters detected: $(join(dup_params, ", "))")) 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 From aa34d5352979399de3f40226bf05b5690b777b9f Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Wed, 31 Jul 2024 21:35:50 -0700 Subject: [PATCH 10/12] start_fabin3: check obs_mean data & meanstructure --- src/additional_functions/start_val/start_fabin3.jl | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/src/additional_functions/start_val/start_fabin3.jl b/src/additional_functions/start_val/start_fabin3.jl index d86b992da..53cf7cff6 100644 --- a/src/additional_functions/start_val/start_fabin3.jl +++ b/src/additional_functions/start_val/start_fabin3.jl @@ -30,13 +30,21 @@ function start_fabin3(observed::SemObservedMissing, imply, optimizer, args...; k return start_fabin3(imply.ram_matrices, observed.em_model.Σ, observed.em_model.μ) end -function start_fabin3(ram_matrices::RAMMatrices, Σ, μ) +function start_fabin3( + ram_matrices::RAMMatrices, + Σ::AbstractMatrix, + μ::Union{AbstractVector, Nothing}, +) A, S, F, M, n_par = ram_matrices.A, ram_matrices.S, ram_matrices.F, ram_matrices.M, nparams(ram_matrices) + if !isnothing(M) && isnothing(μ) + throw(ArgumentError("RAM has meanstructure, but no observed means provided.")) + end + start_val = zeros(n_par) F_var2obs = Dict( i => F.rowval[F.colptr[i]] for i in axes(F, 2) if isobserved_var(ram_matrices, i) From 13aacd0cf87c1f8946f97bc27eb2cfe57f2e21cd Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Fri, 22 Nov 2024 11:28:11 -0800 Subject: [PATCH 11/12] params/vars API tweaks and tests --- .../start_val/start_partable.jl | 31 +++++++------------ src/frontend/specification/ParameterTable.jl | 4 +-- src/frontend/specification/RAMMatrices.jl | 2 +- src/frontend/specification/StenoGraphs.jl | 2 +- test/examples/helper.jl | 2 ++ test/examples/multigroup/build_models.jl | 2 ++ .../political_democracy/constructor.jl | 1 + .../political_democracy.jl | 12 +++++-- 8 files changed, 29 insertions(+), 27 deletions(-) diff --git a/src/additional_functions/start_val/start_partable.jl b/src/additional_functions/start_val/start_partable.jl index 6fb15e365..15f863f5b 100644 --- a/src/additional_functions/start_val/start_partable.jl +++ b/src/additional_functions/start_val/start_partable.jl @@ -21,28 +21,19 @@ function start_parameter_table(observed, imply, optimizer, args...; kwargs...) return start_parameter_table(ram_matrices(imply); kwargs...) end -function start_parameter_table( - ram_matrices::RAMMatrices; - parameter_table::ParameterTable, - kwargs..., -) +function start_parameter_table(ram::RAMMatrices; partable::ParameterTable, kwargs...) start_val = zeros(0) - for param in ram_matrices.params - found = false - for (i, param_table) in enumerate(parameter_table.params) - if param == param_table - push!(start_val, parameter_table.start[i]) - found = true - break - end - end - if !found - throw( - ErrorException( - "At least one parameter could not be found in the parameter table.", - ), - ) + param_indices = Dict(param => i for (i, param) in enumerate(params(ram))) + start_col = partable.columns[:start] + + for (i, param) in enumerate(partable.columns[:param]) + par_ind = get(param_indices, param, nothing) + if !isnothing(par_ind) + par_start = start_col[i] + isfinite(par_start) && (start_val[i] = par_start) + else + throw(ErrorException("Parameter $(param) is not in the parameter table.")) end end diff --git a/src/frontend/specification/ParameterTable.jl b/src/frontend/specification/ParameterTable.jl index 05350fb12..8b7cc0973 100644 --- a/src/frontend/specification/ParameterTable.jl +++ b/src/frontend/specification/ParameterTable.jl @@ -54,8 +54,8 @@ function ParameterTable( 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), + observed_vars = copy(observed_vars(partable)), + latent_vars = copy(latent_vars(partable)), params = params, ) end diff --git a/src/frontend/specification/RAMMatrices.jl b/src/frontend/specification/RAMMatrices.jl index 451f5fd69..0c5722f57 100644 --- a/src/frontend/specification/RAMMatrices.jl +++ b/src/frontend/specification/RAMMatrices.jl @@ -110,7 +110,7 @@ function RAMMatrices(; if any(!isone, spF.nzval) throw(ArgumentError("F should contain only 0s and 1s")) end - return RAMMatrices(A, S, F, M, params, vars) + return RAMMatrices(A, S, F, M, copy(params), vars) end ############################################################################################ diff --git a/src/frontend/specification/StenoGraphs.jl b/src/frontend/specification/StenoGraphs.jl index 64a33f13e..5cf87c07a 100644 --- a/src/frontend/specification/StenoGraphs.jl +++ b/src/frontend/specification/StenoGraphs.jl @@ -42,7 +42,7 @@ function ParameterTable( latent_vars::AbstractVector{Symbol}, params::Union{AbstractVector{Symbol}, Nothing} = nothing, group::Union{Integer, Nothing} = nothing, - param_prefix = :θ, + param_prefix::Symbol = :θ, ) graph = unique(graph) n = length(graph) diff --git a/test/examples/helper.jl b/test/examples/helper.jl index d4c140d67..042f7005f 100644 --- a/test/examples/helper.jl +++ b/test/examples/helper.jl @@ -1,6 +1,8 @@ using LinearAlgebra: norm function test_gradient(model, params; rtol = 1e-10, atol = 0) + @test nparams(model) == length(params) + true_grad = FiniteDiff.finite_difference_gradient(Base.Fix1(objective!, model), params) gradient = similar(params) diff --git a/test/examples/multigroup/build_models.jl b/test/examples/multigroup/build_models.jl index 3f29a6898..6991dd479 100644 --- a/test/examples/multigroup/build_models.jl +++ b/test/examples/multigroup/build_models.jl @@ -8,6 +8,8 @@ model_g1 = Sem(specification = specification_g1, data = dat_g1, imply = RAMSymbo model_g2 = Sem(specification = specification_g2, data = dat_g2, imply = RAM) +@test SEM.params(model_g1.imply.ram_matrices) == SEM.params(model_g2.imply.ram_matrices) + # test the different constructors model_ml_multigroup = SemEnsemble(model_g1, model_g2; optimizer = semoptimizer) model_ml_multigroup2 = SemEnsemble( diff --git a/test/examples/political_democracy/constructor.jl b/test/examples/political_democracy/constructor.jl index 99ef06b3a..bebabf6e0 100644 --- a/test/examples/political_democracy/constructor.jl +++ b/test/examples/political_democracy/constructor.jl @@ -6,6 +6,7 @@ using Random ############################################################################################ model_ml = Sem(specification = spec, data = dat, optimizer = semoptimizer) +@test SEM.params(model_ml.imply.ram_matrices) == SEM.params(spec) model_ml_cov = Sem( specification = spec, diff --git a/test/examples/political_democracy/political_democracy.jl b/test/examples/political_democracy/political_democracy.jl index 2f570302a..2265e2a59 100644 --- a/test/examples/political_democracy/political_democracy.jl +++ b/test/examples/political_democracy/political_democracy.jl @@ -1,5 +1,7 @@ using StructuralEquationModels, Test, FiniteDiff +SEM = StructuralEquationModels + include( joinpath( chop(dirname(pathof(StructuralEquationModels)), tail = 3), @@ -96,9 +98,9 @@ spec = RAMMatrices(; partable = ParameterTable(spec) -# w. meanstructure ------------------------------------------------------------------------- +@test SEM.params(spec) == SEM.params(partable) -x = Symbol.("x" .* string.(1:38)) +# w. meanstructure ------------------------------------------------------------------------- M = [:x32; :x33; :x34; :x35; :x36; :x37; :x38; :x35; :x36; :x37; :x38; 0.0; 0.0; 0.0] @@ -107,7 +109,7 @@ spec_mean = RAMMatrices(; S = S, F = F, M = M, - params = x, + params = [SEM.params(spec); Symbol.("x", string.(32:38))], vars = [ :x1, :x2, @@ -128,6 +130,8 @@ spec_mean = RAMMatrices(; partable_mean = ParameterTable(spec_mean) +@test SEM.params(partable_mean) == SEM.params(spec_mean) + start_test = [fill(1.0, 11); fill(0.05, 3); fill(0.05, 6); fill(0.5, 8); fill(0.05, 3)] start_test_mean = [fill(1.0, 11); fill(0.05, 3); fill(0.05, 6); fill(0.5, 8); fill(0.05, 3); fill(0.1, 7)] @@ -164,6 +168,8 @@ end spec = ParameterTable(spec) spec_mean = ParameterTable(spec_mean) +@test SEM.params(spec) == SEM.params(partable) + partable = spec partable_mean = spec_mean From 8c26c351c5f440df66bc6cb78a92b43007b66c00 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Fri, 22 Mar 2024 15:10:43 -0700 Subject: [PATCH 12/12] generic imply: keep F sparse --- src/imply/RAM/generic.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/imply/RAM/generic.jl b/src/imply/RAM/generic.jl index a16aac179..850934a9c 100644 --- a/src/imply/RAM/generic.jl +++ b/src/imply/RAM/generic.jl @@ -129,7 +129,7 @@ function RAM(; nan_params = fill(NaN, n_par) A_pre = materialize(ram_matrices.A, nan_params) S_pre = materialize(ram_matrices.S, nan_params) - F = Matrix(ram_matrices.F) + F = copy(ram_matrices.F) A_pre = check_acyclic(A_pre, ram_matrices.A)