Multiple refactorings reformatted

Project.toml

@@ -14,9 +14,11 @@ 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"
 ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
+Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
@@ -51,9 +53,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/performance/mixed_differentiation.md

@@ -10,12 +10,9 @@ ml_term = SemML(
     RAMSymbolic(partable)
 )
 
-ridge_term = SemRidge(
-    α_ridge = 0.01,
-    which_ridge = params(ml_term)
-)
+ridge_penalty = SemRidge(partable)
 
-model_ml_ridge = Sem(ml_term, ridge_term)
+model_ml_ridge = Sem(ml_term, ridge_penalty => 0.01)
 
 model_ml_ridge_fit = fit(model_ml_ridge)
 ```

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

ext/SEMBlackBoxOptimExt/DiffEvoFactory.jl

@@ -0,0 +1,196 @@
+"""
+Base class for factories of optimizers for a specific problem.
+"""
+abstract type OptimizerFactory{P <: OptimizationProblem} end
+
+problem(factory::OptimizerFactory) = factory.problem
+
+const OptController_DefaultParameters = ParamsDict(
+    :MaxTime => 60.0,
+    :MaxSteps => 10^8,
+    :TraceMode => :compact,
+    :TraceInterval => 5.0,
+    :RecoverResults => false,
+    :SaveTrace => false,
+)
+
+function generate_opt_controller(alg::Optimizer, optim_factory::OptimizerFactory, params)
+    return BlackBoxOptim.OptController(
+        alg,
+        problem(optim_factory),
+        BlackBoxOptim.chain(
+            BlackBoxOptim.DefaultParameters,
+            OptController_DefaultParameters,
+            params,
+        ),
+    )
+end
+
+function check_population(
+    factory::OptimizerFactory,
+    popmatrix::BlackBoxOptim.PopulationMatrix,
+)
+    ssp = factory |> problem |> search_space
+    for i in 1:popsize(popmatrix)
+        @assert popmatrix[:, i] ∈ ssp "Individual $i is out of space: $(popmatrix[:,i])" # fitness: $(fitness(population, i))"
+    end
+end
+
+initial_search_space(factory::OptimizerFactory, id::Int) = search_space(factory.problem)
+
+function initial_population_matrix(factory::OptimizerFactory, id::Int)
+    #@info "Standard initial_population_matrix()"
+    ini_ss = initial_search_space(factory, id)
+    if !isempty(factory.initial_population)
+        numdims(factory.initial_population) == numdims(factory.problem) || throw(
+            DimensionMismatch(
+                "Dimensions of :Population ($(numdims(factory.initial_population))) " *
+                "are different from the problem dimensions ($(numdims(factory.problem)))",
+            ),
+        )
+        res = factory.initial_population[
+            :,
+            StatsBase.sample(
+                1:popsize(factory.initial_population),
+                factory.population_size,
+            ),
+        ]
+    else
+        res = rand_individuals(ini_ss, factory.population_size, method = :latin_hypercube)
+    end
+    prj = RandomBound(ini_ss)
+    if size(res, 2) > 1
+        apply!(prj, view(res, :, 1), SEM.start_fabin3(factory.problem.model))
+    end
+    if size(res, 2) > 2
+        apply!(prj, view(res, :, 2), SEM.start_simple(factory.problem.model))
+    end
+    return res
+end
+
+# convert individuals in the archive into population matrix
+population_matrix(archive::Any) = population_matrix!(
+    Matrix{Float64}(undef, length(BlackBoxOptim.params(first(archive))), length(archive)),
+    archive,
+)
+
+function population_matrix!(pop::AbstractMatrix{<:Real}, archive::Any)
+    npars = length(BlackBoxOptim.params(first(archive)))
+    size(pop, 1) == npars || throw(
+        DimensionMismatch(
+            "Matrix rows count ($(size(pop, 1))) doesn't match the number of problem dimensions ($(npars))",
+        ),
+    )
+    @inbounds for (i, indi) in enumerate(archive)
+        (i <= size(pop, 2)) || break
+        pop[:, i] .= BlackBoxOptim.params(indi)
+    end
+    if size(pop, 2) > length(archive)
+        @warn "Matrix columns count ($(size(pop, 2))) is bigger than population size ($(length(archive))), last columns not set"
+    end
+    return pop
+end
+
+generate_embedder(factory::OptimizerFactory, id::Int, problem::OptimizationProblem) =
+    RandomBound(search_space(problem))
+
+abstract type DiffEvoFactory{P <: OptimizationProblem} <: OptimizerFactory{P} end
+
+generate_selector(
+    factory::DiffEvoFactory,
+    id::Int,
+    problem::OptimizationProblem,
+    population,
+) = RadiusLimitedSelector(get(factory.params, :selector_radius, popsize(population) ÷ 5))
+
+function generate_modifier(factory::DiffEvoFactory, id::Int, problem::OptimizationProblem)
+    ops = GeneticOperator[
+        MutationClock(UniformMutation(search_space(problem)), 1 / numdims(problem)),
+        BlackBoxOptim.AdaptiveDiffEvoRandBin1(
+            BlackBoxOptim.AdaptiveDiffEvoParameters(
+                factory.params[:fdistr],
+                factory.params[:crdistr],
+            ),
+        ),
+        SimplexCrossover{3}(1.05),
+        SimplexCrossover{2}(1.1),
+        #SimulatedBinaryCrossover(0.05, 16.0),
+        #SimulatedBinaryCrossover(0.05, 3.0),
+        #SimulatedBinaryCrossover(0.1, 5.0),
+        #SimulatedBinaryCrossover(0.2, 16.0),
+        UnimodalNormalDistributionCrossover{2}(
+            chain(BlackBoxOptim.UNDX_DefaultOptions, factory.params),
+        ),
+        UnimodalNormalDistributionCrossover{3}(
+            chain(BlackBoxOptim.UNDX_DefaultOptions, factory.params),
+        ),
+        ParentCentricCrossover{2}(chain(BlackBoxOptim.PCX_DefaultOptions, factory.params)),
+        ParentCentricCrossover{3}(chain(BlackBoxOptim.PCX_DefaultOptions, factory.params)),
+    ]
+    if problem isa SemModelBlackBoxOptimProblem
+        push!(
+            ops,
+            AdamMutation(problem.model, chain(AdamMutation_DefaultOptions, factory.params)),
+        )
+    end
+    FAGeneticOperatorsMixture(ops)
+end
+
+function generate_optimizer(
+    factory::DiffEvoFactory,
+    id::Int,
+    problem::OptimizationProblem,
+    popmatrix,
+)
+    population = FitPopulation(popmatrix, nafitness(fitness_scheme(problem)))
+    BlackBoxOptim.DiffEvoOpt(
+        "AdaptiveDE/rand/1/bin/gradient",
+        population,
+        generate_selector(factory, id, problem, population),
+        generate_modifier(factory, id, problem),
+        generate_embedder(factory, id, problem),
+    )
+end
+
+const Population_DefaultParameters = ParamsDict(
+    :Population => BlackBoxOptim.PopulationMatrix(undef, 0, 0),
+    :PopulationSize => 100,
+)
+
+const DE_DefaultParameters = chain(
+    ParamsDict(
+        :SelectorRadius => 0,
+        :fdistr =>
+            BlackBoxOptim.BimodalCauchy(0.65, 0.1, 1.0, 0.1, clampBelow0 = false),
+        :crdistr =>
+            BlackBoxOptim.BimodalCauchy(0.1, 0.1, 0.95, 0.1, clampBelow0 = false),
+    ),
+    Population_DefaultParameters,
+)
+
+struct DefaultDiffEvoFactory{P <: OptimizationProblem} <: DiffEvoFactory{P}
+    problem::P
+    initial_population::BlackBoxOptim.PopulationMatrix
+    population_size::Int
+    params::ParamsDictChain
+end
+
+DefaultDiffEvoFactory(problem::OptimizationProblem; kwargs...) =
+    DefaultDiffEvoFactory(problem, BlackBoxOptim.kwargs2dict(kwargs))
+
+function DefaultDiffEvoFactory(problem::OptimizationProblem, params::AbstractDict)
+    params = chain(DE_DefaultParameters, params)
+    DefaultDiffEvoFactory{typeof(problem)}(
+        problem,
+        params[:Population],
+        params[:PopulationSize],
+        params,
+    )
+end
+
+function BlackBoxOptim.bbsetup(factory::OptimizerFactory; kwargs...)
+    popmatrix = initial_population_matrix(factory, 1)
+    check_population(factory, popmatrix)
+    alg = generate_optimizer(factory, 1, problem(factory), popmatrix)
+    return generate_opt_controller(alg, factory, BlackBoxOptim.kwargs2dict(kwargs))
+end

ext/SEMBlackBoxOptimExt/SEMBlackBoxOptimExt.jl

@@ -0,0 +1,13 @@
+module SEMBlackBoxOptimExt
+
+using StructuralEquationModels, BlackBoxOptim, Optimisers
+
+SEM = StructuralEquationModels
+
+export SemOptimizerBlackBoxOptim
+
+include("AdamMutation.jl")
+include("DiffEvoFactory.jl")
+include("SemOptimizerBlackBoxOptim.jl")
+
+end