Multiple refactorings reformatted

Project.toml

@@ -13,8 +13,10 @@ LineSearches = "d3d80556-e9d4-5f37-9878-2ab0fcc64255"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 NLSolversBase = "d41bc354-129a-5804-8e4c-c37616107c6c"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
+PackageExtensionCompat = "65ce6f38-6b18-4e1d-a461-8949797d7930"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
+Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
@@ -46,9 +48,12 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 test = ["Test"]
 
 [weakdeps]
+BlackBoxOptim = "a134a8b2-14d6-55f6-9291-3336d3ab0209"
 NLopt = "76087f3c-5699-56af-9a33-bf431cd00edd"
+Optimisers = "3bd65402-5787-11e9-1adc-39752487f4e2"
 ProximalAlgorithms = "140ffc9f-1907-541a-a177-7475e0a401e9"
 
 [extensions]
 SEMNLOptExt = "NLopt"
 SEMProximalOptExt = "ProximalAlgorithms"
+SEMBlackBoxOptimExt = ["BlackBoxOptim", "Optimisers"]

docs/src/tutorials/concept.md

@@ -49,8 +49,8 @@ Available loss functions are
 - [`SemRidge`](@ref): ridge regularization
 
 ## The optimizer part aka `SemOptimizer`
-The optimizer part of a model connects to the numerical optimization backend used to fit the model. 
-It can be used to control options like the optimization algorithm, linesearch, stopping criteria, etc. 
+The optimizer part of a model connects to the numerical optimization backend used to fit the model.
+It can be used to control options like the optimization algorithm, linesearch, stopping criteria, etc.
 There are currently two available backends, [`SemOptimizerOptim`](@ref) connecting to the [Optim.jl](https://github.com/JuliaNLSolvers/Optim.jl) backend, and [`SemOptimizerNLopt`](@ref) connecting to the [NLopt.jl](https://github.com/JuliaOpt/NLopt.jl) backend.
 For more information about the available options see also the tutorials about [Using Optim.jl](@ref) and [Using NLopt.jl](@ref), as well as [Constrained optimization](@ref).
 

docs/src/tutorials/construction/build_by_parts.md

@@ -40,7 +40,7 @@ end
 
 partable = ParameterTable(
     graph,
-    latent_vars = latent_vars, 
+    latent_vars = latent_vars,
     observed_vars = observed_vars)
 ```
 

docs/src/tutorials/fitting/fitting.md

@@ -7,19 +7,19 @@ model_fit = sem_fit(model)
 
 # output
 
-Fitted Structural Equation Model 
-=============================================== 
---------------------- Model ------------------- 
+Fitted Structural Equation Model
+===============================================
+--------------------- Model -------------------
 
-Structural Equation Model 
-- Loss Functions 
+Structural Equation Model
+- Loss Functions
    SemML
-- Fields 
-   observed:  SemObservedData 
-   implied:   RAM 
-   optimizer: SemOptimizerOptim 
+- Fields
+   observed:  SemObservedData
+   implied:   RAM
+   optimizer: SemOptimizerOptim
 
-------------- Optimization result ------------- 
+------------- Optimization result -------------
 
  * Status: success
 

docs/src/tutorials/meanstructure.md

@@ -40,7 +40,7 @@ end
 
 partable = ParameterTable(
     graph,
-    latent_vars = latent_vars, 
+    latent_vars = latent_vars,
     observed_vars = observed_vars)
 ```
 
@@ -78,7 +78,7 @@ end
 
 partable = ParameterTable(
     graph,
-    latent_vars = latent_vars, 
+    latent_vars = latent_vars,
     observed_vars = observed_vars)
 ```
 

docs/src/tutorials/regularization/regularization.md

@@ -2,7 +2,7 @@
 
 ## Setup
 
-For ridge regularization, you can simply use `SemRidge` as an additional loss function 
+For ridge regularization, you can simply use `SemRidge` as an additional loss function
 (for example, a model with the loss functions `SemML` and `SemRidge` corresponds to ridge-regularized maximum likelihood estimation).
 
 For lasso, elastic net and (far) beyond, we provide the `ProximalSEM` package. You can install it and load it alongside `StructuralEquationModels`:
@@ -39,7 +39,7 @@ using ProximalOperators
 ## `SemOptimizerProximal`
 
 `ProximalSEM` provides a new "building block" for the optimizer part of a model, called `SemOptimizerProximal`.
-It connects our package to the [`ProximalAlgorithms.jl`](https://github.com/JuliaFirstOrder/ProximalAlgorithms.jl) optimization backend, providing so-called proximal optimization algorithms. 
+It connects our package to the [`ProximalAlgorithms.jl`](https://github.com/JuliaFirstOrder/ProximalAlgorithms.jl) optimization backend, providing so-called proximal optimization algorithms.
 Those can handle, amongst other things, various forms of regularization.
 
 It can be used as
@@ -87,7 +87,7 @@ end
 
 partable = ParameterTable(
     graph,
-    latent_vars = latent_vars, 
+    latent_vars = latent_vars,
     observed_vars = observed_vars
 )
 

docs/src/tutorials/specification/graph_interface.md

@@ -1,6 +1,6 @@
 # Graph interface
 
-## Workflow 
+## Workflow
 As discussed before, when using the graph interface, you can specify your model as a graph
 
 ```julia
@@ -17,7 +17,7 @@ latent_vars   = ...
 
 partable = ParameterTable(
     graph,
-    latent_vars = latent_vars, 
+    latent_vars = latent_vars,
     observed_vars = observed_vars)
 
 model = Sem(
@@ -32,7 +32,7 @@ In general, there are two different types of parameters: **directed** and **indi
 We allow multiple variables on both sides of an arrow, for example `x → [y z]` or `[a b] → [c d]`. The later specifies element wise edges; that is its the same as `a → c; b → d`. If you want edges corresponding to the cross-product, we have the double lined arrow `[a b] ⇒ [c d]`, corresponding to `a → c; a → d; b → c; b → d`. The undirected arrows ↔ (element-wise) and ⇔ (crossproduct) behave the same way.
 
 !!! note "Unicode symbols in julia"
-    The `→` symbol is a unicode symbol allowed in julia (among many others; see this [list](https://docs.julialang.org/en/v1/manual/unicode-input/)). You can enter it in the julia REPL or the vscode IDE by typing `\to` followed by hitting `tab`. Similarly, 
+    The `→` symbol is a unicode symbol allowed in julia (among many others; see this [list](https://docs.julialang.org/en/v1/manual/unicode-input/)). You can enter it in the julia REPL or the vscode IDE by typing `\to` followed by hitting `tab`. Similarly,
     - `←` = `\leftarrow`,
     - `↔` = `\leftrightarrow`,
     - `⇒` = `\Rightarrow`,
@@ -54,7 +54,7 @@ graph = @StenoGraph begin
     ξ₃ ↔ fixed(1.0)*ξ₃
 end
 ```
-would 
+would
 - fix the directed effects from `ξ₁` to `x1` and from `ξ₂` to `x2` to `1`
 - leave the directed effect from `ξ₃` to `x7` free but instead restrict the variance of `ξ₃` to `1`
 - give the effect from `ξ₁` to `x3` the label `:a` (which can be convenient later if you want to retrieve information from your model about that specific parameter)
@@ -66,7 +66,7 @@ As you saw above and in the [A first model](@ref) example, the graph object need
 ```julia
 partable = ParameterTable(
     graph,
-    latent_vars = latent_vars, 
+    latent_vars = latent_vars,
     observed_vars = observed_vars)
 ```
 
@@ -85,7 +85,7 @@ The variable names (`:x1`) have to be symbols, the syntax `:something` creates a
     _(latent_vars) ⇔ _(latent_vars)
 end
 ```
-creates undirected effects coresponding to 
+creates undirected effects coresponding to
 1. the variances of all observed variables and
 2. the variances plus covariances of all latent variables
 So if you want to work with a subset of variables, simply specify a vector of symbols `somevars = [...]`, and inside the graph specification, refer to them as `_(somevars)`.

docs/src/tutorials/specification/parameter_table.md

@@ -5,5 +5,5 @@ As lavaan also uses parameter tables to store model specifications, we are worki
 
 ## Convert from and to RAMMatrices
 
-To convert a RAMMatrices object to a ParameterTable, simply use `partable = ParameterTable(rammatrices)`. 
+To convert a RAMMatrices object to a ParameterTable, simply use `partable = ParameterTable(rammatrices)`.
 To convert an object of type `ParameterTable` to RAMMatrices, you can use `ram_matrices = RAMMatrices(partable)`.

docs/src/tutorials/specification/ram_matrices.md

@@ -1,6 +1,6 @@
 # RAMMatrices interface
 
-Models can also be specified by an object of type `RAMMatrices`. 
+Models can also be specified by an object of type `RAMMatrices`.
 The RAM (reticular action model) specification corresponds to three matrices; the `A` matrix containing all directed parameters, the `S` matrix containing all undirected parameters, and the `F` matrix filtering out latent variables from the model implied covariance.
 
 The model implied covariance matrix for the observed variables of a SEM is then computed as
@@ -56,9 +56,9 @@ A =[0  0  0  0  0  0  0  0  0  0  0     1.0   0     0
 θ = Symbol.(:θ, 1:31)
 
 spec = RAMMatrices(;
-    A = A, 
-    S = S, 
-    F = F, 
+    A = A,
+    S = S,
+    F = F,
     params = θ,
     colnames = [:x1, :x2, :x3, :y1, :y2, :y3, :y4, :y5, :y6, :y7, :y8, :ind60, :dem60, :dem65]
 )
@@ -71,9 +71,9 @@ model = Sem(
 
 Let's look at this step by step:
 
-First, we specify the `A`, `S` and `F`-Matrices. 
-For a free parameter, we write a `Symbol` like `:θ1` (or any other symbol we like) to the corresponding place in the respective matrix, to constrain parameters to be equal we just use the same `Symbol` in the respective entries. 
-To fix a parameter (as in the `A`-Matrix above), we just write down the number we want to fix it to. 
+First, we specify the `A`, `S` and `F`-Matrices.
+For a free parameter, we write a `Symbol` like `:θ1` (or any other symbol we like) to the corresponding place in the respective matrix, to constrain parameters to be equal we just use the same `Symbol` in the respective entries.
+To fix a parameter (as in the `A`-Matrix above), we just write down the number we want to fix it to.
 All other entries are 0.
 
 Second, we specify a vector of symbols containing our parameters:
@@ -82,14 +82,14 @@ Second, we specify a vector of symbols containing our parameters:
 θ = Symbol.(:θ, 1:31)
 ```
 
-Third, we construct an object of type `RAMMatrices`, passing our matrices and parameters, as well as the column names of our matrices. 
+Third, we construct an object of type `RAMMatrices`, passing our matrices and parameters, as well as the column names of our matrices.
 Those are quite important, as they will be used to rearrange your data to match it to your `RAMMatrices` specification.
 
 ```julia
 spec = RAMMatrices(;
-    A = A, 
-    S = S, 
-    F = F, 
+    A = A,
+    S = S,
+    F = F,
     params = θ,
     colnames = [:x1, :x2, :x3, :y1, :y2, :y3, :y4, :y5, :y6, :y7, :y8, :ind60, :dem60, :dem65]
 )
@@ -109,7 +109,7 @@ According to the RAM, model implied mean values of the observed variables are co
 ```math
 \mu = F(I-A)^{-1}M
 ```
-where `M` is a vector of mean parameters. To estimate the means of the observed variables in our example (and set the latent means to `0`), we would specify the model just as before but add 
+where `M` is a vector of mean parameters. To estimate the means of the observed variables in our example (and set the latent means to `0`), we would specify the model just as before but add
 
 ```julia
 ...
@@ -128,5 +128,5 @@ spec = RAMMatrices(;
 
 ## Convert from and to ParameterTables
 
-To convert a RAMMatrices object to a ParameterTable, simply use `partable = ParameterTable(ram_matrices)`. 
+To convert a RAMMatrices object to a ParameterTable, simply use `partable = ParameterTable(ram_matrices)`.
 To convert an object of type `ParameterTable` to RAMMatrices, you can use `ram_matrices = RAMMatrices(partable)`.

ext/SEMBlackBoxOptimExt/AdamMutation.jl

@@ -0,0 +1,49 @@
+# mutate by moving in the gradient direction
+mutable struct AdamMutation{M <: AbstractSem, O, S} <: MutationOperator
+    model::M
+    optim::O
+    opt_state::S
+    params_fraction::Float64
+
+    function AdamMutation(model::AbstractSem, params::AbstractDict)
+        optim = RAdam(params[:AdamMutation_eta], params[:AdamMutation_beta])
+        params_fraction = params[:AdamMutation_params_fraction]
+        opt_state = Optimisers.init(optim, Vector{Float64}(undef, nparams(model)))
+
+        new{typeof(model), typeof(optim), typeof(opt_state)}(
+            model,
+            optim,
+            opt_state,
+            params_fraction,
+        )
+    end
+end
+
+Base.show(io::IO, op::AdamMutation) =
+    print(io, "AdamMutation(", op.optim, " state[3]=", op.opt_state[3], ")")
+
+"""
+Default parameters for `AdamMutation`.
+"""
+const AdamMutation_DefaultOptions = ParamsDict(
+    :AdamMutation_eta => 1E-1,
+    :AdamMutation_beta => (0.99, 0.999),
+    :AdamMutation_params_fraction => 0.25,
+)
+
+function BlackBoxOptim.apply!(m::AdamMutation, v::AbstractVector{<:Real}, target_index::Int)
+    grad = similar(v)
+    obj = SEM.evaluate!(0.0, grad, nothing, m.model, v)
+    @inbounds for i in eachindex(grad)
+        (rand() > m.params_fraction) && (grad[i] = 0.0)
+    end
+
+    m.opt_state, dv = Optimisers.apply!(m.optim, m.opt_state, v, grad)
+    if (m.opt_state[3][1] <= 1E-20) || !isfinite(obj) || any(!isfinite, dv)
+        m.opt_state = Optimisers.init(m.optim, v)
+    else
+        v .-= dv
+    end
+
+    return v
+end

ext/SEMBlackBoxOptimExt/BlackBoxOptim.jl

@@ -0,0 +1,89 @@
+############################################################################################
+### connect to BlackBoxOptim.jl as backend
+############################################################################################
+
+"""
+"""
+struct SemOptimizerBlackBoxOptim <: SemOptimizer{:BlackBoxOptim}
+    lower_bound::Float64 # default lower bound
+    variance_lower_bound::Float64 # default variance lower bound
+    lower_bounds::Union{Dict{Symbol, Float64}, Nothing}
+
+    upper_bound::Float64 # default upper bound
+    upper_bounds::Union{Dict{Symbol, Float64}, Nothing}
+end
+
+function SemOptimizerBlackBoxOptim(;
+    lower_bound::Float64 = -1000.0,
+    lower_bounds::Union{AbstractDict{Symbol, Float64}, Nothing} = nothing,
+    variance_lower_bound::Float64 = 0.001,
+    upper_bound::Float64 = 1000.0,
+    upper_bounds::Union{AbstractDict{Symbol, Float64}, Nothing} = nothing,
+    kwargs...,
+)
+    if variance_lower_bound < 0.0
+        throw(ArgumentError("variance_lower_bound must be non-negative"))
+    end
+    return SemOptimizerBlackBoxOptim(
+        lower_bound,
+        variance_lower_bound,
+        lower_bounds,
+        upper_bound,
+        upper_bounds,
+    )
+end
+
+SEM.SemOptimizer{:BlackBoxOptim}(args...; kwargs...) =
+    SemOptimizerBlackBoxOptim(args...; kwargs...)
+
+SEM.algorithm(optimizer::SemOptimizerBlackBoxOptim) = optimizer.algorithm
+SEM.options(optimizer::SemOptimizerBlackBoxOptim) = optimizer.options
+
+struct SemModelBlackBoxOptimProblem{M <: AbstractSem} <:
+       OptimizationProblem{ScalarFitnessScheme{true}}
+    model::M
+    fitness_scheme::ScalarFitnessScheme{true}
+    search_space::ContinuousRectSearchSpace
+end
+
+function BlackBoxOptim.search_space(model::AbstractSem)
+    optim = model.optimizer::SemOptimizerBlackBoxOptim
+    varparams = Set(SEM.variance_params(model.implied.ram_matrices))
+    return ContinuousRectSearchSpace(
+        [
+            begin
+                def = in(p, varparams) ? optim.variance_lower_bound : optim.lower_bound
+                isnothing(optim.lower_bounds) ? def : get(optim.lower_bounds, p, def)
+            end for p in SEM.params(model)
+        ],
+        [
+            begin
+                def = optim.upper_bound
+                isnothing(optim.upper_bounds) ? def : get(optim.upper_bounds, p, def)
+            end for p in SEM.params(model)
+        ],
+    )
+end
+
+function SemModelBlackBoxOptimProblem(
+    model::AbstractSem,
+    optimizer::SemOptimizerBlackBoxOptim,
+)
+    SemModelBlackBoxOptimProblem(model, ScalarFitnessScheme{true}(), search_space(model))
+end
+
+BlackBoxOptim.fitness(params::AbstractVector, wrapper::SemModelBlackBoxOptimProblem) =
+    return SEM.evaluate!(0.0, nothing, nothing, wrapper.model, params)
+
+# sem_fit method
+function SEM.sem_fit(
+    optimizer::SemOptimizerBlackBoxOptim,
+    model::AbstractSem,
+    start_params::AbstractVector;
+    MaxSteps::Integer = 50000,
+    kwargs...,
+)
+    problem = SemModelBlackBoxOptimProblem(model, optimizer)
+    res = bboptimize(problem; MaxSteps, kwargs...)
+    return SemFit(best_fitness(res), best_candidate(res), nothing, model, res)
+end