text:
using NLSolvers
using LinearAlgebra
using SparseDiffTools
using SparseArrays
using IterativeSolvers
using ForwardDiff
using Test

@testset "optimization interface" begin
# TODO
# Make a more efficient MeritObjective that returns something that acts as the actual thing if requested (mostly for debug)
# but can also be efficiently used to get cauchy and newton
#
# Make DOGLEG work also with BFGS (why is convergence so slow?)
# # Look into what caches are created
# # Why does SR1 give inf above?
#
# Mutated first r last (easy to make fallback for nonmutating)
# NelderMead
# time limit not enforced in @show solve(NelderMead)
# no convergence crit either
#
# APSO
# no @show solve
# wrong return type (no ConvergenceInfo)
#
# PureRandom search. wrong return type and move bounds to problem
#
# Does not have a ! method, this should be documented. Maybe add it for consistency?
# If the sampler is empty and there are bounds, draw uniformly there in stead of specifying lb and ub in PureRandomSearch
#
# really need a QNmodel for model vars that creates nothing or don't populate fields of a named tuple for Newton for example
# LineObjective and  LineObjective! should just dispatch on the caceh being nothing or not
#
# ConjgtaeGraduent with HZAW fails because it overwrites Py into P∇fz which seems to alias P∇fz. That alias needs to be checked
# and a CGModelVars type should allocate Py where appropriate - could y be overwitten with Py and then recalcualte y afterwards?
#
#
# ADAM needs @show solve and AdaMax
#
# TODO: LineObjetive doesn't need ! when we have problem in there and mstyle

#### OPTIMIZATION
f = OPT_PROBS["himmelblau"]["array"]["mutating"]
x0 = copy(OPT_PROBS["himmelblau"]["array"]["x0"])
prob = OptimizationProblem(f)
prob_oop = OptimizationProblem(f; inplace=false)
prob_bounds = OptimizationProblem(obj=f, bounds=([-5.0,-9.0],[13.0,4.0]))
prob_bounds_oop = OptimizationProblem(obj=f, bounds=([-5.0,-9.0],[13.0,4.0]); inplace=false)
prob_on_bounds = OptimizationProblem(obj=f, bounds=([3.5,-9.0],[13.0,4.0]))
prob_on_bounds_oop = OptimizationProblem(obj=f, bounds=([3.5,-9.0],[13.0,4.0]); inplace=false)

res = solve(prob, x0, NelderMead(), OptimizationOptions())
@test all(x0 .== [3.0,2.0])
@test res.info.minimum == 0.0

x0 = copy(OPT_PROBS["himmelblau"]["array"]["x0"])
res = solve(prob_oop, x0, NelderMead(), OptimizationOptions())
@test all(x0 .== [3.0,1.0])
@test res.info.minimum == 0.0

x0 = copy(OPT_PROBS["himmelblau"]["array"]["x0"])
res = solve(prob_bounds, x0, APSO(), OptimizationOptions())
@test_broken all(x0 .== [3.0,2.0])
@test res.info.minimum == 0.0

x0 = copy(OPT_PROBS["himmelblau"]["array"]["x0"]).+1
res = solve(prob_on_bounds_oop, x0, ActiveBox(), OptimizationOptions())
@test_broken all(x0 .== [3.0,1.0])
xbounds = [ 3.5, 1.6165968467447174]
@test res.info.minimum == NLSolvers.value(prob_on_bounds, xbounds)

x0 = copy(OPT_PROBS["himmelblau"]["array"]["x0"])
res = solve(prob_bounds, x0, SimulatedAnnealing(), OptimizationOptions())
@test_broken all(x0 .== [3.0,2.0])
@test res.info.minimum < 1e-1

solve(prob, PureRandomSearch(lb=[0.0,0.0], ub=[4.0,4.0]), OptimizationOptions())
solve(prob_oop, PureRandomSearch(lb=[0.0,0.0], ub=[4.0,4.0]), OptimizationOptions())

solve(prob, [0.0,0.0], SimulatedAnnealing(), OptimizationOptions())
solve(prob_oop, [0.0,0.0], SimulatedAnnealing(), OptimizationOptions())

#x0 = copy(OPT_PROBS["himmelblau"]["array"]["x0"])
#@show solve(prob_bounds, x0, SIMAN(), OptimizationOptions())
#@test all(x0 .== [3.0,1.0])

x0 = copy(OPT_PROBS["himmelblau"]["array"]["x0"])
res = solve(prob, x0, LineSearch(), OptimizationOptions())
#@test all(x0 .== [3.0,2.0])
@test res.info.minimum < 1e-12

x0 = copy(OPT_PROBS["himmelblau"]["array"]["x0"])
res = solve(prob, x0, LineSearch(SR1()), OptimizationOptions())
@test res.info.minimum < 1e-12

x0 = copy(OPT_PROBS["himmelblau"]["array"]["x0"])
res = solve(prob, x0, LineSearch(DFP()), OptimizationOptions())
@test res.info.minimum < 1e-12

x0 = copy(OPT_PROBS["himmelblau"]["array"]["x0"])
res = solve(prob, x0, LineSearch(BFGS()), OptimizationOptions())
@test res.info.minimum < 1e-12

x0 = copy(OPT_PROBS["himmelblau"]["array"]["x0"])
res = solve(prob, x0, LineSearch(LBFGS()), OptimizationOptions())
@test res.info.minimum < 1e-12

x0 = copy(OPT_PROBS["himmelblau"]["array"]["x0"])
res = solve(prob, x0, LineSearch(LBFGS(), HZAW()), OptimizationOptions())
@test res.info.minimum < 1e-12

x0 = copy(OPT_PROBS["himmelblau"]["array"]["x0"])
res = solve(prob, x0, LineSearch(LBFGS(), Backtracking()), OptimizationOptions())
@test res.info.minimum < 1e-12

#@test all(x0 .== [3.0,2.0])
x0 = copy(OPT_PROBS["himmelblau"]["array"]["x0"])
res = solve(prob, x0, LineSearch(ConjugateGradient(), HZAW()), OptimizationOptions())
@test res.info.minimum < 1e-12

x0 = copy(OPT_PROBS["himmelblau"]["array"]["x0"])
res = solve(prob_oop, x0, LineSearch(ConjugateGradient(), HZAW()), OptimizationOptions())
@test res.info.minimum < 1e-12

x0 = copy(OPT_PROBS["himmelblau"]["array"]["x0"])
res = solve(prob, x0, LineSearch(ConjugateGradient(update=HS()), HZAW()), OptimizationOptions())
@test res.info.minimum < 1e-12

x0 = copy(OPT_PROBS["himmelblau"]["array"]["x0"])
res = solve(prob_oop, x0, LineSearch(ConjugateGradient(update=HS()), HZAW()), OptimizationOptions())
@test res.info.minimum < 1e-12

# Stalls at [3, 1] with default @show solve
x0 = 1.0.+copy(OPT_PROBS["himmelblau"]["array"]["x0"])
res = solve(prob, x0, LineSearch(Newton()), OptimizationOptions())
@test res.info.minimum < 1e-12

#@test all(x0 .== [3.0,2.0])
x0 = 1.0.+copy(OPT_PROBS["himmelblau"]["array"]["x0"])
res = solve(prob, x0, LineSearch(Newton(), HZAW()), OptimizationOptions())
@test res.info.minimum < 1e-12
#@test all(x0 .== [3.0,2.0])

x0 = copy(OPT_PROBS["himmelblau"]["array"]["x0"])
res = solve(prob, x0, TrustRegion(), OptimizationOptions())
@test_broken all(x0 .== [3.0,2.0])
@test res.info.minimum < 1e-16

x0 = copy(OPT_PROBS["himmelblau"]["array"]["x0"])
res = solve(prob, x0, TrustRegion(DBFGS(), Dogleg()), OptimizationOptions())
@test res.info.minimum < 1e-16

x0 = copy(OPT_PROBS["himmelblau"]["array"]["x0"])
res = solve(prob, x0, TrustRegion(BFGS(), Dogleg()), OptimizationOptions())
@test res.info.minimum < 1e-16

x0 = copy(OPT_PROBS["himmelblau"]["array"]["x0"])
res = solve(prob, x0, TrustRegion(SR1(), NTR()), OptimizationOptions())
@test res.info.minimum < 1e-16

# not PSD
#x0 = copy(OPT_PROBS["himmelblau"]["array"]["x0"])
#res = solve(prob, x0, TrustRegion(Newton(), Dogleg()), OptimizationOptions())
#@test res.info.minimum < 1e-16

x0 = copy(OPT_PROBS["himmelblau"]["array"]["x0"])
res = solve(prob, x0, TrustRegion(BFGS(), Dogleg()), OptimizationOptions())
@test res.info.minimum < 1e-16

x0 = copy(OPT_PROBS["himmelblau"]["array"]["x0"])
res = solve(prob, x0, TrustRegion(DBFGS(), Dogleg()), OptimizationOptions())
@test res.info.minimum < 1e-16

x0 = copy(OPT_PROBS["himmelblau"]["array"]["x0"])
res = solve(prob, x0, Adam(), OptimizationOptions(maxiter=20000))
@test res.info.minimum < 1e-16

## Notice that prob is only used for value so this should be extremely generic! It does need a comparison though.
res = solve(prob, PureRandomSearch(lb=[0.0,0.0], ub=[4.0,4.0]), OptimizationOptions())

f = OPT_PROBS["exponential"]["array"]["mutating"]
x0 = copy(OPT_PROBS["exponential"]["array"]["x0"])
prob = OptimizationProblem(f)
prob_bounds = OptimizationProblem(obj=f, bounds=([-5.0,-9.0],[13.0,4.0]))
prob_on_bounds = OptimizationProblem(obj=f, bounds=([3.5,-9.0],[13.0,4.0]))

x0 = copy(OPT_PROBS["exponential"]["array"]["x0"])
res = solve(prob, x0, TrustRegion(DBFGS(), Dogleg()), OptimizationOptions())
@test res.info.minimum == 2.0

x0 = copy(OPT_PROBS["exponential"]["array"]["x0"])
res = solve(prob, x0, TrustRegion(BFGS(), Dogleg()), OptimizationOptions())
@test res.info.minimum == 2.0

x0 = copy(OPT_PROBS["exponential"]["array"]["x0"])
res = solve(prob, x0, TrustRegion(Newton(), Dogleg()), OptimizationOptions())
@test res.info.minimum == 2.0

x0 = copy(OPT_PROBS["exponential"]["array"]["x0"])
res = solve(prob, x0, TrustRegion(Newton(), NTR()), OptimizationOptions())
@test res.info.minimum == 2.0

x0 = copy(OPT_PROBS["exponential"]["array"]["x0"])
res = solve(prob, x0, TrustRegion(Newton(), NWI()), OptimizationOptions())
@test res.info.minimum == 2.0

x0 = copy(OPT_PROBS["exponential"]["array"]["x0"])
res = solve(prob, x0, TrustRegion(SR1(), NTR()), OptimizationOptions())
@test res.info.minimum == 2.0

x0 = copy(OPT_PROBS["exponential"]["array"]["x0"])
res = solve(prob, x0, TrustRegion(SR1(Inverse()), NTR()), OptimizationOptions())
@test res.info.minimum == 2.0
end

const statictest_s0 = OPT_PROBS["himmelblau"]["staticarray"]["state0"]
const statictest_prob = OptimizationProblem(OPT_PROBS["himmelblau"]["staticarray"]["static"]; inplace=false)
@testset "staticopt" begin
    res = solve(statictest_prob, statictest_s0, LineSearch(Newton()), OptimizationOptions())
    @allocated solve(statictest_prob, statictest_s0, LineSearch(Newton()), OptimizationOptions())
    alloc = @allocated solve(statictest_prob, statictest_s0, LineSearch(Newton()), OptimizationOptions())
    @test alloc == 0

    _res = solve(statictest_prob, statictest_s0, LineSearch(Newton()), OptimizationOptions())
    _alloc = @allocated solve(statictest_prob, statictest_s0, LineSearch(Newton()), OptimizationOptions())
    @test _alloc == 0
    @test norm(_res.info.∇fz, Inf) < 1e-8

    _res = solve(statictest_prob, statictest_s0, LineSearch(Newton(), Backtracking()), OptimizationOptions())
    _alloc = @allocated solve(statictest_prob, statictest_s0, LineSearch(Newton(), Backtracking()), OptimizationOptions())
    @test _alloc == 0
    @test norm(_res.info.∇fz, Inf) < 1e-8
end

@testset "newton" begin
    test_x0 = [2.0, 2.0]
    test_prob = OptimizationProblem(OPT_PROBS["himmelblau"]["array"]["mutating"]; inplace=true)
    res = solve(test_prob, copy(test_x0), LineSearch(Newton()), OptimizationOptions())
    @test norm(res.info.∇fz, Inf) < 1e-8

    test_x0 = [2.0, 2.0]
    test_prob = OptimizationProblem(OPT_PROBS["himmelblau"]["array"]["mutating"]; inplace=false)
    res = solve(test_prob, test_x0, LineSearch(Newton()), OptimizationOptions())
    @test norm(res.info.∇fz, Inf) < 1e-8
end
@testset "Newton linsolve" begin
    test_prob = OptimizationProblem(OPT_PROBS["himmelblau"]["array"]["mutating"]; inplace=true)
    res_lu = solve(test_prob, (copy([2.0,2.0]), [0.0 0.0;0.0 0.0]), LineSearch(Newton(; linsolve=(d, B, g)->ldiv!(d, lu(B), g))), OptimizationOptions())
    @test norm(res_lu.info.∇fz, Inf) < 1e-8
    res_qr = solve(test_prob, (copy([2.0,2.0]), [0.0 0.0;0.0 0.0]), LineSearch(Newton(; linsolve=(d, B, g)->ldiv!(d, qr(B), g))), OptimizationOptions())
    @test norm(res_qr.info.∇fz, Inf) < 1e-8

    test_prob = OptimizationProblem(OPT_PROBS["himmelblau"]["array"]["mutating"]; inplace=false)
    res_qr = solve(test_prob, (copy([2.0,2.0]), [0.0 0.0;0.0 0.0]), LineSearch(Newton(; linsolve=(B, g)->qr(B)\g)), OptimizationOptions())
    @test norm(res_qr.info.∇fz, Inf) < 1e-8
    test_prob = OptimizationProblem(OPT_PROBS["himmelblau"]["array"]["mutating"]; inplace=false)
    res_lu = solve(test_prob, (copy([2.0,2.0]), [0.0 0.0;0.0 0.0]), LineSearch(Newton(; linsolve=(B, g)->lu(B)\g)), OptimizationOptions())
    @test norm(res_lu.info.∇fz, Inf) < 1e-8
