text:
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!(∇f, x)
    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

function f∇f(∇f, x)
    if !(∇f == nothing)
        gx = similar(x)
        return f∇f!(gx, x)
    else
        return f∇f!(∇f, x)
    end
end
function f∇fs(∇f, x)
    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]
        return f(x), ∇f
    else
        return f(x)
    end
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)

shortname(::GradientDescent) = "GD  "
shortname(::Inverse) = "(inverse)"
shortname(::Direct) = " (direct)"

function printed_minimize(f∇f, x0, method)
    res = minimize(f∇f, x0, method, OptOptions())
    print("NN  $(shortname(method)) $(shortname(method.approx)): ")
    @printf("%2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
end
function printed_minimize!(obj_inplace, x0, method)
    res = minimize!(obj_inplace, x0, method, OptOptions())
    print("NN! $(shortname(method)) $(shortname(method.approx)): ")
    @printf("%2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
end

printed_minimize(obj_inplace, x0, GradientDescent(Inverse()))
printed_minimize!(obj_inplace, copy(x0), GradientDescent(Inverse()))
res = minimize(obj_static, x0s, GradientDescent(Inverse()))
@printf("NN  GD(S) (inverse): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])

printed_minimize(obj_inplace, x0, GradientDescent(Direct()))
printed_minimize!(obj_inplace, copy(x0), GradientDescent(Direct()))
res = minimize(obj_static, x0s, GradientDescent(Direct()))
@printf("NN  GD(S)  (direct): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
println()

res = minimize(obj_inplace, x0, BFGS(Inverse()), OptOptions())
@test res[4] == 30
@printf("NN  BFGS    (inverse): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
res = minimize!(obj_inplace, copy(x0), (BFGS(Inverse()), Backtracking()), OptOptions())
@test res[4] == 30
@printf("NN! BFGS    (inverse): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
res = minimize!(obj_inplace, copy(x0), (BFGS(Inverse()), Backtracking(interp=FFQuadInterp())), OptOptions())
@test res[4] == 32
@printf("NN! BFGS    (inverse, quad): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
res = minimize(obj_static, x0s, BFGS(Inverse()))
@test res[4] == 30
@printf("NN  BFGS(S) (inverse): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
res = minimize(obj_static, x0s, (BFGS(Inverse()), Backtracking(interp=FFQuadInterp())))
@test res[4] == 32
@printf("NN  BFGS(S) (inverse, quad): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])

res = minimize(obj_inplace, x0, BFGS(Direct()), OptOptions())
@printf("NN  BFGS    (direct): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
res = minimize!(obj_inplace, copy(x0), BFGS(Direct()), OptOptions())
@printf("NN! BFGS    (direct): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
res = minimize(obj_static, x0s, BFGS(Direct()))
@printf("NN  BFGS(S) (direct): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
println()
# res = optimize(obj_inplace, copy(x0), Optim.BFGS(linesearch=Backtracking()))
# @time optimize(obj_inplace, copy(x0), Optim.BFGS(linesearch=Backtracking()))
# @printf("OT! BFGS (inverse): %2.2e  %2.2e\n", norm(Optim.minimizer(res)-xopt,Inf), Optim.g_residual(res))

res = minimize(obj_inplace, x0, DFP(Inverse()))
@printf("NN  DFP    (inverse): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
res = minimize!(obj_inplace, copy(x0), DFP(Inverse()))
@printf("NN! DFP    (inverse): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
res = minimize(obj_static, x0s, DFP(Inverse()))
@printf("NN  DFP(S) (inverse): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
res = minimize(obj_inplace, x0, DFP(Direct()))
@printf("NN  DFP    (direct): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
res = minimize!(obj_inplace, copy(x0), DFP(Direct()))
@printf("NN! DFP    (direct): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
res = minimize(obj_static, x0s, DFP(Direct()))
@printf("NN  DFP(S) (direct): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
println()

res = minimize(obj_inplace, x0, SR1(Inverse()))
@printf("NN  SR1     (inverse): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
res = minimize!(obj_inplace, copy(x0), SR1(Inverse()))
@printf("NN! SR1     (inverse): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
res = minimize(obj_static, x0s, SR1(Inverse()))
@printf("NN  SR1(S)  (inverse): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
res = minimize(obj_inplace, x0, SR1(Direct()))
@printf("NN  SR1     (direct): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
res = minimize!(obj_inplace, copy(x0), SR1(Direct()))
@printf("NN! SR1     (direct): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
res = minimize(obj_static, x0s, SR1(Direct()))
@printf("NN  SR1(S)  (direct): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
res = minimize!(obj_inplace, copy(x0), (SR1(Direct()), NWI()))
@printf("NN! SR1   (direct, NW): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])

xrand = rand(3)
xrands = SVector{3}(xrand)
println("\nFrom a random point: ", xrand)
res = minimize(obj_inplace, xrand, GradientDescent(Inverse()))
@printf("NN  GD   (inverse): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
res = minimize!(obj_inplace, copy(xrand), GradientDescent(Inverse()))
@printf("NN! GD   (inverse): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
res = minimize(obj_inplace, xrand, GradientDescent(Direct()))
@printf("NN  GD    (direct): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
res = minimize!(obj_inplace, copy(xrand), GradientDescent(Direct()))
@printf("NN! GD    (direct): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
println()
res = minimize(obj_inplace, xrand, BFGS(Inverse()))
@printf("NN  BFGS (inverse): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
res = minimize!(obj_inplace, copy(xrand), (BFGS(Inverse()), Backtracking()))
@printf("NN! BFGS (inverse): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
res = minimize!(obj_inplace, copy(xrand), (BFGS(Inverse()), HZAW()))
@printf("NN! BFGS (inverse) HZ: %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
res = minimize(obj_inplace, xrand, BFGS(Direct()))
@printf("NN  BFGS  (direct): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
res = minimize!(obj_inplace, copy(xrand), BFGS(Direct()))
@printf("NN! BFGS  (direct): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
res = minimize(obj_inplace, xrand, DBFGS(Direct()))
@printf("NN  DBFGS  (direct): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
res = minimize!(obj_inplace, copy(xrand), DBFGS(Direct()))
@printf("NN! DBFGS  (direct): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
res = minimize!(obj_inplace, copy(xrand), (DBFGS(Direct()), NWI()))
@printf("NN! DBFGS  (direct, NW): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
res = minimize!(obj_inplace, copy(xrand), (DBFGS(Direct()), Dogleg()))
@printf("NN! DBFGS  (direct, DL): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
println()

res = minimize(obj_inplace, xrand, DFP(Inverse()))
@printf("NN  DFP  (inverse): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
res = minimize!(obj_inplace, copy(xrand), DFP(Inverse()))
@printf("NN! DFP  (inverse): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
res = minimize(obj_inplace, xrand, DFP(Direct()))
@printf("NN  DFP   (direct): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
res = minimize!(obj_inplace, copy(xrand), DFP(Direct()))
@printf("NN! DFP   (direct): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
println()

res = minimize(obj_inplace, xrand, SR1(Inverse()))
@printf("NN  SR1  (inverse): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
res = minimize!(obj_inplace, copy(xrand), SR1(Inverse()))
@printf("NN! SR1  (inverse): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
res = minimize(obj_static, xrand, SR1(Inverse()))
@printf("NN  SR1(S)   (inverse): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
res = minimize(obj_inplace, xrand, SR1(Direct()))
@printf("NN  SR1   (direct): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
res = minimize!(obj_inplace, copy(xrand), SR1(Direct()))
@printf("NN! SR1   (direct): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
res = minimize(obj_static, xrand, SR1(Direct()))
@printf("NN  SR1(S)   (direct): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])

# res = minimize(obj_inplace, xrand, (SR1(Direct()), NWI()))
# @printf("NN  SR1   (direct, tr): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
res = minimize!(obj_inplace, copy(xrand), (SR1(Direct()), NWI()))
@printf("NN! SR1   (direct, NW): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
# res = minimize(obj_static, xrand, (SR1(Direct()), NWI()))
# @printf("NN  SR1(S)   (direct, tr): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])

println()

function himmelblau!(∇f, x)
    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
    return ∇f == nothing ? fx : (fx, ∇f)
end

function himmelblaus(∇f, x)
    f = (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]
        return f, @SVector([∇f1, ∇f2])
    else
        return f
    end
end

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

    return himmelblau!(g, x)
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)
res = minimize(him_outofplace, x0, GradientDescent(Inverse()))
@printf("NN  GD    (inverse): %2.2e  %d\n", norm(res[3], Inf), res[4])
res = minimize!(him_inplace, copy(x0), GradientDescent(Inverse()))
@printf("NN! GD    (inverse): %2.2e  %d\n", norm(res[3], Inf), res[4])
res = minimize(him_static, x0s, GradientDescent(Inverse()))
@printf("NN  GD(S) (inverse): %2.2e  %d\n", norm(res[3], Inf), res[4])
println()
res = minimize(him_outofplace, x0, GradientDescent(Direct()))
@printf("NN  GD    (direct): %2.2e  %d\n", norm(res[3], Inf), res[4])
res = minimize!(him_inplace, copy(x0), GradientDescent(Direct()))
@printf("NN! GD    (direct): %2.2e  %d\n", norm(res[3], Inf), res[4])
res = minimize(him_static, x0s, GradientDescent(Direct()))
@printf("NN  GD(S) (direct): %2.2e  %d\n", norm(res[3], Inf), res[4])
println()

res = minimize(him_outofplace, x0, BFGS(Inverse()))
@printf("NN  BFGS    (inverse): %2.2e  %d\n", norm(res[3], Inf), res[4])
res = minimize!(him_inplace, copy(x0), BFGS(Inverse()))
@printf("NN! BFGS    (inverse): %2.2e  %d\n", norm(res[3], Inf), res[4])
res = minimize(him_outofplace, x0s, BFGS(Inverse()))
@printf("NN  BFGS(S) (inverse): %2.2e  %d\n", norm(res[3], Inf), res[4])
println()
res = minimize(him_outofplace, x0, BFGS(Direct()))
@printf("NN  BFGS    (direct): %2.2e  %d\n", norm(res[3], Inf), res[4])
res = minimize!(him_inplace, copy(x0), BFGS(Direct()))
@printf("NN! BFGS    (direct): %2.2e  %d\n", norm(res[3], Inf), res[4])
res = minimize(him_outofplace, x0s, BFGS(Direct()))
@printf("NN  BFGS(S) (direct): %2.2e  %d\n", norm(res[3], Inf), res[4])
res = minimize!(him_inplace, copy(x0), (DBFGS(Direct()), NWI()))
@printf("NN! DBFGS  (direct, NW): %2.2e  %d\n", norm(res[3], Inf), res[4])
res = minimize!(him_inplace, copy(x0), (DBFGS(Direct()), Dogleg()))
@printf("NN! DBFGS  (direct, DL): %2.2e  %d\n", norm(res[3], Inf), res[4])

println()

res = minimize(him_outofplace, x0, DFP(Inverse()))
@printf("NN  DFP  (inverse): %2.2e  %d\n", norm(res[3], Inf), res[4])
res = minimize!(him_inplace, copy(x0), DFP(Inverse()))
@printf("NN! DFP  (inverse): %2.2e  %d\n", norm(res[3], Inf), res[4])
res = minimize(him_outofplace, x0s, DFP(Inverse()))
@printf("NN  DFP(S)  (inverse): %2.2e  %d\n", norm(res[3], Inf), res[4])
res = minimize(him_outofplace, x0, DFP(Direct()))
println()
@printf("NN  DFP   (direct): %2.2e  %d\n", norm(res[3], Inf), res[4])
res = minimize!(him_inplace, copy(x0), DFP(Direct()))
@printf("NN! DFP   (direct): %2.2e  %d\n", norm(res[3], Inf), res[4])
res = minimize(him_outofplace, x0s, DFP(Direct()))
@printf("NN  DFP(S)   (direct): %2.2e  %d\n", norm(res[3], Inf), res[4])
println()

res = minimize(him_outofplace, x0, SR1(Inverse()))
@printf("NN  SR1  (inverse): %2.2e  %d\n", norm(res[3], Inf), res[4])
res = minimize!(him_inplace, copy(x0), SR1(Inverse()))
@printf("NN! SR1  (inverse): %2.2e  %d\n", norm(res[3], Inf), res[4])
res = minimize(him_outofplace, x0s, SR1(Inverse()))
@printf("NN  SR1(S)  (inverse): %2.2e  %d\n", norm(res[3], Inf), res[4])
res = minimize(him_outofplace, x0, SR1(Direct()))
@printf("NN  SR1   (direct): %2.2e  %d\n", norm(res[3], Inf), res[4])
res = minimize!(him_inplace, copy(x0), SR1(Direct()))
@printf("NN! SR1   (direct): %2.2e  %d\n", norm(res[3], Inf), res[4])
res = minimize(him_outofplace, x0s, SR1(Direct()))
@printf("NN  SR1(S)   (direct): %2.2e  %d\n", norm(res[3], Inf), res[4])
println()
xrand = rand(2)
xrands = SVector{2}(xrand)
println("\nFrom a random point: ", xrand)

res = minimize(him_outofplace, xrand, GradientDescent(Inverse()))
@printf("NN  GD   (inverse): %2.2e  %d\n", norm(res[3], Inf), res[4])
res = minimize!(him_inplace, copy(xrand), GradientDescent(Inverse()))
@printf("NN! GD   (inverse): %2.2e  %d\n", norm(res[3], Inf), res[4])
res = minimize(him_outofplace, xrand, GradientDescent(Direct()))
@printf("NN  GD    (direct): %2.2e  %d\n", norm(res[3], Inf), res[4])
res = minimize!(him_inplace, copy(xrand), GradientDescent(Direct()))
@printf("NN! GD    (direct): %2.2e  %d\n", norm(res[3], Inf), res[4])

res = minimize(him_outofplace, xrand, BFGS(Inverse()))
@printf("NN   BFGS (inverse): %2.2e  %d\n", norm(res[3], Inf), res[4])
res = minimize!(him_inplace, copy(xrand), BFGS(Inverse()))
@printf("NN!  BFGS (inverse): %2.2e  %d\n", norm(res[3], Inf), res[4])
res = minimize(him_static, xrand, BFGS(Direct()))
@printf("NN(S) BFGS (direct): %2.2e  %d\n", norm(res[3], Inf), res[4])
res = minimize(him_outofplace, xrand, BFGS(Direct()))
@printf("NN  BFGS  (direct): %2.2e  %d\n", norm(res[3], Inf), res[4])
res = minimize!(him_inplace, copy(xrand), BFGS(Direct()))
@printf("NN! BFGS  (direct): %2.2e  %d\n", norm(res[3], Inf), res[4])

res = minimize(him_outofplace, xrand, DFP(Inverse()))
@printf("NN  DFP  (inverse): %2.2e  %d\n", norm(res[3], Inf), res[4])
res = minimize!(him_inplace,  copy(xrand), DFP(Inverse()))
@printf("NN! DFP  (inverse): %2.2e  %d\n", norm(res[3], Inf), res[4])
res = minimize(him_outofplace, xrand, DFP(Direct()))
@printf("NN  DFP   (direct): %2.2e  %d\n", norm(res[3], Inf), res[4])
res = minimize!(him_inplace,  copy(xrand), DFP(Direct()))
@printf("NN! DFP   (direct): %2.2e  %d\n", norm(res[3], Inf), res[4])

res = minimize(him_outofplace, xrand, SR1(Inverse()))
@printf("NN  SR1  (inverse): %2.2e  %d\n", norm(res[3], Inf), res[4])
res = minimize!(him_inplace,  copy(xrand), SR1(Inverse()))
@printf("NN! SR1  (inverse): %2.2e  %d\n", norm(res[3], Inf), res[4])
res = minimize(him_outofplace, xrand, SR1(Direct()))
@printf("NN  SR1   (direct): %2.2e  %d\n", norm(res[3], Inf), res[4])
res = minimize!(him_inplace,  copy(xrand), SR1(Direct()))
@printf("NN! SR1   (direct): %2.2e  %d\n", norm(res[3], Inf), res[4])

end


text:
@testset "Sign" begin
    @testset for T in (Posit8, Posit16, Posit32,
                    # pX1_mul is not supported yet by C library
                    # Posit8_1, Posit16_1, Posit24_1,
                    Posit8_2, Posit16_2, Posit24_2)

        @test one(T) == sign(one(T))
        @test zero(T) == sign(zero(T))
        @test minusone(T) == sign(minusone(T))
        @test zero(T) == sign(notareal(T))
    end
end

ext: jl


text:
# From https://en.wikipedia.org/wiki/Fowler%E2%80%93Noll%E2%80%93Vo_hash_function
#algorithm fnv-1 is
# hash := FNV_offset_basis
#    for each byte_of_data to be hashed do
#        hash := hash × FNV_prime
#        hash := hash XOR byte_of_data
#    return hash
#and fnv-1a is the same but with switched order of statements
#in loop body.

const FNV_Primes = Dict(
    32 => UInt32(0x01000193),
    64 => UInt64(0x00000100000001B3),
    128 => UInt128(0x0000000001000000000000000000013B)
)
#const FNV_Prime_256 = parse(BigInt, "374144419156711147060143317175368453031918731002211")

const FNV_Offsets = Dict(
    32 => UInt32(0x811c9dc5),
    64 => UInt64(0xcbf29ce484222325),
    128 => UInt128(0x6c62272e07bb014262b821756295c58d)
)

abstract type AbstractHashFunction{B} end

@inline unsafe_calcwithoffset(hf::AbstractHashFunction, unp::Ptr, len::Integer, offset::Integer = 0) =
    unsafe_calc(hf, reinterpret(Ptr{UInt8}, unp) + offset, len)

@inline calc(hf::AbstractHashFunction, bytes::Vector{UInt8}) =
    unsafe_calc(hf, pointer(bytes), length(bytes))

@inline calc(hf::AbstractHashFunction, s::String) =
    unsafe_calc(hf, pointer(s), length(s))

@inline unsafe_calc(hf::AbstractHashFunction, s::String, offset::Integer, len::Integer) =
    unsafe_calcwithoffset(hf, pointer(s), len, offset)

@inline function calc(hf::AbstractHashFunction, s::String, offset::Integer, len::Integer)
    @assert (offset + len) <= length(s)
    unsafe_calc(hf, s, offset, len)
end

const HashLenToType = Dict(
    32 => UInt32,
    64 => UInt64,
    128 => UInt128
)

struct FNV1A{B,T} <: AbstractHashFunction{B}
    offset::T
    prime::T
end

function unsafe_calc(hf::FNV1A, pnt::Ptr{UInt8}, len::Integer, hash = hf.offset)
    for i in 1:len
        hash = xor(hash, unsafe_load(pnt))
        hash = hash * hf.prime
        pnt += 1
    end
    hash
end

struct FNV1{B,T} <: AbstractHashFunction{B}
    offset::T
    prime::T
end

function unsafe_calc(hf::FNV1, pnt::Ptr{UInt8}, len::Integer, hash = hf.offset)
    for i in 1:len
        hash = hash * hf.prime
        hash = xor(hash, unsafe_load(pnt))
        pnt += 1
    end
    hash
end

function FNVHash(HT, hlen::Integer)
    @assert in(hlen, keys(HashLenToType))
    T = HashLenToType[hlen]
    HT{hlen, T}(FNV_Offsets[hlen], FNV_Primes[hlen])
end

FNV1A(hlen::Integer) = FNVHash(FNV1A, hlen)
FNV1(hlen::Integer) = FNVHash(FNV1, hlen)


text:
using QuasiCopula, LinearAlgebra, Random, GLM
using DataFrames, Statistics
using BenchmarkTools, Test

BLAS.set_num_threads(1)
Threads.nthreads()

p = 3    # number of fixed effects, including intercept
m = 1    # number of variance components
# true parameter values
# βtrue = [0.2; 0.1; -0.05; 0.2; 0.1; -0.1]
Random.seed!(123345)
βtrue = rand(Uniform(-0.2, 0.2), 2 * p)
θtrue = [0.1]

# generate data
intervals = zeros(2 * p + m, 2) #hold intervals
curcoverage = zeros(2 * p + m) #hold current coverage resutls
trueparams = [βtrue; θtrue] #hold true parameters

#simulation parameters
samplesize = 10000
ni = 2

currentind = 1
d1 = Poisson()
d2 = Bernoulli()
vecdist = [d1, d2]

link1 = LogLink()
link2 = LogitLink()
veclink = [link1, link2]

T = Float64
VD = typeof(vecdist)
VL = typeof(veclink)

Random.seed!(12345)
Xstack = [ones(samplesize) randn(samplesize, p - 1)]
gcs = Vector{Poisson_Bernoulli_VCObs{T, VD, VL}}(undef, samplesize)
Γ = θtrue[1] * ones(ni, ni)
vecd = Vector{DiscreteUnivariateDistribution}(undef, ni)

for i in 1:samplesize
    xi = Xstack[i, :]
    X = [transpose(xi) zeros(size(transpose(xi))); zeros(size(transpose(xi))) transpose(xi)]
    η = X * βtrue
    μ = zeros(ni)
    for j in 1:ni
        μ[j] = GLM.linkinv(veclink[j], η[j])
        vecd[j] = typeof(vecdist[j])(μ[j])
    end
    nonmixed_multivariate_dist = NonMixedMultivariateDistribution(vecd, Γ)
    # simuate single vector y
    y = Vector{Float64}(undef, ni)
    res = Vector{Float64}(undef, ni)
    rand(nonmixed_multivariate_dist, y, res)
    V = [ones(ni, ni)]
    gcs[i] = Poisson_Bernoulli_VCObs(y, xi, V, vecdist, veclink)
end

gcm = Poisson_Bernoulli_VCModel(gcs)
# precompile
println("precompiling Bivariate Mixed VCM fit")
gcm2 = deepcopy(gcm);
QuasiCopula.fit!(gcm2, IpoptSolver(print_level = 0, max_iter = 20));

fittime = @elapsed QuasiCopula.fit!(gcm, IpoptSolver(print_level = 5, max_iter = 100, tol = 10^-8, limited_memory_max_history = 16, hessian_approximation = "limited-memory"))
@show fittime
@show gcm.β
@show gcm.θ
@show gcm.∇β
@show gcm.∇θ

loglikelihood!(gcm, true, true)
vcov!(gcm)
@show QuasiCopula.confint(gcm)
# mse under our model
mseβ, mseθ = MSE(gcm, βtrue, θtrue)
@show mseβ
@show mseθ

@test mseβ < 0.01
@test mseθ < 0.01


text:
"""
Package: Forecast

Store results from the functions `acf`, `ccf` and `pacf`

# Arguments
`ccf`    An array with results from ccf, acf and pacf
`N`      Length of ccf
`type    Type of CCF
`lag`    Maximum number of lags
`alph`   CI thresholds
`ci`     CI for the alpha
`auto`   Auto-correlation
`call`   Method called to generate ccf
"""
mutable struct CCF
    ccf::AbstractArray{<:Real}
    N::Integer
    type::String
    lag::Integer
    alpha::Tuple{Real,Real}
    ci::Tuple{Real,Real}
    auto::Bool
    call::String
end


text:
using DataDeps

register(DataDep(
  "communities-and-crime",
  "https://archive.ics.uci.edu/ml/datasets/Communities+and+Crime",
  "https://archive.ics.uci.edu/ml/machine-learning-databases/communities/communities.data",
  "09e0b5c07eae24c1efab19b2edee05e160e7f5743b6f31e31eec3d73624da2ea",
  post_fetch_method=(path -> begin
    UCIData.process_dataset(path,
        target_index=128,
        feature_indices=6:127,
        id_indices=4,
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
module TestCompositionModelsInspection

using Test
using MLJBase
using Tables
import MLJBase
using ..Models
using CategoricalArrays
using OrderedCollections
using Statistics

X = (x1=rand(3), x2=rand(3), x3=rand(3))
y = float.([1, 2, 3])

@load KNNRegressor

@testset "tuple_keyed_on_model_names" begin

    model1 = KNNRegressor(K=1)
    model2 = KNNRegressor(K=2)
    model3 = KNNRegressor(K=3)
    model4 = KNNRegressor(K=4)

    _models = deepcopy.((model4, model3, model1))

    a1 = machine(model1, X, y)
    a2 = machine(model1, X, y)
    b = machine(model2, X, y)
    c = machine(model3, X, y)
    d = machine(model4, X, y)

    fit!.([a1, a2, b, c, d])

    # mutating the models should not effect the result:
    model1.K=0
    model2.K=0
    model3.K=0
    model4.K=0

    item_given_machine =
        LittleDict(a1=>"a1", a2=>"a2", b=>"b", c=>"c", d=>"d")

    _names = (:knn4, :knn3, :knn1)

    nt = MLJBase.tuple_keyed_on_model_names(item_given_machine,
                                            _models,
                                            _names)

    @test nt == (knn4="d", knn3="c", knn1=["a1", "a2"])
end

mutable struct Bar <: DeterministicComposite
    scale::Float64
    rgs
    input_stand
    target_stand
end

function MLJBase.fit(model::Bar, verbosity, X, y)
    Xs = source(X)
    ys = source(y)
    y1 = model.scale*ys
    mach1 = machine(model.input_stand, Xs)
    X1 = transform(mach1, Xs)
    mach2 = machine(model.target_stand, y1)
    z = transform(mach2, ys)
    mach3 = machine(model.rgs, X1, z)
    zhat = predict(mach3, X1)
    yhat = inverse_transform(mach2, zhat)
    mach = machine(Deterministic(), Xs, ys; predict=yhat)
    fit!(mach)
    return mach()
end

scale=0.97
rgs = KNNRegressor()
input_stand = Standardizer()
target_stand = UnivariateStandardizer()
model = Bar(scale, rgs, input_stand, target_stand)

mach = machine(model, X, y)
fit!(mach)

@testset "models_and_names(::Machine{<:Composite})" begin
    _models, _names = MLJBase.models_and_names(mach)
    @test _models == (rgs, input_stand, target_stand)
    @test _names == (:rgs, :input_stand, :target_stand)
end

@testset "user-friendly inspection of reports and fitted params" begin

    # mutating the models should not effect result:
    model.scale = 42.3
    model.rgs.K = 42
    model.input_stand.features=[:x1,]

    r = report(mach)
    @test isempty(r.rgs)
    @test Set(r.input_stand.features_fit) == Set([:x1, :x2, :x3])
    @test isempty(r.target_stand)
    @test r.report_given_machine isa AbstractDict

    fp = fitted_params(mach)
    @test fp.rgs isa NamedTuple{(:tree,)}
    @test fp.input_stand.mean_and_std_given_feature[:x1] |> collect ≈
        [mean(X.x1), std(X.x1)]
    @test fp.target_stand.fitresult |> collect ≈
        [mean(0.97*y), std(0.97*y)]  # scale = 0.97 at fit! call
end

end

true


text:
module SimpleSolvers

    using Distances
    using ForwardDiff
    using LinearAlgebra
    using Printf

    import Base.minimum
    import Base.Callable

    include("utils.jl")

    export solve!
    export config, result, state, status
    export algorithm, objective
    export solution, minimizer, minimum

    export Newton, DFP, BFGS

    include("base/methods.jl")

    export UnivariateObjective,
           MultivariateObjective

    export value, value!, value!!,
           derivative, derivative!, derivative!!,
           gradient, gradient!, gradient!!,
           hessian, hessian!, hessian!!,
           d_calls, f_calls, g_calls, h_calls

    include("base/objectives.jl")

    export Options

    include("base/options.jl")

    export GradientParameters,
           GradientParametersAD,
           GradientParametersFD,
           GradientParametersUser

    export compute_gradient,
           compute_gradient!,
           compute_gradient_ad!,
           compute_gradient_fd!

    export check_gradient,
           print_gradient

    include("base/gradient.jl")

    export Hessian,
           HessianAD,
           HessianUser

    export compute_hessian,
           compute_hessian!,
           compute_hessian_ad!

    export check_hessian,
           print_hessian

    include("base/hessian.jl")

    export JacobianParameters,
           JacobianParametersAD,
           JacobianParametersFD,
           JacobianParametersUser

    export compute_jacobian!,
           compute_jacobian_ad!,
           compute_jacobian_fd!

    export check_jacobian,
           print_jacobian

    include("base/jacobian.jl")

    export LinearSolver, LUSolver, LUSolverLAPACK,
           factorize!

    include("linear/linear_solvers.jl")
    include("linear/lu_solver.jl")
    include("linear/lu_solver_lapack.jl")

    export bracket_minimum

    include("bracketing/bracketing.jl")
    include("bracketing/bracket_minimum.jl")

    export Linesearch, Static
    export Backtracking, backtracking,
           Bisection, bisection,
           Quadratic, quadratic

    include("linesearch/linesearch.jl")
    include("linesearch/static.jl")
    include("linesearch/backtracking.jl")
    include("linesearch/bisection.jl")
    include("linesearch/quadratic.jl")

    export NonlinearSolver, NonlinearSolverException,
           AbstractNewtonSolver, NLsolveNewton, NewtonSolver, QuasiNewtonSolver,
           residual_initial!, residual_absolute!, residual_relative!,
           assess_convergence, assess_convergence!,
           print_status, check_solver_status,
           get_solver_status, get_solver_status!,
           solve!

    include("nonlinear/nonlinear_solver_status.jl")
    include("nonlinear/nonlinear_solver.jl")
    include("nonlinear/abstract_newton_solver.jl")
    include("nonlinear/newton_solver.jl")
    include("nonlinear/quasi_newton_solver.jl")
    include("nonlinear/nlsolve_newton.jl")

    export Optimizer,
           NewtonOptimizer,
           BFGSOptimizer,
           DFPOptimizer,
           HessianBFGS,
           HessianDFP

    include("optimization/optimizer_status.jl")
    include("optimization/optimizer_result.jl")
    include("optimization/optimizer.jl")
    include("optimization/hessian_bfgs.jl")
    include("optimization/hessian_dfp.jl")
    include("optimization/newton_optimizer.jl")

end


text:
struct CoNLL2000{S}
    filepaths::Vector{S}
    year::Int
    trainpath::String
    testpath::String
end

function init_datadeps(::Type{CoNLL2000})
    for (file, checksum) in [("train", "bcbbe17c487d0939d48c2d694622303edb3637ca9c4944776628cd1815c5cb34"),
	                ("test", "2695b931a505b65b1d3fe2a5aa0d749edcd53d786278da6d7726714e1426fbe8")]
        register(DataDep("CoNLL 2000 $file",
	    """
	    Website: https://www.clips.uantwerpen.be/conll2000/chunking/$file.txt.gz
	    Data Website:https://www.clips.uantwerpen.be/conll2000/chunking/

	    The Fourth Conference on Natural Language Learning (CoNLL-2003) was held on September 13 and September 14, 2000 in conjunction with ICGI-2000 and LLL-2000 at the Instituto Superior Técnico in Lisbon, Portugal.

	    The details of the task are in the following paper, please cite the following publication if you are using the corpora.

	    Erik F. Tjong Kim Sang, Sabine Buchholz. "Introduction to the CoNLL-2000 Shared Task: Chunking." In: Proceedings of CoNLL-2000 and LLL-2000, pages 127-132, Lisbon, Portugal, 2000.
	    https://www.clips.uantwerpen.be/conll2000/pdf/12732tjo.pdf
	    """,
	    "https://www.clips.uantwerpen.be/conll2000/chunking/$file.txt.gz",
	    checksum;
	    post_fetch_method = unpack
        ))
    end
end

function CoNLL2000(path1, path2, year=2000)
    @assert(isdir(path1), path1)
    @assert(isdir(path2), path2)

    files = Dict("train" => joinpath(path1, "train.txt"),
             "test" => joinpath(path2, "test.txt"))

    return CoNLL2000(collect(values(files)), year, files["train"],
                  files["test"])
end

CoNLL2000() = CoNLL2000(datadep"CoNLL 2000 train", datadep"CoNLL 2000 test")

MultiResolutionIterators.levelname_map(::Type{CoNLL2000}) = [
    :sent=>1, :sentence=>1,
    :word=>2, :token=>2,
    :char=>3, :character=>3
    ]

function parse_conll2000_tagged_word(line::AbstractString)
    tokens_tags = split(line)
    length(tokens_tags) != 3 && throw("Error parsing line: \"$line\". Invalid Format.")
    return POSTaggedWord(tokens_tags[2], tokens_tags[3], tokens_tags[1])
end

function parse_conll2000file(filename)
    local sent
    sents = @NestedVector(POSTaggedWord, 2)()
    context = Document(intern(basename(filename)), sents)

    function new_sentence()
        sent = @NestedVector(POSTaggedWord,1)()
        push!(sents, sent)
    end

    # words
    get_tagged(line) = push!(sent, parse_conll2000_tagged_word(line))

    new_sentence()
    # parse
    for line in eachline(filename)
        if length(line) == 0 && length(sent) > 0
            new_sentence()
        else
            get_tagged(line)
        end
    end
    isempty(sents[end]) && deleteat!(sents, lastindex(sents))

    return context
end

function load(corpus::CoNLL2000, file="train")
    file == "train" && return parse_conll2000file(corpus.trainpath)
    file == "test" && return parse_conll2000file(corpus.testpath)
    file == "dev" && throw("Dev set not available in $year.")
    throw("Invalid filename! Available datasets are `train` & `test`")
end


text:
include("../src/interface.jl")
include("../src/henon_heiles.jl")
using LinearAlgebra
using Plots
using Plots.PlotMeasures
function getmindif(tab::Vector{Vector{BigFloat}})
    nmmin = Inf
    ret = (0, 0)
    for i = 1:size(tab, 1)
        for j = (i+1):size(tab, 1)
            nm = norm(tab[i] - tab[j], Inf)
            if nm < nmmin
                nmmin = nm
                ret = i, j
            end
        end
    end
    return ret
end
function fctmain(n_tau, prec)
    setprecision(prec)
    u0 = BigFloat.([90, -44, 83, 13] // 100)
    t_max = big"1.0"
    epsilon = big"0.0017"
    println("epsilon=$epsilon")
    nbmaxtest = 12
    ordmax = 17
    debord = 3
    pasord = 2
    y = ones(Float64, nbmaxtest, div(ordmax - debord, pasord) + 1)
    res_gen = Array{Array{BigFloat,1},2}(undef, nbmaxtest, div(ordmax - debord, pasord) + 1)
    x = zeros(Float64, nbmaxtest)
    A = [0 0 1 0; 0 0 0 0; -1 0 0 0; 0 0 0 0]
    prob = HiOscODEProblem(henon_heiles, u0, (big"0.0", t_max), missing, A, epsilon)
    tabsol = Array{Array{BigFloat,1},1}(undef, 0)
    println("ordmax=$ordmax")
    indref = 1
    ind = 1
    for order = debord:pasord:ordmax
        println("eps=$epsilon order=$order")
        nb = 100
        indc = 1
        labels = Array{String,2}(undef, 1, div(order - debord, pasord) + 1)
        resnorm = 0
        resnormprec = 1
        sol = undef
        println("preparation ordre $order + 2")
        par_u0 = missing
        borindc = (order < ordmax) ? nbmaxtest : nbmaxtest + 1
        while indc <= borindc
            if indc > nbmaxtest
                nb = 51 * 2^nbmaxtest
            end
            res = solve(
                prob,
                nb_tau = n_tau,
                order = order,
                nb_t = nb,
                par_u0 = par_u0,
                dense = false,
            )
            par_u0 = res.par_u0
            sol = res[end]
            push!(tabsol, sol)
            if indc <= nbmaxtest
                res_gen[indc, ind] = sol
                x[indc] = 1.0 / nb
            end
            (a, b) = getmindif(tabsol)
            if a != 0
                indref = a
                solref = tabsol[b]
                for i = 1:ind
                    borne = (i < ind) ? size(res_gen, 1) : min(indc, size(res_gen, 1))
                    for j = 1:borne
                        nm = norm(res_gen[j, i] - solref, Inf)
                        if nm == 0
                            nm = norm(res_gen[j, i] - tabsol[a], Inf)
                        end
                        y[j, i] = (nm < 1) ? nm : NaN
                    end
                end
            end
            println("result=$y")
            println("log2(y)=$(log2.(y))")
            nb *= 2
            indc += 1
        end
        for i = debord:pasord:order
            labels[1, (i-debord)÷pasord+1] = " order = $i "
        end
        p = Plots.plot(
            x,
            view(y, :, 1:ind),
            xlabel = "Δt",
            xaxis = :log,
            ylabel = "error",
            yaxis = :log,
            legend = :bottomright,
            label = labels,
            marker = 2,
        )
        pp = Plots.plot(
            x,
            view(y, :, 1:ind),
            xlabel = "Δt",
            xaxis = :log,
            ylabel = "error",
            yaxis = :log,
            legend = :bottomright,
            label = labels,
            marker = 2,
            bottom_margin = 30px,
            left_margin = 60px,
        )
        prec_v = precision(BigFloat)
        eps_v = convert(Float32, epsilon)
        Plots.savefig(p, "out/res4p_$(prec_v)_$(eps_v)_$(order)_$(n_tau)_henon_heiles.pdf")
        Plots.savefig(p, "out/res4p_$(prec_v)_$(eps_v)_$(order)_$(n_tau)_henon_heiles.png")
        Plots.savefig(
            pp,
            "out/res4pp_$(prec_v)_$(eps_v)_$(order)_$(n_tau)_henon_heiles.pdf",
        )
        Plots.savefig(
            pp,
            "out/res4pp_$(prec_v)_$(eps_v)_$(order)_$(n_tau)_henon_heiles.png",
        )
        ind += 1
    end
end

# testODESolver()

# for i=3:9
#     fctMain(2^i)
# end
fctmain(32, 512)


text:
include("utils.jl")
@info "Manifolds.jl Test settings:\n\n" *
      "Testing Float32:  $(TEST_FLOAT32)\n" *
      "Testing Double64: $(TEST_DOUBLE64)\n" *
      "Testing Static:   $(TEST_STATIC_SIZED)\n\n" *
      "Check test/utils.jl if you wish to change these settings."

@testset "Manifolds.jl" begin
    include_test("differentiation.jl")

    @testset "Ambiguities" begin
        # TODO: reduce the number of ambiguities
        if VERSION.prerelease == () #
            @test length(Test.detect_ambiguities(ManifoldsBase)) <= 17
            @test length(Test.detect_ambiguities(Manifolds)) == 0
            @test length(our_base_ambiguities()) <= 24
        else
            @info "Skipping Ambiguity tests for pre-release versions"
        end
    end

    @testset "utils test" begin
        @test Manifolds.usinc_from_cos(-1) == 0
        @test Manifolds.usinc_from_cos(-1.0) == 0.0
    end

    include_test("groups/group_utils.jl")
    include_test("notation.jl")
    # starting with tests of simple manifolds
    include_test("centered_matrices.jl")
    include_test("circle.jl")
    include_test("cholesky_space.jl")
    include_test("elliptope.jl")
    include_test("euclidean.jl")
    include_test("fixed_rank.jl")
    include_test("generalized_grassmann.jl")
    include_test("generalized_stiefel.jl")
    include_test("grassmann.jl")
    include_test("hyperbolic.jl")
    include_test("multinomial_doubly_stochastic.jl")
    include_test("multinomial_symmetric.jl")
    include_test("positive_numbers.jl")
    include_test("probability_simplex.jl")
    include_test("projective_space.jl")
    include_test("rotations.jl")
    include_test("skewsymmetric.jl")
    include_test("spectrahedron.jl")
    include_test("sphere.jl")
    include_test("sphere_symmetric_matrices.jl")
    include_test("stiefel.jl")
    include_test("symmetric.jl")
    include_test("symmetric_positive_definite.jl")
    include_test("symmetric_positive_semidefinite_fixed_rank.jl")

    include_test("essential_manifold.jl")
    include_test("multinomial_matrices.jl")
    include_test("oblique.jl")
    include_test("torus.jl")

    #meta manifolds
    include_test("product_manifold.jl")
    include_test("power_manifold.jl")
    include_test("vector_bundle.jl")
    include_test("graph.jl")

    include_test("metric.jl")
    include_test("statistics.jl")

    # Lie groups and actions
    include_test("groups/groups_general.jl")
    include_test("groups/array_manifold.jl")
    include_test("groups/circle_group.jl")
    include_test("groups/translation_group.jl")
    include_test("groups/special_orthogonal.jl")
    include_test("groups/product_group.jl")
    include_test("groups/semidirect_product_group.jl")
    include_test("groups/special_euclidean.jl")
    include_test("groups/group_operation_action.jl")
    include_test("groups/rotation_action.jl")
    include_test("groups/translation_action.jl")
    include_test("groups/metric.jl")

    include_test("recipes.jl")
end


text:
# TODO: Document additional trick functionality
"""
    gbtranspose!(C::GBMatrix, A::GBMatrix; kwargs...)::Nothing

Eagerly evaluated matrix transpose, storing the output in `C`.

# Arguments
- `C::GBMatrix`: output matrix.
- `A::GBMatrix`: input matrix.

# Keywords
- `mask::Union{Nothing, GBMatrix} = nothing`: optional mask.
- `accum::Union{Nothing, AbstractBinaryOp} = nothing`: binary accumulator operation
    where `C[i,j] = accum(C[i,j], T[i,j])` where T is the result of this function before accum is applied.
- `desc::Union{Nothing, Descriptor} = DEFAULTDESC`
"""
function gbtranspose!(
    C::GBVecOrMat, A::GBArray;
    mask = nothing, accum = nothing, desc = nothing
)
    mask, accum = _handlenothings(mask, accum)
    desc = _handledescriptor(desc; in1=A)
    accum = getaccum(accum, eltype(C))
    libgb.GrB_transpose(C, mask, accum, parent(A), desc)
    return C
end

"""
    gbtranspose(A::GBMatrix; kwargs...)::GBMatrix

Eagerly evaluated matrix transpose which returns the transposed matrix.

# Keywords
- `mask::Union{Nothing, GBMatrix} = nothing`: optional mask.
- `accum::Union{Nothing, AbstractBinaryOp} = nothing`: binary accumulator operation
    where `C[i,j] = accum(C[i,j], T[i,j])` where T is the result of this function before accum is applied.
- `desc::Union{Nothing, Descriptor} = nothing`

# Returns
- `C::GBMatrix`: output matrix.
"""
function gbtranspose(
    A::GBArray;
    mask = nothing, accum = nothing, desc = nothing
)
    C = similar(A, size(A,2), size(A, 1))
    gbtranspose!(C, A; mask, accum, desc)
    return C
end

function LinearAlgebra.transpose(A::GBArray)
    return Transpose(A)
end

function Base.copy!(
    C::GBMatrix, A::LinearAlgebra.Transpose{<:Any, <:GBArray};
    mask = nothing, accum = nothing, desc = nothing
)
    return gbtranspose!(C, A.parent; mask, accum, desc)
end

function Base.copy(
    A::LinearAlgebra.Transpose{<:Any, <:GBArray};
    mask = nothing, accum = nothing, desc = nothing
)
    return gbtranspose(parent(A); mask, accum, desc)
end

#This is ok per the GraphBLAS Slack channel. Should change its effect on Complex input.
LinearAlgebra.adjoint(A::GBMatrix) = transpose(A)

LinearAlgebra.adjoint(v::GBVector) = transpose(v)

#arrrrgh, type piracy.
LinearAlgebra.transpose(::Nothing) = nothing


text:
# Generic equation of state types

export EOS, istempdependent

# EOS type

abstract type EOS end

# Include EOS subsections

include("eos_functions.jl")
include("eos_tables.jl")
include("eos_combined.jl")

# Testing for inclusion

Base.in(xy, eos::EOS) = in(xy..., eos)

# Marking temperature-dependent EOSes

"Does a given EOS have a temperature dependent component?"
function istempdependent(eos::EOS)
    # by default, this throws an error - later types should override
    # in future versions of Julia, this might be easier to implement as a trait
    throw(NotImplementedError())
end


text:
"""
Module providing parallelised versions of [`QuantumOptics.timeevolution.mcwf`](@ref).

`ParallelMCWF` exports functions:
[`pmcwf`](@ref), [`load_trajs`](@ref), [`save_trajs`](@ref),
[`kets_to_dm`](@ref), [`kets_to_obs`](@ref)
"""
module ParallelMCWF

using Distributed, Base.Threads
if myid() == 1
    @info "ParallelMCWF loaded with $(nthreads()) threads per process."
    @info "ParallelMCWF loaded with $(nprocs()) workers."
    #@info "Processes must be added BEFORE using `using ParallelMCWF`, see Julia issue #3674."
end
using ProgressMeter, JLD2
import OrdinaryDiffEq
using QuantumOptics.bases, QuantumOptics.states, QuantumOptics.operators
using QuantumOptics.operators_dense, QuantumOptics.operators_sparse
using QuantumOptics.timeevolution
using QuantumOptics.operators_lazysum, QuantumOptics.operators_lazytensor, QuantumOptics.operators_lazyproduct
@everywhere using QuantumOptics.timeevolution.timeevolution_mcwf, QuantumOptics.timeevolution
const DecayRates = Union{Vector{Float64}, Matrix{Float64}, Nothing}
Base.@pure pure_inference(fout,T) = Core.Compiler.return_type(fout, T)

include("trajs_IO.jl")
export load_trajs, save_trajs
include("trajs_processing.jl")
export kets_to_dm, kets_to_obs
include("pmcwf.jl")
export pmcwf

end # module

using Distributed
@everywhere using QuantumOptics # TO DO: Can be avoided


text:
######################################################################
# KirUtil unit tests
# -----
# Licensed under MIT License.
module Helpers
    using KirUtil

    struct Immutable
        mutable::Mutable{Int}
        flag::Bool
    end
    Immutable(mutable, flag) = Immutable(Mutable(mutable), flag)
    Immutable(mutable) = Immutable(mutable, false)
    Base.:(==)(lhs::Immutable, rhs::Immutable) = lhs.mutable == rhs.mutable && lhs.flag == rhs.flag

    function testmultiply(base::Integer, factor::Number; truncate::Bool = false)
        if truncate
            floor(base*factor)
        else
            base*factor
        end
    end
end

include("./Test.Functions.jl")
include("./Test.Functionals.jl")
include("./Test.Macros.jl")
include("./Test.Types.jl")
include("./Test.XCopy.jl")
include("./Test.CancellableTasks.jl")
include("./Test.Events.jl")


text:
## TODO Update Refactor

using Plots
using PortfolioOpt
using PortfolioOpt.TestUtils

include("./examples/literature_examples/example_utils.jl")

############ Read Prices #############
Prices = get_test_data()
numD, numA = size(Prices) # A: Assets    D: Days

############ Efficient frontier ###################
# Parameters
t = 100
k_back = 60
numS = 10000
α = 0.95

# Statistics
Σ, r̄ = mean_variance(returns[(t - k_back):(t - 1), :])
r, P = returns_montecarlo(Σ, r̄[:, 1], numS)

# test opt functions
R = 0.003
model, w = base_model(numA)
min_cvar_noRf!(model, w, r̄, R, r, P, α);
x, e, Cvar, q1_α = compute_solution_stoc(model, w)
model, w = base_model(numA)
max_return_lim_cvar_noRf!(model, w, r̄, Cvar, r, P, α)
w, e, q1_α = compute_solution_stoc_2(model, w)

# Loop
range = 0:0.0002274:0.005;
len = size(range, 1);
E_noR = zeros(len);
σ_noR = zeros(len);
E_Cvar = zeros(len);
σ_Cvar = zeros(len);
E_limitCvar = zeros(len);
σ_limitCvar = zeros(len);

x = zeros(len, size(r̄, 1));
iter = 0;
for R in range
    global iter += 1

    model, w = base_model(numA)
    min_cvar_noRf!(model, w, r̄, R, r, P, α)
    x[iter, :], E_Cvar[iter], Cvar, q1_α = compute_solution_stoc(model, w)
    model, w = base_model(numA)
    max_return_lim_cvar_noRf!(model, w, r̄, Cvar, r, P, α)
    x2, E_limitCvar[iter] = compute_solution_stoc_2(model, w)

    σ_Cvar[iter] = sqrt(sum(x[iter, :]'Σ * x[iter, :]))
    σ_limitCvar[iter] = sqrt(sum(x2'Σ * x2))

    model, w = base_model(numA)
    po_minvar_limitmean_noRf!(model, w, Σ, r̄, R)
    x_quad, obj, E_noR[iter] = compute_solution(model, w)
    σ_noR[iter] = sqrt(obj)
end

# plot frontier
plt = plot(
    σ_limitCvar,
    E_Cvar;
    title="Efficient Frontier Cvar", # Cvar vs Mean-Variance
    xlabel="σ",
    ylabel="r",
    label="Limit Cvar",
    legend=:topleft,
);
plot!(plt, σ_Cvar, E_Cvar; label="Min Cvar");
plt


text:
# Adapted from arviz's arviz-darkgrid, arviz-whitegrid, and arviz-white styles, which are in
# turn based on Seaborn styles
const ARVIZ_BASE_THEME = PlotThemes.PlotTheme(;
    size=(720, 480),
    dpi=100,
    widen=true,
    ## Background
    background_color=:white,
    ## Legend
    legend=false,
    background_color_legend=false,
    foreground_color_legend=false,
    ## Font styles
    ## these font sizes are technically part of the theme, but they look terrible here
    # legendfontsize=15,
    # titlefontsize=16,
    # guidefontsize=15,
    # tickfontsize=14,
    plot_titlefontcolor=RGB(0.15, 0.15, 0.15),
    legendfontcolor=RGB(0.15, 0.15, 0.15),
    legendtitlefontcolor=RGB(0.15, 0.15, 0.15),
    titlefontcolor=RGB(0.15, 0.15, 0.15),
    foreground_color_text=RGB(0.15, 0.15, 0.15),
    guidefontcolor=RGB(0.15, 0.15, 0.15),
    tickfontcolor=RGB(0.15, 0.15, 0.15),
    ## Colors
    color_palette=[
        colorant"#2a2eec",
        colorant"#fa7c17",
        colorant"#328c06",
        colorant"#c10c90",
        colorant"#933708",
        colorant"#65e5f3",
        colorant"#e6e135",
        colorant"#1ccd6a",
        colorant"#bd8ad5",
        colorant"#b16b57",
    ],
    colorgradient=:viridis,
    ## Marker
    markerstrokewidth=0,
    markersize=7,
)

const THEMES = Dict{Symbol,PlotThemes.PlotTheme}(
    :arviz_darkgrid => PlotThemes.PlotTheme(
        ARVIZ_BASE_THEME;
        ## Background
        background_color_inside=colorant"#eeeeee",
        ## Axes / Ticks
        showaxis=true,
        minorticks=false,
        foreground_color_axis=:white,
        foreground_color_border=:white,
        tick_direction=:out,
        ## Grid
        grid=:true,
        gridstyle=:solid,
        gridalpha=1,
        foreground_color_grid=:white,
    ),
    :arviz_whitegrid => PlotThemes.PlotTheme(
        ARVIZ_BASE_THEME;
        ## Background
        background_color_inside=:white,
        ## Axes / Ticks
        framestyle=:box,
        showaxis=true,
        minorticks=false,
        foreground_color_axis=RGB(0.8, 0.8, 0.8),
        foreground_color_border=RGB(0.8, 0.8, 0.8),
        tick_direction=:out,
        ## Grid
        grid=:true,
        gridstyle=:solid,
        gridalpha=1,
        foreground_color_grid=RGB(0.8, 0.8, 0.8),
    ),
    :arviz_white => PlotThemes.PlotTheme(
        ARVIZ_BASE_THEME;
        ## Background
        background_color_inside=:white,
        ## Axes / Ticks
        showaxis=true,
        minorticks=false,
        foreground_color_axis=:white,
        foreground_color_border=:black,
        tick_direction=:out,
        ## Grid
        grid=:false,
        gridalpha=1,
        foreground_color_grid=:black,
    ),
)


text:
"The effective low frequency bulk modulus and density of a material filled with particles"
function effective_medium(medium::Acoustic{T,Dim}, species::Species{T,Dim};
        numberofparticles::Number = Inf
    ) where {T<:AbstractFloat,Dim}

    scale_number_density = one(T) - one(T) / numberofparticles
    φ = scale_number_density * sum(volume_fraction.(species))

    # calculate effective properties
    β = medium.ρ*medium.c^2 # medium bulk modulus
    if abs(β) == zero(T) || abs(prod(s -> s.particle.medium.ρ * s.particle.medium.c^2, species)) == zero(T)
        β_eff = zero(T)
    else
        β_eff = ((one(T)-φ) / β + scale_number_density * sum(s.volume_fraction / (s.particle.medium.ρ * s.particle.medium.c^2) for s in species))^(-one(T))
    end

    ρ_frac = scale_number_density * sum(
        map(species) do s
            if s.particle.medium.ρ == Inf
                - s.volume_fraction / T(Dim - 1)
            else
                s.volume_fraction * (medium.ρ - s.particle.medium.ρ)/(medium.ρ + T(Dim - 1) * s.particle.medium.ρ)
            end
        end
    )

    ρ_eff = medium.ρ*(1 - ρ_frac)/(1 + T(Dim - 1) * ρ_frac)

    return Acoustic(Dim; ρ = ρ_eff, c = sqrt(β_eff/ρ_eff))
end

"""
reflection_coefficient(PlaneSource, Acoustic[, Halfspace = Halfspace(-psource.direction)])

calculates the reflection coefficient from a homogenious halfspace (assumed to direct incidence if not given), which is also the low frequency reflection from a particulate material when using the effective_medium.
"""
reflection_coefficient(ω::T, source::PlaneSource{T}, reflect_medium::Acoustic{T}, halfspace::Halfspace{T} = Halfspace(-psource.direction)) where T<:AbstractFloat = reflection_transmission_coefficients(ω, psource, reflect_medium, halfspace)[1]

"""
reflection_transmission_coefficients(PlaneSource, Acoustic[, Halfspace = Halfspace(-psource.direction)])

calculates the reflection and transmission coefficient from a homogenious halfspace (assumed to direct incidence if not given), which is also the low frequency reflection from a particulate material when using the effective_medium.
"""
function reflection_transmission_coefficients(ω::T, psource::PlaneSource{T,Dim,1,Acoustic{T,Dim}}, reflect_medium::Acoustic{T,Dim}, halfspace::Halfspace{T,Dim} = Halfspace(-psource.direction)) where {T<:AbstractFloat,Dim}

    k_in = ω / psource.medium.c
    k_r = if abs(reflect_medium.c) == zero(T)
        T(Inf) + zero(T) * im
    else
        ω / reflect_medium.c
    end

    # make sure the normal outward pointing
    n = halfspace.normal / norm(halfspace.normal)
    if real(dot(n,psource.direction)) > 0
        n = -n
    end

    θin = transmission_angle(psource.direction, n)

    vtran = transmission_direction(k_r, k_in * psource.direction, n)
    θtran = transmission_angle(vtran, n)

    q = reflect_medium.c*reflect_medium.ρ / (psource.medium.c*psource.medium.ρ)

    Z = dot(n,halfspace.origin)

    R1 = (exp(2im * k_in * Z * cos(θin)) * field(psource)(zeros(T,Dim),ω)) * (q*cos(θin) - cos(θtran)) / (q*cos(θin) + cos(θtran))
    R1 = R1 * exp(im * k_in * Z * cos(θin)) # assuming a reflected wave  = R * exp(i kR.(X-X1)) where kR.X1 =  k_in * Z * cos(θin)

    T1 = exp(im * Z * (k_in * cos(θin) - k_r * cos(θtran))) * field(psource)(zeros(T,Dim),ω) * 2*q*cos(θin) / (q*cos(θin) + cos(θtran))
    T1 = T1 * exp(im * k_r * Z * cos(θtran)) # assuming a transmitter wave  = T * exp(i k1.(X-X1)) where k1.X1 =  k_r * Z * cos(θtran)

    return [R1,T1]
end

"""
    planewave_amplitudes(psource::PlaneSource, reflect_medium::Acoustic, plate::Plate)

Calculates the coefficients, or amplitudes, of the plane waves scattering from and transmitted inside the plate, for the plane wave source `psource` and the `plate`. The function returns `[R, T, P1, P2]` where `R` is the reflection coefficient, `T` is the coefficient of the transmitted wave, and `P1` (`P2`) are the amplitudes of the wave travling forward (backward) inside the plate.
"""
function planewave_amplitudes(ω::T, psource::PlaneSource{T,Dim,1,Acoustic{T,Dim}}, plate_medium::Acoustic{T,Dim}, plate::Plate{T,Dim}) where {T<:AbstractFloat,Dim}

    k0 = ω / psource.medium.c
    k1 = if abs(plate_medium.c) == zero(T)
        T(Inf) + zero(T) * im
    else
        ω / plate_medium.c
    end

    # make the normal outward pointing
    n = plate.normal / norm(plate.normal)
    if real(dot(n,psource.direction)) > 0
        n = -n
    end

    θ0 = transmission_angle(psource, plate)

    v1 = transmission_direction(k1, k0 .* psource.direction, n)
    θ1 = transmission_angle(v1, n)

    C0 = k0 * cos(θ0) / psource.medium.ρ
    C1 = k1 * cos(θ1) / plate_medium.ρ

    Z0 = dot(-n,plate.origin)
    Z1 = Z0 - plate.width / T(2)
    Z2 = Z0 + plate.width / T(2)

    denom = ((C0 + C1)^2 * exp(2im*k1*Z1*cos(θ1)) - (C0 - C1)^2 * exp(2im*k1*Z2*cos(θ1)))

    UR = (C0 - C1)*(C0 + C1) * exp(2im*k0*Z1*cos(θ0)) * (exp(2im*k1*Z1*cos(θ1)) - exp(2im*k1*Z2*cos(θ1))) / denom
    UR = UR * exp(im*k0*Z1*cos(θ0)) # we assume the reflected wave = UR * exp(im*KR.(X-X1)), where KR.X1 = im*k0*Z1*cos(θ0) = k0 .* [sin(θ0),0,-cos(θ0)], and X1 = [0.0,0.0,Z1], so that KR.X1 = im*k0*Z1*cos(θ0).

    UT = 4*C0*C1*exp(im*(k0*(Z1 - Z2)*cos(θ0) + k1*(Z1 + Z2)*cos(θ1))) / denom
    UT = UT * exp(im*k0*Z2*cos(θ0))

    UP1 = 2C0*(C0 + C1) * exp(im*Z1*(k0*cos(θ0) + k1*cos(θ1))) / denom
    UP1 = UP1 * exp(im*k1*Z1*cos(θ1))

    UP2 = - 2*C0*(C0 - C1)*exp(im*(k0*Z1*cos(θ0) + k1*(Z1 + 2*Z2)*cos(θ1))) / denom
    UP2 = UP2 * exp(im*k1*Z2*cos(θ1))

    return [UR, UT, UP1, UP2] * field(psource)(zeros(T,3),ω)
end

max_stars_repo_name: arturgower/EffectiveWaves


text:
let sys = DiagonalSystem(test_params_diag...)

    @test sys.energy == test_energy_diag
    @test sys.freq == test_freq
    @test sys.lin == test_lin_diag
    @test sys.quad == test_quad_diag

    @test isapprox(sys.deltas, [-4.91666666e-6, -8.66666666e-6, -1.2525e-5])
    @test isapprox(sys.ds, [-1.5e-3 -1.66666666e-3 -1.75e-3; -1.33333333e-3 -1.5e-3 -1.6e-3])

    @test isdiag(sys)

    @test diag(sys) == sys

    @test simplify(sys) == DiagonalSystem(test_energy_diag, test_freq, test_lin_diag, zero(test_quad_diag))

    @test DenseSystem(sys) isa DenseSystem
    @test DenseSystem(sys) == sys
end

@test_throws DomainError DiagonalSystem(test_freq, test_energy_diag, test_lin_diag, test_quad_diag)
@test_throws SurfaceCouplingException DiagonalSystem(test_params...)


text:
export treduce;

"""
    treduce(op::Function, itr::AbstractArray; init)

Multi-threaded version of [reduce(op, itr; init)](https://docs.julialang.org/en/v1.2/base/collections/#Base.reduce-Tuple{Any,Any}). Note that function `op` must **not** change type!
"""
function treduce(op::Function, itr::T; init) where T<:AbstractArray
	ensureThreaded();
	if isempty(itr)
		return init;
	end
	n = Threads.nthreads();
	tmp = Vector{Any}(undef, n);
	ind = eachindex(itr);
	len = length(ind);
	if len > n
		Threads.@threads for i in 1:n
			start = fld(len * (i - 1), n) + firstindex(ind);
			stop = fld(len * i, n) + firstindex(ind) - 1;
			if i == 1
				localInit = reduce(op, itr[ind[start:start]]; init = init);
			else
				localInit = reduce(op, itr[ind[start:start]]);
			end
			tmp[i] = reduce(op, itr[ind[nextind(ind, start):stop]]; init = localInit);
		end
	else
		tmp = itr;
		tmp[firstindex(tmp)] = reduce(op, tmp[firstindex(tmp):firstindex(tmp)]; init = init);
		n = len;
	end
	tmp2 = Vector{Any}(undef, cld(n, 2));
	while n > 1
		nn = cld(n, 2);
		Threads.@threads for j in 1:nn
			stop = min(2j, n);
			localInit = reduce(op, tmp[2j - 1]);
			tmp2[j] = reduce(op, tmp[2j:stop]; init = localInit);
		end
		n = nn;
		tmp = tmp2
	end
	return tmp[1];
end

"""
    treduce(op::Function, itr::AbstractArray)

Multi-threaded version of [reduce(op, itr)](https://docs.julialang.org/en/v1.2/base/collections/#Base.reduce-Tuple{Any,Any}). Note that function `op` must **not** change type!
"""
function treduce(op::Function, itr::T) where T<:AbstractArray
	ensureThreaded();
	if isempty(itr)
		throw(ArgumentError("reducing over an empty collection is not allowed"))
	end
	n = Threads.nthreads();
	tmp = Vector{Any}(undef, n);
	ind = eachindex(itr);
	len = length(ind);
	if len > n
		Threads.@threads for i in 1:n
			start = fld(len * (i - 1), n) + firstindex(ind);
			stop = fld(len * i, n) + firstindex(ind) - 1;
			tmp[i] = reduce(op, itr[ind[start:stop]]);
		end
	else
		tmp = itr;
		n = len;
	end
	tmp2 = Vector{Any}(undef, cld(n, 2));
	while n > 1
		nn = cld(n, 2);
		Threads.@threads for j in 1:nn
			stop = min(2j, n);
			tmp2[j] = reduce(op, tmp[2j - 1:stop]);
		end
		n = nn;
		tmp = tmp2
	end
	return tmp[1];
end


text:
@discord_object struct AuditLogChange
    new_value::Any
    old_value::Any
    key::String
end


text:
module PDBTools

using Parameters
using Printf
using InteractiveUtils
using Formatting
using StaticArrays
using LinearAlgebra: norm

export readPDB, writePDB, getseq, wget, edit!, oneletter, threeletter, residuename
export Atom, printatom, index, index_pdb, name, beta, occup
export Residue, eachresidue, resname, residue, resnum, chain, model, segname
export coor, maxmin, distance, closest
export atomic_number, element, mass, element_name
export formula, stoichiometry
export Sequence

#
# Data structures
#
include("./Atom.jl")
include("./Indexes_mmCIF_fields.jl")
include("./Residue.jl")

#
# Input and output functions
#
include("./all.jl")
include("./oneletter.jl")
include("./threeletter.jl")
include("./residuename.jl")
include("./formula.jl")
include("./check_mmCIF.jl")
include("./parse_int.jl")
include("./read_atom.jl")
include("./write_atom.jl")
include("./readPDB.jl")
include("./edit.jl")
include("./writePDB.jl")
include("./getseq.jl")
include("./coor.jl")
include("./distance.jl")
include("./maxmin.jl")
include("./same_residue.jl")
include("./wget.jl")

# Element properties
include("./elements.jl")
include("./atom_properties.jl")

# For selections
include("./select/select_includes.jl")
include("./alternate_conformation.jl")

end

max_stars_repo_name: mcubeg/PDBTools

max_issues_repo_name: mcubeg/PDBTools

max_forks_repo_name: mcubeg/PDBTools


text:
module StateFulls

using Flux: @functor
import Flux: trainable, reset!

export StateFull

"""
    StateFull{C,S,I}

Wrapper that turns a recurrent cell into a stateful unit, similar to Flux.Recur.
"""
mutable struct StateFull{C,S}
    cell::C
    state::S
    init::S
end # mutable struct StateFull

# convenient constructors
StateFull(c, s, i) = StateFull{typeof(c), typeof(s), typeof(i)}(c, s, i)
StateFull(c, i) = StateFull(c, i, i)

# application: handle state internally
function (sf::StateFull)(x)
    sf.state, y = sf.cell(sf.state, x)
    return y
end

# declarations to turn it into a usable Flux layer
@functor StateFull
trainable(sf::StateFull) = (sf.cell,)
function reset!(sf::StateFull)
    sf.state = sf.init
end # reset

end # module StateFulls


text:
using LinearAlgebra, Test
#JSO packages
using ADNLPModels, NLPModels, NLPModelsIpopt
# , NLPModelsKnitro
#This package
using Stopping, StoppingInterface

@testset "Rosenbrock with ∑x = 1" begin
  nlp = ADNLPModel(
    x -> (x[1] - 1.0)^2 + 100 * (x[2] - x[1]^2)^2,
    [-1.2; 1.0],
    x -> [sum(x) - 1],
    [0.0],
    [0.0],
  )
  sol = [-1.612771347383541; 2.612771347383541]
  stp = NLPStopping(nlp)
  @test status_stopping_to_stats(stp) == :unknown
  fill_in!(stp, sol)
  @test stop!(stp)
  status_stopping_to_stats(stp) == :first_order
  reinit!(stp, rstate = true, x = zeros(2))
  stp.meta.nb_of_stop = stp.meta.max_iter + 1
  fill_in!(stp, stp.current_state.x)
  @test stop!(stp)
  @test status_stopping_to_stats(stp) == :max_iter

  # reinit!(stp)
  # stp = knitro(stp)
  # @show status(stp), stp.current_state.x
  reinit!(stp)
  stp = StoppingInterface.ipopt(stp)
  @test status(stp) == :Optimal
end

for solver in (:lbfgs, :tron, :trunk, :ipopt), T in (Float32, Float64)
  @testset "Rosenbrock with $solver for T=$T" begin
    if solver == :ipopt && T != Float64
      continue
    end
    nlp = ADNLPModel(x -> (x[1] - 1)^2 + 100 * (x[2] - x[1]^2)^2, T[-1.2; 1.0])
    sol = ones(T, 2)
    stp = NLPStopping(nlp)
    @test status_stopping_to_stats(stp) == :unknown
    fill_in!(stp, sol)
    @test stop!(stp)
    status_stopping_to_stats(stp) == :first_order
    reinit!(stp, rstate = true, x = zeros(T, 2))
    stp.meta.nb_of_stop = stp.meta.max_iter + 1
    fill_in!(stp, stp.current_state.x)
    @test stop!(stp)
    @test status_stopping_to_stats(stp) == :max_iter

    # reinit!(stp)
    # stp = knitro(stp)
    # @show status(stp), stp.current_state.x
    reinit!(stp)
    stp.meta.atol = sqrt(eps(T))
    stp = StoppingInterface.eval(solver)(stp)
    @test status(stp) == :Optimal
  end
end

@testset "GSE status to meta" begin
  nlp = ADNLPModel(
    x -> (x[1] - 1.0)^2 + 100 * (x[2] - x[1]^2)^2,
    [-1.2; 1.0],
    x -> [sum(x) - 1],
    [0.0],
    [0.0],
  )
  stp = NLPStopping(nlp)
  StoppingInterface.stats_status_to_meta!(stp, :first_order)
  @test stp.meta.optimal
  reinit!(stp)
  StoppingInterface.stats_status_to_meta!(stp, :acceptable)
  @test stp.meta.suboptimal
  reinit!(stp)
  StoppingInterface.stats_status_to_meta!(stp, :infeasible)
  @test stp.meta.infeasible
  reinit!(stp)
  StoppingInterface.stats_status_to_meta!(stp, :small_step)
  @test stp.meta.stalled
  reinit!(stp)
  StoppingInterface.stats_status_to_meta!(stp, :max_eval)
  @test stp.meta.resources
  reinit!(stp)
  StoppingInterface.stats_status_to_meta!(stp, :max_iter)
  @test stp.meta.iteration_limit
  reinit!(stp)
  StoppingInterface.stats_status_to_meta!(stp, :max_time)
  @test stp.meta.tired
  reinit!(stp)
  StoppingInterface.stats_status_to_meta!(stp, :neg_pred)
  @test stp.meta.fail_sub_pb
  reinit!(stp)
  StoppingInterface.stats_status_to_meta!(stp, :not_desc)
  @test stp.meta.fail_sub_pb
  reinit!(stp)
  StoppingInterface.stats_status_to_meta!(stp, :unbounded)
  @test stp.meta.unbounded
  reinit!(stp)
  StoppingInterface.stats_status_to_meta!(stp, :user)
  @test stp.meta.stopbyuser
  reinit!(stp)
  StoppingInterface.stats_status_to_meta!(stp, :stalled)
  @test stp.meta.stalled
  reinit!(stp)
  StoppingInterface.stats_status_to_meta!(stp, :small_residual)
  @test stp.meta.stalled
  reinit!(stp)
  StoppingInterface.stats_status_to_meta!(stp, :small_step)
  @test stp.meta.stalled
end


text:
module TestView

using Base.Test
using CategoricalArrays
using NullableArrays

# == currently throws an error for Nullables
(==) = isequal

for (A, T) in zip((Array, NullableArray), (CategoricalArray, NullableCategoricalArray))
    for order in (true, false)
        for a in (1:10, 10:-1:1, ["a", "c", "b", "b", "a"])
            for inds in [1:2, :, 1, []]
                x = T(a, ordered=order)
                v = view(x, inds)
                @test levels(v) === levels(x)
                @test unique(v) == (ndims(v) > 0 ? sort(A(unique(a[inds]))) : A([a[inds]]))
                @test isordered(v) === isordered(x)
            end
        end
    end
end

end


text:
#
# User facing struct that unifies a ModelDef and a ModelInstance and delegates
# function calls to one or the other as appropriate.
#
using MacroTools

"""
    modeldef(m)

Return the `ModelDef` contained by Model `m`.
"""
modeldef(m::Model) = m.md

modelinstance(m::Model) = m.mi
modelinstance_def(m::Model) = modeldef(modelinstance(m))

is_built(m::Model) = !(dirty(m.md) || modelinstance(m) === nothing)

is_built(mm::MarginalModel) = (is_built(mm.base) && is_built(mm.modified))

@delegate compinstance(m::Model, name::Symbol) => mi
@delegate has_comp(m::Model, name::Symbol) => md

@delegate number_type(m::Model) => md

@delegate internal_param_conns(m::Model) => md
@delegate external_param_conns(m::Model) => md

@delegate external_params(m::Model) => md
@delegate external_param(m::Model, name::Symbol; missing_ok=false) => md

@delegate connected_params(m::Model) => md
@delegate unconnected_params(m::Model) => md

@delegate add_connector_comps!(m::Model) => md

"""
    connect_param!(m::Model, dst_comp_name::Symbol, dst_par_name::Symbol, src_comp_name::Symbol, src_var_name::Symbol,
    backup::Union{Nothing, Array}=nothing; ignoreunits::Bool=false, backup_offset::Union{Int, Nothing}=nothing)

Bind the parameter `dst_par_name` of one component `dst_comp_name` of model `m`
to a variable `src_var_name` in another component `src_comp_name` of the same model
using `backup` to provide default values and the `ignoreunits` flag to indicate the need
to check match units between the two.  The `backup_offset` argument, which is only valid
when `backup` data has been set, indicates that the backup data should be used for
a specified number of timesteps after the source component begins. ie. the value would be
`1` if the destination componentm parameter should only use the source component
data for the second timestep and beyond.
"""
@delegate connect_param!(m::Model,
                         dst_comp_name::Symbol, dst_par_name::Symbol,
                         src_comp_name::Symbol, src_var_name::Symbol,
                         backup::Union{Nothing, Array}=nothing;
                         ignoreunits::Bool=false, backup_offset::Union{Int, Nothing} = nothing) => md

"""
    connect_param!(m::Model, comp_name::Symbol, param_name::Symbol, ext_param_name::Symbol)

Bind the parameter `param_name` in the component `comp_name` of model `m` to the external parameter
`ext_param_name` already present in the model's list of external parameters.
"""
@delegate connect_param!(m::Model, comp_name::Symbol, param_name::Symbol, ext_param_name::Symbol) => md

"""
    connect_param!(m::Model, dst::Pair{Symbol, Symbol}, src::Pair{Symbol, Symbol}, backup::Array; ignoreunits::Bool=false)

Bind the parameter `dst[2]` of one component `dst[1]` of model `m`
to a variable `src[2]` in another component `src[1]` of the same model
using `backup` to provide default values and the `ignoreunits` flag to indicate the need
to check match units between the two.  The `backup_offset` argument, which is only valid
when `backup` data has been set, indicates that the backup data should be used for
a specified number of timesteps after the source component begins. ie. the value would be
`1` if the destination componentm parameter should only use the source component
data for the second timestep and beyond.

"""
function connect_param!(m::Model, dst::Pair{Symbol, Symbol}, src::Pair{Symbol, Symbol},
                           backup::Union{Nothing, Array}=nothing;
                           ignoreunits::Bool=false, backup_offset::Union{Int, Nothing} = nothing)
    connect_param!(m.md, dst[1], dst[2], src[1], src[2], backup; ignoreunits=ignoreunits, backup_offset=backup_offset)
end

"""
    disconnect_param!(m::Model, comp_name::Symbol, param_name::Symbol)

Remove any parameter connections for a given parameter `param_name` in a given component
`comp_def` in model `m`.
"""
@delegate disconnect_param!(m::Model, comp_name::Symbol, param_name::Symbol) => md

# TBD: these may not be needed as delegators
@delegate set_external_param!(m::Model, name::Symbol, value::ModelParameter) => md

@delegate set_external_param!(m::Model, name::Symbol,
                              value::Union{Number, AbstractArray, AbstractRange, Tuple};
                              param_dims::Union{Nothing,Array{Symbol}} = nothing) => md

@delegate add_internal_param_conn!(m::Model, conn::InternalParameterConnection) => md

# @delegate doesn't handle the 'where T' currently. This is the only instance of it for now...
function set_leftover_params!(m::Model, parameters::Dict{T, Any}) where T
    set_leftover_params!(m.md, parameters)
end

"""
    update_param!(m::Model, name::Symbol, value; update_timesteps = nothing)

Update the `value` of an external model parameter in model `m`, referenced by
`name`. The update_timesteps keyword argument is deprecated, we keep it here
just to provide warnings.
"""
@delegate update_param!(m::Model, name::Symbol, value; update_timesteps = nothing) => md

"""
    update_params!(m::Model, parameters::Dict{T, Any}; update_timesteps = nothing) where T

For each (k, v) in the provided `parameters` dictionary, `update_param!``
is called to update the external parameter by name k to value v.  Each key k
must be a symbol or convert to a symbol matching the name of an external parameter t
hat already exists in the model definition. The update_timesteps keyword argument
is deprecated, but temporarily remains as a dummy argument to allow warning detection.
"""
@delegate update_params!(m::Model, parameters::Dict; update_timesteps = nothing) => md

"""
    add_comp!(
        m::Model, comp_id::ComponentId, comp_name::Symbol=comp_id.comp_name;
        first::NothingInt=nothing,
        last::NothingInt=nothing,
        before::NothingSymbol=nothing,
        after::NothingSymbol=nothing,
        rename::NothingPairList=nothing
    )

Add the component indicated by `comp_id` to the model indicated by `m`. The component is added
at the end of the list unless one of the keywords `before` or `after` is specified. Note
that a copy of `comp_id` is made in the composite and assigned the give name. The optional
argument `rename` can be a list of pairs indicating `original_name => imported_name`. The optional
arguments `first` and `last` indicate the times bounding the run period for the given component,
which must be within the bounds of the model and if explicitly set are fixed.  These default
to flexibly changing with the model's `:time` dimension.
"""
function add_comp!(m::Model, comp_id::ComponentId, comp_name::Symbol=comp_id.comp_name; kwargs...)
    comp_def = add_comp!(m.md, comp_id, comp_name; kwargs...)
    return ComponentReference(m.md, comp_name)
end

"""
    add_comp!(
        m::Model, comp_def::AbstractComponentDef, comp_name::Symbol=comp_id.comp_name;
        first::NothingInt=nothing,
        last::NothingInt=nothing,
        before::NothingSymbol=nothing,
        after::NothingSymbol=nothing,
        rename::NothingPairList=nothing
    )

Add the component `comp_def` to the model indicated by `m`. The component is added at
the end of the list unless one of the keywords, `first`, `last`, `before`, `after`. Note
that a copy of `comp_id` is made in the composite and assigned the give name. The optional
argument `rename` can be a list of pairs indicating `original_name => imported_name`. The optional
arguments `first` and `last` indicate the times bounding the run period for the given component,
which must be within the bounds of the model and if explicitly set are fixed.  These default
to flexibly changing with the model's `:time` dimension.
"""
function add_comp!(m::Model, comp_def::AbstractComponentDef, comp_name::Symbol=comp_def.comp_id.comp_name; kwargs...)
    return add_comp!(m, comp_def.comp_id, comp_name; kwargs...)
end

# DEPRECATION - EVENTUALLY REMOVE
"""
    replace_comp!(
        m::Model, comp_id::ComponentId, comp_name::Symbol=comp_id.comp_name;
        before::NothingSymbol=nothing,
        after::NothingSymbol=nothing,
        reconnect::Bool=true
    )

Deprecated function for replacing the component with name `comp_name` in model `m` with the
new component specified by `comp_id`. Use the following syntax instead:

`replace!(m, comp_name => Mimi.compdef(comp_id); kwargs...)`

See docstring for `replace!` for further description of available functionality.
"""
function replace_comp!(m::Model, comp_id::ComponentId, comp_name::Symbol=comp_id.comp_name; kwargs...)
    error("Function `replace_comp!(m, comp_id, comp_name; kwargs...)` has been deprecated. Use `replace!(m, comp_name => Mimi.compdef(comp_id); kwargs...)` instead.")
end

# DEPRECATION - EVENTUALLY REMOVE
"""
    replace_comp!(
        m::Model, comp_def::ComponentDef, comp_name::Symbol=comp_id.comp_name;
        before::NothingSymbol=nothing,
        after::NothingSymbol=nothing,
        reconnect::Bool=true
    )

Deprecated function for replacing the component with name `comp_name` in model `m` with the
new component specified by `comp_def`. Use the following syntax instead:

`replace!(m, comp_name => comp_def; kwargs...)`

See docstring for `replace!` for further description of available functionality.
"""
function replace_comp!(m::Model, comp_def::ComponentDef, comp_name::Symbol=comp_def.comp_id.comp_name; kwargs...)
    error("Function `replace_comp!(m, comp_def, comp_name; kwargs...)` has been deprecated. Use `replace!(m, comp_name => comp_def; kwargs...)` instead.")
end

"""
    replace!(
        m::Model,
        old_new::Pair{Symbol, ComponentDef},
        before::NothingSymbol=nothing,
        after::NothingSymbol=nothing,
        reconnect::Bool=true
    )

For the pair `comp_name => comp_def` in `old_new`, replace the component with name `comp_name` in
the model `m` with the new component specified by `comp_def`. The new component is added
in the same position as the old component, unless one of the keywords `before` or `after` is
specified for a different position. The optional boolean argument `reconnect` with default value
`true` indicates whether the existing parameter connections should be maintained in the new
component. Returns a ComponentReference for the added component.
"""
function Base.replace!(m::Model, old_new::Pair{Symbol, ComponentDef}; kwargs...)
    comp_name, comp_def = old_new
    _replace!(m.md, comp_name => comp_def.comp_id; kwargs...)
    return ComponentReference(m.md, comp_name)
end

@delegate ComponentReference(m::Model, name::Symbol) => md

"""
    components(m::Model)

Return an iterator on the components in a model's model instance.
"""
@delegate components(m::Model) => mi

@delegate compdefs(m::Model) => md

@delegate compdef(m::Model, comp_name::Symbol) => md

@delegate Base.length(m::Model) => md

@delegate first_and_step(m::Model) => md

@delegate time_labels(m::Model) => md

# Return the number of timesteps a given component in a model will run for.
@delegate getspan(m::Model, comp_name::Symbol) => md

"""
    datumdef(comp_def::ComponentDef, item::Symbol)

Return a DatumDef for `item` in the given component `comp_def`.
"""
function datumdef(comp_def::AbstractComponentDef, item::Symbol)
    if has_variable(comp_def, item)
        return variable(comp_def, item)

    elseif has_parameter(comp_def, item)
        return parameter(comp_def, item)
    else
        error("Cannot access data item; :$item is not a variable or a parameter in component $(comp_def.comp_id).")
    end
end

datumdef(m::Model, comp_name::Symbol, item::Symbol) = datumdef(compdef(m.md, comp_name), item)

"""
    dim_names(m::Model, comp_name::Symbol, datum_name::Symbol)

Return the dimension names for the variable or parameter `datum_name`
in the given component `comp_name` in model `m`.
"""
function dim_names(m::Model, comp_name::Symbol, datum_name::Symbol)
    # the line below would work if the comp_name is in the top level of components in m's component structure
    # return dim_names(compdef(m, comp_name), datum_name)

    paths = _get_all_paths(m)
    comp_path = paths[comp_name]
    comp_def = find_comp(m, comp_path)
    return dim_names(comp_def, datum_name)
end

dim_names(mm::MarginalModel, comp_name::Symbol, datum_name::Symbol) = dim_names(mm.base, comp_name, datum_name)

@delegate dimension(m::Model, dim_name::Symbol) => md

# Allow access of the form my_model[:grosseconomy, :tfp]
@delegate Base.getindex(m::Model, comp_name::Symbol, datum_name::Symbol) => mi

"""
    dim_count(m::Model, dim_name::Symbol)

Return the size of index `dim_name` in model `m`.
"""
@delegate dim_count(m::Model, dim_name::Symbol) => md
@delegate dim_counts(m::Model, dims::Vector{Symbol}) => md
@delegate dim_count_dict(m::Model) => md

"""
    dim_keys(m::Model, dim_name::Symbol)

Return keys for dimension `dim-name` in model `m`.
"""
@delegate dim_keys(m::Model, dim_name::Symbol) => md
"""
     dim_keys(mi::ModelInstance, dim_name::Symbol)

 Return keys for dimension `dim-name` in model instance `mi`.
 """
 @delegate dim_keys(mi::ModelInstance, dim_name::Symbol) => md
"""
    dim_key_dict(m::Model)

Return a dict of dimension keys for all dimensions in model `m`.
"""
@delegate dim_key_dict(m::Model) => md
"""
    dim_values(m::Model, name::Symbol)

Return values for dimension `name` in Model `m`.
"""
@delegate dim_values(m::Model, name::Symbol) => md
"""
    dim_value_dict(m::Model)

Return a dictionary of the values of all dimensions in Model `m`.
"""
@delegate dim_value_dict(m::Model) => md
"""
    set_dimension!(m::Model, name::Symbol, keys::Union{Vector, Tuple, AbstractRange})

Set the values of `m` dimension `name` to integers 1 through `count`, if `keys`` is
an integer; or to the values in the vector or range if `keys`` is either of those types.
"""
@delegate set_dimension!(m::Model, name::Symbol, keys::Union{Int, Vector, Tuple, AbstractRange}) => md

@delegate check_parameter_dimensions(m::Model, value::AbstractArray, dims::Vector, name::Symbol) => md

@delegate parameter_names(m::Model, comp_name::Symbol) => md

@delegate parameter_dimensions(m::Model, comp_name::Symbol, param_name::Symbol) => md

@delegate parameter_unit(m::Model, comp_name::Symbol, param_name::Symbol) => md

parameter(m::Model, comp_def::ComponentDef, param_name::Symbol) = parameter(comp_def, param_name)

parameter(m::Model, comp_name::Symbol, param_name::Symbol) = parameter(m, compdef(m, comp_name), param_name)

"""
    parameters(m::Model, comp_name::Symbol)

Return a list of the parameter definitions for `comp_name` in model `m`.
"""
parameters(m::Model, comp_name::Symbol) = parameters(compdef(m, comp_name))

variable(m::Model, comp_name::Symbol, var_name::Symbol) = variable(compdef(m, comp_name), var_name)

@delegate variable_unit(m::Model, comp_path::ComponentPath, var_name::Symbol) => md

@delegate variable_dimensions(m::Model, comp_path::ComponentPath, var_name::Symbol) => md
@delegate variable_dimensions(m::Model, comp::Symbol, var_name::Symbol) => md
@delegate variable_dimensions(m::Model, comp_path::NTuple{N, Symbol} where N, var_name::Symbol) => md

"""
    variables(m::Model, comp_name::Symbol)

Return an iterator on the variable definitions for `comp_name` in model `m`.
"""
variables(m::Model, comp_name::Symbol) = variables(compdef(m, comp_name))

@delegate variable_names(m::Model, comp_name::Symbol) => md

"""
    set_external_array_param!(m::Model, name::Symbol, value::Union{AbstractArray, TimestepArray}, dims)

Add a one or two dimensional (optionally, time-indexed) array parameter `name`
with value `value` to the model `m`.
"""
@delegate set_external_array_param!(m::Model, name::Symbol, value::Union{AbstractArray, TimestepArray}, dims) => md

"""
    set_external_scalar_param!(m::Model, name::Symbol, value::Any)

Add a scalar type parameter `name` with value `value` to the model `m`.
"""
@delegate set_external_scalar_param!(m::Model, name::Symbol, value::Any) => md

"""
    delete!(m::Model, component::Symbol; deep::Bool=false)

Delete a `component` by name from a model `m`'s ModelDef, and nullify the ModelInstance.
If `deep=true` then any external model parameters connected only to
this component will also be deleted.
"""
@delegate Base.delete!(m::Model, comp_name::Symbol; deep::Bool=false) => md

"""
    delete_param!(m::Model, external_param_name::Symbol)

Delete `external_param_name` from a model `m`'s ModelDef's list of external parameters, and
also remove all external parameters connections that were connected to `external_param_name`.
"""
@delegate delete_param!(m::Model, external_param_name::Symbol) => md

"""
    set_param!(m::Model, comp_name::Symbol, param_name::Symbol, value; dims=nothing)

Set the parameter of a component `comp_name` in a model `m` to a given `value`.
The `value` can by a scalar, an array, or a NamedAray. Optional keyword argument 'dims'
is a list of the dimension names of the provided data, and will be used to check
that they match the model's index labels.
"""
@delegate set_param!(m::Model, comp_name::Symbol, param_name::Symbol, value; dims=nothing) => md

"""
    set_param!(m::Model, comp_name::Symbol, param_name::Symbol, ext_param_name::Symbol, value; dims=nothing)

Set the parameter `param_name` of a component `comp_name` in a model `m` to a given `value`,
storing the value in the model's external parameter list by the provided name `ext_param_name`.
The `value` can by a scalar, an array, or a NamedAray. Optional keyword argument 'dims'
is a list of the dimension names of the provided data, and will be used to check
that they match the model's index labels.
"""
@delegate set_param!(m::Model, comp_name::Symbol, param_name::Symbol, ext_param_name::Symbol, value; dims=nothing) => md

"""
    set_param!(m::Model, param_name::Symbol, value; dims=nothing)

Set the value of a parameter in all components of the model that have a parameter of
the specified name.
"""
@delegate set_param!(m::Model, param_name::Symbol, value; dims=nothing, ignoreunits::Bool=false) => md

@delegate import_params!(m::Model) => md

"""
    run(m::Model)

Run model `m` once.
"""
function Base.run(m::Model; ntimesteps::Int=typemax(Int), rebuild::Bool=false,
                  dim_keys::Union{Nothing, Dict{Symbol, Vector{T} where T <: DimensionKeyTypes}}=nothing)
    if length(m) == 0
        error("Cannot run a model with no components.")
    end

    if (rebuild || ! is_built(m))
        build!(m)
    end

    # println("Running model...")
    mi = modelinstance(m)
    run(mi, ntimesteps, dim_keys)
    nothing
end


text:
using MetabolicFluxAnalysis

reload("MetabolicFluxAnalysis")

#s_data = open("D:\\Thomas\\Downloads\\IsoDesign-v1.2.1\\influx_s-v2.11.1\\test\\e_coli.ftbl")
s_data = open("testdata/e_coli.ftbl")

ftbl_data = MetabolicFluxAnalysis.parse_ftbl_data(s_data)

fluxi_from_ftbl = MetabolicFluxAnalysis.ftbl_to_fluxi(ftbl_data)

MetabolicFluxAnalysis.get_wxflux_net(fluxi_from_ftbl.fluxidata.net,fluxi_from_ftbl.wx_transform,"tk1")
MetabolicFluxAnalysis.get_wxflux_xch(fluxi_from_ftbl.fluxidata.net,fluxi_from_ftbl.wx_transform,"tk1")
MetabolicFluxAnalysis.find_wx_flux(fluxi_from_ftbl.fluxidata.net,fluxi_from_ftbl.wx_transform,"tk1")

MetabolicFluxAnalysis.create_all_wx_fluxnames(fluxi_from_ftbl.fluxidata.net,fluxi_from_ftbl.wx_transform)[16]

all_wx_fluxes = MetabolicFluxAnalysis.create_all_wx_fluxnames(fluxi_from_ftbl.fluxidata.net,fluxi_from_ftbl.wx_transform)

# create wx parametrization-based flux polytope..
wx_constraints = [ fluxi_from_ftbl.wx_constraints_fluxes ; fluxi_from_ftbl.wx_constraints_eq_and_ineq ; fluxi_from_ftbl.wx_constraints_flux_measurements ]

for wci in wx_constraints
    println(wci)
end

#fp_wx_constraints = MetabolicFluxAnalysis.lc2polyhedron(fluxi_from_ftbl.fluxidata.net,fluxi_from_ftbl.wx_transform, wx_constraints[1:10] )
#fp_wx_bounded = MetabolicFluxAnalysis.Polyhedron( fp_wx_constraints , MetabolicFluxAnalysis.bounding_box(zeros(170),2*ones(170)) )
#fp_wx_rp = MetabolicFluxAnalysis.get_rounded_fulldimensional_parametrization( fp_wx_bounded ; remove_red_constr=1)

fp_wx = MetabolicFluxAnalysis.get_flux_polyhedron( fluxi_from_ftbl.fluxidata.net , fluxi_from_ftbl.wx_transform , wx_constraints ; max_val=20. )

for zi=1:132
    vi = 1.0*(1:132 .==zi )  #rand(132)
    result_i = MetabolicFluxAnalysis.solve_lp_minmax(fp_wx,vi./norm(vi))
    print("$(zi) : $(all_wx_fluxes[zi]) : $(result_i[1]) - $(result_i[2]) \n")
end

# TODO: make this work reliably, i.e. figure out what the problems are during parametrization..
if(false)
    fp_wx_rp = MetabolicFluxAnalysis.get_rounded_fulldimensional_parametrization( MetabolicFluxAnalysis.Polyhedron(fp_wx,MetabolicFluxAnalysis.bounding_box( -100*ones(132) , 100*ones(132) )  ) ; remove_red_constr=0)

    vix = rand(132,5000)
    vix = broadcast( / , vix , sqrt.(sum(vix,1)) )
    results_lp = MetabolicFluxAnalysis.solve_lp_minmax.([fp_wx],  [ slicedim( vix,2,i) for i in 1:1000 ] )
    minimum( map( ri -> ri[2] , results_lp) )

    for zi=1:10000
        vi = rand(132)
        result_i = MetabolicFluxAnalysis.solve_lp_minmax(fp_wx,vi./norm(vi))
        print("$(vi[1]) - $(vi[2]) \n")
    end
end

## analyze fluxes..

##

MetabolicFluxAnalysis.heuristically_find_infeasible_constraints( fp_wx ; relaxation_step = 0.1 , max_tries=500 )

find( fp_wx.A[35,:].!=0)
find( fp_wx.A[134-size(fp_wx.A,1),:].!=0)
find( fp_wx.A[141-size(fp_wx.A,1),:].!=0)

fp_wx.G[1,14]
fp_wx.h[1]

find( fp_wx.G[108,:].!=0 )

find( x -> x>2 , sum( [fp_wx.G;fp_wx.A] .!=0 , 1 ) )

find( x -> x>2 , sum( [fp_wx.G] .!=0 , 1 ) )


text:
module Term
# general utils
include("__text_utils.jl")
include("_ansi.jl")
include("_utils.jl")

# don't import other modules
include("measure.jl")
include("color.jl")
include("theme.jl")
include("highlight.jl")

# rely on other modules
include("markup.jl")
include("style.jl")
include("segment.jl")
include("macros.jl")

# renderables, rely heavily on other modules
include("box.jl")
include("console.jl")
include("renderables.jl")
include("layout.jl")
include("panel.jl")
include("errors.jl")
include("logging.jl")
include("tprint.jl")
include("progress.jl")
include("tree.jl")
include("logo.jl")
include("inspect.jl")

export RenderableText, Panel, TextBox
export Spacer, vLine, hLine
export theme, highlight
export inspect, typestree
export @red, @black, @green, @yellow, @blue, @magenta, @cyan, @white, @default
export @bold, @dim, @italic, @underline, @style
export tprint, tprintln
export install_stacktrace
export install_term_logger, uninstall_term_logger
export track
export Tree

# ----------------------------------- base ----------------------------------- #
using .measure: measure
using .measure: Measure

# ----------------------------------- style ---------------------------------- #
using .markup: extract_markup, MarkupTag, pairup_tags

using .color: NamedColor, BitColor, RGBColor, get_color

using .style: MarkupStyle, extract_style, apply_style

using .segment: Segment

"""
    Measure(seg::Segment)

gives the measure of a segment
"""
measure.Measure(seg::Segment) = seg.measure

"""
    Measure(segments::Vector{Segment})

gives the measure of a vector of segments
"""
function measure.Measure(segments::Vector{Segment})
    return Measure(
        max([seg.measure.w for seg in segments]...),
        sum([seg.measure.h for seg in segments]),
    )
end

# -------------------------------- renderables ------------------------------- #
using .box

using .consoles: Console, console, err_console, console_height, console_width

using .renderables: AbstractRenderable, Renderable, RenderableText

using .layout: Padding, vstack, hstack, Spacer, vLine, hLine, pad

using .panel: Panel, TextBox

# define additional methods for measure functions
measure.width(text::AbstractString) = Measure(text).w
measure.width(seg::Segment) = seg.measure.w
measure.width(ren::AbstractRenderable) = ren.measure.w

measure.height(text::AbstractString) = Measure(text).h
measure.height(seg::Segment) = seg.measure.h
measure.height(ren::AbstractRenderable) = ren.measure.h

# ---------------------------------- others ---------------------------------- #
using .errors: install_stacktrace

using .logging: install_term_logger, uninstall_term_logger, TermLogger

using .Tprint: tprint, tprintln

using .progress: ProgressBar, update, track

using .tree: Tree

using .introspection: inspect, typestree

end


text:
@testset "Constraint Stamp" begin

    # CStamp
    N = 10
    p = 4
    s0 = CStamp(:v, :col, 1, 2, 3)
    @test valid(s0, N, p)
    s0 = CStamp(:v, :col, 1, 1, 3)
    @test !valid(s0, N, p)
    s0 = CStamp(:v, :col, 1, 3, 3)
    @test valid(s0, N, p)
    s0 = CStamp(:v, :col, 1, 5, 3)
    @test !valid(s0, N, p)
    s0 = CStamp(:v, :col, 0, 3, 3)
    @test !valid(s0, N, p)
    s0 = CStamp(:v, :col, 1, 3, 1)
    @test !valid(s0, N, p)
    s0 = CStamp(:v, :col, 3, 1, 3)
    @test !valid(s0, N, p)
    s0 = CStamp(:v, :col, 1, 2, 11)
    @test !valid(s0, N, p)

    s0 = CStamp(:h, :col, 1, 3, 3)
    @test valid(s0, N, p)
    s0 = CStamp(:h, :col, 2, 2, 3)
    @test !valid(s0, N, p)
    s0 = CStamp(:h, :col, 3, 1, 3)
    @test valid(s0, N, p)
    s0 = CStamp(:h, :col, 3, 1, 1)
    @test !valid(s0, N, p)
    s0 = CStamp(:h, :col, 3, 1, 11)
    @test !valid(s0, N, p)
    s0 = CStamp(:h, :col, 5, 1, 10)
    @test !valid(s0, N, p)
    s0 = CStamp(:h, :col, 4, 1, 10)
    @test valid(s0, N, p)

    s1 = CStamp()
    @test s1 != s0
    @test !(s1 == s0)

    s2 = stampify(:v,:col,1,2,10)
    stampify!(s1,:v,:col,1,2,10)
    @test s2 == s1
    @test isequal(s2, s1)
    @test hash(s2) == hash(s1)

end