# ] add IntervalArithmetic#fastpowers
# ] add IntervalConstraintProgramming#separators_from_contractors  # ?
# ] add ModelingToolkit#master

using IntervalArithmetic, IntervalOptimisation
using IntervalConstraintProgramming, IntervalRootFinding
using StaticArrays, ModelingToolkit, LinearAlgebra

include("constraints.jl"); include("unify_solutions.jl")

G(X) = 1 + sum(X.^2) / 4000 - prod( cos(X[i] / √i) for i in 1:length(X) )   # Griewank
h(x) = G(x - 1.5) * G(x-1) * G(x+1) * G(x + 1.5)

using Plots; gr()

p = plot(h, -50:0.001:50, leg=false)

@format midpoint 3

using ForwardDiff

@time rts = IntervalRootFinding.roots(x->ForwardDiff.derivative(h, x), -50..50);  rts

mids = mid.(interval.(rts))
second_deriv(f, x) = ForwardDiff.derivative(x->ForwardDiff.derivative(f, x), x)

scatter!(mids, h.(mids), marker_z = second_deriv.(h, mids))

@time globalmin, minimizers = IntervalOptimisation.minimize(h, -50..50, tol=1e-5);

globalmin

sort(unify_solutions(minimizers))

minimizers2 = mid.(minimizers);

scatter!(minimizers2, h.(minimizers2), m=:x, ms=3)

xlims!(-3, 3)
ylims!(-1, 1)
plot!()

V(r) = 4 * (r^-12 - r^-6)

using Plots
plot(0.01:0.01:3, V, ylim=(-2, 5))

all(isunique, rts)

second_deriv(G, x) = ForwardDiff.derivative(x->ForwardDiff.derivative(G, x), x)

p = plot(h, -50:0.001:50)
scatter!(minimizers, [m for i in 1:length(minimizers)])

mids = mid.(interval.(rts))
scatter!(mids, h.(mids), ms=2, marker_z=second_deriv.(G, mids))

G2(X) = 1 + sum(X.^2) / 4000 - prod( cos(X[i] / √i) for i in 1:length(X) )
G2(x, y) = G2(SVector(x, y))

rts = roots(∇(G), IntervalBox(-10..10, 2), Krawczyk)

@time rts = roots(∇(G), IntervalBox(-600..600, 2), Krawczyk, 1e-8)

stability(f, X) = count(eigvals(ForwardDiff.hessian(f, X)) .> 0)  # can make rigorous

r = -10:0.1:10
surface(r, r, G, alpha=0.8)

mids = mid.(interval.(rts))

scatter!(first.(mids), last.(mids), G2.(mids), marker_z=stability.(G2, mids), alpha=0.7)

using IntervalArithmetic

@format standard
X = -1..2   

X^2 

X = -1..2;  Y = 3..4;  (X, Y)

X + Y, X - Y, X * Y, X / Y

X ∩ Y, X ∪ Y

X - X

using Interact, Plots
gr()

function mince(x::Interval, n) 
    nodes = range(x.lo, x.hi, length=n+1)
    
    return [Interval(nodes[i], nodes[i+1]) for i in 1:length(nodes) - 1]
end

function bound_function(f, X, N, ylims=(-2, 2))  
    @manipulate for N in slider(1:200, value=1)
        intervals = mince(X, N)    
        images = f.(intervals)

        boxes = [IntervalBox(intervals[i], images[i]) for i in 1:length(intervals)]
        contains_zero = [0 ∈ images[i] for i in 1:length(images)]

        if length(boxes[contains_zero]) > 0
            plot(boxes[contains_zero], ylim=ylims)   # uses "plot recipe" for IntervalBoxes
        end

        if length(boxes[.!(contains_zero)]) > 0
            plot!(boxes[.!(contains_zero)], ylim=ylims)  
        end
        
        plot!(X.lo:0.0001:X.hi, f, lw=3, leg=false)
    end
end

bound_function(x -> x^2 - 2x, -1..3, 50, (-10, 15))

bound_function(x -> cos(x) + 0.5*sin(2*x), -5..5,  50, (-2,2))

bound_function(x -> (1/50) * log(abs(3 *  (1 - x) + 1)) + x^2 + 1, 0.8..2.0, 100, (1, 4))

bound_function(x -> sin(1 / x), 0.01..0.2, (-1, 1))

using IntervalArithmetic
@format standard

X = 3..4;  (X, typeof(X))  # make an interval

X^2

f(x) = x^2 - 2;  f(X)   

0 ∈ f(X) 

