Merge pull request #200 from StructuralEquationModels/devel

Release 0.2.4

Author: Maximilian-Stefan-Ernst

.github/workflows/CI.yml

@@ -3,6 +3,8 @@ on:
   push:
     branches: '*'
     tags: '*'
+  pull_request:
+    branches: '*'
   workflow_dispatch:
 jobs:
   test:

.github/workflows/FormatCheck.yml

@@ -1,6 +1,6 @@
 name: Format suggestions
 on:
-  pull_request:
+  pull_request_target:
     # this argument is not required if you don't use the `suggestion-label` input
     types: [ opened, reopened, synchronize, labeled, unlabeled ]
 jobs:

Project.toml

@@ -1,7 +1,7 @@
 name = "StructuralEquationModels"
 uuid = "383ca8c5-e4ff-4104-b0a9-f7b279deed53"
 authors = ["Maximilian Ernst", "Aaron Peikert"]
-version = "0.2.3"
+version = "0.2.4"
 
 [deps]
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"

src/StructuralEquationModels.jl

@@ -4,6 +4,7 @@ using LinearAlgebra,
     Optim,
     NLSolversBase,
     Statistics,
+    StatsBase,
     SparseArrays,
     Symbolics,
     NLopt,
@@ -18,6 +19,8 @@ using LinearAlgebra,
 import DataFrames: DataFrame
 export StenoGraphs, @StenoGraph, meld
 
+const SEM = StructuralEquationModels
+
 # type hierarchy
 include("types.jl")
 include("objective_gradient_hessian.jl")
@@ -139,6 +142,7 @@ export AbstractSem,
     update_partable!,
     update_estimate!,
     update_start!,
+    update_se_hessian!,
     Fixed,
     fixed,
     Start,

src/additional_functions/helper.jl

@@ -1,16 +1,20 @@
-function neumann_series(mat::SparseMatrixCSC)
+# Neumann series representation of (I - mat)⁻¹
+function neumann_series(mat::SparseMatrixCSC; maxiter::Integer = size(mat, 1))
     inverse = I + mat
     next_term = mat^2
 
+    n = 1
     while nnz(next_term) != 0
+        (n <= maxiter) || error("Neumann series did not converge in $maxiter steps")
         inverse += next_term
         next_term *= mat
+        n += 1
     end
 
     return inverse
 end
 
-#= 
+#=
 function make_onelement_array(A)
     isa(A, Array) ? nothing : (A = [A])
     return A
