feat: Add slower but more accurate "optimal" path

This uses the pairwise greedy optimisation algorithm.

src/opt.jl

@@ -67,9 +67,9 @@ SemioticOpt.hooks(a::AnalyticOpt) = a.hooks
 Perform the analytic optimisation.
 """
 function SemioticOpt.iteration(f::Function, a::AnalyticOpt)
-    ixs = f(nothing)
-    Ω = a.Ω
-    ψ = a.ψ
+    ixs = f(SemioticOpt.x(a))
+    Ω = a.Ω[ixs]
+    ψ = a.ψ[ixs]
     σ = a.σ
     v = dual(Ω, ψ, σ)
     x = primal(Ω, ψ, v)
@@ -94,7 +94,7 @@ julia> AllocationOpt.optimizeanalytic(Ω, ψ, σ)
 ```
 """
 function optimizeanalytic(Ω, ψ, σ)
-    f(::Nothing) = 1:length(ψ)  # This doesn't matter; we're not using it
+    f(::Any) = 1:length(ψ)
     alg = AllocationOpt.AnalyticOpt(;
         x=zero(ψ), Ω=Ω, ψ=ψ, σ=σ, hooks=[StopWhen((a; kws...) -> kws[:i] > 1)]
     )
@@ -242,6 +242,67 @@ function optimizek(xopt, Ω, ψ, σ, k, Φ, Ψ)
     return floor.(SemioticOpt.x(sol); digits=1)
 end
 
+"""
+    optimize(::Val{:optimal}, Ω, ψ, σ, K, Φ, Ψ, g)
+
+Find the optimal solution vector given allocations of other indexers `Ω`, signals
+`ψ`, available stake `σ`, new tokens issued `Φ`, total signal `Ψ`, and gas in grt `g`.
+
+# Example
+```julia
+julia> using AllocationOpt
+julia> Ω = [1.0, 1.0]
+julia> ψ = [10.0, 10.0]
+julia> σ = 5.0
+julia> K = 2
+julia> Φ = 1.0
+julia> Ψ = 20.0
+julia> g = 0.01
+julia> xs, nonzeros, profits = AllocationOpt.optimize(
+           Val(:optimal), Ω, ψ, σ, K, Φ, Ψ, g
+       )
+([2.5; 2.5;;], Int32[2], [0.34714285714285714; 0.34714285714285714;;])
+```
+"""
+function optimize(::Val{:optimal}, Ω, ψ, σ, K, Φ, Ψ, g)
+    # Helper function to compute profit
+    obj = x -> profit.(indexingreward.(x, Ω, ψ, Φ, Ψ), g) |> sum
+
+    # Preallocate solution vectors for in-place operations
+    _x = zeros(length(ψ), 1)
+    profits = Matrix{Float64}(undef, length(ψ), 1)
+    nonzeros = Vector{Int32}(undef, 1)
+
+    f(x, ixs) = ixs
+
+    # Set up optimizer
+    function makeanalytic(x)
+        return AllocationOpt.AnalyticOpt(;
+            x=x, Ω=Ω, ψ=ψ, σ=σ, hooks=[StopWhen((a; kws...) -> kws[:i] > 1)]
+        )
+    end
+    alg = PairwiseGreedyOpt(;
+        kmax=K,
+        x=zeros(length(ψ)),
+        xinit=zeros(length(ψ)),
+        f=f,
+        a=makeanalytic,
+        hooks=[
+            StopWhen((a; kws...) -> kws[:f](kws[:z]) ≤ kws[:f](SemioticOpt.x(a)))
+            StopWhen(
+                (a; kws...) -> length(kws[:z]) == length(SemioticOpt.nonzeroixs(kws[:z]))
+            )
+        ],
+    )
+    sol = minimize!(obj, alg)
+
+    _x[:, 1] .= SemioticOpt.x(sol)
+    nonzeros[1] = _x[:, 1] |> nonzero |> length
+    profits[:, 1] .= profit.(indexingreward.(_x[:, 1], Ω, ψ, Φ, Ψ), g)
+
+    return _x, nonzeros, profits
+end
+
 """
     lipschitzconstant(ψ, Ω)
 

test/opt.jl

@@ -7,7 +7,7 @@
         Ω = [1.0, 1.0]
         ψ = [10.0, 10.0]
         σ = 5.0
-        f(::Nothing) = 1:length(ψ)
+        f(::Any) = 1:length(ψ)
         alg = AllocationOpt.AnalyticOpt(;
             x=x, Ω=Ω, ψ=ψ, σ=σ, hooks=[StopWhen((a; kws...) -> kws[:i] > 1)]
         )
@@ -104,6 +104,23 @@
             @test nonzeros == [1, 2]
             @test isapprox(profits, [[0.41 0.35]; [0.0 0.35]]; atol=0.1)
         end
+
+        @testset "optimal" begin
+            Ω = [1.0, 1.0]
+            ψ = [10.0, 10.0]
+            σ = 5.0
+            K = 2
+            Φ = 1.0
+            Ψ = 20.0
+            g = 0.01
+            xs, nonzeros, profits = AllocationOpt.optimize(
+                Val(:optimal), Ω, ψ, σ, K, Φ, Ψ, g
+            )
+            @test xs == [2.5; 2.5;;]
+            @test nonzeros == [2]
+            @test isapprox(profits, [0.35; 0.35]; atol=0.1)
+        end
+
         @testset "dispatch" begin
             config = Dict("opt_mode" => "fast")
             Ω = [1.0, 1.0]
@@ -114,6 +131,19 @@
             Ψ = 20.0
             g = 0.01
             xs, nonzeros, profits = AllocationOpt.optimize(Ω, ψ, σ, K, Φ, Ψ, g, config)
+
+            config = Dict("opt_mode" => "optimal")
+            Ω = [1.0, 1.0]
+            ψ = [10.0, 10.0]
+            σ = 5.0
+            K = 2
+            Φ = 1.0
+            Ψ = 20.0
+            g = 0.01
+            xs, nonzeros, profits = AllocationOpt.optimize(Ω, ψ, σ, K, Φ, Ψ, g, config)
+            @test xs == [2.5; 2.5;;]
+            @test nonzeros == [2]
+            @test isapprox(profits, [0.35; 0.35]; atol=0.1)
         end
     end
 end