Multiple Refactorings

Project.toml

@@ -12,8 +12,8 @@ LazyArtifacts = "4af54fe1-eca0-43a8-85a7-787d91b784e3"
 LineSearches = "d3d80556-e9d4-5f37-9878-2ab0fcc64255"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 NLSolversBase = "d41bc354-129a-5804-8e4c-c37616107c6c"
-NLopt = "76087f3c-5699-56af-9a33-bf431cd00edd"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
+PackageExtensionCompat = "65ce6f38-6b18-4e1d-a461-8949797d7930"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
@@ -42,3 +42,16 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
 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"
+ProximalCore = "dc4f5ac2-75d1-4f31-931e-60435d74994b"
+ProximalOperators = "f3b72e0c-5f3e-4b3e-8f3e-3f4f3e3e3e3e"
+
+[extensions]
+SEMNLOptExt = "NLopt"
+SEMProximalOptExt = ["ProximalCore", "ProximalAlgorithms", "ProximalOperators"]
+SEMBlackBoxOptimExt = ["BlackBoxOptim", "Optimisers"]

ext/SEMBlackBoxOptimExt/AdamMutation.jl

@@ -0,0 +1,44 @@
+# 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,75 @@
+############################################################################################
+### 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.imply.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,138 @@
+"""
+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

ext/SEMBlackBoxOptimExt/SemOptimizerBlackBoxOptim.jl

@@ -0,0 +1,76 @@
+############################################################################################
+### 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
+    return ContinuousRectSearchSpace(
+        SEM.lower_bounds(optim.lower_bounds, model, default=optim.lower_bound, variance_default=optim.variance_lower_bound),
+        SEM.upper_bounds(optim.upper_bounds, model, default=optim.upper_bound))
+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;
+    Method::Symbol = :adaptive_de_rand_1_bin_with_gradient,
+    MaxSteps::Integer = 50000,
+    kwargs...)
+
+    problem = SemModelBlackBoxOptimProblem(model, optimizer)
+    if Method == :adaptive_de_rand_1_bin_with_gradient
+        # custom adaptive differential evolution with mutation that moves along the gradient
+        bbopt_factory = DefaultDiffEvoFactory(problem; kwargs...)
+        bbopt = bbsetup(bbopt_factory; MaxSteps, kwargs...)
+    else
+        bbopt = bbsetup(problem; Method, MaxSteps, kwargs...)
+    end
+    res = bboptimize(bbopt)
+    return SemFit(best_fitness(res), best_candidate(res),
+                  nothing, model, res)
+end

ext/SEMNLOptExt.jl

@@ -0,0 +1,12 @@
+module SEMNLOptExt
+
+using StructuralEquationModels, NLopt
+
+SEM = StructuralEquationModels
+
+export SemOptimizerNLopt, NLoptConstraint
+
+include("diff/NLopt.jl")
+include("optimizer/NLopt.jl")
+
+end