
using Convex
using SCS # I think I fixed it
#using ECOS
using Plots
using LinearAlgebra


include("plotregion.jl")

Main.PlotRegion


A1 = [-2.0 1; 
      -1 2; 
       1 0]
b = [2.0; 7; 3]

A = [A1 Matrix{Float64}(I,3,3)] # form the problem with slacks.

PlotRegion.plotregion(A,b)


# Convert the problem into standard form.
cs = [-1 -2 0 0 0]'
AS = A

3×5 Matrix{Float64}:
 -2.0  1.0  1.0  0.0  0.0
 -1.0  2.0  0.0  1.0  0.0
  1.0  0.0  0.0  0.0  1.0


""" Solve the central-path problem for interior point methods. """
function ip_central(c,A,b,tau)
    x = Variable(length(c))
    constraints = Constraint[
        A*x == b
    ]
    p = minimize(c'*x - tau*sum(log(x)), constraints)
    
    solve!(p, SCS.Optimizer; silent=true)
    #solve!(p, ECOS.Optimizer; silent = true)
    return x.value, p
end
ip_central(cs,AS,b,10.0)[1]

5×1 Matrix{Float64}:
 2.1944265108421326
 2.1971437784015126
 4.191709445576553
 4.800138904383412
 0.8055738457513615


ip_central(cs,AS,b,1e-7)[1]

5×1 Matrix{Float64}:
 2.9999995053211297
 4.999999384661257
 2.9999996153139215
 7.2995288733118e-7
 4.927979362203644e-7


taus = vec([10 7.5 5 3.5 2 1 0.75 0.5 0.35 0.20 10.0.^(range(-1,stop=-7,length=10))'])

p = PlotRegion.plotregion(AS,b)
for tau in taus
    x = ip_central(cs, AS, b, tau)[1]
    scatter!([x[1]],[x[2]],label="", marker_z=log10(tau))
end
p


taus = vec([10 7.5 5 3.5 2 1 0.75 0.5 0.35 0.20 10.0.^(range(-1,stop=-7,length=10))'])

p = PlotRegion.plotregion(AS,b)
for tau in taus
    x = ip_central([-1 0 0 0 0.0]', AS, b, tau)[1]
    scatter!([x[1]],[x[2]],label="", marker_z=log10(tau))
end
p


x = ip_central([-1,-2.0,0,0,0], AS, b, 0.0002)
@show x[1]

x[1] = [3.00000036770714; 4.999971118925543; 3.000029604328589; 5.83496815133341e-5; -2.567714595444165e-7;;]

5×1 Matrix{Float64}:
  3.00000036770714
  4.999971118925543
  3.000029604328589
  5.83496815133341e-5
 -2.567714595444165e-7


A1 = [-2.0 1; 
      -1 2; 
       1 0]
b = [2.0; 7; 3]

AS = [A1 Matrix{Float64}(I,3,3)] # form the problem with slacks.
cs = [1 0 0 0 0]'
#cs = [-1 -2 0 0 0]'

taus = vec([10 7.5 5 3.5 2 1 0.75 0.5 0.35 0.20 10.0.^(range(-1,stop=-7,length=10))'])

p = PlotRegion.plotregion(AS,b)
for tau in taus
    x = ip_central(cs, AS, b, tau)[1]
    scatter!([x[1]],[x[2]],label="", color="red")
end
p


# Convert the problem into standard form.
cs = [-1 -2 0 0 0]'
AS = A

tau = 1.0
x0,prob = ip_central(cs,AS,b,tau)
x = copy(x0)
x[1] += 0.5
x[2] += 0.5

lam = vec(prob.constraints[1].dual)
@show lam

# show the region and the starting point
plt = PlotRegion.plotregion(AS,b)
scatter!([x[1]],[x[2]],label="", color="red")

# compute the steps 
n = length(cs)
m = size(AS,1)

s = tau./x
J = [zeros(n,n) AS'         Matrix{Float64}(I,n,n); 
     AS        zeros(m,m)   zeros(m,n);
    Diagonal(vec(s))  zeros(n,m) Diagonal(vec(x))]

mu = dot(x,s)/n
sigma = 0.5

F = [s + AS'*lam - cs; AS*x - b; x.*s]

Fc = [s + AS'*lam .- cs; AS*x .- b; x.*s .- sigma*mu ]

p = J\-F
pc = J\-Fc

plot!([x[1];x[1] + p[1]], [x[2];x[2] + p[2]], label="Affine")
plot!([x[1];x[1] + pc[1]], [x[2];x[2] + pc[2]], label = "Centered")

plt


xf = [x; lam; s];

[xf+p xf+pc]


taus = vec([10 7.5 5 3.5 2 1 0.75 0.5 0.35 0.20 10.0.^(range(-1,stop=-7,length=10))'])

p = PlotRegion.plotregion(AS,b)
for tau in taus
    x = ip_central(cs, AS, b, tau)[1]
    @show tau
    @show x
    scatter!([x[1]],[x[2]],label="", color="red")
end
p


x = ip_central([-1,-2.0,0,0,0], AS, b, 0.0002)
@show x[1]

Show a centered vs. uncentered step¶