end

const static_x0 = OPT_PROBS["fletcher_powell"]["staticarray"]["x0"][1]
const static_prob_qn = OPT_PROBS["fletcher_powell"]["staticarray"]["static"]
@testset "no alloc static" begin

    @testset "no alloc" begin
        @allocated solve(static_prob_qn, static_x0, LineSearch(BFGS(Inverse()), Backtracking()), OptimizationOptions())
        _alloc = @allocated solve(static_prob_qn, static_x0, LineSearch(BFGS(Inverse()), Backtracking()), OptimizationOptions())
        @test _alloc == 0

        solve(static_prob_qn, static_x0, LineSearch(BFGS(Direct()), Backtracking()), OptimizationOptions())
        _alloc = @allocated solve(static_prob_qn, static_x0, LineSearch(BFGS(Direct()), Backtracking()), OptimizationOptions())
        @test _alloc == 0

        solve(static_prob_qn, static_x0, LineSearch(DFP(Inverse()), Backtracking()), OptimizationOptions())
        _alloc = @allocated solve(static_prob_qn, static_x0, LineSearch(DFP(Inverse()), Backtracking()), OptimizationOptions())
        @test _alloc == 0

        solve(static_prob_qn, static_x0, LineSearch(DFP(Direct()), Backtracking()), OptimizationOptions())
        _alloc = @allocated solve(static_prob_qn, static_x0, LineSearch(DFP(Direct()), Backtracking()), OptimizationOptions())
        @test _alloc == 0

        solve(static_prob_qn, static_x0, LineSearch(SR1(Inverse()), Backtracking()), OptimizationOptions())
        _alloc = @allocated solve(static_prob_qn, static_x0, LineSearch(SR1(Inverse()), Backtracking()), OptimizationOptions())
        @test _alloc == 0

        solve(static_prob_qn, static_x0, LineSearch(SR1(Direct()), Backtracking()), OptimizationOptions())
        _alloc = @allocated solve(static_prob_qn, static_x0, LineSearch(SR1(Direct()), Backtracking()), OptimizationOptions())
        @test _alloc == 0
    end
end

Random.seed!(4568532)
solve(static_prob_qn, rand(3), Adam(), OptimizationOptions(maxiter=1000))
solve(static_prob_qn, rand(3), AdaMax(), OptimizationOptions(maxiter=1000))

@testset "bound newton" begin
    f = OPT_PROBS["himmelblau"]["array"]["mutating"]
    x0 = copy(OPT_PROBS["himmelblau"]["array"]["x0"])
    prob = OptimizationProblem(f)
    prob_oop = OptimizationProblem(f; inplace=false)
    prob_bounds = OptimizationProblem(obj=f, bounds=([-5.0,-9.0],[13.0,4.0]))
    prob_bounds_oop = OptimizationProblem(obj=f, bounds=([-5.0,-9.0],[13.0,4.0]); inplace=false)
    prob_on_bounds = OptimizationProblem(obj=f, bounds=([3.5,-9.0],[13.0,4.0]))
    prob_on_bounds_oop = OptimizationProblem(obj=f, bounds=([3.5,-9.0],[13.0,4.0]); inplace=false)

    start = [3.7,2.0]

    res_unc = solve(prob_bounds, copy(start), LineSearch(Newton(), Backtracking()), OptimizationOptions())
    @test res_unc.info.minimizer ≈ [3.0, 2.0]
    res_con = solve(prob_bounds, copy(start), ActiveBox(), OptimizationOptions())
    @test res_con.info.minimizer ≈ [3.0, 2.0]
    res_unc = solve(prob_on_bounds, copy(start), LineSearch(Newton(), Backtracking()), OptimizationOptions())
    @test res_unc.info.minimizer ≈ [3.0, 2.0]
    res_con = solve(prob_on_bounds, copy(start), ActiveBox(), OptimizationOptions())
    @test res_con.info.minimizer ≈ [3.5, 1.6165968467448326]
end

function fourth_f(x)
    fx = x^4 + sin(x)
    return fx
end
function fourth_fg(∇f, x)
    ∇f = 4x^3 + cos(x)

    fx = x^4 + sin(x)
    return fx, ∇f
end

function fourth_fgh(∇f, ∇²fx, x)
    ∇²f = 12x^2 - sin(x)
    ∇f = 4x^3 + cos(x)

    fx = x^4 + sin(x)
    return fx, ∇f, ∇²f
end

const scalar_prob_oop = OptimizationProblem(ScalarObjective(fourth_f, nothing, fourth_fg, fourth_fgh, nothing, nothing, nothing, nothing); inplace=false)
@testset "scalar no-alloc" begin
    _alloc = @allocated solve(scalar_prob_oop, 4.0, LineSearch(SR1(Direct())), OptimizationOptions())
    _alloc = @allocated solve(scalar_prob_oop, 4.0, LineSearch(SR1(Direct())), OptimizationOptions())
    @test _alloc == 0

    _alloc = @allocated solve(scalar_prob_oop, 4.0, LineSearch(BFGS(Direct())), OptimizationOptions())
    _alloc = @allocated solve(scalar_prob_oop, 4.0, LineSearch(BFGS(Direct())), OptimizationOptions())
    @test _alloc == 0

    _alloc = @allocated solve(scalar_prob_oop, 4.0, LineSearch(DFP(Direct())), OptimizationOptions())
    _alloc = @allocated solve(scalar_prob_oop, 4.0, LineSearch(DFP(Direct())), OptimizationOptions())
    @test _alloc == 0

    _alloc = @allocated solve(scalar_prob_oop, 4.0, LineSearch(Newton()), OptimizationOptions())
    _alloc = @allocated solve(scalar_prob_oop, 4.0, LineSearch(Newton()), OptimizationOptions())
    @test _alloc == 0

    _alloc = @allocated solve(scalar_prob_oop, 4.0, LineSearch(Newton()), OptimizationOptions())
    _alloc = @allocated solve(scalar_prob_oop, 4.0, LineSearch(Newton()), OptimizationOptions())
    @test _alloc == 0
end

using DoubleFloats
@testset "Test double floats" begin
	function fdouble(x)
	    fx = (1 - x[1])^2 + 100 * (x[2] - x[1]^2)^2
	    return fx
	end
	function fgdouble(G, x)
	    fx = (1 - x[1])^2 + 100 * (x[2] - x[1]^2)^2

        G1 = -2 * (1 - x[1]) - 400 * (x[2] - x[1]^2) * x[1]
        G2 = 200 * (x[2] - x[1]^2)
        G = [G1,G2]
	    return fx, G
	end
	f_obj = OptimizationProblem(ScalarObjective(fdouble, nothing, fgdouble, nothing, nothing, nothing, nothing, nothing))
	res = res = solve(f_obj, Double64.([1,2]), LineSearch(GradientDescent(Inverse())), OptimizationOptions(;g_abstol=1e-32, maxiter=100000))
    @test res.info.minimum < 1e-45
    res = res = solve(f_obj, Double64.([1,2]), LineSearch(BFGS(Inverse())), OptimizationOptions(;g_abstol=1e-32))
    @test res.info.minimum < 1e-45
	res = res = solve(f_obj, Double64.([1,2]), LineSearch(DFP(Inverse())), OptimizationOptions(;g_abstol=1e-32))
    @test res.info.minimum < 1e-45
	res = res = solve(f_obj, Double64.([1,2]), LineSearch(SR1(Inverse())), OptimizationOptions(;g_abstol=1e-32))
    @test res.info.minimum < 1e-45
end

function myfun(x::T) where T
    fx = T(x^4 + sin(x))
    return fx
end
function myfun(∇f, x::T) where T
    ∇f = T(4*x^3 + cos(x))
    fx = myfun(x)
    fx, ∇f
end
function myfun(∇f, ∇²f, x::T) where T<:Real
    ∇²f = T(12*x^2 - sin(x))
    fx, ∇f = myfun(∇f, x)
    T(fx), ∇f, ∇²f
end
@testset "scalar return types" begin
    for T in (Float16, Float32, Float64, Rational{BigInt}, Double32, Double64)
        if T == Rational{BigInt}
            options = OptimizationOptions()
        else
            options = OptimizationOptions(g_abstol=eps(T), g_reltol=T(0))
        end
        for M in (SR1, BFGS, DFP, Newton)
            if M == Newton
                obj = OptimizationProblem(ScalarObjective(myfun, nothing, myfun, myfun, nothing, nothing, nothing, nothing); inplace=false)
                res = solve(obj, T(3.1), LineSearch(M()), options)
                @test all(isa.([res.info.minimum, res.info.∇fz, res.info.minimizer], T))
            else
                obj = OptimizationProblem(ScalarObjective(myfun, nothing, myfun, myfun, nothing, nothing, nothing, nothing); inplace=false)
                res = solve(obj, T(3.1), LineSearch(M(Direct())), options)
                @test all(isa.([res.info.minimum, res.info.∇fz, res.info.minimizer], T))
                res = solve(obj, T(3.1), LineSearch(M(Inverse())), options)
                @test all(isa.([res.info.minimum, res.info.∇fz, res.info.minimizer], T))
            end
        end
    end
end

