Refactor objective/gradient/hessian evaluation

src/StructuralEquationModels.jl

@@ -95,6 +95,12 @@ export AbstractSem,
     Sem,
     SemFiniteDiff,
     SemEnsemble,
+    MeanStruct,
+    NoMeanStruct,
+    HasMeanStruct,
+    HessianEval,
+    ExactHessian,
+    ApproxHessian,
     SemImply,
     RAMSymbolic,
     RAM,

src/additional_functions/helper.jl

@@ -33,15 +33,7 @@
     return sem_vec
 end
 
-function get_observed(rowind, data, semobserved; args = (), kwargs = NamedTuple())
-    observed_vec = Vector{semobserved}(undef, length(rowind))
-    for i in 1:length(rowind)
-        observed_vec[i] = semobserved(args...; data = Matrix(data[rowind[i], :]), kwargs...)
-    end
-    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)

src/frontend/fit/standard_errors/hessian.jl

@@ -1,58 +1,54 @@
 """
-    se_hessian(semfit::SemFit; hessian = :finitediff)
+    se_hessian(fit::SemFit; method = :finitediff)
 
-Return hessian based standard errors.
+Return hessian-based standard errors.
 
 # Arguments
-- `hessian`: how to compute the hessian. Options are
+- `method`: how to compute the hessian. Options are
     - `:analytic`: (only if an analytic hessian for the model can be computed)
     - `:finitediff`: for finite difference approximation
 """
-function se_hessian(sem_fit::SemFit; hessian = :finitediff)
-    c = H_scaling(sem_fit.model)
-
-    if hessian == :analytic
-        par = solution(sem_fit)
-        H = zeros(eltype(par), length(par), length(par))
-        hessian!(H, sem_fit.model, sem_fit.solution)
-    elseif hessian == :finitediff
-        H = FiniteDiff.finite_difference_hessian(
-            Base.Fix1(objective!, sem_fit.model),
-            sem_fit.solution,
-        )
-    elseif hessian == :optimizer
-        throw(
-            ArgumentError(
-                "standard errors from the optimizer hessian are not implemented yet",
-            ),
-        )
-    elseif hessian == :expected
-        throw(
-            ArgumentError(
-                "standard errors based on the expected hessian are not implemented yet",
-            ),
+function se_hessian(fit::SemFit; method = :finitediff)
+    c = H_scaling(fit.model)
+    params = solution(fit)
+    H = similar(params, (length(params), length(params)))
+
+    if method == :analytic
+        evaluate!(nothing, nothing, H, fit.model, params)
+    elseif method == :finitediff
+        FiniteDiff.finite_difference_hessian!(
+            H,
+            p -> evaluate!(zero(eltype(H)), nothing, nothing, fit.model, p),
+            params,
         )
+    elseif method == :optimizer
+        error("Standard errors from the optimizer hessian are not implemented yet")
+    elseif method == :expected
+        error("Standard errors based on the expected hessian are not implemented yet")
     else
-        throw(ArgumentError("I don't know how to compute `$hessian` standard-errors"))
+        throw(ArgumentError("Unsupported hessian calculation method :$method"))
     end
 
-    invH = c * inv(H)
-    se = sqrt.(diag(invH))
-
-    return se
+    H_chol = cholesky!(Symmetric(H))
+    H_inv = LinearAlgebra.inv!(H_chol)
+    return [sqrt(c * H_inv[i]) for i in diagind(H_inv)]
 end
 
 # Addition functions -------------------------------------------------------------
-H_scaling(model::AbstractSemSingle) =
-    H_scaling(model, model.observed, model.imply, model.optimizer, model.loss.functions...)
+function H_scaling(model::AbstractSemSingle)
+    if length(model.loss.functions) > 1
+        @warn "Hessian scaling for multiple loss functions is not implemented yet"
+    end
+    return H_scaling(model.loss.functions[1], model)
+end
 
-H_scaling(model, obs, imp, optimizer, lossfun::SemML) = 2 / (nsamples(model) - 1)
+H_scaling(lossfun::SemML, model::AbstractSemSingle) = 2 / (nsamples(model) - 1)
 
-function H_scaling(model, obs, imp, optimizer, lossfun::SemWLS)
+function H_scaling(lossfun::SemWLS, model::AbstractSemSingle)
     @warn "Standard errors for WLS are only correct if a GLS weight matrix (the default) is used."
     return 2 / (nsamples(model) - 1)
 end
 
-H_scaling(model, obs, imp, optimizer, lossfun::SemFIML) = 2 / nsamples(model)
+H_scaling(lossfun::SemFIML, model::AbstractSemSingle) = 2 / nsamples(model)
 
 H_scaling(model::SemEnsemble) = 2 / nsamples(model)

src/frontend/specification/ParameterTable.jl

@@ -362,22 +362,22 @@ end
     update_se_hessian!(
         partable::AbstractParameterTable,
         fit::SemFit;
-        hessian = :finitediff)
+        method = :finitediff)
 
 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.
+- `method::Symbol`: how to compute the hessian, see [se_hessian](@ref) for more information.
 
 # Examples
 
 """
 function update_se_hessian!(
     partable::AbstractParameterTable,
     fit::SemFit;
-    hessian = :finitediff,
+    method = :finitediff,
 )
-    se = se_hessian(fit; hessian = hessian)
+    se = se_hessian(fit; method)
     return update_partable!(partable, :se, params(fit), se)
 end
 

src/frontend/specification/StenoGraphs.jl

@@ -73,12 +73,9 @@ function ParameterTable(
             )
         end
         if element isa ModifiedEdge
-            if any(Base.Fix2(isa, Fixed), values(element.modifiers)) & any(Base.Fix2(isa, Label), values(element.modifiers))
-                throw(
-                    ArgumentError(
-                        "It is not allowed to label fixed parameters."
-                    )
-                )
+            if any(Base.Fix2(isa, Fixed), values(element.modifiers)) &&
+               any(Base.Fix2(isa, Label), values(element.modifiers))
+                throw(ArgumentError("It is not allowed to label fixed parameters."))
             end
             for modifier in values(element.modifiers)
                 if isnothing(group) &&