Huge Make-Over, Version 0.0.2

Project.toml

@@ -1,9 +1,10 @@
 name = "Compromise"
 uuid = "254bc946-86ae-484f-a9da-8147cb79ba93"
 authors = ["Manuel Berkemeier <[email protected]> and contributors"]
-version = "0.0.1"
+version = "0.0.2"
 
 [deps]
+ConcurrentUtils = "3df5f688-6c4c-4767-8685-17f5ad261477"
 Dictionaries = "85a47980-9c8c-11e8-2b9f-f7ca1fa99fb4"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 ElasticArrays = "fdbdab4c-e67f-52f5-8c3f-e7b388dad3d4"
@@ -14,9 +15,11 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 NLopt = "76087f3c-5699-56af-9a33-bf431cd00edd"
 Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
+Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+StridedViews = "4db3bf67-4bd7-4b4e-b153-31dc3fb37143"
 StructHelpers = "4093c41a-2008-41fd-82b8-e3f9d02b504f"
 
 [weakdeps]
@@ -30,7 +33,8 @@ julia = "1"
 
 [extras]
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test", "ForwardDiff"]
+test = ["SafeTestsets", "Test", "ForwardDiff"]

README.md

@@ -2,8 +2,6 @@
 
 [![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://manuelbb-upb.github.io/Compromise.jl/stable/)
 [![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://manuelbb-upb.github.io/Compromise.jl/dev/)
-[![Build Status](https://github.com/manuelbb-upb/Compromise.jl/actions/workflows/CI.yml/badge.svg?branch=main)](https://github.com/manuelbb-upb/Compromise.jl/actions/workflows/CI.yml?query=branch%3Amain)
-[![Coverage](https://codecov.io/gh/manuelbb-upb/Compromise.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/manuelbb-upb/Compromise.jl)
 
 ## About “CoMPrOMISE”
 
@@ -18,11 +16,68 @@ Box constraints and linear constraints are supported and passed down to an inner
 Nonlinear constraint functions can be modelled and are dealt with by incorporating a filter.
 They are *relaxable*, i.e., all other functions must be computable even when the
 constraints are violated.
+### Release Notes
+
+I don't really keep up a consistent versioning scheme.
+But the changes in this section have been significant enough to warrant some comments.
+
+#### Version 0.0.2
+
+##### RBF Surrogate Models
+The internals of how RBF Surrogates are constructed has been redone.
+As before, the construction is based off “Derivative-Free Optimization Algorithms For Computationally Expensive Functions,” (Wild, 2009).
+In the old version, I did not really care for the matrix factorizations.
+Finding a poised set of points for fully-linear interpolation needs repeated QR factorizations
+of the point matrix.
+Afterwards, additional points are found by updating the Cholesky factorization of
+some symmetric matrix product involving the RBF Gram matrix.
+
+* To save allocations, the QR factorizations now make use of structures similar to
+those in `FastLapackInterface`.
+  Once I manage to [make a pull request](https://github.com/DynareJulia/FastLapackInterface.jl/issues/40)
+  to avoid even more allocations, we can also make `FastLapackInterface` a dependency.
+* The Cholesky factorization is now updated by using the formulas from the reference,
+  and the factors are used to compute the surrogate coefficients.
+* In both cases, buffers are pre-allocated to support a maximum number of interpolation
+  points, and we work with views mainly.
+  Such temporary buffers are stored in `RBFTrainingBuffers`.
+* An `RBFModel` now only needs `RBFParameters` for successful training and
+  reproducible evaluation.
+  Most importantly, evaluation is decoupled from the `RBFDatabase`!!
+  In older versions, we would view into the database to query interpolation points.
+  These are now copied instead, so that changes to the database don't invalidate a model.
+* With [commit cc709fa](https://github.com/manuelbb-upb/Compromise.jl/commit/cc709fa0391d4a796b543a0733c31fb6f2e6ad46)
+  we can thus share a database in multiple optimization runs.
+* As an aside, there is a whole new backend for RBFs, which can be used in a standalone manner, too.
+
+For most of the RBF related changes, [commit ab5cba8](https://github.com/manuelbb-upb/Compromise.jl/commit/ab5cba8f3d4ba39cd4bf2757a072bb655b4f0cc2)
+is most relevant.
+
+##### Other changes
+
+* There likely was a bug in how descent directions were computed.
+  In old versions, I tried to avoid the computation of an initial steplength by making it part
+  of the descent direction sub-problem, but accounting for the change in criticality measure
+  did not work out well.
+  [Commit f1386c2](https://github.com/manuelbb-upb/Compromise.jl/commit/f1386c29e19b7d5c3bad28fc03af8024015666c5)
+  makes things look a bit more elegant.
+* At some point in time, I messed up affine scaling. Should work now, and there is tests for it.
+* Threaded parallel execution is now supported internally (but improvised).
+* Lots of tests.
+* Changes to `AbstractNonlinearOperator` interface.
+  A new `AbstractNonlinearOperatorWrapper <: AbstractNonlinearOperator`.
+* New defaults in `AlgorithmOptions`. Stopping based mainly on minimum radius.
+* A new return type (`ReturnObject`).
 
 ## Example
 
 ### Constrained Two-Parabolas Problem
 
+This first example uses Taylor Polynomials to approximate the function locally.
+For that we need a gradient function.
+But we also show, how to add functions with derivative-free surrogates -- in this case
+nonlinear constraints.
+
 First we load the optimizer package, “Compromise.jl”:
 
 ````julia
@@ -46,7 +101,7 @@ nonlinear inequality constraints (`nl_ineq_constraints`)
 and nonlinear equality constraints (`nl_eq_constraints`).
 By default, these fields default to `nothing`.
 Alternatively, they need objects of type `Compromise.NonlinearFunction`.
-We have helpers to support normal julia functions.
+We have helpers to support normal Julia functions.
 For example, consider this vector-to-vector function:
 
 ````julia
@@ -75,9 +130,8 @@ To add these functions to `mop`, we call the helper method
 `add_objectives` and also specify the model class to be used.
 There are shorthand symbols, for example `:exact` or `taylor1`,
 for objectives with known gradients.
-Alternatively, give some `Compromise.A_bstractSurrogateModelConfig`.
 We also have to tell the optimizer about the function signature.
-`func_iip=true` would imply an in-place objective with sigunature
+`func_iip=true` would imply an in-place objective with signature
 `objective_function!(fx, x)`.
 `dim_out` is a mandatory argument.
 
@@ -108,24 +162,6 @@ for derivative calculations.
 Instead, you can also use derivative-free models, such as
 radial basis function (RBF) models.
 
-### Excursion: RBF Kernels
-By default, a cubic kernel is used, if we use the `:rbf`
-option with `add_nl_ineq_constraints!`.
-To use the Gaussian kernel ``φ_ε(r) = \\exp(-(εr)^2)``
-with fixed shape paramater `10`, do
-```julia
-rbf_config = RBFConfig(; kernel=GaussianKernel(10))
-```
-There is other cool features. If you want to vary the shape parameter
-with the trust-region radius, you can give a function ``Δ ↦ ε(Δ)``
-to `GaussianKernel` (as well as `InverseMultiQuadricKernel`).
-
-````julia
-gk = GaussianKernel(del -> min(1/del, 1_000))
-Compromise.RBFModels.shape_paramater(gk, 1) # equals 1.0
-Compromise.RBFModels.shape_paramater(gk, 0.0005) ## equals 1000
-````
-
 For now, we stick with the fixed shape parameter and finalize
 our problem:
 
@@ -139,25 +175,152 @@ The `MutableMOP` is turned into a `TypedMOP` during initialization.
 We can thus simply pass it to `optimize`:
 
 ````julia
-final_vals, ret = optimize(mop, [-2.0, 0.5])
+ret = optimize(mop, [-2.0, 0.5])
+````
+
+`ret` is the return object.
+You can query it using functions like `opt_vars` etc.
+Final argument vector:
+
+````julia
+opt_vars(ret)
+````
+
+Final value vector:
+
+````julia
+opt_objectives(ret)
+````
+
+Final constraint vector:
+
+````julia
+opt_nl_ineq_constraints(ret)
+````
+
+### More RBF Options
+Instead of passing `:rbf`, you can also pass an `RBFConfig`.
+To use the Gaussian kernel:
+
+````julia
+cfg = RBFConfig(; kernel=GaussianKernel())
+````
+
+Or the inverse multiquadric:
+
+````julia
+cfg = RBFConfig(; kernel=InverseMultiQuadricKernel())
+````
+
+Then:
+
+````julia
+mop = MutableMOP(; num_vars=2)
+add_objectives!(
+    mop, objective_function, cfg; dim_out=2, func_iip=false,
+)
+ret = optimize(mop, [-2.0, 0.5])
+````
+
+See the docstring for more options.
+
+#### Sharing an `RBFDatabase`
+Normally, each optimization run initializes a new database.
+But a database is only ever referenced.
+We can thus pre-initialize a database and use it in multiple runs:
+
+````julia
+rbf_db = Compromise.RBFModels.init_rbf_database(2, 2, nothing, nothing, Float64)
+cfg = RBFConfig(; database=rbf_db)
+````
+
+Set up problem:
+
+````julia
+mop = MutableMOP(; num_vars=2)
+objf_counter = Ref(0)
+function counted_objf(x)
+    global objf_counter[] += 1
+    return objective_function(x)
+end
+add_objectives!(
+    mop, counted_objf, cfg; dim_out=2, func_iip=false,
+)
+````
+
+First run:
+
+````julia
+ret = optimize(mop, [-2.0, 0.5])
+objf_counter[]
 ````
 
-`ret` is the return code. If no budget or time based stopping criterion
-was used, then `CRITICAL` is a nice return value, indicating a first-order
-critical point.
+Second run:
+
+````julia
+ret = optimize(mop, [-2.0, 0.0])
+objf_counter[]
+````
+
+##### Parallelism
+The RBF update algorithm has a lock to access the database in a safe way (?) when
+multiple optimization runs are done concurrently.
+There even is an “algorithm” for this:
+
+````julia
+X0 = [
+    -2.0    -2.0    0.0
+    0.5     0.0     0.0
+]
+opts = Compromise.ThreadedOuterAlgorithmOptions(;
+    inner_opts=AlgorithmOptions(;
+        stop_delta_min = 1e-8,
+    )
+)
+rets = Compromise.optimize_with_algo(mop, opts, X0)
+````
+
+### Stopping based on Number of Function Evaluations
+The restriction of evaluation budget is a property of the evaluators.
+Because of this, it is not configurable with `AlgorithmOptions`.
+You can pass `max_func_calls` as a keyword argument to `add_objectives!` and similar functions.
+Likewise, `max_grad_calls` restricts the number of gradient calls,
+`max_hess_calls` limits Hessian computations.
+
+For historic reasons, the count is kept between runs.
+To reset the count between runs (sequential or parallel), indicate it when setting up
+the MOP.
+
+````julia
+mop = MutableMOP(; num_vars=2, reset_call_counters=false)   # default
+add_objectives!(
+    mop, objective_function, :rbf; dim_out=2, func_iip=false, max_func_calls=10
+)
+ret1 = optimize(mop, [-2, .5])
+````
+
+Now, there is no budget left for a second run:
+
+````julia
+ret2 = optimize(mop, [-2, -.5])
+ismissing(opt_vars(ret2))
+````
+
+Here is a remedy:
+
+````julia
+mop.reset_call_counters=true
+ret1 = optimize(mop, [-2, .5])
+````
 
-`final_vals` is a `ValueArrays` object holding values at the final iteration
-site.
-`final_vals.x` is the parameter vector, `final_vals.fx` has the
-objective values.
-The nonlinear equality constraints are `final_vals.hx`, the inequality
-constraints are `final_vals.gx`...
+Now, there is no budget left for a second run:
 
 ````julia
-final_vals.x, final_vals.fx
+ret2 = optimize(mop, [-2, -.5])
+!ismissing(opt_vars(ret2))
 ````
 
-### Automatic Diffentiation
+### Automatic Differentiation
 There is an optional `ForwardDiff` extension.
 To use a derivative-based model without specifying the gradients by-hand,
 first load `ForwardDiff`.