@@ -37,13 +41,8 @@ function get_observed(rowind, data, semobserved; args = (), kwargs = NamedTuple(
     return observed_vec
 end
 
-function skipmissing_mean(mat)
-    means = Vector{Float64}(undef, size(mat, 2))
-    for i in 1:size(mat, 2)
-        @views means[i] = mean(skipmissing(mat[:, i]))
-    end
-    return means
-end
+skipmissing_mean(mat::AbstractMatrix) = 
+    [mean(skipmissing(coldata)) for coldata in eachcol(mat)]
 
 function F_one_person(imp_mean, meandiff, inverse, data, logdet)
     F = logdet
@@ -52,10 +51,10 @@ function F_one_person(imp_mean, meandiff, inverse, data, logdet)
     return F
 end
 
-function remove_all_missing(data)
+function remove_all_missing(data::AbstractMatrix)
     keep = Vector{Int64}()
-    for i in 1:size(data, 1)
-        if any(.!ismissing.(data[i, :]))
+    for (i, coldata) in zip(axes(data, 1), eachrow(data))
+        if any(!ismissing, coldata)
             push!(keep, i)
         end
     end
@@ -108,11 +107,8 @@ function sparse_outer_mul!(C, A, B::Vector, ind) #computes A*S*B -> C, where ind
 end
 
 function cov_and_mean(rows; corrected = false)
-    data = transpose(reduce(hcat, rows))
-    size(rows, 1) > 1 ? obs_cov = Statistics.cov(data; corrected = corrected) :
-    obs_cov = reshape([0.0], 1, 1)
-    obs_mean = vec(Statistics.mean(data, dims = 1))
-    return obs_cov, obs_mean
+    obs_mean, obs_cov = StatsBase.mean_and_cov(reduce(hcat, rows), 2, corrected = corrected)
+    return obs_cov, vec(obs_mean)
 end
 
 function duplication_matrix(nobs)

src/additional_functions/parameters.jl

@@ -1,5 +1,14 @@
-function fill_A_S_M(A, S, M, A_indices, S_indices, M_indices, parameters)
-    for (iA, iS, par) in zip(A_indices, S_indices, parameters)
+# fill A, S, and M matrices with the parameter values according to the parameters map
+function fill_A_S_M!(
+    A::AbstractMatrix,
+    S::AbstractMatrix,
+    M::Union{AbstractVector, Nothing},
+    A_indices::AbstractArrayParamsMap,
+    S_indices::AbstractArrayParamsMap,
+    M_indices::Union{AbstractArrayParamsMap, Nothing},
+    parameters::AbstractVector,
+)
+    @inbounds for (iA, iS, par) in zip(A_indices, S_indices, parameters)
         for index_A in iA
             A[index_A] = par
         end
@@ -10,22 +19,28 @@ function fill_A_S_M(A, S, M, A_indices, S_indices, M_indices, parameters)
     end
 
     if !isnothing(M)
-        for (iM, par) in zip(M_indices, parameters)
+        @inbounds for (iM, par) in zip(M_indices, parameters)
             for index_M in iM
                 M[index_M] = par
             end
         end
     end
 end
 
-function get_parameter_indices(parameters, M; linear = true, kwargs...)
-    M_indices = [findall(x -> (x == par), M) for par in parameters]
-
-    if linear
-        M_indices = cartesian2linear.(M_indices, [M])
+# build the map from the index of the parameter to the linear indices
+# of this parameter occurences in M
+# returns ArrayParamsMap object
+function array_parameters_map(parameters::AbstractVector, M::AbstractArray)
+    params_index = Dict(param => i for (i, param) in enumerate(parameters))
+    T = Base.eltype(eachindex(M))
+    res = [Vector{T}() for _ in eachindex(parameters)]
+    for (i, val) in enumerate(M)
+        par_ind = get(params_index, val, nothing)
+        if !isnothing(par_ind)
+            push!(res[par_ind], i)
+        end
     end
-
-    return M_indices
+    return res
 end
 
 function eachindex_lower(M; linear_indices = false, kwargs...)
@@ -49,9 +64,6 @@ function linear2cartesian(ind_lin, dims)
     return ind_cart
 end
 
-cartesian2linear(ind_cart, A::AbstractArray) = cartesian2linear(ind_cart, size(A))
-linear2cartesian(ind_linear, A::AbstractArray) = linear2cartesian(ind_linear, size(A))
-
 function set_constants!(M, M_pre)
     for index in eachindex(M)
         δ = tryparse(Float64, string(M[index]))
@@ -74,116 +86,52 @@ function check_constants(M)
     return false
 end
 
-function get_matrix_derivative(M_indices, parameters, n_long)
-    ∇M = [
-        sparsevec(M_indices[i], ones(length(M_indices[i])), n_long) for
-        i in 1:length(parameters)
-    ]
-
-    ∇M = reduce(hcat, ∇M)
-
-    return ∇M
+# construct length(M)×length(parameters) sparse matrix of 1s at the positions,
+# where the corresponding parameter occurs in the M matrix
+function matrix_gradient(M_indices::ArrayParamsMap, M_length::Integer)
+    rowval = reduce(vcat, M_indices)
+    colptr =
+        pushfirst!(accumulate((ptr, M_ind) -> ptr + length(M_ind), M_indices, init = 1), 1)
+    return SparseMatrixCSC(
+        M_length,
+        length(M_indices),
+        colptr,
+        rowval,
+        ones(length(rowval)),
+    )
 end
 
-function fill_matrix(M, M_indices, parameters)
+# fill M with parameters
+function fill_matrix!(
+    M::AbstractMatrix,
+    M_indices::AbstractArrayParamsMap,
+    parameters::AbstractVector,
+)
     for (iM, par) in zip(M_indices, parameters)
         for index_M in iM
             M[index_M] = par
         end
     end
+    return M
 end
 
-function get_partition(A_indices, S_indices)
-    n_par = length(A_indices)
-
-    first_A = "a"
-    first_S = "a"
-    last_A = "a"
-    last_S = "a"
-
-    for i in 1:n_par
-        if length(A_indices[i]) != 0
-            first_A = i
-            break
-        end
-    end
-
-    for i in 1:n_par
-        if length(S_indices[i]) != 0
-            first_S = i
-            break
-        end
-    end
-
-    for i in n_par + 1 .- (1:n_par)
-        if length(A_indices[i]) != 0
-            last_A = i
-            break
-        end
-    end
-
-    for i in n_par + 1 .- (1:n_par)
-        if length(S_indices[i]) != 0
-            last_S = i
-            break
-        end
-    end
-
-    for i in first_A:last_A
-        if length(A_indices[i]) == 0
-            throw(
-                ErrorException(
-                    "Your parameter vector is not partitioned into directed and undirected effects",
-                ),
-            )
-            return nothing
-        end
-    end
-
-    for i in first_S:last_S
-        if length(S_indices[i]) == 0
-            throw(
-                ErrorException(
-                    "Your parameter vector is not partitioned into directed and undirected effects",
-                ),
-            )
-            return nothing
-        end
-    end
-
-    return first_A:last_A, first_S:last_S
-end
-
-function get_partition(M_indices)
-    n_par = length(M_indices)
-
-    first_M = "a"
-    last_M = "a"
-
-    for i in 1:n_par
-        if length(M_indices[i]) != 0
-            first_M = i
-            break
-        end
-    end
-
-    for i in n_par + 1 .- (1:n_par)
-        if length(M_indices[i]) != 0
-            last_M = i
-            break
-        end
-    end
-
-    for i in first_M:last_M
-        if length(M_indices[i]) == 0
-            throw(
-                ErrorException(
-                    "Your parameter vector is not partitioned into directed, undirected and mean effects",
-                ),
-            )
-            return nothing
+# range of parameters that are referenced in the matrix
+function param_range(mtx_indices::AbstractArrayParamsMap)
+    first_i = findfirst(!isempty, mtx_indices)
+    last_i = findlast(!isempty, mtx_indices)
+
+    if !isnothing(first_i) && !isnothing(last_i)
+        for i in first_i:last_i
+            if isempty(mtx_indices[i])
+                # TODO show which parameter is missing in which matrix
+                throw(
+                    ErrorException(
+                        "Your parameter vector is not partitioned into directed and undirected effects",
+                    ),
+                )
+            end
         end
     end
 
-    return first_M:last_M
+    return first_i:last_i
 end