Merge pull request #203 from alyst/refactor_semspec

Refactor SEM Specification

Author: Maximilian-Stefan-Ernst

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"

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.

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.

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

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

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,

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]`.

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)
 ```

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:

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
 )