0 ∈ f(2..∞)

XX = (3..4) × (5..6)   # \times <TAB>

XX = IntervalBox(3..4, 5..6)

f(X) = ( (x, y) = X; IntervalBox(x^2 + y^2 - 1, x - y) )
f(XX)

XX = IntervalBox(3..4, 10)

using ForwardDiff, StaticArrays

ForwardDiff.derivative(X -> X^2 - 2X, -2..2)   

ff(XX) = ( (X, Y) = XX; SVector(X^2 + Y^2, X - Y) )
 
J = ForwardDiff.jacobian(ff, SVector(big(3..4), 5..6) )

J^-1  

m = 10.0   # running upper bound on global minimum

f(x) = x^2 - 2

f(5..10) > m  # throw away box 5..10

using Interact, Plots, IntervalArithmetic

function calculate_branch_bound(f, X, N)
    interval_lists = [[X]]
    m = mid(X)
    
    upper_bound = X.hi
    upper_bounds = [upper_bound]
    
    working_list = [X]
    
    for i in 1:N
        X = popfirst!(working_list)
        
        upper_bound = min(upper_bound, f(interval(mid(X))).hi)
        
        if f(X).lo <= upper_bound
            X1, X2 = bisect(X)
            push!(working_list, X1, X2)
        end
        
        push!(interval_lists, copy(working_list))
        push!(upper_bounds, upper_bound)
    end
            
    
    return interval_lists, upper_bounds
end

function visualize_branch_bound(f, X, N)
    interval_lists, upper_bounds = calculate_branch_bound(f, X, N)
    
    @manipulate for i in slider(1:N, value=1)
        Xs = interval_lists[i]
        plot(IntervalBox.(Xs[2:end], f.(Xs[2:end])))
        plot!(IntervalBox(Xs[1], f(Xs[1])), ylim=(-10, 5))
        hline!([upper_bounds[i+1]], ls=:dash, lw=2)
        plot!(X.lo:0.001:X.hi, f, lw=3, ylim=(-12, 12), leg=false, alpha=0.8)
        scatter!([ (mid(Xs[1]), f(mid(Xs[1])) ) ])

    end
end

gr()
visualize_branch_bound(x -> x * sin(x), -10..8, 100)

function simple_minimize(f, X, ϵ=1e-3)
    L = [X]; R = typeof(X)[ ]; m = ∞   

    while !isempty(L)
        X = popfirst!(L)  # first element: queue

        f(X).lo > m && continue 

        diam(X) < ϵ && ( push!(R, X); continue )

        current_max = sup.(f(Interval.(mid(X))))   # evaluate at midpoint
        current_max < m && (m = current_max)

        push!(L, bisect(X)...)
    end

    minimizers = filter(X -> f(X).lo <= m, R)  # remove infeasible boxes
    lower_bound = reduce(∪, f.(minimizers)).lo

    return interval(lower_bound, m), minimizers
end

@format 4

@time m, R = simple_minimize(x->(x^2 - 2)^2 - 2, -10..10)
m, unify_solutions(R)

G(X) = 1 + sum(abs2, X) / 4000 - prod(cos(X[i] / √i) for i in 1:length(X))
G(x, y) = G([x, y])

using Plots
plotlyjs()

r = -10:0.1:10
surface(r, r, G, alpha=0.8)

minimize(G, -∞..∞)

minimize(G, (-∞..∞) × (-∞..∞))

@time minimize(G, IntervalBox(-∞..∞, 10))

@time [minimize(G, IntervalBox(-1..1, n)) for n in 1:100];

using Statistics, IntervalOptimisation, DataStructures

ns = round.(Int, exp10.(0.5:0.1:2));

times = [mean(@elapsed minimize(G, IntervalBox(-Inf..Inf, n), structure=PriorityQueue, tol=1e-3) for i in 1:1) for n in ns]

using Plots

using LaTeXStrings, Interact

gr()
times_plot = plot(ns, times, m=:o, scale=:log10, ylims=(1e-3, 10), ms=3, label="data", lw=3)
plot!(x->2e-4*x^2, ls=:dash, label=L"\sim n^2", leg=:topleft, bottom_margin = 10*Interact.mm, lw=3)
xlabel!("dimension")
ylabel!("CPU time (s)")
xlims!(2.5, 150)
ylims!(0.001, 100)

times_plot

savefig("griewank_times.pdf")

using IntervalRootFinding

@time rts = roots(∇(G), IntervalBox(-10..10, 2), Krawczyk)

@time rts = roots(∇(G), IntervalBox(-600..600, 2), Krawczyk, 1e-8)

