Implement Stableswap for 2 assets

Project.toml

@@ -5,6 +5,7 @@ version = "0.1.0"
 
 [deps]
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LBFGSB = "5be7bae1-8223-5378-bac3-9e7378a2f6e6"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"

examples/arbitrage-stableswap.jl

@@ -0,0 +1,55 @@
+#=
+# Arbitrage
+This example illustrates how to use CFMMRouter.jl to solve the multi-market
+arbitrage problem
+=#
+using Pkg
+Pkg.activate("..")
+Pkg.instantiate()
+using CFMMRouter
+using LinearAlgebra
+
+
+## Create three pools of the same tokens, no fees (γ=1)
+# equal_pool = ProductTwoCoin([1e6, 1e6], 1, [1, 2])
+equal_pool = StableswapTwoCoin([1e6, 1e6], 1, [1, 2], 10.)
+unequal_small_pool = ProductTwoCoin([1e3, 2e3], 1, [1, 2])
+weighted_pool = GeometricMeanTwoCoin([1e4, 2e4], [.4, .6], 1, [1, 2])
+
+## Build a routing problem with price vector = [1.0, 1.0]
+prices = ones(2)
+router = Router(
+    LinearNonnegative(prices),
+    [equal_pool, unequal_small_pool, weighted_pool],
+    2,
+)
+
+## Optimize!
+route!(router)
+
+## Print results
+Ψ = round.(Int, netflows(router))
+println("Net trade: $Ψ")
+println("Profit: $(dot(prices, Ψ))")
+
+#=
+We can also see the list of individual trades with each CFMM:
+=#
+## Print individual trades
+for (i, (Δ, Λ)) in enumerate(zip(router.Δs, router.Λs))
+    tokens = router.cfmms[i].Ai
+    println("CFMM $i:")
+    println("\tTendered basket:")
+    for (ind, δ) in enumerate(Δ)
+        if δ > eps()
+            print("\t  $(tokens[ind]): $(round(Int, δ)), ")
+        end
+    end
+    println("\n\tRecieved basket:")
+    for (ind, λ) in enumerate(Λ)
+        if λ > eps()
+            print("\t  $(tokens[ind]): $(round(Int, λ)), ")
+        end
+    end
+    print("\n")
+end

examples/arbitrage.jl

@@ -6,6 +6,7 @@ arbitrage problem
 using CFMMRouter
 using LinearAlgebra
 
+
 ## Create three pools of the same tokens, no fees (γ=1)
 equal_pool = ProductTwoCoin([1e6, 1e6], 1, [1, 2])
 unequal_small_pool = ProductTwoCoin([1e3, 2e3], 1, [1, 2])

examples/convex-opt-test-jump.jl

@@ -0,0 +1,82 @@
+using JuMP, LinearAlgebra, Ipopt
+# using ForwardDiff
+
+# R = [100.0, 80.0, 50.0]
+# π = [0.9, 1.2, 1.3]
+
+R = [100.0, 90.0]
+π = [1., 0.9]
+
+# A = 10000.0
+A = 10.0
+
+γ = 0.996
+# γ = 1.
+
+ϕ(α, R...) = sum(R) - α / prod(R)
+# ϕ(α, R...) = solve_D(R, A)
+# ϕ(α, R...) = prod(R)
+# ϕ(α, R...) = sum(R)
+
+
+# function ∇ϕ(R, α)
+#     return ForwardDiff.gradient(x -> ϕ(x, α), R)
+# end
+
+# Solve D for given reserves. Taken from the Python mockup:
+# https://github.com/curvefi/curve-contract/blob/b0bbf77f8f93c9c5f4e415bce9cd71f0cdee960e/tests/simulation.py#L30-L52
+function solve_D(R, A, tol=1e-12)
+    Dprev = 0
+    n = length(R)
+    S = sum(R)
+    D = S
+    Ann = A * n^n
+    while abs(D - Dprev) > tol
+        D_P = D
+        for x in R
+            D_P = D_P * D / (n * x)
+        end
+        Dprev = D
+        D = (Ann * S + D_P * n) * D / ((Ann - 1) * D + (n + 1) * D_P)
+    end
+    return D
+end
+
+# Solves the maximum arbitrage problem for the generic Stableswap case.
+function find_arb(R, π, A, γ)
+    n = length(R)
+    D = solve_D(R, A)
+    α = D^(n + 1) / (A * n^(2 * n))
+
+    # @assert A*n^n *sum(R) + D - A*D*n^n - D^(n+1)/n^n/prod(R) <= 1e-10
+
+    model = Model(Ipopt.Optimizer)
+    register(model, :ϕ, n + 1, ϕ; autodiff=true)
+
+    @variable(model, Δ[1:n] >= 0)
+    @variable(model, Λ[1:n] >= 0)
+
+    ex1 = @NLexpression(model, [i = 1:n], R[i] + Δ[i] - Λ[i] / γ)
+    # ex1 = @NLexpression(model, [i = 1:n], R[i] + γ * Δ[i] - Λ[i])
+    # @NLconstraint(model, ϕ(R + Δ - Λ / γ, α) >= ϕ(R, α))
+    @NLconstraint(model, ϕ(α, ex1...) >= ϕ(α, R...))
+    @NLconstraint(model, [i=1:n], Δ[i] * Λ[i] == 0)
+    @NLconstraint(model, [i=1:n], ex1[i] >= 0)
+
+    @objective(model, Max, sum(π[i] * (Λ[i] - Δ[i]) for i = 1:n))
+    optimize!(model)
+
+    print(model)
+    println("Objective value: ", objective_value(model))
+    println("Δ = ", value.(Δ))
+    println("Λ = ", value.(Λ))
+
+    R_new = R + value.(Δ) - value.(Λ) / γ
+
+    println("R                = ", R)
+    println("R + Δ - Λ / γ    = ", R_new)
+    println("ϕ(R)             = ", ϕ(α, R...))
+    println("ϕ(R + Δ - Λ / γ) = ", ϕ(α, R_new...))
+end
+
+find_arb(R, π, A, γ)

examples/liquidate.jl

@@ -21,10 +21,10 @@ router = Router(
     cfmms,
     maximum([maximum(cfmm.Ai) for cfmm in cfmms]),
 )
-    
+
 ## Optimize!
 route!(router)
-    
+
 ## Print results
 Ψ = round.(Int, netflows(router))
 println("Input Basket: $(round.(Int, Δin))")