Merge pull request #198 from alyst/fixes_2

Misc.Fixes 2

Author: Maximilian-Stefan-Ernst

src/StructuralEquationModels.jl

@@ -142,6 +142,7 @@ export AbstractSem,
     update_partable!,
     update_estimate!,
     update_start!,
+    update_se_hessian!,
     Fixed,
     fixed,
     Start,

src/additional_functions/helper.jl

@@ -1,10 +1,14 @@
-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
@@ -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

src/additional_functions/parameters.jl

@@ -86,15 +86,19 @@ 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
 
 # fill M with parameters
@@ -111,97 +115,23 @@ function fill_matrix!(
     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