using ForwardDiff
stability(f, X) = count(eigvals(ForwardDiff.hessian(f, X)) .> 0)  # can make rigorous

r = -10:0.1:10
surface(r, r, G, alpha=0.8)

mids = mid.(interval.(rts))
scatter!(first.(mids), last.(mids), G.(mids), marker_z=stability.(G, mids), alpha=0.7)

@time m, R = simple_minimize(x -> (x-2) * (x+2), -10..10, 1e-7);  length(R)

using ForwardDiff, LinearAlgebra

ForwardDiff.gradient(f, X::IntervalBox) = ForwardDiff.gradient(f, X.v)

function mean_value_form(f, X)
    m = IntervalBox(mid.(X))
    
    return f(m) + ForwardDiff.gradient(f, X) ⋅ (X - m)
end

mean_value_form(f) = X -> mean_value_form(f, X)

f(x, y) = x^2 - 2x*y + y^2 + 2x
f(X) = f(X[1], X[2])

ϵ = 1e-2
@time m, R1 = minimize(f, IntervalBox(-10..10, 2), tol=ϵ); @show length(R1)
@time m, R2 = minimize(mean_value_form(f), IntervalBox(-10..10, 2), tol=ϵ); @show length(R2);

using TreeView

@tree x^2 + y^2

using IntervalArithmetic 
   
x = y = -∞..∞  # intervals

a = x^2    
b = y^2

c = a + b 

c = c ∩ (-∞..1)   # c = x^2 + y^2 <= 1

c

# c = a + b  ∴  a = c - b

a,  (c - b)  # currently inconsistent

a = a ∩ (c - b)  # make consistent

# a = x^2
x = √a ∪ (-√a)  

using IntervalConstraintProgramming

# macro syntax:

C = @contractor x^2 + y^2

 using ModelingToolkit
 
@variables x y      # master branch

f(x, y ) = x^2 + y^2

f(x, y)

C = Contractor(f(x, y))

using MacroTools

MacroTools.striplines(C.forward.code)

MacroTools.striplines(C.backward.code)  # using reverse (inverse) functions

using Interact
gr()

using IntervalConstraintProgramming 

C = @constraint 2 <= x^3 + y^2 <= 5

X = (-5..5) × (-5..5)  # \times<TAB>

@manipulate for n in slider(-0.0:0.5:8, value=-0.0)
    ϵ = 2.0^(-n)
    paving = pave(C, X, ϵ)
    
    plot(paving.inner, aspect_ratio=:equal, leg=false)
    plot!(paving.boundary, title="area: $(Vol(paving).bounds)")
end

using IntervalArithmetic, IntervalConstraintProgramming, ModelingToolkit
include("constraints.jl")

function unconstrained_minimize(C::Contractor, X, ϵ=1e-3)
    
    L = [X]; R = typeof(X)[ ]; m = ∞   

    while !isempty(L)
        X = popfirst!(L)  # first element: queue

        C(X).lo > m && continue  # C(X) is the usual forward function

        diam(X) < ϵ && ( push!(R, X); continue )

        current_max = sup.(C(IntervalBox(mid(X))))
        current_max < m && (m = current_max)
        
        X = C(interval(-∞, m), X)  # C(constraint, X) performs constraint propagation

        !isempty(X) && push!(L, bisect(X)...)
            
    end

    minimizers = filter(X -> C(X).lo <= m, R)
    lower_bound = reduce(∪, C.(minimizers)).lo

    return interval(lower_bound, m), minimizers

end

f(x, y) = (1.5 - x * (1 - y))^2 + (2.25 - x * (1 - y^2))^2 + (2.625 - x * (1 - y^3))^2
f(X) = f(X[1], X[2])

C = Contractor(f)

using IntervalOptimisation, DataStructures

@time m, R = minimize(f, IntervalBox(-1e2..1e2, 2), structure=PriorityQueue);

# Particularly favourable case:

@time unconstrained_minimize(C, IntervalBox(-1e2..1e2, 2), 1e-8)

using IntervalArithmetic, IntervalConstraintProgramming, ModelingToolkit
include("constraints.jl")