using GeometryTypes
@testset "GeometryTypes" begin
    function fu(G, x)
        fx = (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2

        if !(G == nothing)
            G1 = -2.0 * (1.0 - x[1]) - 400.0 * (x[2] - x[1]^2) * x[1]
            G2 = 200.0 * (x[2] - x[1]^2)
            gx = Point(G1,G2)

            return fx, gx
        else
            return fx
        end
    end
    fu(x) =  (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
    f_obj = OptimizationProblem(ScalarObjective(fu, nothing, fu, nothing, nothing, nothing, nothing, nothing); inplace=false)
    res = solve(f_obj, Point(1.3,1.3), LineSearch(GradientDescent(Inverse())), OptimizationOptions())
    @test res.info.minimum < 1e-9
    res = solve(f_obj, Point(1.3,1.3), LineSearch(BFGS(Inverse())), OptimizationOptions())
    @test res.info.minimum < 1e-10
    res = solve(f_obj, Point(1.3,1.3), LineSearch(DFP(Inverse())), OptimizationOptions())
    @test res.info.minimum < 1e-10
    res = solve(f_obj, Point(1.3,1.3), LineSearch(SR1(Inverse())), OptimizationOptions())
    @test res.info.minimum < 1e-10

    res = solve(f_obj, Point(1.3,1.3), LineSearch(GradientDescent(Direct())), OptimizationOptions())
    @test res.info.minimum < 1e-9
    res = solve(f_obj, Point(1.3,1.3), LineSearch(BFGS(Direct())), OptimizationOptions())
    @test res.info.minimum < 1e-10
    res = solve(f_obj, Point(1.3,1.3), LineSearch(DFP(Direct())), OptimizationOptions())
    @test res.info.minimum < 1e-10
    res = solve(f_obj, Point(1.3,1.3), LineSearch(SR1(Direct())), OptimizationOptions())
    @test res.info.minimum < 1e-10
end

@testset "quadratics" begin
    A = rand(2,2)
    A = abs.(A)
    A = Symmetric(A*A')
    x = rand(2)
    b = A*x

    quadf(x) = -dot(b, x) + dot(x, A*x)/2
    function quadfg(G, x)
        G.=A*x-b
        quadf(x), G
    end
    function quadfgh(G, H, x)
        H .= A
        f, G = quadfg(G, x)
        f, G, H
    end

    quadprob = OptimizationProblem(ScalarObjective(quadf, nothing, quadfg, quadfgh, nothing, nothing, nothing, nothing); inplace=true)

    for approx in (GradientDescent(), BFGS(Inverse()), BFGS(Direct()), DBFGS(), SR1(Inverse()), SR1(Direct()), DFP(), Newton(), BB(), LBFGS()) # CBB
        lsres =  solve(quadprob, zeros(2), LineSearch(approx, Backtracking()), OptimizationOptions(maxiter=20000))
        println(rpad(summary(approx), 40), "   ||   $(rpad(lsres.info.iter, 5))   ||   $(lsres.info.∇fz)")
    end
end
@testset "batched" begin

end

@testset "MArray" begin
    f = OPT_PROBS["himmelblau"]["array"]["mutating"]
    x0 = copy(OPT_PROBS["himmelblau"]["array"]["x0"])
    prob = OptimizationProblem(f)

    x0m = @MVector [-1.0, 0.0, 0.0]
    x0 = [-1.0, 0.0, 0.0]
    @time res = solve(prob, copy(x0m), ConjugateGradient(update=VPRP()), OptimizationOptions());
    @time res = solve(prob, copy(x0), ConjugateGradient(update=VPRP()), OptimizationOptions());
    # workaround for https://github.com/JuliaArrays/StaticArrays.jl/issues/828
    @time res = solve(prob, (copy(x0m), MArray(I+x0m*x0m')), LineSearch(BFGS()), OptimizationOptions());
    @time res = solve(prob, copy(x0), LineSearch(BFGS()), OptimizationOptions());
    @time res = solve(prob, (copy(x0m), MArray(I+x0m*x0m')), LineSearch(SR1()), OptimizationOptions());
    @time res = solve(prob, copy(x0), LineSearch(SR1()), OptimizationOptions());
#    @time res = solve(prob, (copy(x0m), MArray(I+x0m*x0m')), TrustRegion(DBFGS(), Dogleg()), OptimizationOptions());
    @time res = solve(prob, copy(x0), TrustRegion(DBFGS(), Dogleg()), OptimizationOptions());
end

# using NLSolvers, StaticArrays, Test
# @testset "mixed optimization problems" begin
# function theta(x)
#    if x[1] > 0
#        return atan(x[2] / x[1]) / (2.0 * pi)
#    else
#        return (pi + atan(x[2] / x[1])) / (2.0 * pi)
#    end
# end
# f(x) = 100.0 * ((x[3] - 10.0 * theta(x))^2 + (sqrt(x[1]^2 + x[2]^2) - 1.0)^2) + x[3]^2

# function f∇f!(x, ∇f)
#     if !(∇f==nothing)
#         if ( x[1]^2 + x[2]^2 == 0 )
#             dtdx1 = 0;
#             dtdx2 = 0;
#         else
#             dtdx1 = - x[2] / ( 2 * pi * ( x[1]^2 + x[2]^2 ) );
#             dtdx2 =   x[1] / ( 2 * pi * ( x[1]^2 + x[2]^2 ) );
#         end
#         ∇f[1] = -2000.0*(x[3]-10.0*theta(x))*dtdx1 +
#             200.0*(sqrt(x[1]^2+x[2]^2)-1)*x[1]/sqrt( x[1]^2+x[2]^2 );
#         ∇f[2] = -2000.0*(x[3]-10.0*theta(x))*dtdx2 +
#             200.0*(sqrt(x[1]^2+x[2]^2)-1)*x[2]/sqrt( x[1]^2+x[2]^2 );
#         ∇f[3] =  200.0*(x[3]-10.0*theta(x)) + 2.0*x[3];
#     end

#     fx = f(x)
#     objective_return(fx, ∇f)
# end

# function f∇f(x, ∇f)
#     if !(∇f == nothing)
#         ∇f = similar(x)
#     end
#     fx, ∇f = f∇f!(gx, x)
#     objective_return(fx, ∇f)
# end
# function f∇fs(x, ∇f)
#     if !(∇f == nothing)
#         if ( x[1]^2 + x[2]^2 == 0 )
#             dtdx1 = 0;
#             dtdx2 = 0;
#         else
#             dtdx1 = - x[2] / ( 2 * pi * ( x[1]^2 + x[2]^2 ) )
#             dtdx2 =   x[1] / ( 2 * pi * ( x[1]^2 + x[2]^2 ) )
#         end

#         s1 = -2000.0*(x[3]-10.0*theta(x))*dtdx1 +
#             200.0*(sqrt(x[1]^2+x[2]^2)-1)*x[1]/sqrt( x[1]^2+x[2]^2 )
#         s2 = -2000.0*(x[3]-10.0*theta(x))*dtdx2 +
#             200.0*(sqrt(x[1]^2+x[2]^2)-1)*x[2]/sqrt( x[1]^2+x[2]^2 )
#         s3 = 200.0*(x[3]-10.0*theta(x)) + 2.0*x[3]
#         ∇f = @SVector [s1, s2, s3]
#     end
#     objective_return(f(x), ∇f)
# end
# obj_inplace = OnceDiffed(f∇f!)
# obj_outofplace = OnceDiffed(f∇f)
# obj_static = OnceDiffed(f∇fs)

# x0 = [-1.0, 0.0, 0.0]
# xopt = [1.0, 0.0, 0.0]
# x0s = @SVector [-1.0, 0.0, 0.0]

# println("Starting  from: ", x0)
# println("Targeting     : ", xopt)

# res = minimize(obj_inplace, copy(x0), NelderMead(), OptimizationOptions())
# print("NN  $(summary(NelderMead()))         ")
# @printf("%2.2e  %2.2e %2.2e %d\n",  norm(res.info.minimizer-xopt,Inf), norm(res.info.nm_obj, Inf), res.info.minimum, res.info.iter)
# res = minimize(obj_inplace, copy(x0), SimulatedAnnealing(), OptimizationOptions())
# print("NN  $(summary(SimulatedAnnealing()))         ")
# @printf("%2.2e  %2.2e %d\n",  norm(res.info.minimizer-xopt,Inf),  res.info.minimum, res.info.iter)
# res = minimize(obj_inplace, copy(x0), ConjugateGradient(update=HZ()), OptimizationOptions())
# print("NN  $(summary(ConjugateGradient(update=HZ())))         ")
# @printf("%2.2e  %2.2e %d\n",  norm(res.info.minimizer-xopt,Inf),  res.info.minimum, res.info.iter)
# res = minimize(obj_inplace, copy(x0), ConjugateGradient(update=CD()), OptimizationOptions())
# print("NN  $(summary(ConjugateGradient(update=CD())))         ")
# @printf("%2.2e  %2.2e %d\n",  norm(res.info.minimizer-xopt,Inf),  res.info.minimum, res.info.iter)
# res = minimize(obj_inplace, copy(x0), ConjugateGradient(update=HS()), OptimizationOptions())
# print("NN  $(summary(ConjugateGradient(update=HS())))         ")
# @printf("%2.2e  %2.2e %d\n",  norm(res.info.minimizer-xopt,Inf),  res.info.minimum, res.info.iter)
# res = minimize(obj_inplace, copy(x0), ConjugateGradient(update=FR()), OptimizationOptions())
# print("NN  $(summary(ConjugateGradient(update=FR())))         ")
# @printf("%2.2e  %2.2e %d\n",  norm(res.info.minimizer-xopt,Inf),  res.info.minimum, res.info.iter)
# res = minimize(obj_inplace, copy(x0), ConjugateGradient(update=PRP()), OptimizationOptions())
# print("NN  $(summary(ConjugateGradient(update=PRP(;plus=false))))         ")
# @printf("%2.2e  %2.2e %d\n",  norm(res.info.minimizer-xopt,Inf),  res.info.minimum, res.info.iter)
# res = minimize(obj_inplace, copy(x0), ConjugateGradient(update=PRP(plus=true)), OptimizationOptions())
# print("NN  $(summary(ConjugateGradient(update=PRP(;plus=true))))         ")
# @printf("%2.2e  %2.2e %d\n",  norm(res.info.minimizer-xopt,Inf),  res.info.minimum, res.info.iter)
# res = minimize(obj_inplace, copy(x0), ConjugateGradient(update=VPRP()), OptimizationOptions())
# print("NN  $(summary(ConjugateGradient(update=VPRP())))         ")
# @printf("%2.2e  %2.2e %d\n",  norm(res.info.minimizer-xopt,Inf),  res.info.minimum, res.info.iter)
# res = minimize(obj_inplace, copy(x0), ConjugateGradient(update=LS()), OptimizationOptions())
# print("NN  $(summary(ConjugateGradient(update=LS())))         ")
# @printf("%2.2e  %2.2e %d\n",  norm(res.info.minimizer-xopt,Inf),  res.info.minimum, res.info.iter)
# res = minimize(obj_inplace, copy(x0), ConjugateGradient(update=DY()), OptimizationOptions())
# print("NN  $(summary(ConjugateGradient(update=DY())))         ")
# @printf("%2.2e  %2.2e %d\n",  norm(res.info.minimizer-xopt,Inf),  res.info.minimum, res.info.iter)

# for _method in (GradientDescent, LBFGS, BFGS, DBFGS, DFP, SR1)
# # for _method in (GradientDescent, BFGS, DBFGS, DFP, SR1)
#     methodtxt = summary(_method())
#     for m in (Inverse(), Direct())
#         mtxt = m isa Inverse ? "(inverse): " : "(direct):  "
#         if _method == LBFGS && m isa Inverse
#             res = minimize!(obj_inplace, copy(x0), LineSearch(_method(m)), OptimizationOptions())
#             print("NN! $_method    $mtxt")
#             @printf("%2.2e  %2.2e %2.2e  %d\n", norm(res.info.minimizer-xopt,Inf), norm(res.info.∇fz, Inf), res.info.minimum, res.info.iter)
#         elseif _method !== LBFGS
#             res = minimize(obj_inplace, copy(x0), LineSearch(_method(m)), OptimizationOptions())
#             print("NN  $_method    $mtxt")
#             @printf("%2.2e  %2.2e %2.2e  %d\n", norm(res.info.minimizer-xopt,Inf), norm(res.info.∇fz, Inf), res.info.minimum, res.info.iter)
#             res = minimize!(obj_inplace, copy(x0), LineSearch(_method(m)), OptimizationOptions())
#             print("NN! $_method    $mtxt")
#             @printf("%2.2e  %2.2e %2.2e  %d\n", norm(res.info.minimizer-xopt,Inf), norm(res.info.∇fz, Inf), res.info.minimum, res.info.iter)
#             res = minimize(obj_static, x0s, LineSearch(_method(m)), OptimizationOptions())
#             print("NN  $_method(S) $mtxt")
#             @printf("%2.2e  %2.2e %2.2e  %d\n", norm(res.info.minimizer-xopt,Inf), norm(res.info.∇fz, Inf), res.info.minimum, res.info.iter)
#             if m isa Direct && !(_method == GradientDescent)
#                 println("Trust region: NWI")
#                 res = minimize!(obj_inplace, copy(x0), TrustRegion(_method(Direct()), NWI()), OptimizationOptions())
#                 print("NN! $_method    $mtxt")
#                 @printf("%2.2e  %2.2e %2.2e  %d\n", norm(res.info.minimizer-xopt,Inf), norm(res.info.∇fz, Inf), res.info.minimum, res.info.iter)
#                 println("Trust region: NTR")
#                 res = minimize!(obj_inplace, copy(x0), TrustRegion(_method(Direct()), NTR()), OptimizationOptions())
#                 print("NN! $_method    $mtxt")
#                 @printf("%2.2e  %2.2e %2.2e  %d\n", norm(res.info.minimizer-xopt,Inf), norm(res.info.∇fz, Inf), res.info.minimum, res.info.iter)
#             end
#         end
#     end
# end
# println()

# res = minimize(obj_inplace, x0, LineSearch(BFGS(Inverse())), OptimizationOptions())
# @test res.info.iter == 30
# @printf("NN  BFGS    (inverse): %2.2e  %2.2e %2.2e %d\n", norm(res.info.minimizer-xopt,Inf), norm(res.info.∇fz, Inf), res.info.minimum, res.info.iter)
# res = minimize!(obj_inplace, copy(x0), LineSearch(BFGS(Inverse()), Backtracking()), OptimizationOptions())
# @test res.info.iter == 30
# @printf("NN! BFGS    (inverse): %2.2e  %2.2e %2.2e %d\n",  norm(res.info.minimizer-xopt,Inf), norm(res.info.∇fz, Inf), res.info.minimum, res.info.iter)
# res = minimize!(obj_inplace, copy(x0), LineSearch(BFGS(Inverse()), Backtracking(interp=FFQuadInterp())), OptimizationOptions())
# @test res.info.iter == 30
# @printf("NN! BFGS    (inverse, quad): %2.2e %2.2e %2.2e %d\n",  norm(res.info.minimizer-xopt,Inf), norm(res.info.∇fz, Inf), res.info.minimum, res.info.iter)
# res = minimize(obj_static, x0s, LineSearch(BFGS(Inverse())), OptimizationOptions())
# @test res.info.iter == 30
# @printf("NN  BFGS(S) (inverse): %2.2e %2.2e %2.2e %d\n",  norm(res.info.minimizer-xopt,Inf), norm(res.info.∇fz, Inf), res.info.minimum, res.info.iter)
# res = minimize(obj_static, x0s, LineSearch(BFGS(Inverse()), Backtracking(interp=FFQuadInterp())), OptimizationOptions())
# @test res.info.iter == 30
# @printf("NN  BFGS(S) (inverse, quad): %2.2e %2.2e %2.2e %d\n",  norm(res.info.minimizer-xopt,Inf), norm(res.info.∇fz, Inf), res.info.minimum, res.info.iter)

# x0 = rand(3)
# x0s = SVector{3}(x0)
# println("\nFrom a random point: ", x0)
# res = minimize(obj_inplace, copy(x0), NelderMead(), OptimizationOptions())
# print("NN  $(summary(NelderMead()))         ")
# @printf("%2.2e  %2.2e %2.2e %d\n",  norm(res.info.minimizer-xopt,Inf), norm(res.info.nm_obj, Inf), res.info.minimum, res.info.iter)
# res = minimize(obj_inplace, copy(x0), SimulatedAnnealing(), OptimizationOptions())
# print("NN  $(summary(SimulatedAnnealing()))         ")
# @printf("%2.2e  %2.2e %d\n",  norm(res.info.minimizer-xopt,Inf),  res.info.minimum, res.info.iter)
# res = minimize(obj_inplace, copy(x0), ConjugateGradient(), OptimizationOptions())
# print("NN  $(summary(ConjugateGradient()))         ")
# @printf("%2.2e  %2.2e %d\n",  norm(res.info.minimizer-xopt,Inf),  res.info.minimum, res.info.iter)

# for _method in (GradientDescent, LBFGS, BFGS, DBFGS, DFP, SR1)
#     methodtxt = summary(_method())
#     for m in (Inverse(), Direct())
#         mtxt = m isa Inverse ? "(inverse): " : "(direct):  "
#         if _method == LBFGS && m isa Inverse
#             res = minimize!(obj_inplace, copy(x0), LineSearch(_method(m)), OptimizationOptions())
#             print("NN! $_method    $mtxt")
#             @printf("%2.2e  %2.2e %2.2e  %d\n", norm(res.info.minimizer-xopt,Inf), norm(res.info.∇fz, Inf), res.info.minimum, res.info.iter)
#         elseif _method !== LBFGS
#             res = minimize(obj_inplace, copy(x0), LineSearch(_method(m)), OptimizationOptions())
#             print("NN  $_method    $mtxt")
#             @printf("%2.2e  %2.2e %2.2e %d\n",  norm(res.info.minimizer-xopt,Inf), norm(res.info.∇fz, Inf), res.info.minimum, res.info.iter)
#             res = minimize!(obj_inplace, copy(x0), LineSearch(_method(m)), OptimizationOptions())
#             print("NN! $_method    $mtxt")
#             @printf("%2.2e  %2.2e %2.2e %d\n",  norm(res.info.minimizer-xopt,Inf), norm(res.info.∇fz, Inf), res.info.minimum, res.info.iter)
#             res = minimize(obj_static, x0s, LineSearch(_method(m)), OptimizationOptions())
#             print("NN  $_method(S) $mtxt")
#             @printf("%2.2e  %2.2e %2.2e %d\n",  norm(res.info.minimizer-xopt,Inf), norm(res.info.∇fz, Inf), res.info.minimum, res.info.iter)
#             if m isa Direct && !(_method == GradientDescent)
#                 println("Trust region: NWI")
#                 res = minimize!(obj_inplace, copy(x0), TrustRegion(_method(Direct()), NWI()), OptimizationOptions())
#                 print("NN! $_method    $mtxt")
#                 @printf("%2.2e  %2.2e %2.2e %d\n",  norm(res.info.minimizer-xopt,Inf), norm(res.info.∇fz, Inf), res.info.minimum, res.info.iter)
#                 println("Trust region: NTR")
#                 res = minimize!(obj_inplace, copy(x0), TrustRegion(_method(Direct()), NTR()), OptimizationOptions())
#                 print("NN! $_method    $mtxt")
#                 @printf("%2.2e  %2.2e %2.2e %d\n",  norm(res.info.minimizer-xopt,Inf), norm(res.info.∇fz, Inf), res.info.minimum, res.info.iter)
#             end
#         end
#     end
# end
# println()

# function himmelblau!(x, ∇f)
#     if !(∇f == nothing)
#         ∇f[1] = 4.0 * x[1]^3 + 4.0 * x[1] * x[2] -
#             44.0 * x[1] + 2.0 * x[1] + 2.0 * x[2]^2 - 14.0
#         ∇f[2] = 2.0 * x[1]^2 + 2.0 * x[2] - 22.0 +
#             4.0 * x[1] * x[2] + 4.0 * x[2]^3 - 28.0 * x[2]
#     end

#     fx = (x[1]^2 + x[2] - 11)^2 + (x[1] + x[2]^2 - 7)^2
#     objective_return(fx, ∇f)
# end

# function himmelblaus(x, ∇f)
#     fx = (x[1]^2 + x[2] - 11)^2 + (x[1] + x[2]^2 - 7)^2
#     if !(∇f == nothing)
#         ∇f1 = 4.0 * x[1]^3 + 4.0 * x[1] * x[2] -
#             44.0 * x[1] + 2.0 * x[1] + 2.0 * x[2]^2 - 14.0
#         ∇f2 = 2.0 * x[1]^2 + 2.0 * x[2] - 22.0 +
#             4.0 * x[1] * x[2] + 4.0 * x[2]^3 - 28.0 * x[2]
#         ∇f = @SVector([∇f1, ∇f2])
#     end
#     objective_return(fx, ∇f)
# end

# function himmelblau(x, ∇f)
#     g = ∇f == nothing ? ∇f : similar(x)

#     return himmelblau!(x, g)
# end

# him_inplace = OnceDiffed(himmelblau!)
# him_static = OnceDiffed(himmelblaus)
# him_outofplace = OnceDiffed(himmelblau)

# println("\nHimmelblau function")
# x0 = [3.0, 1.0]
# x0s = SVector{2}(x0)
# minimizers = [[3.0,2.0],[-2.805118,3.131312],[-3.779310,-3.283186],[3.584428,-1.848126]]
# for _method in (GradientDescent, LBFGS, BFGS, DBFGS, DFP, SR1)
# # for _method in (GradientDescent, BFGS, DBFGS, DFP, SR1)
#     methodtxt = summary(_method())
#     for m in (Inverse(), Direct())
#         mtxt = m isa Inverse ? "(inverse): " : "(direct):  "
#         if _method == LBFGS && m isa Inverse
#             res = minimize!(him_inplace, copy(x0), LineSearch(_method(m)), OptimizationOptions())
#             print("NN! $_method    $mtxt")
#             @printf("%2.2e  %2.2e %2.2e  %d\n", minimum([norm(res.info.minimizer-xopt,Inf) for xopt in minimizers]), norm(res.info.∇fz, Inf), res.info.minimum, res.info.iter)
#         elseif _method !== LBFGS
#             res = minimize(him_outofplace, copy(x0), LineSearch(_method(m)), OptimizationOptions())
#             print("NN  $_method    $mtxt")
#             @printf("%2.2e  %2.2e %2.2e  %d\n", minimum([norm(res.info.minimizer-xopt,Inf) for xopt in minimizers]), norm(res.info.∇fz, Inf), res.info.minimum, res.info.iter)
#             res = minimize!(him_inplace, copy(x0), LineSearch(_method(m)), OptimizationOptions())
#             print("NN! $_method    $mtxt")
#             @printf("%2.2e  %2.2e %2.2e  %d\n", minimum([norm(res.info.minimizer-xopt,Inf) for xopt in minimizers]), norm(res.info.∇fz, Inf), res.info.minimum, res.info.iter)
#             res = minimize(him_static, x0s, LineSearch(_method(m)), OptimizationOptions())
#             print("NN  $_method(S) $mtxt")
#             @printf("%2.2e  %2.2e %2.2e  %d\n", minimum([norm(res.info.minimizer-xopt,Inf) for xopt in minimizers]), norm(res.info.∇fz, Inf), res.info.minimum, res.info.iter)
#             if m isa Direct && !(_method == GradientDescent)
#                 println("Trust region: NWI")
#                 res = minimize!(him_inplace, copy(x0), TrustRegion(_method(Direct()), NWI()), OptimizationOptions())
#                 print("NN! $_method    $mtxt")
#                 @printf("%2.2e  %2.2e %2.2e  %d\n", minimum([norm(res.info.minimizer-xopt,Inf) for xopt in minimizers]), norm(res.info.∇fz, Inf), res.info.minimum, res.info.iter)
#                 println("Trust region: NTR")
#                 res = minimize!(him_inplace, copy(x0), TrustRegion(_method(Direct()), NTR()), OptimizationOptions())
#                 print("NN! $_method    $mtxt")
#                 @printf("%2.2e  %2.2e %2.2e  %d\n", minimum([norm(res.info.minimizer-xopt,Inf) for xopt in minimizers]), norm(res.info.∇fz, Inf), res.info.minimum, res.info.iter)
#             end
#         end
#     end
# end
# println()
# println()
# xrand = rand(2)
# xrands = SVector{2}(xrand)
# println("\nFrom a random point: ", xrand)

# for _method in (GradientDescent, LBFGS, BFGS, DBFGS, SR1)
# # for _method in (GradientDescent, BFGS, DBFGS, DFP, SR1)
#     methodtxt = summary(_method())
#     for m in (Inverse(), Direct())
#         mtxt = m isa Inverse ? "(inverse): " : "(direct):  "
#         if _method == LBFGS && m isa Inverse
#             res = minimize!(him_inplace, copy(x0), LineSearch(_method(m)), OptimizationOptions())
#             print("NN! $_method    $mtxt")
#             @printf("%2.2e  %2.2e %2.2e  %d\n", minimum([norm(res.info.minimizer-xopt,Inf) for xopt in minimizers]), norm(res.info.∇fz, Inf), res.info.minimum, res.info.iter)
#         elseif _method !== LBFGS
#             res = minimize(him_outofplace, copy(xrand), LineSearch(_method(m)), OptimizationOptions())
#             print("NN  $_method    $mtxt")
#             @printf("%2.2e  %2.2e %2.2e  %d\n", minimum([norm(res.info.minimizer-xopt,Inf) for xopt in minimizers]), norm(res.info.∇fz, Inf), res.info.minimum, res.info.iter)
#             res = minimize!(him_inplace, copy(xrand), LineSearch(_method(m)), OptimizationOptions())
#             print("NN! $_method    $mtxt")
#             @printf("%2.2e  %2.2e %2.2e  %d\n", minimum([norm(res.info.minimizer-xopt,Inf) for xopt in minimizers]), norm(res.info.∇fz, Inf), res.info.minimum, res.info.iter)
#             res = minimize(him_static, xrands, LineSearch(_method(m)), OptimizationOptions())
#             print("NN  $_method(S) $mtxt")
#             @printf("%2.2e  %2.2e %2.2e  %d\n", minimum([norm(res.info.minimizer-xopt,Inf) for xopt in minimizers]), norm(res.info.∇fz, Inf), res.info.minimum, res.info.iter)
#             if m isa Direct && !(_method == GradientDescent)
#                 println("Trust region: NWI")
#                 res = minimize!(him_inplace, copy(xrand), TrustRegion(_method(Direct()), NWI()), OptimizationOptions())
#                 print("NN! $_method    $mtxt")
#                 @printf("%2.2e  %2.2e %2.2e  %d\n", minimum([norm(res.info.minimizer-xopt,Inf) for xopt in minimizers]), norm(res.info.∇fz, Inf), res.info.minimum, res.info.iter)
#                 println("Trust region: NTR")
#                 res = minimize!(him_inplace, copy(xrand), TrustRegion(_method(Direct()), NTR()), OptimizationOptions())
#                 print("NN! $_method    $mtxt")
#                 @printf("%2.2e  %2.2e %2.2e  %d\n", minimum([norm(res.info.minimizer-xopt,Inf) for xopt in minimizers]), norm(res.info.∇fz, Inf), res.info.minimum, res.info.iter)
#             end
#         end
#     end
# end
# println()

# end

#=
using NLSolvers
function f∇f!(x, ∇f)
   if !(∇f==nothing)
       if ( x[1]^2 + x[2]^2 == 0 )
           dtdx1 = 0;
           dtdx2 = 0;
       else
           dtdx1 = - x[2] / ( 2 * pi * ( x[1]^2 + x[2]^2 ) );
           dtdx2 =   x[1] / ( 2 * pi * ( x[1]^2 + x[2]^2 ) );
       end
       ∇f[1] = -2000.0*(x[3]-10.0*theta(x))*dtdx1 +
           200.0*(sqrt(x[1]^2+x[2]^2)-1)*x[1]/sqrt( x[1]^2+x[2]^2 );
       ∇f[2] = -2000.0*(x[3]-10.0*theta(x))*dtdx2 +
           200.0*(sqrt(x[1]^2+x[2]^2)-1)*x[2]/sqrt( x[1]^2+x[2]^2 );
       ∇f[3] =  200.0*(x[3]-10.0*theta(x)) + 2.0*x[3];
   end

   fx = f(x)
   return ∇f==nothing ? fx : (fx, ∇f)
end

f(x) = 100.0 * ((x[3] - 10.0 * theta(x))^2 + (sqrt(x[1]^2 + x[2]^2) - 1.0)^2) + x[3]^2
function theta(x)
  if x[1] > 0
      return atan(x[2] / x[1]) / (2.0 * pi)
  else
      return (pi + atan(x[2] / x[1])) / (2.0 * pi)
  end
end

#= test solve interface =#
obj! = OnceDiffed(f∇f!)
nm_prob! = MinProblem(obj=obj!)
solve(nm_prob!, rand(3), NelderMead(), OptimizationOptions())

solve(obj!, -rand(3)*9 .- 3, NLSolvers.NelderMead(), OptimizationOptions())
V = [[1.0,1.0,1.0], [0.0,1.0,1.0],[0.40,0.0,0.0],[-1.0,2.0,.03]]
F = obj!.(V)

splx = NLSolvers.ValuedSimplex(V, F)

solve(obj!, splx, NLSolvers.NelderMead(), OptimizationOptions())
solve(nm_prob!, splx, NelderMead(), OptimizationOptions())

function powell(x, ∇f)
    fx = (x[1] + 10.0 * x[2])^2 + 5.0 * (x[3] - x[4])^2 +
        (x[2] - 2.0 * x[3])^4 + 10.0 * (x[1] - x[4])^4

    if !isa(∇f, Nothing)
        ∇f[1] = 2.0 * (x[1] + 10.0 * x[2]) + 40.0 * (x[1] - x[4])^3
        ∇f[2] = 20.0 * (x[1] + 10.0 * x[2]) + 4.0 * (x[2] - 2.0 * x[3])^3
        ∇f[3] = 10.0 * (x[3] - x[4]) - 8.0 * (x[2] - 2.0 * x[3])^3
        ∇f[4] = -10.0 * (x[3] - x[4]) - 40.0 * (x[1] - x[4])^3
        return ∇f, fx
    end
    return fx
end

obj_powell = OnceDiffed(powell)
# Define vertices in V
x0 = [1.0,1.0,1.0,1.0]
V = [copy(x0)]
for i = 1:4
    push!(V, x0+39*rand(4))
end
V
F = obj_powell.(V)
splx = NLSolvers.ValuedSimplex(V, F)
@time solve(obj_powell, splx, NLSolvers.NelderMead(), OptimizationOptions(maxiter=3000))
@time solve(obj_powell, copy(x0), NLSolvers.NelderMead(), OptimizationOptions(maxiter=3000))
@allocated solve(obj_powell, splx, NLSolvers.NelderMead(), OptimizationOptions(maxiter=3000))
solve(obj_powell, copy(x0), NLSolvers.NelderMead(), OptimizationOptions(maxiter=3000))

function extros!(x, storage)
   n = length(x)
   jodd = 1:2:n-1
   jeven = 2:2:n
   xt = similar(x)
   @. xt[jodd] = 10.0 * (x[jeven] - x[jodd]^2)
   @. xt[jeven] = 1.0 - x[jodd]

   if !isa(storage, Nothing)
       @. storage[jodd] = -20.0 * x[jodd] * xt[jodd] - xt[jeven]
       @. storage[jeven] = 10.0 * xt[jodd]
       return 0.5*sum(abs2, xt), storage
   end
   return 0.5*sum(abs2, xt)
end
x0=rand(300)
V = [copy(x0)]
for i = 1:300
    push!(V, x0+39*rand(length(x0)))
end
V
extros_obj! = OnceDiffed(extros!)
F = extros_obj!.(V)
splx = NLSolvers.ValuedSimplex(V, F)
@time solve(extros_obj!, splx, NLSolvers.NelderMead(), OptimizationOptions(maxiter=3000))
=#


text:
import Base: randn,transpose,transpose!,range

export
    randn,
    randu,
    constant,
    lookup,
    transpose,
    transpose!,
	range

immutable Create <: AFImpl
    randn::Ptr{Void}
    randu::Ptr{Void}
    constant::Ptr{Void}
    constantLong::Ptr{Void}
    constantULong::Ptr{Void}
    lookup::Ptr{Void}
    transpose::Ptr{Void}
    transposeInPlace::Ptr{Void}
	range::Ptr{Void}

    function Create(ptr)
        new(
            Libdl.dlsym(ptr, :af_randn),
            Libdl.dlsym(ptr, :af_randu),
            Libdl.dlsym(ptr, :af_constant),
            Libdl.dlsym(ptr, :af_constant_long),
            Libdl.dlsym(ptr, :af_constant_ulong),
            Libdl.dlsym(ptr, :af_lookup),
            Libdl.dlsym(ptr, :af_transpose),
            Libdl.dlsym(ptr, :af_transpose_inplace),
			Libdl.dlsym(ptr, :af_range)
        )
    end
end

function randn{B, T}(af::ArrayFire{B}, ::Type{T}, dims::DimT...)
    ptr = af.results.ptr
    dims2 = collectDims(af.results.dims, dims)
    err = ccall(af.create.randn,
        Cint, (Ptr{Ptr{Void}}, Cuint, Ptr{DimT}, DType),
        ptr, length(dims2), pointer(dims2), asDType(T))
    assertErr(err)
    AFArray{B}(af, ptr[])
end

function randu{B, T}(af::ArrayFire{B}, ::Type{T}, dims::DimT...)
    ptr = af.results.ptr
    dims2 = collectDims(af.results.dims, dims)
    err = ccall(af.create.randu,
        Cint, (Ptr{Ptr{Void}}, Cuint, Ptr{DimT}, DType),
        ptr, length(dims2), pointer(dims2), asDType(T))
    assertErr(err)
    AFArray{B}(af, ptr[])
end

function constant{B, T<:Real}(af::ArrayFire{B}, value::T, dims::DimT...)
    ptr = af.results.ptr
    dims2 = collectDims(af.results.dims, dims)
    err = ccall(af.create.constant,
        Cint, (Ptr{Ptr{Void}}, Float64, Cuint, Ptr{DimT}, DType),
        ptr, Float64(value), length(dims2), pointer(dims2), asDType(T))
    assertErr(err)
    AFArray{B}(af, ptr[])
end

function constant{B}(af::ArrayFire{B}, value::Int64, dims::DimT...)
    ptr = af.results.ptr
    dims2 = collectDims(af.results.dims, dims)
    err = ccall(af.create.constantLong,
        Cint, (Ptr{Ptr{Void}}, Int64, Cuint, Ptr{DimT}, DType),
        ptr, value, length(dims2), pointer(dims2), asDType(Int64))
    assertErr(err)
    AFArray{B}(af, ptr[])
end

function constant{B}(af::ArrayFire{B}, value::UInt64, dims::DimT...)
    ptr = af.results.ptr
    dims2 = collectDims(af.results.dims, dims)
    err = ccall(af.create.constantLong,
        Cint, (Ptr{Ptr{Void}}, UInt64, Cuint, Ptr{DimT}, DType),
        ptr, value, length(dims2), pointer(dims2), asDType(UInt64))
    assertErr(err)
    AFArray{B}(af, ptr[])
end

@afCall_Arr_Arr_Arr_Unsigned(lookup, create, lookup)

@afCall_Arr_Arr_Bool(transpose, create, transpose, false)

function transpose!{B}(arr::AFArray{B}, conjugate::Bool = false)
    verifyAccess(arr)
    af = arr.af
    err = ccall(af.create.transposeInPlace,
        Cint, (Ptr{Void}, Bool),
        arr.ptr, conjugate)
    assertErr(err)
    arr
end

function range{B, T}(af::ArrayFire{B}, ::Type{T}, dims::Dim4, seqDim = 0)
    ptr = af.results.ptr
    err = ccall(af.create.range,
        Cint, (Ptr{Ptr{Void}}, Cuint, Ptr{DimT}, Cint, DType),
        ptr, length(dims), pointer(dims), seqDim, asDType(T))
    assertErr(err)
    AFArray{B}(af, ptr[])
end

ext: jl

max_stars_repo_name: unbornchikken/julia-ml-proto

max_issues_repo_name: unbornchikken/julia-ml-proto

max_forks_repo_name: unbornchikken/julia-ml-proto


text:
module test_bukdu_plugs_conn

using Test
import HTTP
import Dates: DateTime

@testset "cookies" begin
    req = HTTP.Messages.Request()
    z = HTTP.Cookies.Cookie("special-2", "z")
    HTTP.Messages.setheader(req, "Cookie" => String(z))
    @test HTTP.Messages.header(req.headers, "Cookie") == "special-2=z"
    cookies = HTTP.Cookies.cookies(req)
    @test cookies isa Vector{HTTP.Cookies.Cookie}
    @test first(cookies) == z
    cookies = HTTP.Cookies.readcookies(req.headers, "special-2")
    @test first(cookies) == z

    expires_at = DateTime(2021,06,09, 10,18,14)
    resp_cookie = HTTP.Cookies.Cookie("theme", "light", expires=expires_at)
    HTTP.Messages.setheader(req.response, "Set-Cookie" => String(resp_cookie, false))
    @test HTTP.Messages.header(req.response.headers, "Set-Cookie") == "theme=light; Expires=Wed, 09 Jun 2021 10:18:14 GMT"
end

end # module test_bukdu_plugs_conn

max_stars_repo_name: wookay/TestBukdu

max_issues_repo_name: wookay/TestBukdu

max_forks_repo_name: wookay/TestBukdu


text:
module JTFunctions

    ######
    using TaylorSeries, TaylorIntegration, Distributions, Plots, DataFrames
    pyplot()
    import TaylorSeries: NumberNotSeries

    export evaluate_neighborhood, evaluate_distribution,
           circle2, square2, ξmax, ξmax_lastT, anal_vs_taylor2D,
                area_of_polygon, separation_rate, remoteness,
           grid_ξmax, grid_FTLE, grid_seprate, vectorField_plot,
           harmonic_oscillator!, simple_pendulum!, artificial_ode!,
           stability_matrix, simplecticity, hp_filter,
           myfonts
    ######

    ##### Plotting attributes #####
    fontXSmall = Plots.font("Helvetica", 9) #,hcenter, :vcenter, 0.0, RGB{N0f8})
    fontSmall = Plots.font("Helvetica", 11)
    fontMed = Plots.font("Helvetica", 12)
    fontBig = Plots.font("Helvetica", 14)
    fontXBig =Plots.font("Helvetica", 16)
    myfonts = Dict(:guidefont=>fontBig, :xtickfont=>fontMed, :ytickfont=>fontMed, :legendfont=>fontSmall)

    ##### Functions & Indicators ####
    # Neighborhood Evaluation
    function evaluate_neighborhood{T<:NumberNotSeries}(δU::Function,ϕ::Taylor1{T},num_vals=100)
        τrange = linspace(0,1,num_vals)
        length(δU(τrange[1])) != 1 ? error("δU is not of the same dimension of the phase space.") : nothing

        return evaluate.(ϕ, δU.(τrange))
    end

    function evaluate_neighborhood{T<:NumberNotSeries}(δU::Function,ϕ::TaylorN{T},num_vals::Integer=100,ξ::Real=0.1)
        τrange = linspace(0,1,num_vals)
        length(δU(τrange[1])) != get_numvars() ? error("δU is not of the same dimension of the phase space.") : nothing
        !(δU(τrange[1]) ≈ δU(τrange[end])) ? error("δU is not a closed surface in the range given.") : nothing

        return evaluate.(ϕ, δU.(τrange,r=ξ))
        end

    function evaluate_neighborhood{T<:AbstractSeries}(δU::Function,ϕ::Array{T,1};num_vals::Integer=100,ξ::Real=0.1)
        return evaluate_neighborhood.(δU,ϕ,num_vals,ξ)
    end

    function evaluate_distribution{T<:NumberNotSeries}(ϕ::TaylorN{T};U::Array{T}=zeros(10,2) )
        distributed_ϕ = Vector{eltype(ϕ)}()

        for i in eachindex(U[:,1])
            push!(distributed_ϕ, ϕ(U[i,:]))
        end

        return distributed_ϕ
    end

        dof = get_numvars()

    function evaluate_distribution{T<:NumberNotSeries}(ϕ::Vector{TaylorN{T}},U)
        return evaluate_distribution.(ϕ,U=U)
    end

    function evaluate_distribution{T<:NumberNotSeries}(ϕ::Array{TaylorN{T},2},dist::Distribution=Normal(0,0.01);num_vals::Integer=200)
        dof = get_numvars()
        U = rand(dist, num_vals, dof)

        return evaluate_distribution(ϕ[:,1],U), evaluate_distribution(ϕ[:,2],U), U
    end

    # Some Neighborhood Parametrization (this should deserve a special type Neighborhood)
    circle2(τ;r=0.1) = [r*cos(2π*τ),r*sin(2π*τ)]
    square2(τ) = nothing

    # Area of Neighborhood
    function area_of_polygon(pjet,qjet)
        #points_per_polygon = length(pjet[1])
        num_time_steps = length(pjet)

        area_per_polygon = zeros(eltype(pjet[1]), num_time_steps)
        for i in 1:num_time_steps

            p_aux = copy(pjet[i])

            p_end = shift!(p_aux)
            push!(p_aux,p_end)

            A_1 = sum(p_aux.*qjet[i])

            q_aux = copy(qjet[i])

            q_end = shift!(q_aux)
            push!(q_aux,q_end)

            A_2 = sum(q_aux.*pjet[i])

            area_per_polygon[i] = (A_1 - A_2)/2.0
        end
        return area_per_polygon
    end

    # Maximum Box Size
    function ξmax(ϕ::Vector{T};ϵ=1e-5) where T<:TaylorN
        dof = length(ϕ)
        ord = ϕ[1].order
        ξ = Inf

        for m in 1:dof
            a_k = ϕ[m].coeffs[ord+1].coeffs

            #In case the last term is zero.
            k = 1
            while a_k == zeros(a_k)
                if k > ord
                    break
                end
                a_k = ϕ[m].coeffs[ord+1-k].coeffs
                k += 1
            end

            for a_jk in a_k
                ξ = min(ξ,abs(ϵ/a_jk)^(1./ord) )
            end
        end

        return ξ
    end

    function ξmax(ϕ::Vector{T};ϵ=1e-5) where T<:Taylor1
        dof = length(ϕ)
        ord = ϕ[1].order
        ξ = Inf

        for m in 1:dof
            a_k = ϕ[m].coeffs[ord+1]

            #In case the last term is zero.
            k = 1
            while a_k == zero(a_k)
                if k > ord
                    break
                end
                a_k = ϕ[m].coeffs[ord+1-k]
                k += 1
            end

            ξ = min(ξ,abs(ϵ/a_k)^(1./ord) )
        end

        return ξ
    end

    function ξmax(ϕ::Matrix{T};ϵ=1e-5) where T<:AbstractSeries

        total_timesteps = size(ϕ)[1]
        ξmax_T = zeros(total_timesteps)

        for t in 1:total_timesteps
            ξmax_T[t] = ξmax(ϕ[t,:];ϵ=ϵ)
        end

        return ξmax_T
    end

    ξmax_lastT(ϕ;ϵ=1e-5) = ξmax(ϕ[end,:],ϵ=ϵ)

    remoteness(ϕμ, Δμ) = sqrt.( (ϕμ[:,1](Δμ) - ϕμ[1,1](Δμ)).^2 + (ϕμ[:,2](Δμ) - ϕμ[1,2](Δμ)).^2 )
    function separation_rate(ϕN;neighborhood_vals::Integer=100)

        ξ_max = ξmax_lastT(ϕN) #should ξ_max be an argument?
        qjet = evaluate_neighborhood(circle2,ϕN[end,1],neighborhood_vals,ξ_max) .- evaluate(ϕN[end,1])
        pjet = evaluate_neighborhood(circle2,ϕN[end,2],neighborhood_vals,ξ_max) .- evaluate(ϕN[end,2])

        P_n = hcat(qjet,pjet)

        P_max = -Inf
        P_min =  Inf

        τ = linspace(0,1,neighborhood_vals)
        θ_max = 0.0
        θ_min = 0.0

        #i_min = 0
        #i_max = 0

        for i in 1:length(qjet)
            P_n_norm = norm(P_n[i,:])

            if P_n_norm > P_max
                P_max = P_n_norm
                θ_max = τ[i]*2π
                #i_max = i
            end

            if P_n_norm < P_min
                P_min = P_n_norm
                θ_min = τ[i]*2π
                #i_min = i
            end
        end

        return P_max/ξ_max, P_min/ξ_max, θ_max, θ_min #, i_max, i_min
    end

    # Compare taylorinteg solution `x,y` vs analytical ones `xa`
    anal_vs_taylor2D(x,y,xa) = sqrt.((x - xa[:,1]).^2 + (y - xa[:,2]).^2)
    # # Compare 2 solutions
    # function JT_accuracy(ϕ,ϕ_JT)
    #     dims = size(ϕ)[2]
    #     diffs = zero(ϕ[:,1])
    #     for dim in 1:dims
    #         diffs += (ϕ[:,dim] - ϕ_JT[:,dim]).^2
    #     end
    #     return sqrt.(diffs)
    # end

    #### PLOTTING & GRIDS ####
    #No parallelization is done yet.

    #ξmax grid without the plot.. should the plot be included?
    #### This should be eliminated after PR is merged
    import TaylorSeries: get_variables
    get_variables(order) = [TaylorN(i,order=order) for i in 1:get_numvars()]
    ####

    function grid_ξmax{T<:Real}(eqs_diff!::Function,t0::T,tmax::T,μ::T;
                        xlim::Tuple=(-1,1),ylim::Tuple=(-1,1),num_points::Integer=30,
                        order_jet::Integer=3,order_taylor::Integer=25,abstol=1e-10)
        @assert get_numvars() == 2 "ξmax must be 2-dimensional."
        δx,δy = get_variables(order_jet)

        xgrid= linspace(xlim...,num_points)
        ygrid= linspace(ylim...,num_points)
        heatgrid = zeros(length(xgrid),length(ygrid))

        L3_ =  -(1 + (5μ/12))
        L2_ = 1 + (μ/3)^(1/3)
        ε₁,ε₂ = (-μ - L3_)/1.5, ( L2_ - (1-μ) )/1.5

        #Parallelize here!
        for (i,x) in enumerate(xgrid)
            for (j,y) in enumerate(ygrid)

                if (norm( [x - (-μ), y] ) > ε₁) && (norm( [x - (1-μ), y] ) > ε₂)
                q0TN = [x+δx,y+δy,0.0,0.0]
                _,ϕN = taylorinteg(eqs_diff!,q0TN,t0,tmax,order_taylor,abstol,maxsteps=800);
                heatgrid[i,j] = ξmax_lastT(ϕN)
                end

            end

        end
        return xgrid, ygrid, heatgrid'
    end

    #FTLE scalar field using TaylorIntegration lyapunov expectrum.
    #TODO: FTLEs could also be made with MY method... (JT of order = 1 jets)
    function grid_FTLE{T<:Real}(eqs_diff!::Function,t0::T,tmax::T,μ::T;
                            xlim::Tuple=(-1,1),ylim::Tuple=(-1,1),num_points::Integer=30,
                            order_taylor::Integer=25,abstol=1e-10)

        qgrid = linspace(xlim...,num_points)
        pgrid = linspace(ylim...,num_points)
        FTLE = zeros(num_points,num_points)

        L3_ =  -(1 + (5μ/12))
        L2_ = 1 + (μ/3)^(1/3)
        ε₁,ε₂ = (-μ - L3_)/1.5, ( L2_ - (1-μ) )/1.5

        #Parallelize here!
        for (i,q) in enumerate(qgrid)
            for (j,p) in enumerate(pgrid)

                if (norm( [x - (-μ), y] ) > ε₁) && (norm( [x - (1-μ), y] ) > ε₂)
                x0 = [q,p]
                _,__,λ = liap_taylorinteg(eqs_diff!,x0,t0,tmax,order_taylor,abstol)
                FTLE[i,j] = maximum(λ[end,:])
                end

            end
        end

        return qgrid, pgrid, FTLE'
    end

    function grid_seprate{T<:Real}(eqs_diff!::Function,t0::T,tmax::T,μ::T;
                        xlim::Tuple=(-1,1),ylim::Tuple=(-1,1),num_points::Integer=30,
                        order_jet::Integer=3,order_taylor::Integer=25,abstol=1e-10,
                        neighborhood_vals::Integer=100)
        @assert get_numvars() == 2 "ξmax must be 2-dimensional."
        δx,δy = get_variables(order_jet)

        #allocation
        #grid for heatmap allocation
        xgrid = linspace(xlim...,num_points)
        ygrid = linspace(ylim...,num_points)
        #matrix allocation for ξmax
        heatgrid_ξmax = zeros(length(xgrid),length(ygrid))
        #matrix allocation for separation rates
        heatgrid_sepRate_max = zeros(length(xgrid),length(ygrid))
        heatgrid_sepRate_min = zeros(length(xgrid),length(ygrid))

        k  = 0
        Δr = max(xgrid.step.hi,ygrid.step.hi)
        #grid allocation for separation rates vector field
        xy_grid = Vector{Tuple{Float64,Float64}}(num_points^2)
        vecField_sepRate_max = Vector{Tuple{Float64,Float64}}(num_points^2)
        vecField_sepRate_min = Vector{Tuple{Float64,Float64}}(num_points^2)

        L3_ =  -(1 + (5μ/12))
        L2_ = 1 + (μ/3)^(1/3)
        ε₁,ε₂ = (-μ - L3_)/1.5, ( L2_ - (1-μ) )/1.5

        #FTLE is missing...

        #Parallelize here!
        for (i,x) in enumerate(xgrid)
            for (j,y) in enumerate(ygrid)

                if (norm( [x - (-μ), y] ) > ε₁) && (norm( [x - (1-μ), y] ) > ε₂)
                    #initial TaylorN condition
                    q0TN = [x+δx,y+δy,0,0]
                    #Jet Trasnport Integration
                    _,ϕN = taylorinteg(eqs_diff!,q0TN,t0,tmax,order_taylor,abstol);
                    #indicators' evaluations
                    ξ_max = ξmax_lastT(ϕN)
                    P_max,P_min,θ_max,θ_min = separation_rate(ϕN,neighborhood_vals=neighborhood_vals)

                    #vector field grids as series of tuples
                    k = num_points*(i-1)+j
                    xy_grid[k] = (x,y)
                    vecField_sepRate_max[k] = (0.45 * Δr) .* ( cos(θ_max), sin(θ_max) )
                    vecField_sepRate_min[k] = (0.45 * Δr) .* ( cos(θ_min), sin(θ_min) )

                    #heatmap matrices
                    heatgrid_ξmax[i,j] = ξ_max
                    heatgrid_sepRate_max[i,j] = P_max
                    heatgrid_sepRate_min[i,j] = P_min
                else
                    k = num_points*(i-1)+j
                    xy_grid[k] = (x,y)
                    vecField_sepRate_max[k] = (0.0, 0.0 )
                    vecField_sepRate_min[k] = (0.0, 0.0 )
                end

            end

        end
        return xgrid, ygrid, heatgrid_ξmax', heatgrid_sepRate_max', heatgrid_sepRate_min', xy_grid, vecField_sepRate_max, vecField_sepRate_min
    end

    #2D vector Field plotting for an TaylorIntegration defined eqs_motion!
    #See some options for beautifying (optional)
    function vectorField_plot(vecField!::Function,t0=0.0;xlim::Tuple=(-1,1),ylim::Tuple=(-1,1),
                              num_points::Int=15,arrowsize::Real=1.2,gridlog::Bool=false)
        X = linspace(xlim...,num_points)
        Y = linspace(ylim...,num_points)

        Δr = max(X.step.hi,Y.step.hi)
        sizeGrid = num_points^2

        #allocation
        xy_grid = Vector{Tuple{Float64,Float64}}(sizeGrid)
        vecField_normed = Vector{Tuple{Float64,Float64}}(sizeGrid)
        normFactor = -Inf
        dx = zeros(Float64,2)

        for (i,x_) in enumerate(X)
            for (j,y_) in enumerate(Y)

                k = num_points*(i-1)+j
                xy_grid[k] = (x_,y_)
                vecField!(t0,xy_grid[k],dx)
                normFactor = max(normFactor, norm(dx))
                vecField_normed[k] = Tuple(arrowsize^2 * Δr * dx)

                if gridlog #sign(f) is because log(-a) not in Real
                    vecField_normed[k] = sign.(vecField_normed[k]).*log.(norm.(vecField_normed[k]).+1)
                end

            end
        end

        if gridlog
            normFactor = log(normFactor+1)
        end

        vecField_normed = [vecField_normed[i]./normFactor for i in eachindex(vecField_normed)]

        #Plotting
        quiver(xy_grid, quiver=vecField_normed)
        xlims!(xlim...)
        ylims!(ylim...)
    end

    function stability_matrix(eqs_diff!::Function,δξ,x0,t0=0.0)
        x0 = x0 .+ δξ
        dx = zeros(x0)
        eqs_diff!(t0,x0,dx)

        dof = length(x0)
        A = zeros(typeof(x0[1]), dof, dof)

        for i in 1:dof
            for j in 1:dof
                A[i,j] = derivative(dx[i],j)
            end
        end

        return A
    end

    simplecticity(ϕ::Vector{T}, ω) where T<:Number = jacobian(ϕ)' * ω * jacobian(ϕ) - ω
    function simplecticity(ϕ::Matrix{T}) where T<:Number

        time_units, D = size(ϕ)
        z_matrix = zeros(D÷2,D÷2)
        ω = hcat(z_matrix, I)
        tmp = hcat(-I, z_matrix)
        ω = vcat(ω, tmp)

        symp = zeros(time_units)
        for t in 1:time_units
            symp[t] = sum(abs,simplecticity(ϕ[t,:],ω))
        end
        return symp
    end

    #### ODE systems ####
    function harmonic_oscillator!(t,x,dx)
        dx[1] = -x[2]
        dx[2] = x[1]
        nothing
    end

    function simple_pendulum!(t,x,dx)
        θ,p = (x[1],x[2])
        dx[1] = p/(m*l^2)
        dx[2] = -m*g*l*sin(θ)
        nothing
    end
    m,g,l = 1.0,1.0,1.0

    function artificial_ode!(t,x,dx)
        dx[1] = 2.0*x[1]*x[2]
        dx[2] = -x[2]^2 + x[1]
        nothing
    end

    ### Other functions... ###
    ### Taken from Sébastien Villemot at http://www.econforge.org/posts/2014/juil./28/cef2014-julia/ ###
    function hp_filter(y::Vector{Float64}, lambda::Int64)
        n = length(y)
        @assert n >= 4

        diag2 = lambda*ones(n-2)
        diag1 = [ -2lambda; -4lambda*ones(n-3); -2lambda ]
        diag0 = [ 1+lambda; 1+5lambda; (1+6lambda)*ones(n-4); 1+5lambda; 1+lambda ]

        D = spdiagm((diag2, diag1, diag0, diag1, diag2), (-2,-1,0,1,2))

        D\y, y - D\y
    end

    function hp_filter(df::T, lambda::Int64) where T <: DataFrame
        num_cols = size(df)[2]
        cycle = Vector{Any}(num_cols) # Vector{Float64}
        trend = Vector{Any}(num_cols)

        i = 1
        for ts in df.columns
            tmp = hp_filter(ts,lambda)
            trend[i] = tmp[1]
            cycle[i] = tmp[2]
            i+=1
        end

        trend = DataFrame(trend,df.colindex)
        cycle = DataFrame(cycle,df.colindex)
        return trend, cycle
    end

end #module

max_stars_repo_name: blas-ko/tesis

max_issues_repo_name: blas-ko/tesis

max_forks_repo_name: blas-ko/tesis


text:
using Rust
using Base.Test
using Base.Libdl

# write your own tests here
@testset "Compile code string" begin
    rustlib = compile_code("""
        extern crate libc;
        use libc::uint32_t;

        #[no_mangle]
        pub extern fn addition(a: uint32_t, b: uint32_t) -> uint32_t {
           a + b
        }
    """)

    addition = dlsym(rustlib, :addition)
    @test addition != C_NULL
    @test UInt32(UInt32(2) + UInt32(10)) === ccall(addition, UInt32, (UInt32, UInt32), 2, 10)
end

@testset "Compile Julia function" begin
    @testset "Emission" begin
        function addition(a::UInt32, b::UInt32)
            a + b
        end

        name, emission = Rust.AST.generate_rust(addition, (UInt32, UInt32))
        @test isa(emission, String)
        @test name == "addition"
    end

    @testset "Compilation" begin
        function addition(a::UInt32, b::UInt32)
            a + b
        end

        addition = compile_function(addition, (UInt32, UInt32))

        @test addition != C_NULL
        @test UInt32(UInt32(2) + UInt32(10)) === ccall(addition, UInt32, (UInt32, UInt32), 2, 10)
    end

    @testset "Every function" begin
        function everything(a::UInt32, b::UInt32)
            (b >> a * b + div(b, a) * (a % b) - a | b) + (a & b + (~b $ a))
        end

        everything_rust = compile_function(everything, (UInt32, UInt32))

        @test everything_rust != C_NULL
        @test everything(UInt32(2), UInt32(10)) === ccall(everything_rust, UInt32, (UInt32, UInt32), 2, 10)
    end
end


text:
using Distributions

# walker type, with a position and if it is a wanderer
#  mutable struct walker
#      x::Float64
#      y::Float64
#      loner::Bool
#  end

"""
    distance(walker1,walker2)

Computes de distance between two walkers.
"""
function distance(w1::walker,w2::walker)
    on_x = (w1.x - w2.x)^2
    on_y = (w1.y - w2.y)^2
    sqrt(on_x+on_y)
end

"""
    freeze(walker1,walker2,tol)

Changes the loner state of a walker if it's distance
to another one is < tol
"""
function freeze(w1::walker,w2::walker,tol::Float64)
    if w1.loner == false || w2.loner == false
        if distance(w1,w2) <= tol
            return w1.loner=w2.loner=0
        end
    end
end

"""
    look_around(element,W,tol)

Computes the distance of a walker to all the other walkers
in an array of them and freezes them if their distance to
the nearest one is < tol
"""
function look_around(element::Int64,W::Array,tol::Float64)
    not_me = W[1:end .!= element]
    map(x -> freeze(the_walkers[element],x,tol),not_me)
end

"""
    random_walk(w,min,max)

Makes a walker do a random walk inside a restriction box with
length Δ
"""
function random_walk(w::walker,min::Number,max::Number,Δ=0.1)
    Δx = rand([-Δ,0,Δ])
    Δy = rand([-Δ,0,Δ])
    if w.x + Δx > min && w.x + Δx < max
        w.x = w.x + Δx;
    end
    if w.y + Δy > min && w.y + Δy < max
        w.y = w.y + Δy;
    end
end

"""
    walker_set(count,min_box,max_box)

Gives a set of walkers  bounded in a box
"""
function walker_set(count::Int64,min_box,max_box)
    the_walkers = [walker(round(rand(Uniform(min_box,max_box));digits=1),
                          round(rand(Uniform(min_box,max_box));digits=1),
                          1) for i in 1:count]
    return the_walkers
end

"""
dla(walker_array,tol,min_box,max_box;max_no_mov)

Performs a DLA on a set of walkers in a bounding box,
makes them stop moving if their distance < tol to a
non-moving walker. Stops when you get the maximum of
non-moving ones of desire
"""
function dla(walker_array::Array,tol,min_box,max_box;max_no_mov=0)
    while sum([x.loner for x ∈ walker_array]) > max_no_mov
        global the_movers = findall(x -> x.loner == true, walker_array);
        for i ∈ the_movers
            look_around(i,walker_array,tol)
        end
        the_movers = findall(x -> x.loner == true, the_walkers);
        map(x -> random_walk(walker_array[x],min_box,max_box),the_movers);
    end
    walker_array
end

## Main program
#  the_walkers = [walker(round(rand(Uniform(-5,5));digits=1),
#                        round(rand(Uniform(-5,5));digits=1),
#                        1) for i in 1:5000]
#  unshift!(the_walkers,walker(0.,0.,0))
#
#  while sum([x.loner for x ∈ the_walkers]) > 0
#      global the_movers = findall(x -> x.loner == true, the_walkers);;
#      for i ∈ the_movers
#          look_around(i,the_walkers,0.1)
#      end
#      the_movers = findall(x -> x.loner == true, the_walkers);
#      map(x -> random_walk(the_walkers[x],-5,5),the_movers);
#  end
#
#  ## Plotting
#  the_x = [w.x for w in the_walkers]
#  the_y = [w.y for w in the_walkers]
#  p = scatter(the_x,the_y,markersize=3.,
#              leg=false,ticks=nothing,
#              markerstrokewidth=0.,border=false)
#  png(p,"/figs/initial")


text:
# Load DataFrame package
require("DataFrames")
using DataFrames

# Load a CSV file into a DataFrame
df = read_table(Pkg.dir("DataFrames", "demo", "toy_example.csv"))

# Look at the data structure
# dump() is like R's str()
dump(df)

# Look at the internal data structure
idump(df)

# Print out a summary of each column
describe(df)

# head/tail of a DataFrame
head(df)
head(df, 3)
tail(df, 2)

# Select all rows where column A is greater than 4.0
# Element-wise operators in Julia usually have a leading "."
df[:( A .> 4.0 ), :]

# Make a new column using within
df2 = within(df, :( D = A + C ))

# This is similar, but now changes apply directly to df
within!(df, quote
    D = A + C
    E = A + sum(C)
end)

dump(df)

# Create a new DataFrame based on operations on another DataFrame
# This is similar to plyr's summarise()
df3 = based_on(df, quote
    ct = cut(replaceNA(A, 0.0), 3)
    sum_A = sum(A)
end)

# cut makes a PooledDataVector that is like R's factors, but
# PooledDataVectors can contain more than just strings. Here's
# the internal structure of a PooledDataVector:
idump(df3["ct"])

# In DataFrame, copies of data are minimized, especially for column
# operations
# These are both the same entity; change one, and you change the other:
df2 = df

colA = df2["A"]
colA[1] = 99.0   # This changes df and df2, too.
df


text:
include("./src/chiEFTint.jl")
using .chiEFTint

#make_chiEFTint(;is_plot=true)
make_chiEFTint()

## for profiling
#using StatProfilerHTML
#@profilehtml make_chiEFTint()


text:
using Test
using CryptoGroups
import CryptoGroups: F2PB, F2GNB, FP, Reducer, Weierstrass, AffinePoint, double, Reducer, BinaryCurve, oncurve, specialize

eq = Weierstrass{1, 1}

x1 = FP{23}(3)
y1 = FP{23}(10)

p1 = AffinePoint{eq}(x1, y1)

x2 = FP{23}(9)
y2 = FP{23}(7)

p2 = AffinePoint{eq}(x2, y2)

@test p1 + p2 == AffinePoint{eq, FP{23}}(17, 20)
@test double(p1) == AffinePoint{eq, FP{23}}(7, 12)
@test oncurve(p1) == true
@test oncurve(p2) == true

f = reverse(bin"10011")
R = Reducer(f)

α = F2PB{R}(bin"0100")

a = α^4
b = one(α)

#eq = BEQ(a, b)
eq = specialize(BinaryCurve, a, b)

p1 = AffinePoint{eq}(α^6, α^8)
p2 = AffinePoint{eq}(α^3, α^13)

@test p1 + p2 == AffinePoint{eq}(one(α), α^13)
@test double(p1) == AffinePoint{eq}(α^10, α^8)

@test oncurve(p1) == true
@test oncurve(p2) == true

α = F2GNB{4, 1}(bin"1100")

#eq = BEQ(zero(α), α^3)

eq = specialize(BinaryCurve, zero(α), α^3)

#error("here")

G = AffinePoint{eq}(α^3, α^5)

@test double(G) == AffinePoint{eq}(α^4, α^3)
@test double(G) + G == AffinePoint{eq}(α^13, α^2)

@test oncurve(G) == true
@test oncurve(G*3) == true

### Experimenting with 𝔽₂ field in Weierstrass equation. The doubling operation though does not work and thus is type constrained.

α = F2GNB{4, 1}(bin"1100")

#eq = WeierstrassEQ(one(α), one(α))
eq = specialize(Weierstrass, one(α), one(α))

p1 = AffinePoint{eq}(α^5, α^10)
p2 = AffinePoint{eq}(α^14, α^5)

@test oncurve(p1) == true
@test oncurve(p2) == true

@test oncurve(p1 + p2) == true

# Testing point multiplication
eq = Weierstrass{1, 1}

p = AffinePoint{eq, FP{23}}(3, 10)

@test oncurve(3*p) == true
@test double(p) == 2p
@test double(p) + p == 3p

@test double(p) + p + double(p) == 5p
@test double(double(p)) + p == 5p

############################### Real curve examples #####################

##################### P_192 ####
let

    p = 6277101735386680763835789423207666416083908700390324961279
    n = 6277101735386680763835789423176059013767194773182842284081

    a = p - 3
    b = parse(BigInt, "64210519 e59c80e7 0fa7e9ab 72243049 feb8deec c146b9b1", base=16)

    Gx = parse(BigInt, "188da80e b03090f6 7cbf20eb 43a18800 f4ff0afd 82ff1012", base=16)
    Gy = parse(BigInt, "07192b95 ffc8da78 631011ed 6b24cdd5 73f977a1 1e794811", base=16)

    #point = AffinePoint{WeierstrassEQ(a, b), (FP/p)}(Gx, Gy)
    point = AffinePoint{specialize(Weierstrass, a, b), specialize(FP, p)}(Gx, Gy)

    @test oncurve(point) == true
    @test point * n == zero(point)

end

######################  K_163 #################

import CryptoGroups: hex2bits, <|

let

    a = 1
    b = 1

    ########### Polynomial basis

    #R = Reducer([163, 7, 6, 3, 0])
    #F = F2PB{R}

    F = specialize(F2PB, [163, 7, 6, 3, 0])

    EQ = specialize(BinaryCurve, F <| a, F <| b)

    Gx = "2 fe13c053 7bbc11ac aa07d793 de4e6d5e 5c94eee8"
    Gy = "2 89070fb0 5d38ff58 321f2e80 0536d538 ccdaa3d9"

    point = AffinePoint{EQ, F}(Gx, Gy)
    @test oncurve(point)

    ######### Normal Basis

    T = 4

    Gx = "0 5679b353 caa46825 fea2d371 3ba450da 0c2a4541"
    Gy = "2 35b7c671 00506899 06bac3d9 dec76a83 5591edb2"

    F = F2GNB{163, 4}
    eq = specialize(BinaryCurve, F <| a, F <| b)

    point = AffinePoint{eq, F}(Gx, Gy)

    @test oncurve(point)
    @test oncurve(point*3)

end

max_stars_repo_name: PeaceFounder/DiffieHellmanGroups

max_issues_repo_name: PeaceFounder/DiffieHellmanGroups

max_forks_repo_name: PeaceFounder/DiffieHellmanGroups


text:
__precompile__(true)

module MapImages

export MapImage

import Parameters: @with_kw, @unpack
import Sario
import Sario: DemRsc
import Base: length, size, similar, step, parent, getindex, setindex!

# TODO: tests

# Note: modelling off OffsetArrays
# https://github.com/JuliaArrays/OffsetArrays.jl/blob/master/src/OffsetArrays.jl
# mutable struct MapImage{T, N, AA<:AbstractArray} <: AbstractArray{T, N}
@with_kw mutable struct MapImage{T, N} <: AbstractArray{T, N}
    image::AbstractArray{T, N}
    demrsc::DemRsc
    # filename::Union{Nothing, AbstractString}  # Do I care to store this??
end
MapImage(A::AbstractArray{T,N}, demrsc::DemRsc) where {T,N} = MapImage{T,N}(A, demrsc)

# Extra functions manipulating MapImages
include("./latlon.jl")
include("./gps.jl")
include("./geotiff.jl")

# TODO: maybe move this to Sario
function _get_dem(filename)
    if Sario._is_h5(filename) && "dem_rsc" in names(filename)
        demrsc = Sario.load_dem_from_h5(filename)
    else
        demrscfile = Sario.find_rsc_file(filename)
        demrsc = Sario.load(demrscfile)
    end
    return demrsc
end

# TODO: add permutes to pass through to load
function MapImage(filename::AbstractString; dset_name::AbstractString="")
    demrsc = _get_dem(filename)
    return MapImage(Sario.load(filename, dset_name=dset_name), demrsc)
end

Base.eachindex(::IndexCartesian, A::MapImage) = CartesianIndices(axes(A))

Base.parent(A::MapImage) = A.image

# Base.parent(A::MapImage) = A.image  # TODO: do we need?
Base.size(A::MapImage) = size(A.image)
Base.size(A::MapImage, d) = size(A.image, d)

# Taken from https://docs.julialang.org/en/v1/manual/interfaces/ and OffsetArray examples
Base.similar(A::MapImage) = MapImage(similar(A.image), copy(A.demrsc))
Base.similar(A::MapImage, ::Type{S}) where S = MapImage(similar(A.image, S, size(A)), copy(A.demrsc))
Base.similar(A::MapImage, dims::Dims) = MapImage(similar(A.image, eltype(A.image), size(A)),copy(A.demrsc))

function Base.similar(A::MapImage, ::Type{T}, dims::Dims) where T
    B = similar(A.image, T, dims)
    return MapImage(B, copy(A.demrsc))
end

# Allow broadasting, while keeping the output as a MapImage (instead of Array)
# Note: Example taken directly from here:
# https://docs.julialang.org/en/v1/manual/interfaces/#man-interfaces-broadcasting-1
Base.BroadcastStyle(::Type{<:MapImage}) = Broadcast.ArrayStyle{MapImage}()

function Base.similar(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{MapImage}}, ::Type{ElType}) where ElType
    # Scan the inputs for the MapImage:
    A = find_obj(bc)
    # Use the DemRsc field of A to create the output
    MapImage(similar(Array{ElType}, axes(bc)), copy(A.demrsc))
end

"`A = find_obj(As)` returns the first MapImage among the arguments."
find_obj(bc::Base.Broadcast.Broadcasted) = find_obj(bc.args)
find_obj(args::Tuple) = find_obj(find_obj(args[1]), Base.tail(args))
find_obj(x) = x
find_obj(a::MapImage, rest) = a
find_obj(::Any, rest) = find_obj(rest)

# Catchall: just getindex on image, ignore demrsc
Base.getindex(A::MapImage{T,N}, I::Vararg{Any,N}) where {T,N} = A.image[I...]

ColonOrRange = Union{Colon, <:AbstractRange{Int}}
# For DemRsc adjustment purposes
Base.step(x::Colon) = 1
start(x::ColonOrRange) = typeof(x) <: AbstractRange ? x.start : 1
Base.length(x::ColonOrRange) = typeof(x) <: AbstractRange ? Base.length(x) : nothing

# 2D array subarray: handles ranges and colons
function Base.getindex(A::MapImage{T,2}, I::Vararg{ColonOrRange, 2}) where {T}
    subimg = A.image[I...]

    # DemRsc adjusting:
    rowrange, colrange = I
    newdemrsc = _new_dem_data(A.demrsc, rowrange, colrange)

    return MapImage(subimg, newdemrsc)
end

# TODO: 3D

# 2D array single floats + colons (for lat/lon)
# First: A[lat, lon]
function Base.getindex(A::MapImage{T,2}, I::Vararg{<:Real, 2}) where {T}
    row = _get_row(A.demrsc, I[1])
    col = _get_col(A.demrsc, I[2])

    subimg = A.image[row, col]

    # We assume this is no longer a value image since it's not a range:
    # just return image data
    return subimg
end

# 2D array subarray: For float ranges, we want a tuple of floats
# This is because doing something like 30.1:33.4 won't tell you what the
# final number was.
# The StepRangeLen gives the start, len, offset, and step, so you will lose what
# the end of the range was
# E.G.
#   A[(30.1, 32.2), (-104.1, :)]
function Base.getindex(A::MapImage{T,2}, I::Vararg{Tuple{<:Real, <:Real}, 2}) where {T}
    lats, lons = I

    rowrange, colrange = latlon_to_rowcol(A.demrsc, lats, lons)
    subimg = A.image[rowrange, colrange]

    # DemRsc adjusting:
    newdemrsc = _new_dem_data(A.demrsc, rowrange, colrange)

    return MapImage(subimg, newdemrsc)
end

# Making changes to the DemRsc: prep incoming data
function _new_dem_data(d::DemRsc, rowrange::ColonOrRange, colrange::ColonOrRange)
    new_rows = isnothing(length(rowrange)) ? d.rows : length(rowrange)
    new_cols = isnothing(length(colrange)) ? d.cols : length(colrange)
    row_step = step(rowrange)
    col_step = step(colrange)
    row_start = start(rowrange)
    col_start = start(colrange)
    return crop_rsc_data(d, new_rows, new_cols,
                         row_start, col_start,
                         row_step, col_step)
end

"""Adjusts the old demrsc of a MapImage after slicing

Takes the 'rows' and 'cols' from demrsc
and adjusts for the smaller size with a new dict

The row/col_start are the indices where we have sliced:
    e.g. img[2:3:11, :] would have
    row_start = 2, row_step = 3, new_rows = 4
"""
function crop_rsc_data(demrsc::DemRsc,
                       new_rows=demrsc.rows, new_cols=demrsc.cols,
                       row_start=1, col_start=1,
                       row_step=1, col_step=1)
    @unpack rows, cols, y_step, x_step, y_first, x_first = demrsc

    # Move forward the starting row/col from where it used to be
    new_x_first = x_first + x_step * (col_start - 1)
    new_y_first = y_first + y_step * (row_start - 1)

    # now adjust step sizes (say, if we take only every other row)
    new_x_step = x_step * col_step
    new_y_step = y_step * row_step

    return DemRsc(demrsc, width=new_cols, cols=new_cols,
                  file_length=new_rows, rows=new_rows,
                  x_first=new_x_first, y_first=new_y_first,
                  x_step=new_x_step, y_step=new_y_step)
end

# function Base.setindex!(A::MapImage{T,N}, val, I::Vararg{Int,N}) where {T,N}
function Base.setindex!(A::MapImage{T,N}, val, I::Vararg{Any,N}) where {T,N}
    # setindex!(A, X, inds...)  corresponds to A[inds...] = X
    A.image[I...] = val
    return val
end

#TODO:
#https://docs.julialang.org/en/v1/base/io-network/index.html#Base.show-Tuple{Any,Any,Any}
# MIME"image/png"
# show(io, ::MIME"image/png", x::MapImage) = .

end # module


text:
using Manopt, Manifolds, Test

@testset "Bezier Tests" begin
    @testset "General Bezier Tests" begin
        repr(BezierSegment([[0.0, 0.0], [0.0, 0.0]])) ==
        "BezierSegment([[0.0, 0.0], [0.0, 0.0]])"
    end
    @testset "Spherical Test" begin
        M = Sphere(2)
        pC = [0.0, 1.0, 0.0]
        pT = exp(M, pC, [0.0, 0.0, 0.7])
        pB = exp(M, pC, [0.0, 0.0, -0.7])
        B = [
            BezierSegment(shortest_geodesic(M, pT, pC, [0.0, 1 / 3, 2 / 3, 1.0])),
            BezierSegment(shortest_geodesic(M, pC, pB, [0.0, 1 / 3, 2 / 3, 1.0])),
        ]
        # this is equispaced, so the pure cost is zero and the gradient is a zero-vector
        t = collect(range(0.0, 1.0; length=5))
        pts = shortest_geodesic(M, pT, pB, t)
        pts2 = de_casteljau(M, B, 2 .* t)
        @test sum(distance.(Ref(M), pts, pts2)) ≈ 0
        aX = log(M, pT, pC)
        aT1 = adjoint_differential_bezier_control(M, BezierSegment([pT, pC]), 0.5, aX).pts
        aT2 = [
            adjoint_differential_geodesic_startpoint(M, pT, pC, 0.5, aX),
            adjoint_differential_geodesic_endpoint(M, pT, pC, 0.5, aX),
        ]
        @test aT1 ≈ aT2
        #
        @test sum(norm.(
            ∇acceleration_bezier(M, B[1], collect(range(0.0, 1.0; length=20))).pts .-
            [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]],
        )) ≈ 0 atol = 10^(-12)

        # cost and gradient
        T = collect(range(0.0, 2.0; length=51))
        degrees = get_bezier_degrees(M, B)
        Bvec = get_bezier_points(M, B, :differentiable)
        Mp = PowerManifold(M, NestedPowerRepresentation(), length(Bvec))
        @test cost_acceleration_bezier(M, Bvec, degrees, T) ≈ 0 atol = 10^-10
        z = zero_tangent_vector(Mp, Bvec)
        distance(Mp, ∇acceleration_bezier(M, Bvec, degrees, T), z)
        @test norm(Mp, Bvec, ∇acceleration_bezier(M, Bvec, degrees, T) - z) ≈ 0 atol =
            10^(-12)

        d = [pT, exp(M, pC, [0.3, 0.0, 0.0]), pB]
        λ = 3.0

        # cost and gradient with data term
        @test cost_L2_acceleration_bezier(M, Bvec, degrees, T, λ, [pT, pC, pB]) ≈ 0 atol =
            10^(-10)
        @test cost_L2_acceleration_bezier(M, Bvec, degrees, T, λ, d) ≈
              λ / 2 * distance(M, d[2], pC) .^ 2
        # when the data are the junctions
        @test norm(
            Mp, Bvec, ∇L2_acceleration_bezier(M, Bvec, degrees, T, λ, [pT, pC, pB]) - z
        ) ≈ 0 atol = 10^(-12)
        z[4][1] = -0.9
        @test norm(Mp, Bvec, ∇L2_acceleration_bezier(M, Bvec, degrees, T, λ, d) - z) ≈ 0 atol =
            10^(-12)
        # when the data is weighted with zero
        @test cost_L2_acceleration_bezier(M, Bvec, degrees, T, 0.0, d) ≈ 0 atol = 10^(-10)
        z[4][1] = 0.0
        @test norm(Mp, Bvec, ∇L2_acceleration_bezier(M, Bvec, degrees, T, 0.0, d) - z) ≈ 0 atol =
            10^(-12)
    end
    @testset "de Casteljau variants" begin
        M = Sphere(2)
        B = artificial_S2_composite_bezier_curve()
        b = B[2]
        b2s = BezierSegment([b.pts[1], b.pts[end]])
        # (a) 2 points -> geo special case
        f1 = de_casteljau(M, b2s) # fct -> recursive
        pts1 = f1.([0.0, 0.5, 1.0])
        pts2 = de_casteljau(M, b2s, [0.0, 0.5, 1.0])
        @test pts1 ≈ pts2
        # (b) one segment
        f2 = de_casteljau(M, b) # fct -> recursive
        pts3 = f2.([0.0, 0.5, 1.0])
        pts4 = de_casteljau(M, b, [0.0, 0.5, 1.0])
        @test pts3 ≈ pts4
        # (c) whole composites
        f3 = de_casteljau(M, B) # fct -> recursive
        pts5 = f3.([0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0])
        pts6 = de_casteljau(M, B, [0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0])
        @test pts5 ≈ pts6
        @test_throws DomainError f3(-0.1)
        @test_throws DomainError de_casteljau(M, B, -0.1)
        @test_throws DomainError f3(3.5)
        @test_throws DomainError de_casteljau(M, B, 3.5)
    end
    @testset "Spherical Data" begin
        M = Sphere(2)
        B = artificial_S2_composite_bezier_curve()
        @test de_casteljau(M, B, [0.0, 1.0, 2.0, 3.0]) ≈
              [B[1].pts[1], B[2].pts[1], B[3].pts[1], B[3].pts[4]]
        @test get_bezier_junction_tangent_vectors(M, B) ≈ [
            log(M, B[1].pts[1], B[1].pts[2]),
            log(M, B[1].pts[4], B[1].pts[3]),
            log(M, B[2].pts[1], B[2].pts[2]),
            log(M, B[2].pts[4], B[2].pts[3]),
            log(M, B[3].pts[1], B[3].pts[2]),
            log(M, B[3].pts[4], B[3].pts[3]),
        ]
        @test get_bezier_junction_tangent_vectors(M, B[1]) ≈
              [log(M, B[1].pts[1], B[1].pts[2]), log(M, B[1].pts[4], B[1].pts[3])]
        @test get_bezier_junctions(M, B[1]) == B[1].pts[[1, end]]
        @test get_bezier_inner_points(M, B) ==
              [B[1].pts[2], B[1].pts[3], B[2].pts[2], B[2].pts[3], B[3].pts[2], B[3].pts[3]]
        @test get_bezier_inner_points(M, B[1]) == [B[1].pts[2], B[1].pts[3]]

        @test get_bezier_points(M, B) == cat([[b.pts...] for b in B]...; dims=1)
        @test get_bezier_points(M, B, :continuous) ==
              cat([[b.pts[[1:3]...]...] for b in B]..., [B[3].pts[4]]; dims=1)
        @test get_bezier_points(M, B, :differentiable) ==
              cat([B[1].pts[[1, 2]]...], [b.pts[[3, 4]] for b in B]...; dims=1)
        @test get_bezier_points(M, B[1]) == B[1].pts
        # for segments just check that they
        d = get_bezier_degrees(M, B)
        A = get_bezier_segments(M, get_bezier_points(M, B), d)
        @test [A[i].pts for i in 1:3] == [B[i].pts for i in 1:3]
        A = get_bezier_segments(M, get_bezier_points(M, B, :continuous), d, :continuous)
        @test [A[i].pts for i in 1:3] == [B[i].pts for i in 1:3]
        A = get_bezier_segments(
            M, get_bezier_points(M, B, :differentiable), d, :differentiable
        )
        @test [A[i].pts for i in 1:3] == [B[i].pts for i in 1:3]

        # out of range
        @test_throws ErrorException adjoint_differential_bezier_control(
            M, B, 7.0, zero_tangent_vector(M, B[1].pts[1])
        )
        # a shortcut to evaluate the adjoint at several points is equal to seperate evals
        b = B[2]
        @test isapprox(
            adjoint_differential_bezier_control(
                M, b, [0.0, 1.0], [log(M, b.pts[1], b.pts[2]), -log(M, b.pts[4], b.pts[3])]
            ).pts,
            adjoint_differential_bezier_control(M, b, 0.0, log(M, b.pts[1], b.pts[2])).pts +
            adjoint_differential_bezier_control(M, b, 1.0, -log(M, b.pts[4], b.pts[3])).pts,
        )
        # differential
        X = BezierSegment([
            log(M, b.pts[1], b.pts[2]), [zero_tangent_vector(M, b.pts[i]) for i in 2:4]...
        ])
        @test differential_bezier_control(M, b, 0.0, X) ≈ X.pts[1]
        dT1 = differential_bezier_control.(Ref(M), Ref(b), [0.0, 1.0], Ref(X))
        dT2 = differential_bezier_control(M, b, [0.0, 1.0], X)
        @test dT1 ≈ dT2
        X2 = [
            BezierSegment([[0.0, 0.0, 0.0] for i in 1:4]),
            X,
            BezierSegment([[0.0, 0.0, 0.0] for i in 1:4]),
        ]
        @test_throws DomainError differential_bezier_control(M, B, 20.0, X2)
        dbT2a = differential_bezier_control(M, B, 1.0, X2)
        @test dbT2a ≈ X.pts[1]
        dbT2 = differential_bezier_control(M, B, [1.0, 2.0], X2)
        dbT1 = differential_bezier_control.(Ref(M), Ref(B), [1.0, 2.0], Ref(X2))
        @test dT1 ≈ dbT1
        @test dbT2 ≈ dbT1
    end
end


text:
######################################
## define getEntries, double index ##
######################################

function getEntries(td::TimeData.AbstractTimedata,
                    rowInds::Array{Int}, colInds::Array{Int})

    nObs = length(rowInds)
    nams = names(td)

    ## preallocation
    idxs = Array(typeof(td.idx[1]), nObs)
    variables = Array(Symbol, nObs)
    values = Array(Any, nObs)

    for ii=1:nObs
        values[ii] = get(td, rowInds[ii], colInds[ii])
        idxs[ii] = TimeData.idx(td)[rowInds[ii]]
        variables[ii] = nams[colInds[ii]]
    end
    df = DataFrame(variable = variables, value = values)
    return TimeData.Timedata(df, idxs)
end

######################################
## define getEntries, single index ##
######################################

function getEntries(td::TimeData.AbstractTimedata,
                    singleInd::Array{Int})
    (rowInds, colInds) = ind2sub(size(td), singleInd)
    return TimeData.getEntries(td, rowInds, colInds)
end

###################################
## getEntries, logical indexing ##
###################################

function getEntries(td::TimeData.AbstractTimedata,
                    logInds::Array{Bool, 2},
                    sortDates::Bool=false)

    if size(td) != size(logInds)
        error("TimeData and logical index array must be of equal
dimensions for logical indexing")
    end

    if sortDates
        singleInds = find(logInds')
        (colInds, rowInds) = ind2sub(size(logInds'), singleInds)
    else
        singleInds = find(logInds)
        (rowInds, colInds) = ind2sub(size(td), singleInds)
    end
    return getEntries(td, rowInds, colInds)
end

#######################################################
## getEntries, logical indexing with Timedata object ##
#######################################################

function getEntries(td::TimeData.AbstractTimedata,
                    td2::TimeData.AbstractTimedata;
                    getNAs::Bool=false,
                    sortDates::Bool=false)

    if size(td) != size(td2)
        error("TimeData objects must be of equal dimensions for
logical indexing")
    end

    logInds = asArr(td2, Bool, getNAs)
    return getEntries(td, logInds, sortDates)
end


text:
include("../utils/functions.jl")
include("../optimizers/activeset.jl")

using Gadfly

## Toy example

A = [1 1   1 0;
     2 1/2 0 1]
b = Vector{Float64}([5, 8])
c = Vector{Float64}([-4, -2, 0, 0])

#assts = [true,false,true,false]
assts = [false,false,true,true]

#x = [4.,0.,1.,0.]
x = [0.,0.,5.,8.]
run_simplex(assts,x,A,b,c,2)
#simplex_iteration(assts, x, A, b, c)

# phase I of simplex
x0 = [0.,0.,0.,0.,5.,8.] # last 2 entries are z1, z2
A0 = [1 1   1 0    1 0;
     2 1/2 0 1    0 1]
c0 = [0., 0., 0., 0., 1., 1.]

assts0 = [false,false,false,false,true,true]

assts_start,x_start = run_simplex(assts0,x0,A0,b,c0,5)

function run_program(m, num_iters=100)
    n = 3*m
    # generate random values for A
    A = randn(m,n)
    b = randn(m) + 1
    c = randn(n) * 2 + 2

    return run_two_stage_simplex(A, b, c, num_iters)
end

function run_experiment(num_rounds, num_iters)
    times = []
    for m in range(10,10,num_rounds)
        x, assts, time = run_program(m, num_iters)
        push!(times,time)
    end
    return times
end

times = run_experiment(5,100)

num_rounds = 20

length(times)

#Gadfly.plot(x=range(10,10,num_rounds), y=times, Guide.XLabel("m"),
#    Guide.YLabel("time (s)"),)