From 4edbbd5a96f40bc0ce18325367f80cf910e829f3 Mon Sep 17 00:00:00 2001 From: James Rising Date: Wed, 30 Nov 2022 16:04:22 -0500 Subject: [PATCH 1/4] add paramtrans --- Project.toml | 1 - src/OptiMimi.jl | 71 +++++++++++++++++++++++++++------------------ src/linproghouse.jl | 1 - 3 files changed, 43 insertions(+), 30 deletions(-) diff --git a/Project.toml b/Project.toml index 22040ac..087e824 100644 --- a/Project.toml +++ b/Project.toml @@ -3,7 +3,6 @@ uuid = "61d5e40c-f051-5c04-90a0-b0f47eb95ad5" [deps] BlackBoxOptim = "a134a8b2-14d6-55f6-9291-3336d3ab0209" -Clp = "e2554f3b-3117-50c0-817c-e040a3ddf72d" DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" DiffResults = "163ba53b-c6d8-5494-b064-1a9d43ac40c5" DiffRules = "b552c78f-8df3-52c6-915a-8e097449b14b" diff --git a/src/OptiMimi.jl b/src/OptiMimi.jl index 761e1f2..105526d 100644 --- a/src/OptiMimi.jl +++ b/src/OptiMimi.jl @@ -24,6 +24,7 @@ mutable struct OptimizationProblem names::Vector{Symbol} opt::Opt constraints::Vector{Function} + paramtrans::Union{Vector{Function}, Nothing} end mutable struct BlackBoxOptimizationProblem{T} @@ -34,6 +35,7 @@ mutable struct BlackBoxOptimizationProblem{T} objective::Function lowers::Vector{T} uppers::Vector{T} + paramtrans::Union{Vector{Function}, Nothing} end mutable struct LinprogOptimizationProblem{T} @@ -64,24 +66,30 @@ end """ Initialize a model with the results from an optimization. """ -function setparameters(model::Model, components::Vector{Symbol}, names::Vector{Symbol}, xx::Vector) - startindex = 1 - for (ii, len, isscalar) in Channel((chnl) -> nameindexes(model, components, names, chnl)) - if isscalar - set_or_update_param!(model, components[ii], names[ii], xx[startindex]) - else - shape = getdims(model, components[ii], names[ii]) - reshaped = reshape(collect(model.md.number_type, xx[startindex:(startindex+len - 1)]), tuple(shape...)) - set_or_update_param!(model, components[ii], names[ii], reshaped) +function setparameters(model::Model, components::Vector{Symbol}, names::Vector{Symbol}, xx::Vector, paramtrans::Union{Vector{Function}, Nothing}) + if isnothing(paramtrans) + startindex = 1 + for (ii, len, isscalar) in Channel((chnl) -> nameindexes(model, components, names, chnl)) + if isscalar + set_or_update_param!(model, components[ii], names[ii], xx[startindex]) + else + shape = getdims(model, components[ii], names[ii]) + reshaped = reshape(collect(model.md.number_type, xx[startindex:(startindex+len - 1)]), tuple(shape...)) + set_or_update_param!(model, components[ii], names[ii], reshaped) + end + startindex += len + end + else + for (ii, len, isscalar) in Channel((chnl) -> nameindexes(model, components, names, chnl)) + set_or_update_param!(model, components[ii], names[ii], paramtrans[ii](xx)) end - startindex += len end end -function sensitivity(model::Model, component::Symbol, parameter::Symbol, objective::Function, points::Vector{Float64}) +function sensitivity(model::Model, component::Symbol, parameter::Symbol, objective::Function, points::Vector{Float64}, paramtrans::Union{Vector{Function}, Nothing}=nothing) results = [] for point in points - setparameters(model, [component], [parameter], [point]) + setparameters(model, [component], [parameter], [point], paramtrans) run(model) push!(results, objective(model)) end @@ -90,7 +98,7 @@ function sensitivity(model::Model, component::Symbol, parameter::Symbol, objecti end """Generate the form of objective function used by the optimization, taking parameters rather than a model.""" -function unaryobjective(model::Model, components::Vector{Symbol}, names::Vector{Symbol}, objective::Function) +function unaryobjective(model::Model, components::Vector{Symbol}, names::Vector{Symbol}, objective::Function, paramtrans::Union{Vector{Function}, Nothing}) function my_objective(xx::Vector) if allverbose println(xx) @@ -99,7 +107,7 @@ function unaryobjective(model::Model, components::Vector{Symbol}, names::Vector{ global objevals objevals += 1 - setparameters(model, components, names, xx) + setparameters(model, components, names, xx, paramtrans) run(model) objective(model) end @@ -108,8 +116,8 @@ function unaryobjective(model::Model, components::Vector{Symbol}, names::Vector{ end """Create an NLopt-style objective function which does not use its grad argument.""" -function gradfreeobjective(model::Model, components::Vector{Symbol}, names::Vector{Symbol}, objective::Function) - myunaryobjective = unaryobjective(model, components, names, objective) +function gradfreeobjective(model::Model, components::Vector{Symbol}, names::Vector{Symbol}, objective::Function, paramtrans::Union{Vector{Function}, Nothing}) + myunaryobjective = unaryobjective(model, components, names, objective, paramtrans) function myobjective(xx::Vector, grad::Vector) myunaryobjective(xx) end @@ -118,8 +126,8 @@ function gradfreeobjective(model::Model, components::Vector{Symbol}, names::Vect end """Create an NLopt-style objective function which computes an autodiff gradient.""" -function autodiffobjective(model::Model, components::Vector{Symbol}, names::Vector{Symbol}, objective::Function) - myunaryobjective = unaryobjective(model, components, names, objective) +function autodiffobjective(model::Model, components::Vector{Symbol}, names::Vector{Symbol}, objective::Function, paramtrans::Union{Vector{Function}, Nothing}) + myunaryobjective = unaryobjective(model, components, names, objective, paramtrans) function myobjective(xx::Vector, grad::Vector) out = DiffResults.GradientResult(xx) ForwardDiff.gradient!(out, myunaryobjective, xx) @@ -160,8 +168,14 @@ function expandlimits(model::Model, components::Vector{Symbol}, names::Vector{Sy end """Setup an optimization problem.""" -function problem(model::Model, components::Vector{Symbol}, names::Vector{Symbol}, lowers::Vector{T}, uppers::Vector{T}, objective::Function; constraints::Vector{Function}=Function[], algorithm::Symbol=:LN_COBYLA_OR_LD_MMA) where T <: Real - my_lowers, my_uppers, totalvars = expandlimits(model, components, names, lowers, uppers) +function problem(model::Model, components::Vector{Symbol}, names::Vector{Symbol}, lowers::Vector{T}, uppers::Vector{T}, objective::Function; constraints::Vector{Function}=Function[], algorithm::Symbol=:LN_COBYLA_OR_LD_MMA, paramtrans::Union{Vector{Function}, Nothing}=nothing) where T <: Real + if isnothing(paramtrans) + my_lowers, my_uppers, totalvars = expandlimits(model, components, names, lowers, uppers) + else + my_lowers = lowers + my_uppers = uppers + totalvars = length(lowers) + end if algorithm == :GUROBI_LINPROG # Make no changes to objective! @@ -171,10 +185,10 @@ function problem(model::Model, components::Vector{Symbol}, names::Vector{Symbol} end if string(algorithm)[2] == 'N' warn("Model is autodifferentiable, but optimizing using a derivative-free algorithm.") - myobjective = gradfreeobjective(model, components, names, objective) + myobjective = gradfreeobjective(model, components, names, objective, paramtrans) else println("Using AutoDiff objective.") - myobjective = autodiffobjective(model, components, names, objective) + myobjective = autodiffobjective(model, components, names, objective, paramtrans) end else if algorithm == :LN_COBYLA_OR_LD_MMA @@ -184,7 +198,7 @@ function problem(model::Model, components::Vector{Symbol}, names::Vector{Symbol} algorithm = :LN_COBYLA end - myobjective = gradfreeobjective(model, components, names, objective) + myobjective = gradfreeobjective(model, components, names, objective, paramtrans) end if algorithm == :GUROBI_LINPROG @@ -195,7 +209,7 @@ function problem(model::Model, components::Vector{Symbol}, names::Vector{Symbol} warn("Functional constraints not supported for BBO algorithms.") end - BlackBoxOptimizationProblem(algorithm, model, components, names, myobjective, my_lowers, my_uppers) + BlackBoxOptimizationProblem(algorithm, model, components, names, myobjective, my_lowers, my_uppers, paramtrans) else opt = Opt(algorithm, totalvars) lower_bounds!(opt, my_lowers) @@ -207,7 +221,7 @@ function problem(model::Model, components::Vector{Symbol}, names::Vector{Symbol} for constraint in constraints let this_constraint = constraint function my_constraint(xx::Vector, grad::Vector) - setparameters(model, components, names, xx) + setparameters(model, components, names, xx, paramtrans) run(model) this_constraint(model) end @@ -216,7 +230,7 @@ function problem(model::Model, components::Vector{Symbol}, names::Vector{Symbol} end end - OptimizationProblem(model, components, names, opt, constraints) + OptimizationProblem(model, components, names, opt, constraints, paramtrans) end end end @@ -236,7 +250,7 @@ function solution(optprob::OptimizationProblem, generator::Function; maxiter=Inf while attempts < maxiter initial = generator() - setparameters(optprob.model, optprob.components, optprob.names, initial) + setparameters(optprob.model, optprob.components, optprob.names, initial, optprob.paramtrans) run(optprob.model) valid = true @@ -269,6 +283,7 @@ function solution(optprob::OptimizationProblem, generator::Function; maxiter=Inf end (minf,minx,ret) = optimize(optprob.opt, initial) + print(ret) (minf, minx) end @@ -304,7 +319,7 @@ function solution(optprob::LinprogOptimizationProblem, verbose=false) end if optprob.model.md.number_type == Number - myobjective = unaryobjective(optprob.model, optprob.components, optprob.names, optprob.objective) + myobjective = unaryobjective(optprob.model, optprob.components, optprob.names, optprob.objective, nothing) f, b, A = lpconstraints(optprob.model, optprob.components, optprob.names, myobjective, objectiveconstraints) else f, b, A = lpconstraints(optprob.model, optprob.components, optprob.names, optprob.objective, optprob.objectiveconstraints) diff --git a/src/linproghouse.jl b/src/linproghouse.jl index 5b78827..7689f54 100644 --- a/src/linproghouse.jl +++ b/src/linproghouse.jl @@ -1,6 +1,5 @@ using MathProgBase using DataFrames -using Clp using SparseArrays, LinearAlgebra import Base.*, Base.-, Base.+, Base./, Base.max From f7b33005e141a68c85699a98e1a2396ad6494f16 Mon Sep 17 00:00:00 2001 From: James Rising Date: Wed, 30 Nov 2022 16:04:52 -0500 Subject: [PATCH 2/4] smooth optimization --- src/smooth.jl | 86 +++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 86 insertions(+) create mode 100644 src/smooth.jl diff --git a/src/smooth.jl b/src/smooth.jl new file mode 100644 index 0000000..9d81d8c --- /dev/null +++ b/src/smooth.jl @@ -0,0 +1,86 @@ +using Interpolations +using Dierckx +using Distributions + +using NLopt +import OptiMimi.OptimizationProblem, OptiMimi.getdims, OptiMimi.gradfreeobjective +# using Plots + +function piecewisesolution_single(optprob::OptimizationProblem, objective::Function, initgen::Function, maxiter=Inf, verbose=false; reltol=1e-6) + len = getdims(optprob.model, optprob.components[1], optprob.names[1])[1] + numsegs = Int64.([2 .^(0:(floor(log2(len-1))-1)); len-1]) + + last_xx = last_soln = nothing + + for numseg in numsegs + print(numseg) + xx = trunc.(Int, range(1, stop=len, length=numseg+1)) + + keeps = [true for ii in 1:length(xx)] + if numseg > 2 + # Drop xx where [LinearInterp - CubicInterp] < reltol + linears = LinearInterpolation(last_xx, last_soln)(xx) + cubics = evaluate(Spline1D(last_xx, last_soln, k=min(3, length(last_xx) - 1)), xx) + keeps = abs.(linears .- cubics) ./ ((linears .+ cubics) / 2) .> reltol + keeps[1] = true + keeps[length(keeps)] = true + + # Drop knots only if can be linearly guessed + for ii in 2:numseg + if xx[ii] ∈ last_xx + jj = findfirst(last_xx .== xx[ii]) + guess = LinearInterpolation([last_xx[jj-1], last_xx[jj+1]], [last_soln[jj-1], last_soln[jj+1]])(xx[ii]) + if abs(guess - last_soln[jj]) / last_soln[jj] > reltol + keeps[ii] = true + end + end + end + + println(keeps) + end + + paramtrans = convert(Vector{Function}, [yy -> LinearInterpolation(xx[keeps], yy)(1:len)]) + + opt = Opt(optprob.opt.algorithm, sum(keeps)) + lower_bounds!(opt, optprob.opt.lower_bounds[xx[keeps]]) + upper_bounds!(opt, optprob.opt.upper_bounds[xx[keeps]]) + xtol_rel!(opt, reltol) + + myobjective = gradfreeobjective(optprob.model, optprob.components, optprob.names, objective, paramtrans) + + max_objective!(opt, myobjective) + + for constraint in optprob.constraints + let this_constraint = constraint + function my_constraint(xx::Vector, grad::Vector) + setparameters(optprob.model, optprob.components, optprob.names, xx, paramtrans) + run(optprob.model) + this_constraint(optprob.model) + end + + inequality_constraint!(opt, my_constraint) + end + end + + if sum(keeps) == len - 1 + subprob = OptimizationProblem(optprob.model, optprob.components, optprob.names, opt, optprob.constraints, + nothing) + else + subprob = OptimizationProblem(optprob.model, optprob.components, optprob.names, opt, optprob.constraints, + paramtrans) + end + + soln = solution(subprob, initgen(sum(keeps), xx[keeps], last_xx, last_soln)) + + if numseg > 2 + last_soln = LinearInterpolation(last_xx, last_soln)(xx) + last_soln[keeps] = soln[2] + else + last_soln = soln[2] + end + last_xx = xx + println([last_xx, last_soln]) + end + + last_soln +end From f3268f98c06f473e449a1555b7014caa07c5cb15 Mon Sep 17 00:00:00 2001 From: James Rising Date: Tue, 13 Dec 2022 13:45:51 -0500 Subject: [PATCH 3/4] uncertain and smooth --- Project.toml | 2 + src/OptiMimi.jl | 5 +- src/smooth.jl | 23 +++++---- src/uncertainty.jl | 117 ++++++++++++++++++++++++++++----------------- 4 files changed, 92 insertions(+), 55 deletions(-) diff --git a/Project.toml b/Project.toml index 087e824..528255a 100644 --- a/Project.toml +++ b/Project.toml @@ -4,10 +4,12 @@ uuid = "61d5e40c-f051-5c04-90a0-b0f47eb95ad5" [deps] BlackBoxOptim = "a134a8b2-14d6-55f6-9291-3336d3ab0209" DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" +Dierckx = "39dd38d3-220a-591b-8e3c-4c3a8c710a94" DiffResults = "163ba53b-c6d8-5494-b064-1a9d43ac40c5" DiffRules = "b552c78f-8df3-52c6-915a-8e097449b14b" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" +Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MathProgBase = "fdba3010-5040-5b88-9595-932c9decdf73" Mimi = "e4e893b0-ee5e-52ea-8111-44b3bdec128c" diff --git a/src/OptiMimi.jl b/src/OptiMimi.jl index 105526d..5beb446 100644 --- a/src/OptiMimi.jl +++ b/src/OptiMimi.jl @@ -22,6 +22,7 @@ mutable struct OptimizationProblem model::Model components::Vector{Symbol} names::Vector{Symbol} + objective::Function opt::Opt constraints::Vector{Function} paramtrans::Union{Vector{Function}, Nothing} @@ -49,6 +50,8 @@ mutable struct LinprogOptimizationProblem{T} exuppers::Vector{T} end +include("smooth.jl") + BlackBoxAlgorithms = [:separable_nes, :xnes, :dxnes] """Returns (ii, len, isscalar) with the index of each symbol and its length.""" @@ -230,7 +233,7 @@ function problem(model::Model, components::Vector{Symbol}, names::Vector{Symbol} end end - OptimizationProblem(model, components, names, opt, constraints, paramtrans) + OptimizationProblem(model, components, names, objective, opt, constraints, paramtrans) end end end diff --git a/src/smooth.jl b/src/smooth.jl index 9d81d8c..01a575a 100644 --- a/src/smooth.jl +++ b/src/smooth.jl @@ -1,12 +1,12 @@ using Interpolations -using Dierckx +import Dierckx using Distributions using NLopt -import OptiMimi.OptimizationProblem, OptiMimi.getdims, OptiMimi.gradfreeobjective +# import OptiMimi.OptimizationProblem, OptiMimi.getdims, OptiMimi.gradfreeobjective # using Plots -function piecewisesolution_single(optprob::OptimizationProblem, objective::Function, initgen::Function, maxiter=Inf, verbose=false; reltol=1e-6) +function piecewisesolution_single(optprob::OptimizationProblem, initgen::Function, maxiter=Inf, verbose=false; reltol=1e-6) len = getdims(optprob.model, optprob.components[1], optprob.names[1])[1] numsegs = Int64.([2 .^(0:(floor(log2(len-1))-1)); len-1]) @@ -20,7 +20,7 @@ function piecewisesolution_single(optprob::OptimizationProblem, objective::Funct if numseg > 2 # Drop xx where [LinearInterp - CubicInterp] < reltol linears = LinearInterpolation(last_xx, last_soln)(xx) - cubics = evaluate(Spline1D(last_xx, last_soln, k=min(3, length(last_xx) - 1)), xx) + cubics = Dierckx.evaluate(Dierckx.Spline1D(last_xx, last_soln, k=min(3, length(last_xx) - 1)), xx) keeps = abs.(linears .- cubics) ./ ((linears .+ cubics) / 2) .> reltol keeps[1] = true keeps[length(keeps)] = true @@ -46,7 +46,7 @@ function piecewisesolution_single(optprob::OptimizationProblem, objective::Funct upper_bounds!(opt, optprob.opt.upper_bounds[xx[keeps]]) xtol_rel!(opt, reltol) - myobjective = gradfreeobjective(optprob.model, optprob.components, optprob.names, objective, paramtrans) + myobjective = gradfreeobjective(optprob.model, optprob.components, optprob.names, optprob.objective, paramtrans) max_objective!(opt, myobjective) @@ -63,14 +63,17 @@ function piecewisesolution_single(optprob::OptimizationProblem, objective::Funct end if sum(keeps) == len - 1 - subprob = OptimizationProblem(optprob.model, optprob.components, optprob.names, opt, optprob.constraints, - nothing) + subprob = OptimizationProblem(optprob.model, optprob.components, optprob.names, optprob.objective, opt, optprob.constraints, nothing) else - subprob = OptimizationProblem(optprob.model, optprob.components, optprob.names, opt, optprob.constraints, - paramtrans) + subprob = OptimizationProblem(optprob.model, optprob.components, optprob.names, optprob.objective, opt, optprob.constraints, paramtrans) end - soln = solution(subprob, initgen(sum(keeps), xx[keeps], last_xx, last_soln)) + if isnothing(last_xx) + lin_yy = nothing + else + lin_yy = LinearInterpolation(last_xx, last_soln)(xx[keeps]) + end + soln = solution(subprob, initgen(sum(keeps), lin_yy)) if numseg > 2 last_soln = LinearInterpolation(last_xx, last_soln)(xx) diff --git a/src/uncertainty.jl b/src/uncertainty.jl index 2e7613a..05522b5 100644 --- a/src/uncertainty.jl +++ b/src/uncertainty.jl @@ -4,12 +4,14 @@ ## Two algorithms available: ## biological: https://github.com/robertfeldt/BlackBoxOptim.jl -## sampled: Deterministic optimization under different Monte Carlo +## sampled: Deterministic optimization under different Monte Carlos +## montecarlo: Optimization of Monte Carlo mean of objective ## draws. ##### Biological Genetics ##### -using BlackBoxOptim +import BlackBoxOptim +using Random ########## @@ -17,7 +19,7 @@ mutable struct UncertainOptimizationProblem model::Model components::Vector{Symbol} names::Vector{Symbol} - montecarlo::Function + montecarlo::Union{Function, SimulationDef} objective::Function lowers::Vector{Float64} uppers::Vector{Float64} @@ -35,64 +37,91 @@ end """Setup an optimization over Monte Carlo uncertainty.""" function uncertainproblem(model::Model, components::Vector{Symbol}, names::Vector{Symbol}, lowers::Vector{Float64}, uppers::Vector{Float64}, objective::Function, montecarlo::Function=(model) -> nothing) my_lowers, my_uppers, totalvars = expandlimits(model, components, names, lowers, uppers) - my_objective = unaryobjective(model, components, names, objective) + UncertainOptimizationProblem(model, components, names, montecarlo, objective, my_lowers, my_uppers) +end + +function uncertainproblem(model::Model, components::Vector{Symbol}, names::Vector{Symbol}, lowers::Vector{Float64}, uppers::Vector{Float64}, objective::Function, sim::SimulationDef, mcnum=Int64) + function objectivetopayload(sim_inst::SimulationInstance, trialnum::Int, ntimesteps::Int, tup::Nothing) + model = sim_inst.models[1] + results = Mimi.payload(sim_inst) + results[trialnum] = objective(model) + end + + function uncertainobjective(model::Model) + Mimi.set_payload!(sim, zeros(mcnum)) + Random.seed!(0) + mcout = run(sim, model, mcnum, post_trial_func=objectivetopayload) + values = Mimi.payload(mcout) + mean(values[isfinite.(values)]) + end - UncertainOptimizationProblem(model, components, names, montecarlo, my_objective, my_lowers, my_uppers) + problem(model, components, names, lowers, uppers, uncertainobjective) end """Solve an optimization over Monte Carlo uncertainty.""" function solution(optprob::UncertainOptimizationProblem, generator::Function, algorithm::Symbol, mcperlife::Int64, samples::Int64) - function sample_objective(parameters::Vector{Float64}) - total = 0 - for iter in 1:mcperlife - optprob.montecarlo(optprob.model) - total += optprob.objective(parameters) - end + if algorithm ∈ [:biological, :sampled] + my_objective = unaryobjective(model, components, names, objective) - total / mcperlife - end + function sample_objective(parameters::Vector{Float64}) + total = 0 + for iter in 1:mcperlife + optprob.montecarlo(optprob.model) + total += optprob.objective(parameters) + end - if algorithm == :biological - res = bboptimize(p -> -sample_objective(p); SearchRange=collect(zip(optprob.lowers, optprob.uppers)), MaxFuncEvals=samples, Method=:separable_nes) - - fmean = NaN - fserr = NaN - xmean = best_candidate(res) - xserr = repmat([NaN], length(optprob.lowers)) - extra = res - elseif algorithm == :sampled - opt = Opt(:LN_COBYLA, length(optprob.lowers)) - lower_bounds!(opt, optprob.lowers) - upper_bounds!(opt, optprob.uppers) - xtol_rel!(opt, minimum(1e-6 * (optprob.uppers - optprob.lowers))) - - sampleseed = 0 - function hold_objective(parameters::Vector{Float64}, grad::Vector{Float64}) - srand(sampleseed) - sample_objective(parameters) + total / mcperlife end - max_objective!(opt, hold_objective) + if algorithm == :biological + res = BlackBoxOptim.bboptimize(p -> -sample_objective(p); SearchRange=collect(zip(optprob.lowers, optprob.uppers)), MaxFuncEvals=samples, Method=:separable_nes) - sampleprob = OptimizationProblem(optprob.model, optprob.components, optprob.names, opt, Function[]) + fmean = NaN + fserr = NaN + xmean = best_candidate(res) + xserr = repmat([NaN], length(optprob.lowers)) + extra = res + elseif algorithm == :sampled + opt = Opt(:LN_COBYLA, length(optprob.lowers)) + lower_bounds!(opt, optprob.lowers) + upper_bounds!(opt, optprob.uppers) + xtol_rel!(opt, minimum(1e-6 * (optprob.uppers - optprob.lowers))) - allparams = Vector{Float64}[] + sampleseed = 0 + function hold_objective(parameters::Vector{Float64}, grad::Vector{Float64}) + srand(sampleseed) + sample_objective(parameters) + end - for sample in 1:samples - sampleseed = sample + max_objective!(opt, hold_objective) - minf, minx = solution(sampleprob, generator, maxiter=10000) - push!(allparams, Float64[minf; minx]) - end + sampleprob = OptimizationProblem(optprob.model, optprob.components, optprob.names, optprob.objective, opt, Function[]) - allparams = transpose(hcat(allparams...)) + allparams = Vector{Float64}[] - fmean = mean(allparams[:, 1]) - fserr = std(allparams[:, 1]) / sqrt(samples) + for sample in 1:samples + sampleseed = sample + + minf, minx = solution(sampleprob, generator, maxiter=10000) + push!(allparams, Float64[minf; minx]) + end - xmean = vec(mean(allparams[:, 2:end], 1)) - xserr = vec(std(allparams[:, 2:end], 1) / sqrt(samples)) + allparams = transpose(hcat(allparams...)) + fmean = mean(allparams[:, 1]) + fserr = std(allparams[:, 1]) / sqrt(samples) + + xmean = vec(mean(allparams[:, 2:end], 1)) + xserr = vec(std(allparams[:, 2:end], 1) / sqrt(samples)) + + extra = nothing + end + elseif algorithm == :montecarlo + minf, minx = solution(subprob, generator, samples) + fmean = minf + fserr = 0 + xmean = minx + xserr = 0 extra = nothing else error("Unknown algorithm") From 7a8e1758d707af014d6a9a1d3544a215da29a128 Mon Sep 17 00:00:00 2001 From: James Rising Date: Thu, 15 Jun 2023 16:38:58 -0400 Subject: [PATCH 4/4] updates to handle smooth optimization --- src/OptiMimi.jl | 36 ++++++++++++++++++++---------------- 1 file changed, 20 insertions(+), 16 deletions(-) diff --git a/src/OptiMimi.jl b/src/OptiMimi.jl index 105526d..1c5e3c7 100644 --- a/src/OptiMimi.jl +++ b/src/OptiMimi.jl @@ -4,7 +4,7 @@ using NLopt using ForwardDiff, DiffResults using MathProgBase -import Mimi: Model +import Mimi: Model, AbstractModel export problem, solution, unaryobjective, objevals, setparameters, nameindexes, sensitivity, uncertainproblem, setsolution @@ -19,7 +19,7 @@ allverbose = false objevals = 0 mutable struct OptimizationProblem - model::Model + model::AbstractModel components::Vector{Symbol} names::Vector{Symbol} opt::Opt @@ -29,7 +29,7 @@ end mutable struct BlackBoxOptimizationProblem{T} algorithm::Symbol - model::Model + model::AbstractModel components::Vector{Symbol} names::Vector{Symbol} objective::Function @@ -39,7 +39,7 @@ mutable struct BlackBoxOptimizationProblem{T} end mutable struct LinprogOptimizationProblem{T} - model::Model + model::AbstractModel components::Vector{Symbol} names::Vector{Symbol} objective::Function @@ -52,7 +52,7 @@ end BlackBoxAlgorithms = [:separable_nes, :xnes, :dxnes] """Returns (ii, len, isscalar) with the index of each symbol and its length.""" -function nameindexes(model::Model, components::Vector{Symbol}, names::Vector{Symbol}, chnl) +function nameindexes(model::AbstractModel, components::Vector{Symbol}, names::Vector{Symbol}, chnl) for ii in 1:length(components) dims = getdims(model, components[ii], names[ii]) if length(dims) == 0 @@ -66,7 +66,7 @@ end """ Initialize a model with the results from an optimization. """ -function setparameters(model::Model, components::Vector{Symbol}, names::Vector{Symbol}, xx::Vector, paramtrans::Union{Vector{Function}, Nothing}) +function setparameters(model::AbstractModel, components::Vector{Symbol}, names::Vector{Symbol}, xx::Vector, paramtrans::Union{Vector{Function}, Nothing}) if isnothing(paramtrans) startindex = 1 for (ii, len, isscalar) in Channel((chnl) -> nameindexes(model, components, names, chnl)) @@ -86,7 +86,7 @@ function setparameters(model::Model, components::Vector{Symbol}, names::Vector{S end end -function sensitivity(model::Model, component::Symbol, parameter::Symbol, objective::Function, points::Vector{Float64}, paramtrans::Union{Vector{Function}, Nothing}=nothing) +function sensitivity(model::AbstractModel, component::Symbol, parameter::Symbol, objective::Function, points::Vector{Float64}, paramtrans::Union{Vector{Function}, Nothing}=nothing) results = [] for point in points setparameters(model, [component], [parameter], [point], paramtrans) @@ -98,7 +98,7 @@ function sensitivity(model::Model, component::Symbol, parameter::Symbol, objecti end """Generate the form of objective function used by the optimization, taking parameters rather than a model.""" -function unaryobjective(model::Model, components::Vector{Symbol}, names::Vector{Symbol}, objective::Function, paramtrans::Union{Vector{Function}, Nothing}) +function unaryobjective(model::AbstractModel, components::Vector{Symbol}, names::Vector{Symbol}, objective::Function, paramtrans::Union{Vector{Function}, Nothing}) function my_objective(xx::Vector) if allverbose println(xx) @@ -116,7 +116,7 @@ function unaryobjective(model::Model, components::Vector{Symbol}, names::Vector{ end """Create an NLopt-style objective function which does not use its grad argument.""" -function gradfreeobjective(model::Model, components::Vector{Symbol}, names::Vector{Symbol}, objective::Function, paramtrans::Union{Vector{Function}, Nothing}) +function gradfreeobjective(model::AbstractModel, components::Vector{Symbol}, names::Vector{Symbol}, objective::Function, paramtrans::Union{Vector{Function}, Nothing}) myunaryobjective = unaryobjective(model, components, names, objective, paramtrans) function myobjective(xx::Vector, grad::Vector) myunaryobjective(xx) @@ -126,7 +126,7 @@ function gradfreeobjective(model::Model, components::Vector{Symbol}, names::Vect end """Create an NLopt-style objective function which computes an autodiff gradient.""" -function autodiffobjective(model::Model, components::Vector{Symbol}, names::Vector{Symbol}, objective::Function, paramtrans::Union{Vector{Function}, Nothing}) +function autodiffobjective(model::AbstractModel, components::Vector{Symbol}, names::Vector{Symbol}, objective::Function, paramtrans::Union{Vector{Function}, Nothing}) myunaryobjective = unaryobjective(model, components, names, objective, paramtrans) function myobjective(xx::Vector, grad::Vector) out = DiffResults.GradientResult(xx) @@ -142,7 +142,7 @@ function autodiffobjective(model::Model, components::Vector{Symbol}, names::Vect end """Create a 0 point.""" -function make0(model::Model, components::Vector{Symbol}, names::Vector{Symbol}) +function make0(model::AbstractModel, components::Vector{Symbol}, names::Vector{Symbol}) initial = Float64[] for (ii, len, isscalar) in Channel((chnl) -> nameindexes(model, components, names, chnl)) append!(initial, [0. for jj in 1:len]) @@ -152,7 +152,7 @@ function make0(model::Model, components::Vector{Symbol}, names::Vector{Symbol}) end """Expand parameter constraints to full vectors for every numerical parameter.""" -function expandlimits(model::Model, components::Vector{Symbol}, names::Vector{Symbol}, lowers::Vector{T}, uppers::Vector{T}) where T <: Real +function expandlimits(model::AbstractModel, components::Vector{Symbol}, names::Vector{Symbol}, lowers::Vector{T}, uppers::Vector{T}) where T <: Real my_lowers = T[] my_uppers = T[] @@ -168,7 +168,7 @@ function expandlimits(model::Model, components::Vector{Symbol}, names::Vector{Sy end """Setup an optimization problem.""" -function problem(model::Model, components::Vector{Symbol}, names::Vector{Symbol}, lowers::Vector{T}, uppers::Vector{T}, objective::Function; constraints::Vector{Function}=Function[], algorithm::Symbol=:LN_COBYLA_OR_LD_MMA, paramtrans::Union{Vector{Function}, Nothing}=nothing) where T <: Real +function problem(model::AbstractModel, components::Vector{Symbol}, names::Vector{Symbol}, lowers::Vector{T}, uppers::Vector{T}, objective::Function; constraints::Vector{Function}=Function[], algorithm::Symbol=:LN_COBYLA_OR_LD_MMA, paramtrans::Union{Vector{Function}, Nothing}=nothing) where T <: Real if isnothing(paramtrans) my_lowers, my_uppers, totalvars = expandlimits(model, components, names, lowers, uppers) else @@ -179,7 +179,7 @@ function problem(model::Model, components::Vector{Symbol}, names::Vector{Symbol} if algorithm == :GUROBI_LINPROG # Make no changes to objective! - elseif model.md.number_type == Number + elseif hasfield(typeof(model), :md) && model.md.number_type == Number if algorithm == :LN_COBYLA_OR_LD_MMA algorithm = :LD_MMA end @@ -302,7 +302,7 @@ end """Setup an optimization problem.""" -function problem(model::Model, components::Vector{Symbol}, names::Vector{Symbol}, lowers::Vector{T}, uppers::Vector{T}, objective::Function, objectiveconstraints::Vector{Function}, matrixconstraints::Vector{MatrixConstraintSet}) where T <: Real +function problem(model::AbstractModel, components::Vector{Symbol}, names::Vector{Symbol}, lowers::Vector{T}, uppers::Vector{T}, objective::Function, objectiveconstraints::Vector{Function}, matrixconstraints::Vector{MatrixConstraintSet}) where T <: Real my_lowers, my_uppers, totalvars = expandlimits(model, components, names, lowers, uppers) LinprogOptimizationProblem(model, components, names, objective, objectiveconstraints, matrixconstraints, my_lowers, my_uppers) end @@ -318,7 +318,7 @@ function solution(optprob::LinprogOptimizationProblem, verbose=false) println("Optimizing...") end - if optprob.model.md.number_type == Number + if hasfield(typeof(optprob.model), :md) && optprob.model.md.number_type == Number myobjective = unaryobjective(optprob.model, optprob.components, optprob.names, optprob.objective, nothing) f, b, A = lpconstraints(optprob.model, optprob.components, optprob.names, myobjective, objectiveconstraints) else @@ -334,4 +334,8 @@ function solution(optprob::LinprogOptimizationProblem, verbose=false) sol.sol end +function include_smooth() + include(joinpath(dirname(@__FILE__), "smooth.jl")) +end + end # module