function constrained_minimize(C::Contractor, X, constraints::Vector{<:Constraint}, ϵ=1e-3)
    
    L = [X]; R = typeof(X)[ ]; m = ∞   

    while !isempty(L)
        X = popfirst!(L)  # first element: queue

        C(X).lo > m && continue  # C(X) is the usual forward function

        diam(X) < ϵ && ( push!(R, X); continue )
        
        for constraint in constraints
            X = constraint(X)
        end
        
        isempty(X) && continue

        image = C(X)   # can no longer assume we can find feasible point
        image.hi < m && (m = image.hi)
        
        X = C(interval(-∞, m), X)

        isempty(X) && continue
        
        push!(L, bisect(X)...)
    end

    minimizers = filter(X -> C(X).lo <= m, R)
    lower_bound = reduce(∪, C.(minimizers)).lo

    return interval(lower_bound, m), minimizers
end

X = IntervalBox(-1e5..1e5, 2)

@variables x y
@time m, R = constrained_minimize(Contractor((x,y) -> x + y), X, 
                            [Constraint(x^2 + y^2 < 1)], 1e-7);

m

unify_solutions(R)

# Jaulin book, p. 125 / Hansen book section 17.8

c1(p1, p2) = 12p1^2 - 6.3p1^4 + p1^6 + 6p1*p2 + 6p2^2
c1(p) = c1(p[1], p[2])

X = IntervalBox(-2..4, 2)

@variables p1 p2

constraints = Constraint.( [
                1 - 16p1^2 - 25p2^2 < 0,
                13p1^3 - 145p1 + 85p2 - 400 < 0,
                p1*p2 - 4 < 0
                ]  );

C = Contractor(c1)

@time m, R = constrained_minimize(Contractor(c1), X, constraints, 1e-6);
m

include("unify_solutions.jl")
unify_solutions(R)

separators = Separator.(constraints)

paving = pave(reduce(∩, separators), X, 0.01)

gr()

plot(paving.inner, leg=false, aspect_ratio=1, lw=0, alpha=0.5)
plot!(paving.boundary)

m = mid(R[1])
scatter!( [(m[1], m[2]), (-m[1], -m[2])])

using IntervalArithmetic, IntervalRootFinding
using Interact, Plots

plotlyjs()

@manipulate for A in 1:20
    f(x, y) = 2A + x^2 - A*cos(2π*x) + y^2 - A*cos(2π*y)
    f(X) = f(X[1], X[2])

    L = 5.0; X = IntervalBox( -L..(L+1), 2)

    @time rts = roots(∇(f), X, Newton, 1e-10);

    surface(-5:0.1:6, -5:0.1:6, f, alpha=0.8)
    mids = mid.(interval.(rts))
    scatter!(first.(mids), last.(mids), f.(mids), alpha=0.5, ms=1, c=:blue)
end

stability(X) = count(eigvals(ForwardDiff.hessian(f, X)) .> 0)  # can make rigorous

surface(-5:0.1:6, -5:0.1:6, f, alpha=0.8)

mids = mid.(interval.(rts))
stabilities = stability.(mids)
mins = mids[stabilities .== 2]; maxs = mids[stabilities .== 0]; stats = mids[stabilities .== 1]
scatter!(first.(mins), last.(mins), f.(mins), alpha=0.5, ms=2, c=:blue)
scatter!(first.(maxs), last.(maxs), f.(maxs), alpha=0.5, ms=2, c=:red)
scatter!(first.(stats), last.(stats), f.(stats), alpha=0.5, ms=2, c=:green)

struct SimpleInterval{T <: AbstractFloat}
    lo::T
    hi::T
end

Base.:+(x::SimpleInterval, y::SimpleInterval) = @round(x.lo + y.lo, x.hi + y.hi)

Interval( +(x.lo, y.lo, RoundDown), +(x.hi, y.hi, RoundUp) )

a = 3.1..3.2

@macroexpand IntervalArithmetic.@round(a.lo + a.lo, a.hi + a.hi)



bound_function(x->sin(1/x), 0.01..0.2, (-1,1))

struct SimpleInterval   
    lo::Float64
    hi::Float64
end

SimpleInterval(a) = SimpleInterval(a, a)

a = SimpleInterval(0.1)
b = SimpleInterval(0.2)  

import Base: +, show

show(io::IO, x::SimpleInterval) = 
    print(io, "[$(x.lo), $(x.hi)]")

+(x::SimpleInterval, y::SimpleInterval) = 
    SimpleInterval(x.lo + y.lo, x.hi + y.hi) 

a = SimpleInterval(0.1, 0.1) 
b = SimpleInterval(0.2, 0.2)
a + b

@code_native a + b

using ErrorfreeArithmetic

ErrorfreeArithmetic.two_sum(0.1, 0.3)

+(0.1,   0.3, RoundDown), +(0.1, 0.3, RoundUp)

