diff --git a/src/dual_model.jl b/src/dual_model.jl index 4794bfc..6d98c6a 100755 --- a/src/dual_model.jl +++ b/src/dual_model.jl @@ -9,7 +9,8 @@ end The dual objective function of the MMA model, `model`. The approximate objective and constraint values are cached in `approxfg` when the dual objective is called. For every input `λ`, the optimal primal solution is computed and cached in the field `x`. """ -struct MMADualObj{TApprox <: AbstractMMAApprox, Tx <: AbstractVector, Ta <: AbstractVector} <: AbstractFunction +struct MMADualObj{TApprox<:AbstractMMAApprox,Tx<:AbstractVector,Ta<:AbstractVector} <: + AbstractFunction model::MMAApproxModel{TApprox} x::Tx # primal solution approxfg::Ta # approximate [f(x); g(x)] @@ -86,7 +87,7 @@ function optimizeprimal!(f::MMADualObj{TA}, λ::AbstractVector) where {TA} if TA <: XMMAApprox @unpack a, c, d = model.approx_objective_ineq_constraints - map!(view(x, size(p, 2) + 1 : size(p, 2) + length(c)), 1:length(c)) do i + map!(view(x, size(p, 2)+1:size(p, 2)+length(c)), 1:length(c)) do i temp = c[i] - λ[i] if temp >= 0 return zero(T) @@ -155,31 +156,23 @@ getoptimalx(m::MMADualModel) = getoptimalx(m.obj) getparent(m::MMADualModel) = m.parent -@doc docσ -getσ(m::MMADualModel) = getσ(getparent(m)) +@doc docσ getσ(m::MMADualModel) = getσ(getparent(m)) -@doc docσ -setσ!(m::MMADualModel, σ) = setσ!(getparent(m), σ) +@doc docσ setσ!(m::MMADualModel, σ) = setσ!(getparent(m), σ) -@doc docρ -getρ(m::MMADualModel) = getρ(getparent(m)) +@doc docρ getρ(m::MMADualModel) = getρ(getparent(m)) -@doc docρ -setρ!(m::MMADualModel, ρ) = setρ!(getparent(m), ρ) +@doc docρ setρ!(m::MMADualModel, ρ) = setρ!(getparent(m), ρ) -@doc docupdateapprox! -function updateapprox!(m::MMADualModel, args...) +@doc docupdateapprox! function updateapprox!(m::MMADualModel, args...) updateapprox!(getparent(m), args...) return m end getobjective(m::MMADualModel) = m.obj -@doc docxk -getxk(m::MMADualModel) = getxk(getparent(m)) +@doc docxk getxk(m::MMADualModel) = getxk(getparent(m)) -@doc docfk -getfk(m::MMADualModel) = getfk(getparent(m)) +@doc docfk getfk(m::MMADualModel) = getfk(getparent(m)) -@doc doc∇fk -get∇fk(m::MMADualModel) = get∇fk(getparent(m)) +@doc doc∇fk get∇fk(m::MMADualModel) = get∇fk(getparent(m)) diff --git a/src/mma_algorithm.jl b/src/mma_algorithm.jl index 88dd955..3aaff66 100755 --- a/src/mma_algorithm.jl +++ b/src/mma_algorithm.jl @@ -6,7 +6,7 @@ The original method of moving asymptotes (MMA) algorithm from the [1987 paper](https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.1620240207). """ @params struct MMA87 <: AbstractOptimizer - dualoptimizer + dualoptimizer::Any end function MMA87(; dualoptimizer = Optim.GradientDescent()) return MMA87(dualoptimizer) @@ -20,7 +20,7 @@ const MMA = MMA87 The globally convergent method of moving asymptotes (MMA) algorithm from the [2002 paper](https://epubs.siam.org/doi/abs/10.1137/S1052623499362822). """ @params struct MMA02 <: AbstractOptimizer - dualoptimizer + dualoptimizer::Any end function MMA02(; dualoptimizer = Optim.GradientDescent()) return MMA02(dualoptimizer) @@ -44,7 +44,12 @@ A struct that stores all the options of the MMA algorithms. Th following are the - `convcriteria`: an instance of [`ConvergenceCriteria`](@ref) that specifies the convergence criteria of the MMA algorithm. - `verbose`: true/false, when true prints convergence statistics. """ -@with_kw mutable struct MMAOptions{T, Ttol <: Tolerance, TSubOptions <: Optim.Options, TC <: ConvergenceCriteria} +@with_kw mutable struct MMAOptions{ + T, + Ttol<:Tolerance, + TSubOptions<:Optim.Options, + TC<:ConvergenceCriteria, +} maxiter::Int = 1000 outer_maxiter::Int = 10^8 maxinner::Int = 10 @@ -60,7 +65,7 @@ A struct that stores all the options of the MMA algorithms. Th following are the allow_outer_f_increases = true, allow_f_increases = true, iterations = 1000, - outer_iterations=1000, + outer_iterations = 1000, ) convcriteria::TC = KKTCriteria() verbose::Bool = true @@ -105,11 +110,11 @@ A struct that stores all the intermediate states and memory allocations needed f tempx::AbstractVector callback::Function optimizer::AbstractOptimizer - options + options::Any trace::Trace outer_iter::Int iter::Int - fcalls::Int + fcalls::Int end function MMAWorkspace( model::VecModel, @@ -119,7 +124,11 @@ function MMAWorkspace( plot_trace::Bool = false, show_plot::Bool = plot_trace, save_plot = nothing, - callback::Function = plot_trace ? LazyPlottingCallback(; show_plot = show_plot, save_plot = save_plot) : NoCallback(), + callback::Function = plot_trace ? + LazyPlottingCallback(; + show_plot = show_plot, + save_plot = save_plot, + ) : NoCallback(), kwargs..., ) where {T} init!(model) @@ -128,7 +137,14 @@ function MMAWorkspace( # Convergence λ = ones(getdim(getineqconstraints(model))) solution = Solution(dualmodel, λ) - assess_convergence!(solution, model, options.tol, options.convcriteria, options.verbose, 0) + assess_convergence!( + solution, + model, + options.tol, + options.convcriteria, + options.verbose, + 0, + ) correctsolution!(solution, model, options) # Trace @@ -178,7 +194,7 @@ default_options(model::VecModel, alg::MMA02) = MMAOptions() init!(model::VecModel) = model -function Workspace(model::VecModel, optimizer::Union{MMA87, MMA02}, args...; kwargs...) +function Workspace(model::VecModel, optimizer::Union{MMA87,MMA02}, args...; kwargs...) return MMAWorkspace(model, optimizer, args...; kwargs...) end @@ -221,7 +237,9 @@ function optimize!(workspace::MMAWorkspace) # Initialize the workspace's lift ρ, only used in MMA02 initializeρ!(workspace) - while (!hasconverged(solution) || outer_iter == 0) && iter < maxiter && outer_iter < outer_maxiter + while (!hasconverged(solution) || outer_iter == 0) && + iter < maxiter && + outer_iter < outer_maxiter outer_iter += 1 # Adapt workspace's trust region σ using s_incr and s_decr @@ -240,19 +258,23 @@ function optimize!(workspace::MMAWorkspace) # Assume the convex approximation is not an upper bound approximation at x upperbounded = false prev_iter = iter - while (!hasconverged(solution) || outer_iter == 1) && !upperbounded && iter < maxiter && iter < prev_iter + options.maxinner + while (!hasconverged(solution) || outer_iter == 1) && + !upperbounded && + iter < maxiter && + iter < prev_iter + options.maxinner iter += 1 # Solve the dual problem by minimizing negative the dual objective value if length(λ) > 0 λ .= 1.0 - λ .= Optim.optimize( - Optim.only_fg!(optimobj), - λl, - λu, - λ, - Optim.Fminbox(dualoptimizer), - dual_options, - ).minimizer + λ .= + Optim.optimize( + Optim.only_fg!(optimobj), + λl, + λu, + λ, + Optim.Fminbox(dualoptimizer), + dual_options, + ).minimizer end if debugging[] @show λ @@ -270,10 +292,7 @@ function optimize!(workspace::MMAWorkspace) # Evaluates the exact objective and constraints and their gradients at the optimal x optimalx = getoptimalx(dualmodel) - fg, ∇fg = value_jacobian( - approxmodel.objective_ineq_constraints, - optimalx, - ) + fg, ∇fg = value_jacobian(approxmodel.objective_ineq_constraints, optimalx) fcalls += 1 # Scale the objective appropriately @@ -303,7 +322,7 @@ function optimize!(workspace::MMAWorkspace) if options.keep_best best_solution = get_best_solution(solution, best_solution, options) else - best_solution = deepcopy(solution) + best_solution = deepcopy(solution) end hasconverged(best_solution) && break @@ -316,7 +335,7 @@ function optimize!(workspace::MMAWorkspace) # Reset the objective scaling factor to 1 set_objective_multiple!(model, 1) callback(best_solution, update = true) - + results = GenericResult( optimizer, x0, @@ -337,13 +356,13 @@ end Scales the objective of `model` using the ratio of the ∞ norms of the constraint and objective values and gradients. After scaling, the ∞ norms will be the same. """ function scaleobjective!(model::VecModel, fg, ∇fg, approxfg, λ) - @views normratio = norm(∇fg[2:end,:], Inf) / norm(∇fg[1,:], Inf) + @views normratio = norm(∇fg[2:end, :], Inf) / norm(∇fg[1, :], Inf) @views normratio = max(normratio, norm(fg[2:end], Inf) / abs(fg[1])) normratio = isfinite(normratio) ? normratio : one(normratio) set_objective_multiple!(model, get_objective_multiple(model) * normratio) fg[1] = fg[1] * normratio approxfg[1] = approxfg[1] * normratio - ∇fg[1,:] .= ∇fg[1,:] .* normratio + ∇fg[1, :] .= ∇fg[1, :] .* normratio scalequadweight!(model, normratio) λ .*= normratio return model @@ -358,8 +377,8 @@ end function get_best_solution(solution::Solution, best_solution::Solution, options) best_infeas = best_solution.convstate.infeas - if max(solution.convstate.infeas, solution.convstate.kkt_residual) <= - max(best_solution.convstate.infeas, best_solution.convstate.kkt_residual) + if max(solution.convstate.infeas, solution.convstate.kkt_residual) <= + max(best_solution.convstate.infeas, best_solution.convstate.kkt_residual) best_solution = deepcopy(solution) end return best_solution @@ -418,7 +437,7 @@ function updateσ!(workspace::MMAWorkspace) if !isfinite(diff) diff = 1000 * one(T) end - min = diff/100 + min = diff / 100 max = 10diff if d <= min return min @@ -449,14 +468,14 @@ function increaseρ!( @unpack x, prevx = solution T = eltype(x) w = zero(T) - for j in 1:length(x) + for j = 1:length(x) diff2 = (x[j] - prevx[j])^2 w += diff2 / (σ[j]^2 - diff2) end w /= 2 upperbounded = true - for i in 1:length(ρ) + for i = 1:length(ρ) if fg[i] > approxfg[i] upperbounded = false ρ[i] = min(10ρ[i], 1.1 * (ρ[i] + (fg[i] - approxfg[i]) / w)) diff --git a/src/mma_approx.jl b/src/mma_approx.jl index e0165eb..23f84b4 100755 --- a/src/mma_approx.jl +++ b/src/mma_approx.jl @@ -73,7 +73,7 @@ end function MMAApprox( parent::AbstractFunction, x::AbstractVector, - f::Union{Real, AbstractVector}, + f::Union{Real,AbstractVector}, ∇f::AbstractMatrix; σ = fill(10.0, length(x)), ρ = f isa Real ? Ref(zero(f)) : zeros(length(f)), @@ -101,9 +101,9 @@ Returns an approximation of the function `f.parent` at point `x`. See [`MMAAppro """ function (f::MMAApprox{<:Base.RefValue{<:Real}})(x::AbstractVector) @unpack p, q, l, u, r, out = f - @assert all(l[j] < x[j] < u[j] for j in 1:length(x)) + @assert all(l[j] < x[j] < u[j] for j = 1:length(x)) _out = r[] - for j in 1:length(x) + for j = 1:length(x) _out += p[j] / (u[j] - x[j]) + q[j] / (x[j] - l[j]) end out[] = ForwardDiff.value(_out) @@ -111,11 +111,11 @@ function (f::MMAApprox{<:Base.RefValue{<:Real}})(x::AbstractVector) end function (f::MMAApprox{<:AbstractVector})(x::AbstractVector) @unpack p, q, l, u, r, out = f - @assert all(l[j] < x[j] < u[j] for j in 1:length(x)) + @assert all(l[j] < x[j] < u[j] for j = 1:length(x)) if eltype(x) <: ForwardDiff.Dual dual_out = map(1:getdim(f)) do i output = r[i] - for j in 1:length(x) + for j = 1:length(x) output += p[i, j] / (u[j] - x[j]) + q[i, j] / (x[j] - l[j]) end return output @@ -123,9 +123,9 @@ function (f::MMAApprox{<:AbstractVector})(x::AbstractVector) out .= ForwardDiff.value.(dual_out) return dual_out else - for i in 1:getdim(f) + for i = 1:getdim(f) out[i] = r[i] - for j in 1:length(x) + for j = 1:length(x) out[i] += p[i, j] / (u[j] - x[j]) + q[i, j] / (x[j] - l[j]) end end @@ -135,7 +135,7 @@ end function ChainRulesCore.rrule(f::MMAApprox{<:Base.RefValue{<:Real}}, x::AbstractVector) @unpack p, q, l, u, r = f - @assert all(f.l[j] < x[j] < f.u[j] for j in 1:length(x)) + @assert all(f.l[j] < x[j] < f.u[j] for j = 1:length(x)) out = f(x) adjoint = Δ -> begin if p isa AbstractMatrix @@ -149,7 +149,7 @@ function ChainRulesCore.rrule(f::MMAApprox{<:Base.RefValue{<:Real}}, x::Abstract end function ChainRulesCore.rrule(f::MMAApprox{<:AbstractVector}, x::AbstractVector) @unpack p, q, l, u, r = f - @assert all(f.l[j] < x[j] < f.u[j] for j in 1:length(x)) + @assert all(f.l[j] < x[j] < f.u[j] for j = 1:length(x)) out = f(x) adjoint = Δ -> begin (nothing, (p ./ (u' .- x') .^ 2 .- q ./ (x' .- l') .^ 2)' * Δ) @@ -217,17 +217,11 @@ end getparent(approx::MMAApprox) = approx.parent getdim(approx::MMAApprox) = getdim(getparent(approx)) -function updateapprox!( - approx::MMAApprox{<:Base.RefValue{<:Real}}, - x::AbstractVector, -) +function updateapprox!(approx::MMAApprox{<:Base.RefValue{<:Real}}, x::AbstractVector) val, grad = value_gradient(getparent(approx), x) return updateapprox!(approx, x, val, grad) end -function updateapprox!( - approx::MMAApprox{<:AbstractVector}, - x::AbstractVector, -) +function updateapprox!(approx::MMAApprox{<:AbstractVector}, x::AbstractVector) val, jac = value_jacobian(getparent(approx), x) return updateapprox!(approx, x, val, jac) end @@ -235,7 +229,7 @@ function updateapprox!( approx::MMAApprox{<:Base.RefValue{<:Real}}, x::AbstractVector, f::Real, - ∇f::Union{AbstractVector, Adjoint{<:Any, <:AbstractVector}}, + ∇f::Union{AbstractVector,Adjoint{<:Any,<:AbstractVector}}, ) setxk!(approx, x) setfk!(approx, f) @@ -258,12 +252,12 @@ function updateapprox!(approx::MMAApprox{<:Base.RefValue{<:Real}}) T = eltype(∇fk) l .= xk .- σ u .= xk .+ σ - for j in 1:length(p) + for j = 1:length(p) p[j] = σ[j] * σ[j] * max(0, ∇fk[j]) + ρ[] * σ[j] / 4 q[j] = σ[j] * σ[j] * max(0, -∇fk[j]) + ρ[] * σ[j] / 4 end r[] = fk[] - for j in 1:length(p) + for j = 1:length(p) r[] -= (p[j] + q[j]) / σ[j] end return approx @@ -274,15 +268,15 @@ function updateapprox!(approx::MMAApprox{<:AbstractVector}) l .= xk .- σ u .= xk .+ σ # All matrices are stored as Adjoint{<:Real, <:AbstractMatrix} for cache efficiency - for i in 1:size(p, 1) - for j in 1:size(p, 2) + for i = 1:size(p, 1) + for j = 1:size(p, 2) p[i, j] = σ[j] * σ[j] * max(0, ∇fk[i, j]) + ρ[i] * σ[j] / 4 q[i, j] = σ[j] * σ[j] * max(0, -∇fk[i, j]) + ρ[i] * σ[j] / 4 end end - for i in 1:size(p, 1) + for i = 1:size(p, 1) r[i] = fk[i] - for j in 1:size(p, 2) + for j = 1:size(p, 2) r[i] -= (p[i, j] + q[i, j]) / σ[j] end end diff --git a/src/mma_approx_docs.jl b/src/mma_approx_docs.jl index 0266a37..f68ca15 100755 --- a/src/mma_approx_docs.jl +++ b/src/mma_approx_docs.jl @@ -23,13 +23,11 @@ docσ = """ Get or set the move limit `σ` of the MMA approximation `approx`. See [`MMAApprox`](@ref) for an explanation. """ -@doc docσ -function getσ(approx::AbstractMMAApprox) +@doc docσ function getσ(approx::AbstractMMAApprox) throw("`getσ` is not defined.") end -@doc docσ -function setσ!(approx::AbstractMMAApprox, σ) +@doc docσ function setσ!(approx::AbstractMMAApprox, σ) throw("`setσ!` is not defined.") end @@ -40,13 +38,11 @@ docρ = """ Get or set the value of `ρ` in `approx`. See [`MMAApprox`](@ref) for an explanation. """ -@doc docρ -function getρ(approx::AbstractMMAApprox) +@doc docρ function getρ(approx::AbstractMMAApprox) throw("`getρ` is not defined.") end -@doc docρ -function setρ!(approx::AbstractMMAApprox, ρ) +@doc docρ function setρ!(approx::AbstractMMAApprox, ρ) throw("`setρ!` is not defined.") end @@ -57,13 +53,11 @@ docfk = """ Get or set the function value `fk` of the approximated function at the approximation point `xk` of `approx`. """ -@doc docfk -function getfk(approx::AbstractMMAApprox) +@doc docfk function getfk(approx::AbstractMMAApprox) throw("`getfk` is not defined.") end -@doc docfk -function setfk!(approx::AbstractMMAApprox, f) +@doc docfk function setfk!(approx::AbstractMMAApprox, f) throw("`setfk!` is not defined.") end @@ -74,13 +68,11 @@ docxk = """ Get or set the approximation point `xk` of `approx`. """ -@doc docxk -function getxk(approx::AbstractMMAApprox) +@doc docxk function getxk(approx::AbstractMMAApprox) throw("`getfk` is not defined.") end -@doc docxk -function setxk!(approx::AbstractMMAApprox, x::AbstractVector) +@doc docxk function setxk!(approx::AbstractMMAApprox, x::AbstractVector) throw("`setxk!` is not defined.") end @@ -91,13 +83,11 @@ doc∇fk = """ Get or set the value of the gradient/Jacobian, `∇fk`, of the approximated function at the approximation point `xk`. """ -@doc doc∇fk -function get∇fk(approx::AbstractMMAApprox) +@doc doc∇fk function get∇fk(approx::AbstractMMAApprox) throw("`get∇fk` is not defined.") end -@doc doc∇fk -function set∇fk!(approx::AbstractMMAApprox, ∇f::AbstractVecOrMat) +@doc doc∇fk function set∇fk!(approx::AbstractMMAApprox, ∇f::AbstractVecOrMat) throw("`set∇fk!` is not defined.") end @@ -107,8 +97,7 @@ docupdateapprox! = """ Updates the MMA approximation in `f` making it around the point `x`. """ -@doc docupdateapprox! -function updateapprox!(approx::AbstractMMAApprox, x::AbstractVector) +@doc docupdateapprox! function updateapprox!(approx::AbstractMMAApprox, x::AbstractVector) throw("`updateapprox` is not defined.") end diff --git a/src/mma_model.jl b/src/mma_model.jl index d96248f..f372a3f 100755 --- a/src/mma_model.jl +++ b/src/mma_model.jl @@ -16,7 +16,13 @@ An approximate restricted model that uses the `MMAApprox` or `XMMAApprox` approx - `box_min`: the restricted lower bounds on the decision variables. - `box_max`: the restricted upper bounds on the decision variables. """ -struct MMAApproxModel{TApprox <: AbstractMMAApprox, Tp <: VecModel, To <: AbstractFunction, Tmi <: AbstractVector, Tma <: AbstractVector} <: AbstractModel +struct MMAApproxModel{ + TApprox<:AbstractMMAApprox, + Tp<:VecModel, + To<:AbstractFunction, + Tmi<:AbstractVector, + Tma<:AbstractVector, +} <: AbstractModel parent::Tp objective_ineq_constraints::To approx_objective_ineq_constraints::TApprox @@ -29,12 +35,7 @@ end Constructs an MMA approximation of the model `parent` around the point `x`. If `extended` is true, an [`XMMAApprox`](@ref) approximation is used. Otherwise, the default [`MMAApprox`](@ref) approximation is used. `kwargs` can be used to pass additional options to `XMMAApprox`, e.g. setting the coefficients. """ -function MMAApproxModel( - parent::VecModel, - x::AbstractVector; - extended = false, - kwargs..., -) +function MMAApproxModel(parent::VecModel, x::AbstractVector; extended = false, kwargs...) T = eltype(x) σ = map(1:length(x)) do j diff = getmax(parent, j) - getmin(parent, j) @@ -45,12 +46,7 @@ function MMAApproxModel( end ρ = zeros(length(getineqconstraints(parent)) + 1) obj_constr = getobjectiveconstraints(parent) - approx_obj_constr = MMAApprox( - obj_constr, - x; - σ = σ, - ρ = ρ, - ) + approx_obj_constr = MMAApprox(obj_constr, x; σ = σ, ρ = ρ) model = MMAApproxModel( parent, obj_constr, @@ -62,50 +58,40 @@ function MMAApproxModel( return model end -getmin(m::MMAApproxModel)= m.box_min +getmin(m::MMAApproxModel) = m.box_min getmax(m::MMAApproxModel) = m.box_max -@doc docρ -getρ(model::MMAApproxModel) = getρ(model.approx_objective_ineq_constraints) +@doc docρ getρ(model::MMAApproxModel) = getρ(model.approx_objective_ineq_constraints) -@doc docρ -function setρ!(model::MMAApproxModel, ρ) +@doc docρ function setρ!(model::MMAApproxModel, ρ) setρ!(model.approx_objective_ineq_constraints, ρ) return model end -@doc docσ -getσ(model::MMAApproxModel) = getσ(model.approx_objective_ineq_constraints) +@doc docσ getσ(model::MMAApproxModel) = getσ(model.approx_objective_ineq_constraints) -@doc docσ -function setσ!(model::MMAApproxModel, σ) +@doc docσ function setσ!(model::MMAApproxModel, σ) setσ!(model.approx_objective_ineq_constraints, σ) return model end -@doc docxk -getxk(model::MMAApproxModel) = getxk(model.approx_objective_ineq_constraints) +@doc docxk getxk(model::MMAApproxModel) = getxk(model.approx_objective_ineq_constraints) -@doc docxk -function setxk!(model::MMAApproxModel, x::AbstractVector) +@doc docxk function setxk!(model::MMAApproxModel, x::AbstractVector) setxk!(model.approx_objective_ineq_constraints, x) return model end -@doc docfk -getfk(model::MMAApproxModel) = getfk(model.approx_objective_ineq_constraints) +@doc docfk getfk(model::MMAApproxModel) = getfk(model.approx_objective_ineq_constraints) -@doc docfk -function setfk!(model::MMAApproxModel, f) +@doc docfk function setfk!(model::MMAApproxModel, f) setfk!(model.approx_objective_ineq_constraints, f) return model end -@doc doc∇fk -get∇fk(model::MMAApproxModel) = get∇fk(model.approx_objective_ineq_constraints) +@doc doc∇fk get∇fk(model::MMAApproxModel) = get∇fk(model.approx_objective_ineq_constraints) -@doc doc∇fk -function set∇fk!(model::MMAApproxModel, ∇f::AbstractVecOrMat) +@doc doc∇fk function set∇fk!(model::MMAApproxModel, ∇f::AbstractVecOrMat) set∇fk!(model.approx_objective_ineq_constraints, ∇f) return model end diff --git a/src/xmma_approx.jl b/src/xmma_approx.jl index 1382e4f..82e6dd6 100755 --- a/src/xmma_approx.jl +++ b/src/xmma_approx.jl @@ -74,11 +74,13 @@ function (approx::XMMAApprox)(xyz::AbstractVector) @unpack mma, a, c, d = approx T = eltype(xyz) z = xyz[end] - y = xyz[end - length(c) : end - 1] - x = xyz[begin : end - length(c) - 1] + y = xyz[end-length(c):end-1] + x = xyz[begin:end-length(c)-1] approxfg = mma(x) ε = 1e-4 - out = approxfg .+ [a[1] * z + ε * z^2 + sum(c .* y .+ d .* y.^2 ./ 2); .-a[2:end] .* z .- y] + out = + approxfg .+ + [a[1] * z + ε * z^2 + sum(c .* y .+ d .* y .^ 2 ./ 2); .-a[2:end] .* z .- y] saveout!(approx, out) return out end @@ -133,10 +135,7 @@ end getparent(approx::XMMAApprox) = getparent(approx.mma) getdim(approx::XMMAApprox) = getdim(approx.mma) -function updateapprox!( - approx::XMMAApprox, - x::AbstractVector, -) +function updateapprox!(approx::XMMAApprox, x::AbstractVector) updateapprox!(approx.mma, x) return approx end diff --git a/test/approximation.jl b/test/approximation.jl index 8f139fe..d803546 100644 --- a/test/approximation.jl +++ b/test/approximation.jl @@ -2,7 +2,7 @@ using NonconvexMMA, LinearAlgebra, Test, Zygote, FiniteDifferences const FDM = FiniteDifferences f(x::AbstractVector) = sqrt(x[2]) -g(x::AbstractVector, a, b) = (a*x[1] + b)^3 - x[2] +g(x::AbstractVector, a, b) = (a * x[1] + b)^3 - x[2] _m = Model(f) addvar!(_m, [0.0, 0.0], [10.0, 10.0]) @@ -22,10 +22,10 @@ x = fill(0.5, 2) grad5 = FDM.grad(central_fdm(5, 1), approxf, x)[1] @test val1 ≈ val2 ≈ val3 ≈ val4 @test grad1 ≈ grad2 ≈ grad3 ≈ grad4 ≈ grad5 - for i in 1:100 - y = rand(2)*10 + for i = 1:100 + y = rand(2) * 10 H = Zygote.hessian(approxf, y) - for i in 1:2, j in 1:2 + for i = 1:2, j = 1:2 if i == j @test H[i, j] >= 0 else @@ -43,10 +43,10 @@ x = fill(0.5, 2) grad5 = FDM.grad(central_fdm(5, 1), approxf, x)[1] @test val1 ≈ val2 ≈ val3 ≈ val4 @test grad1 ≈ grad2 ≈ grad3 ≈ grad4 ≈ grad5 - for i in 1:100 - y = rand(2)*10 + for i = 1:100 + y = rand(2) * 10 H = Zygote.hessian(approxf, y) - for i in 1:2, j in 1:2 + for i = 1:2, j = 1:2 if i == j @test H[i, j] >= 0 else @@ -59,11 +59,13 @@ x = fill(0.5, 2) approxf = NonconvexMMA.MMAApprox(exactf, x) val1, jac1 = NonconvexCore.value_jacobian(exactf, x) val2 = exactf(x) - jac2 = vcat([Zygote.gradient(x -> exactf(x)[i], x)[1]' for i in 1:length(val2)]...) + jac2 = vcat([Zygote.gradient(x -> exactf(x)[i], x)[1]' for i = 1:length(val2)]...) val3, jac3 = NonconvexCore.value_jacobian(approxf, x) val4 = approxf(x) - jac4 = vcat([Zygote.gradient(x -> approxf(x)[i], x)[1]' for i in 1:length(val2)]...) - jac5 = vcat([FDM.grad(central_fdm(5, 1), x -> approxf(x)[i], x)[1]' for i in 1:length(val2)]...) + jac4 = vcat([Zygote.gradient(x -> approxf(x)[i], x)[1]' for i = 1:length(val2)]...) + jac5 = vcat( + [FDM.grad(central_fdm(5, 1), x -> approxf(x)[i], x)[1]' for i = 1:length(val2)]..., + ) @test val1 ≈ val2 ≈ val3 ≈ val4 @test jac1 ≈ jac2 ≈ jac3 ≈ jac4 ≈ jac5 @@ -106,7 +108,7 @@ end grad2 = FDM.grad(central_fdm(5, 1), x -> dualobj(x, λ), optimalx)[1] @test grad1 ≈ grad2 @test all(1:length(x)) do j - abs(grad1[j]) <= 1e-6 || + abs(grad1[j]) <= 1e-6 || optimalx[j] == 0 && grad1[j] > 0 || optimalx[j] == 10 && grad1[j] < 0 end diff --git a/test/mma.jl b/test/mma.jl index 305ca44..3b2be39 100755 --- a/test/mma.jl +++ b/test/mma.jl @@ -1,7 +1,7 @@ using NonconvexMMA, LinearAlgebra, Test, Zygote f(x::AbstractVector) = x[2] < 0 ? Inf : sqrt(x[2]) -g(x::AbstractVector, a, b) = (a*x[1] + b)^3 - x[2] +g(x::AbstractVector, a, b) = (a * x[1] + b)^3 - x[2] @testset "Simple constraints" begin m = Model(f) @@ -17,8 +17,8 @@ g(x::AbstractVector, a, b) = (a*x[1] + b)^3 - x[2] convcriteria, ) r = NonconvexMMA.optimize(m, alg, [1.234, 2.345], options = options) - @test abs(r.minimum - sqrt(8/27)) < 1e-6 - @test norm(r.minimizer - [1/3, 8/27]) < 1e-6 + @test abs(r.minimum - sqrt(8 / 27)) < 1e-6 + @test norm(r.minimizer - [1 / 3, 8 / 27]) < 1e-6 end end end @@ -36,8 +36,8 @@ end convcriteria, ) r = NonconvexMMA.optimize(m, alg, [1.234, 2.345], options = options) - @test abs(r.minimum - sqrt(8/27)) < 1e-6 - @test norm(r.minimizer - [1/3, 8/27]) < 1e-6 + @test abs(r.minimum - sqrt(8 / 27)) < 1e-6 + @test norm(r.minimizer - [1 / 3, 8 / 27]) < 1e-6 end end end @@ -57,8 +57,8 @@ end convcriteria, ) r = NonconvexMMA.optimize(m, alg, [1.234, 2.345], options = options) - @test abs(r.minimum - sqrt(8/27)) < 1e-6 - @test norm(r.minimizer - [1/3, 8/27]) < 1e-6 + @test abs(r.minimum - sqrt(8 / 27)) < 1e-6 + @test norm(r.minimizer - [1 / 3, 8 / 27]) < 1e-6 end end end @@ -76,8 +76,8 @@ end convcriteria, ) r = NonconvexMMA.optimize(m, alg, [1.234, 2.345], options = options) - @test abs(r.minimum - sqrt(8/27)) < 1e-6 - @test norm(r.minimizer - [1/3, 8/27]) < 1e-6 + @test abs(r.minimum - sqrt(8 / 27)) < 1e-6 + @test norm(r.minimizer - [1 / 3, 8 / 27]) < 1e-6 end end end @@ -95,8 +95,8 @@ end convcriteria, ) r = NonconvexMMA.optimize(m, alg, [1.234, 2.345], options = options) - @test abs(r.minimum - sqrt(8/27)) < 1e-6 - @test norm(r.minimizer - [1/3, 8/27]) < 1e-6 + @test abs(r.minimum - sqrt(8 / 27)) < 1e-6 + @test norm(r.minimizer - [1 / 3, 8 / 27]) < 1e-6 end end end @@ -112,13 +112,10 @@ end add_ineq_constraint!(m, x -> g(x, -1, 1)) for convcriteria in (KKTCriteria(), IpoptCriteria()) - options = MMAOptions(; - tol = Tolerance(kkt = 1e-6, f = 0.0), - s_init = 0.1, - convcriteria, - ) + options = + MMAOptions(; tol = Tolerance(kkt = 1e-6, f = 0.0), s_init = 0.1, convcriteria) r = NonconvexMMA.optimize(m, MMA02(), [0.4, 0.5], options = options) - @test abs(r.minimum - sqrt(8/27)) < 1e-6 - @test norm(r.minimizer - [1/3, 8/27]) < 1e-6 + @test abs(r.minimum - sqrt(8 / 27)) < 1e-6 + @test norm(r.minimizer - [1 / 3, 8 / 27]) < 1e-6 end end diff --git a/test/runtests.jl b/test/runtests.jl index ae74d59..c3bfe2b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,4 +1,8 @@ using SafeTestsets -@safetestset "MMA approximation" begin include("approximation.jl") end -@safetestset "MMA algorithms" begin include("mma.jl") end +@safetestset "MMA approximation" begin + include("approximation.jl") +end +@safetestset "MMA algorithms" begin + include("mma.jl") +end