function mince(x::Interval, n) 
    nodes = linspace(x.lo, x.hi, n+1)
    
    return [Interval(nodes[i], nodes[i+1]) for i in 1:length(nodes) - 1]
end

f(x) = cos(x) + 0.5*sin(2*x)
 
function image(f, X)
    II = f(X)
   
    if II.lo == -Inf
        II = -100..(II.hi)
    end
    
    return II
end

intervals = mince(-5..5, 10)
    
images = image.(f, intervals)

using Interact, Plots, IntervalArithmetic; gr()

function bound_function(f, X, N ,ylims=(-2,2))  
    @manipulate for N in slider(1:200, value=1)
        intervals = mince(X, N)    
        images = image.(f, intervals)

        boxes = [IntervalBox(intervals[i], images[i]) for i in 1:length(intervals)]
        contains_zero = [0 ∈ images[i] for i in 1:length(images)]

        if length(boxes[contains_zero]) > 0
            plot(boxes[contains_zero], ylim=ylims)   # uses "plot recipe" for IntervalBoxes
        end

        if length(boxes[.!(contains_zero)]) > 0
            plot!(boxes[.!(contains_zero)], ylim=ylims)  
        end
        plot!(f)
    end
end

f(x) = cos(x) + 0.5*sin(2*x)

bound_function(f, -5..5,  50, (-2,2))


f(x) = (1/50)*log(abs(3*(1-x) + 1)) + x^2 + 1

bound_function(f, 0.8..2.0, 100, (1,4))


bound_function(x->sin(1/x), 0.01..0.2, (-1,1))

G(X) = 1 + sum(X.^2) / 4000 - prod( cos(X[i] / √i) for i in 1:length(X) )

using IntervalOptimisation

@time m, R = minimize(x->(x^2 - 2)^2 - 2, -10..10)   # gives lower bound for global minimum too


carrom(x, y) = - (exp(2*abs(1 - sqrt(x^2 + y^2) / pi)) * cos(x)^2 * cos(y)^2) / 30

carrom(X) = carrom(X[1], X[2])

using Plots
plotly()
surface(-5:0.1:5, -5:0.1:5, carrom)

using IntervalArithmetic, IntervalOptimisation, Plots

@time globalmin, minimisers = simple_minimize(carrom, 
                                       (-10..10) × (-10..10), 
                                       2e-3)

@time globalmin, minimisers = minimize(carrom, 
                                       (-10..10) × (-10..10), 
                                       1e-3)

gr()
contour(-10:0.05:10, -10:0.05:10, carrom, aspect_ratio=1)
plot!(minimisers)

XX = IntervalBox(3..4, 10)

import ForwardDiff
using StaticArrays, LinearAlgebra

ForwardDiff.gradient(f, X::IntervalBox) = ForwardDiff.gradient(f, X.v)

function mean_value_form(f, X)
    m = IntervalBox(mid.(X))
    
    return f(m) + ForwardDiff.gradient(f, X) ⋅ (X - m)
end

mean_value_form(f) = X -> mean_value_form(f, X)

using StaticArrays, IntervalArithmetic
 
# f(X) = X[1] * (X[1] - X[2]) + X[2] * (-X[1] + X[2]) + 1  # (X[1] - X[2])^2  # infinite number of solutions

# f(X) = (X[1] - X[2])^2 + (X[1] + X[2])^2

f(X) = X[1]^2 - 2X[1]*X[2] + X[2]^2 + X[1]^2 + 2X[1]*X[2] + X[2]^2


f(X) = (X[1]^2 - 2)^2 + (X[2]^2 - 2)^2

@time m, R = minimize(f, IntervalBox(-3..3, 2), structure=PriorityQueue, tol=1e-3);

length(R)

R

g(X) = mean_value_form(f)(X)  ∩ f(X)

@time m, R = minimize(g, IntervalBox(-10..10, 2), structure=PriorityQueue, tol=1e-1);

R

R

@time m, R = minimize(X->mean_value_form(f)(X)  ∩ f(X), IntervalBox(-10..10, 2), structure=PriorityQueue, tol=1e-3);

MM = @SMatrix [1 -1; -1 1]
f(X) = dot(X, MM*X)

@time m, R = minimize(f, IntervalBox(-10..10, 2), structure=PriorityQueue, tol=1e-2);

length(R)

f(x) = 1 - x[1]^4 + x[1]^5

m, R = minimize(f, 0..1)

length(R)

m, R = minimize(X->mean_value_form(f)(X) ∩ f(X), IntervalBox(0..1), tol=1e-5)

length(R)