text:
module NonlinearUnscentedTest

using Test
using LinearAlgebra
using ForneyLab
using ForneyLab: outboundType, isApplicable, sigmaPointsAndWeights, Unscented, requiresBreaker, breakerParameters
using ForneyLab: SPNonlinearUTOutNG, SPNonlinearUTIn1GG, SPNonlinearUTOutNGX, SPNonlinearUTInGX, MNonlinearUTInGX
using ForneyLab: unscentedStatistics, smoothRTS, smoothRTSMessage, collectStatistics, marginalizeGaussianMV, concatenateGaussianMV, split, isMultiIn

f(x) = x

g(x::Float64) = x^2 - 5.0
g(x::Vector{Float64}) = x.^2 .- 5.0
g_inv(y::Float64) = sqrt(y + 5.0)
g_inv(y::Vector{Float64}) = sqrt.(y .+ 5.0)

h(x::Float64, y::Float64) = x^2 - y
h(x::Vector{Float64}, y::Vector{Float64}) = x.^2 .- y
h(x::Float64, y::Vector{Float64}) = x^2 .- y
h_inv_x(z::Float64, y::Float64) = sqrt(z + y)
h_inv_x(z::Vector{Float64}, y::Vector{Float64}) = sqrt.(z .+ y)

#--------
# Helpers
#--------

@testset "sigmaPointsAndWeights" begin
    (sigma_points, weights_m, weights_c) = sigmaPointsAndWeights(0.0, 1.0, alpha=1e-3)
    @test sigma_points == [0.0, 0.0010000000000143778, -0.0010000000000143778]
    @test weights_m == [-999998.9999712444, 499999.9999856222, 499999.9999856222]
    @test weights_c == [-999995.9999722444, 499999.9999856222, 499999.9999856222]

    (sigma_points, weights_m, weights_c) = sigmaPointsAndWeights([0.0], mat(1.0), alpha=1e-3)
    @test sigma_points == [[0.0], [0.0010000000000143778], [-0.0010000000000143778]]
    @test weights_m == [-999998.9999712444, 499999.9999856222, 499999.9999856222]
    @test weights_c == [-999995.9999722444, 499999.9999856222, 499999.9999856222]
end

@testset "unscentedStatistics" begin
    m_x = 0.0
    v_x = 1.0
    m_y = 2.0
    v_y = 3.0

    # Single univariate inbound
    (m_tilde, V_tilde, C_tilde) = unscentedStatistics(m_x, v_x, g)
    @test m_tilde == -4.0
    @test V_tilde == 1.9999999997671694
    @test C_tilde == 0.0

    # Multiple univariate inbounds
    (m_tilde, V_tilde, C_tilde) = unscentedStatistics([m_x, m_y], [v_x, v_y], h)
    @test m_tilde == -1.0000000000582077
    @test V_tilde == 5.0000009994837455
    @test C_tilde == [0.0, -2.9999999999442934]

    # Single multivariate inbound
    (m_tilde, V_tilde, C_tilde) = unscentedStatistics([m_x, m_y], Diagonal([v_x, v_y]), g)
    @test m_tilde == [-4.000000000465661, 1.9999999997380655]
    @test V_tilde == [2.000000997039024 5.999996995087713; 5.999996995087713 66.00000899611041]
    @test C_tilde == [0.0 0.0; 5.547917680814862e-11 12.000000000165528]

    # Multiple multivariate inbounds
    (m_tilde, V_tilde, C_tilde) = unscentedStatistics([[m_x, m_x], [m_y, m_y]], [Diagonal([v_x, v_x]), Diagonal([v_y, v_y])], h)
    @test m_tilde == [-1.000000000174623, -1.000000000174623]
    @test V_tilde == [5.000002998916898 1.9999989990901668; 1.999998999031959 5.000002998946002]
    @test C_tilde == [0.0 0.0; 0.0 0.0; -2.9999999999998863 0.0; 0.0 -2.9999999999998863]

    # Mixed variate inbounds
    (m_tilde, V_tilde, C_tilde) = unscentedStatistics([m_x, [m_y, m_y]], [v_x, Diagonal([v_y, v_y])], h)
    @test m_tilde == [-1.0, -1.0]
    @test V_tilde == [5.000002000015229 2.000002000015229; 2.000002000015229 5.000002000015229]
    @test C_tilde == [0.0 0.0; -3.0000000000001137 0.0; 0.0 -3.0000000000001137]
end

@testset "smoothRTSMessage" begin
    @test smoothRTSMessage(-4.0, 2.0, 1.0, 3.0, 5.0, 4.0, 6.0) == (43.0, 195.0)
    @test smoothRTSMessage([-4.0], mat(2.0), mat(1.0), [3.0], mat(5.0), [4.0], mat(6.0)) == ([43.0], mat(195.0))
end

@testset "smoothRTS" begin
    @test smoothRTS(-4.0, 2.0, 1.0, 3.0, 5.0, 4.0, 6.0) == (4.0, 4.875)
    @test smoothRTS([-4.0], mat(2.0), mat(1.0), [3.0], mat(5.0), [4.0], mat(6.0)) == ([4.0], mat(4.875))
end

@testset "collectStatistics" begin
    @test collectStatistics(Message(Univariate, GaussianMeanVariance, m=0.0, v=1.0), nothing, Message(Univariate, GaussianMeanVariance, m=2.0, v=3.0)) == ([0.0, 2.0], [1.0, 3.0])
    @test collectStatistics(Message(Multivariate, GaussianMeanVariance, m=[0.0], v=mat(1.0)), nothing, Message(Multivariate, GaussianMeanVariance, m=[2.0], v=mat(3.0))) == ([[0.0], [2.0]], [mat(1.0), mat(3.0)])
    @test collectStatistics(Message(Univariate, GaussianMeanVariance, m=0.0, v=1.0), nothing, Message(Multivariate, GaussianMeanVariance, m=[2.0], v=mat(3.0))) == ([0.0, [2.0]], [1.0, mat(3.0)])
end

@testset "marginalizeGaussianMV" begin
    @test marginalizeGaussianMV([0.0, 1.0], [2.0 0.5; 0.5 3.0], [(), ()], 1) == (0.0, 2.0) # Univariate
    @test marginalizeGaussianMV([0.0, 1.0, 2.0], [2.0 0.0 0.5; 0.0 3.0 0.0; 0.5 0.0 4.0], [(), (2,)], 1) == (0.0, 2.0) # Univariate
    @test marginalizeGaussianMV([0.0, 1.0, 2.0], [2.0 0.0 0.5; 0.0 3.0 0.0; 0.5 0.0 4.0], [(), (2,)], 2) == ([1.0, 2.0], [3.0 0.0; 0.0 4.0]) # Multivariate
end

@testset "concatenateGaussianMV" begin
    @test concatenateGaussianMV([1.0, 2.0, 3.0], [4.0, 5.0, 6.0]) == ([1.0, 2.0, 3.0], Diagonal([4.0, 5.0, 6.0]), [(), (), ()])
    @test concatenateGaussianMV([[1.0], [2.0, 3.0]], [mat(4.0), Diagonal([5.0, 6.0])]) == ([1.0, 2.0, 3.0], [4.0 0.0 0.0; 0.0 5.0 0.0; 0.0 0.0 6.0], [(1,), (2,)])
    @test concatenateGaussianMV([1.0, [2.0, 3.0]], [4.0, Diagonal([5.0, 6.0])]) == ([1.0, 2.0, 3.0], [4.0 0.0 0.0; 0.0 5.0 0.0; 0.0 0.0 6.0], [(), (2,)])
end

@testset "split" begin
    @test split([1.0, 2.0], [(), ()]) == [1.0, 2.0]
    @test split([1.0, 2.0, 3.0, 4.0], [(2,), (2,)]) == [[1.0, 2.0], [3.0, 4.0]]
    @test split([1.0, 2.0, 3.0], [(), (2,)]) == [1.0, [2.0, 3.0]]
    @test split([1.0, 2.0, 3.0, 4.0, 5.0, 6.0], [(2,), (2,2)]) == [[1.0, 2.0], [3.0 5.0; 4.0 6.0]]
end

@testset "requiresBreaker and breakerParameters" begin
    # Without given inverse
    fg = FactorGraph()
    x = Variable()
    y = Variable()
    nd = GaussianMeanVariance(x, 0.0, 1.0)
    Nonlinear{Unscented}(y, x, g=g)

    @test requiresBreaker(nd.i[:out])
    @test_throws Exception breakerParameters(nd.i[:out].partner)
    @test breakerParameters(nd.i[:out]) == (Message{GaussianMeanVariance, Univariate}, ())

    # With given inverse
    fg = FactorGraph()
    x = Variable()
    y = Variable()
    nd = GaussianMeanVariance(x, 0.0, 1.0)
    Nonlinear{Unscented}(y, x, g=g, g_inv=g)

    @test !requiresBreaker(nd.i[:out])
    @test_throws Exception breakerParameters(nd.i[:out])
end

@testset "isMultiIn" begin
    fg = FactorGraph()
    x = Variable()
    y = Variable()
    GaussianMeanVariance(x, 0.0, 1.0)
    nd = Nonlinear{Unscented}(y, x, g=g)
    Clamp(y, 1.0)
    @test !isMultiIn(nd)

    fg = FactorGraph()
    x = Variable()
    y = Variable()
    z = Variable()
    GaussianMeanVariance(x, 0.0, 1.0)
    GaussianMeanVariance(z, 0.0, 1.0)
    nd = Nonlinear{Unscented}(y, x, z, g=g)
    Clamp(y, 1.0)
    @test isMultiIn(nd)

    fg = FactorGraph()
    x = Variable()
    y = Variable()
    z = Variable()
    GaussianMeanVariance(x, 0.0, 1.0)
    Clamp(z, 1.0)
    nd = Nonlinear{Unscented}(y, x, z, g=g)
    Clamp(y, 1.0)
    @test !isMultiIn(nd)
end

#-------------
# Update rules
#-------------

@testset "SPNonlinearUTOutNG" begin
    @test SPNonlinearUTOutNG <: SumProductRule{Nonlinear{Unscented}}
    @test outboundType(SPNonlinearUTOutNG) == Message{GaussianMeanVariance}
    @test isApplicable(SPNonlinearUTOutNG, [Nothing, Message{GaussianMeanVariance}])

    @test ruleSPNonlinearUTOutNG(g, nothing, Message(Univariate, GaussianMeanVariance, m=2.0, v=3.0)) == Message(Univariate, GaussianMeanVariance, m=2.0000000001164153, v=66.00000000093132)
    @test ruleSPNonlinearUTOutNG(g, nothing, Message(Univariate, GaussianMeanVariance, m=2.0, v=3.0), alpha=1.0) == Message(Univariate, GaussianMeanVariance, m=2.0, v=66.0)
    @test ruleSPNonlinearUTOutNG(g, nothing, Message(Multivariate, GaussianMeanVariance, m=[2.0], v=mat(3.0))) == Message(Multivariate, GaussianMeanVariance, m=[2.0000000001164153], v=mat(66.00000000093132))
    @test ruleSPNonlinearUTOutNG(g, nothing, Message(Multivariate, GaussianMeanVariance, m=[2.0], v=mat(3.0)), alpha=1.0) == Message(Multivariate, GaussianMeanVariance, m=[2.0], v=mat(66.0))
end

@testset "SPNonlinearUTOutNGX" begin
    @test SPNonlinearUTOutNGX <: SumProductRule{Nonlinear{Unscented}}
    @test outboundType(SPNonlinearUTOutNGX) == Message{GaussianMeanVariance}
    @test !isApplicable(SPNonlinearUTOutNGX, [Nothing, Message{Gaussian}])
    @test isApplicable(SPNonlinearUTOutNGX, [Nothing, Message{Gaussian}, Message{Gaussian}])
    @test !isApplicable(SPNonlinearUTOutNGX, [Message{Gaussian}, Nothing, Message{Gaussian}])

    @test ruleSPNonlinearUTOutNGX(h, nothing, Message(Univariate, GaussianMeanVariance, m=2.0, v=3.0), Message(Univariate, GaussianMeanVariance, m=5.0, v=1.0)) == Message(Univariate, GaussianMeanVariance, m=1.9999999997671694, v=67.00000899657607)
    @test ruleSPNonlinearUTOutNGX(h, nothing, Message(Multivariate, GaussianMeanVariance, m=[2.0], v=mat(3.0)), Message(Multivariate, GaussianMeanVariance, m=[5.0], v=mat(1.0))) == Message(Multivariate, GaussianMeanVariance, m=[1.9999999997671694], v=mat(67.00000899657607))
end

@testset "SPNonlinearUTIn1GG" begin
    @test SPNonlinearUTIn1GG <: SumProductRule{Nonlinear{Unscented}}
    @test outboundType(SPNonlinearUTIn1GG) == Message{GaussianMeanVariance}
    @test isApplicable(SPNonlinearUTIn1GG, [Message{Gaussian}, Nothing])

    # Without given inverse
    @test ruleSPNonlinearUTIn1GG(g, Message(Univariate, GaussianMeanVariance, m=2.0, v=3.0), Message(Univariate, GaussianMeanVariance, m=2.0, v=1.0)) == Message(Univariate, GaussianMeanVariance, m=2.499999999868301, v=0.3125000002253504)
    @test ruleSPNonlinearUTIn1GG(g, Message(Univariate, GaussianMeanVariance, m=2.0, v=3.0), Message(Univariate, GaussianMeanVariance, m=2.0, v=1.0), alpha=1.0) == Message(Univariate, GaussianMeanVariance, m=2.5, v=0.3125)
    @test ruleSPNonlinearUTIn1GG(g, Message(Multivariate, GaussianMeanVariance, m=[2.0], v=mat(3.0)), Message(Multivariate, GaussianMeanVariance, m=[2.0], v=mat(1.0))) == Message(Multivariate, GaussianMeanVariance, m=[2.499999999868301], v=mat(0.31250000021807445))
    @test ruleSPNonlinearUTIn1GG(g, Message(Multivariate, GaussianMeanVariance, m=[2.0], v=mat(3.0)), Message(Multivariate, GaussianMeanVariance, m=[2.0], v=mat(1.0)), alpha=1.0) == Message(Multivariate, GaussianMeanVariance, m=[2.5], v=mat(0.3125))

    # With given inverse
    @test ruleSPNonlinearUTIn1GG(g, g_inv, Message(Univariate, GaussianMeanVariance, m=2.0, v=3.0), nothing) == Message(Univariate, GaussianMeanVariance, m=2.6255032138433307, v=0.10796282966583703)
    @test ruleSPNonlinearUTIn1GG(g, g_inv, Message(Univariate, GaussianMeanVariance, m=2.0, v=3.0), nothing, alpha=1.0) == Message(Univariate, GaussianMeanVariance, m=2.6251028535207217, v=0.10968772603524787)
    @test ruleSPNonlinearUTIn1GG(g, g_inv, Message(Multivariate, GaussianMeanVariance, m=[2.0], v=mat(3.0)), nothing) == Message(Multivariate, GaussianMeanVariance, m=[2.6255032138433307], v=mat(0.10796282966583703))
    @test ruleSPNonlinearUTIn1GG(g, g_inv, Message(Multivariate, GaussianMeanVariance, m=[2.0], v=mat(3.0)), nothing, alpha=1.0) == Message(Multivariate, GaussianMeanVariance, m=[2.6251028535207217], v=mat(0.10968772603524787))
end

@testset "SPNonlinearUTInGX" begin
    @test SPNonlinearUTInGX <: SumProductRule{Nonlinear{Unscented}}
    @test outboundType(SPNonlinearUTInGX) == Message{Gaussian}
    @test !isApplicable(SPNonlinearUTInGX, [Message{Gaussian}, Nothing])
    @test !isApplicable(SPNonlinearUTInGX, [Nothing, Message{Gaussian}, Message{Gaussian}])
    @test isApplicable(SPNonlinearUTInGX, [Message{Gaussian}, Nothing, Message{Gaussian}])

    # Without given inverse
    @test ruleSPNonlinearUTInGX(h, 1, Message(Univariate, GaussianMeanVariance, m=2.0, v=3.0), Message(Univariate, GaussianMeanVariance, m=2.0, v=1.0), Message(Univariate, GaussianMeanVariance, m=5.0, v=1.0)) == Message(Univariate, GaussianWeightedMeanPrecision, xi=6.666665554127903, w=2.6666662216903676)
    @test ruleSPNonlinearUTInGX(h, 1, Message(Multivariate, GaussianMeanVariance, m=[2.0], v=mat(3.0)), Message(Multivariate, GaussianMeanVariance, m=[2.0], v=mat(1.0)), Message(Multivariate, GaussianMeanVariance, m=[5.0], v=mat(1.0))) == Message(Multivariate, GaussianWeightedMeanPrecision, xi=[6.666665554127903], w=mat(2.666666221690368))

    # With given inverse
    @test ruleSPNonlinearUTInGX(h, h_inv_x, Message(Univariate, GaussianMeanVariance, m=2.0, v=3.0), nothing, Message(Univariate, GaussianMeanVariance, m=5.0, v=1.0)) == Message(Univariate, GaussianMeanVariance, m=2.6187538476660848, v=0.14431487274498522)
    @test ruleSPNonlinearUTInGX(h, h_inv_x, Message(Multivariate, GaussianMeanVariance, m=[2.0], v=mat(3.0)), nothing, Message(Multivariate, GaussianMeanVariance, m=[5.0], v=mat(1.0))) == Message(Multivariate, GaussianMeanVariance, m=[2.6187538476660848], v=mat(0.14431487274475785))
end

@testset "MNonlinearUTInGX" begin
    @test MNonlinearUTInGX <: MarginalRule{Nonlinear{Unscented}}
    @test isApplicable(MNonlinearUTInGX, [Nothing, Message{Gaussian}, Message{Gaussian}])

    @test ruleMNonlinearUTInGX(h, Message(Univariate, GaussianMeanVariance, m=2.0, v=3.0), Message(Univariate, GaussianMeanVariance, m=2.0, v=1.0), Message(Univariate, GaussianMeanVariance, m=5.0, v=1.0)) == ProbabilityDistribution(Multivariate, GaussianMeanVariance, m=[2.3636363470614055, 4.9090909132334355], v=[0.2727273058237252 0.1818181735464949; 0.18181817354649488 0.9545454566127697])
    @test ruleMNonlinearUTInGX(h, Message(Multivariate, GaussianMeanVariance, m=[2.0], v=mat(3.0)), Message(Multivariate, GaussianMeanVariance, m=[2.0], v=mat(1.0)), Message(Multivariate, GaussianMeanVariance, m=[5.0], v=mat(1.0))) == ProbabilityDistribution(Multivariate, GaussianMeanVariance, m=[2.3636363470609245, 4.909090913233555], v=[0.2727273058246874 0.18181817354625435; 0.18181817354625435 0.9545454566128299])
end

#------------
# Integration
#------------

@testset "Nonlinear integration via UT with given inverse" begin
    fg = FactorGraph()

    @RV x ~ GaussianMeanVariance(2.0, 1.0)
    @RV y ~ GaussianMeanVariance(2.0, 3.0)
    n = Nonlinear{Unscented}(y, x, g=g, g_inv=g_inv)

    @test isa(n, Nonlinear{Unscented})

    # Forward; g_inv should not be present in call
    pfz = PosteriorFactorization(fg)
    algo = messagePassingAlgorithm(y)
    code = algorithmSourceCode(algo)
    @test occursin("ruleSPNonlinearUTOutNG(g, nothing, messages[2])", code)
    @test !occursin("g_inv", code)

    # Backward; g_inv should be present in call
    pfz = PosteriorFactorization(fg)
    algo = messagePassingAlgorithm(x)
    code = algorithmSourceCode(algo)
    @test occursin("ruleSPNonlinearUTIn1GG(g, g_inv, messages[2], nothing)", code)
end

@testset "Multi-argument nonlinear integration via UT" begin
    fg = FactorGraph()

    @RV x ~ GaussianMeanVariance(2.0, 1.0)
    @RV y ~ GaussianMeanVariance(2.0, 3.0)
    @RV z ~ GaussianMeanVariance(5.0, 1.0)
    n = Nonlinear{Unscented}(y, x, z, g=h, g_inv=[h_inv_x, nothing])

    # Forward; h_inv_x should not be present in call
    pfz = PosteriorFactorization(fg)
    algo = messagePassingAlgorithm(y)
    code = algorithmSourceCode(algo)
    @test occursin("ruleSPNonlinearUTOutNGX(h, nothing, messages[3], messages[1])", code)
    @test !occursin("h_inv_x", code)

    # Backward with given inverse; h_inv_x should be present in call
    pfz = PosteriorFactorization(fg)
    algo = messagePassingAlgorithm(x)
    code = algorithmSourceCode(algo)
    @test occursin("ruleSPNonlinearUTInGX(h, h_inv_x, messages[3], nothing, messages[1])", code)

    # Backward without given inverse
    pfz = PosteriorFactorization(fg)
    algo = messagePassingAlgorithm(z)
    code = algorithmSourceCode(algo)
    @test occursin("ruleSPNonlinearUTInGX(h, 2, messages[3], messages[2], messages[1])", code)
    @test !occursin("h_inv_x", code)
    @test occursin("messages[1] = Message(vague(GaussianMeanVariance))", code)
end

@testset "Nonlinear integration via UT with given alpha" begin
    fg = FactorGraph()

    @RV x ~ GaussianMeanVariance(2.0, 1.0)
    @RV y ~ GaussianMeanVariance(2.0, 3.0)
    n = Nonlinear{Unscented}(y, x, g=g, alpha=1.0)

    # Forward; alpha should be present in call
    pfz = PosteriorFactorization(fg)
    algo = messagePassingAlgorithm(y)
    code = algorithmSourceCode(algo)
    @test occursin("ruleSPNonlinearUTOutNG(g, nothing, messages[1], alpha=1.0)", code)
end

@testset "Nonlinear integration via UT without given inverse" begin
    fg = FactorGraph()

    @RV x ~ GaussianMeanVariance(2.0, 1.0)
    @RV y ~ GaussianMeanVariance(2.0, 3.0)
    n = Nonlinear{Unscented}(y, x, g=g)

    # Forward; g_inv should not be present in call
    pfz = PosteriorFactorization(fg)
    algo = messagePassingAlgorithm(y)
    code = algorithmSourceCode(algo)
    @test occursin("ruleSPNonlinearUTOutNG(g, nothing, messages[1])", code)
    @test !occursin("g_inv", code)

    # Backward; g_inv should not be present in call,
    # both messages should be required, and initialization should take place
    pfz = PosteriorFactorization(fg)
    algo = messagePassingAlgorithm(x)
    code = algorithmSourceCode(algo)
    @test occursin("ruleSPNonlinearUTIn1GG(g, messages[2], messages[1])", code)
    @test !occursin("g_inv", code)
    @test occursin("messages[1] = Message(vague(GaussianMeanVariance))", code)
end

end #module


text:
## Feature extraction
# All features calculated relative to current player to move

"""
Discretize a 2D integer matrix into `maxval` binary layers
Layers correspond to == 1, == 2, ... == (maxval - 1), >= maxval
"""
function discretize(counts::Array{Int, 2}; maxval=8)
    out = BitArray(size(counts)..., maxval)
    for i in 1:(maxval-1)
        out[:, :, i] = counts .== i
    end
    out[:, :, maxval] = counts .>= maxval
    out
end

#
function stone_color(board::Board)
    result = BitArray(N, N, 3)
    cp = current_player(board)
    result[:, :, 1] = board.board .== cp
    result[:, :, 2] = board.board .== -cp
    result[:, :, 3] = board.board .== EMPTY
    result
end

function ones(board::Board)
    fill!(BitArray(N, N), true)
end

function zeros(board::Board)
    fill!(BitArray(N, N), false)
end

function turns_since(board::Board; maxval=8)
    discretize(board.cmove - board.order, maxval=maxval)
end

function liberties(board::Board; maxval=3)
    liberties = Array{Int}(N,N)
    for x in 1:N
        for y in 1:N
            point = Point(y,x)
            group = board.groups[point]
            liberties[point] = length(group.liberties)
        end
    end
    discretize(liberties, maxval=maxval)
end

"Encode how many liberty counts broken down by stone color"
function perplayer_liberties(board::Board, maxval=3)
    cp = current_player(board)
    cur_liberties = Array{Int}(N,N)
    op_liberties = Array{Int}(N,N)
    fill!(cur_liberties, 0)
    fill!(op_liberties, 0)
    for x in 1:N
        for y in 1:N
            point = Point(y,x)
            group = board.groups[point]
            count = length(group.liberties)
            if board[point] == cp
                cur_liberties[point] = count
            elseif board[point] == -cp
                op_liberties[point] = count
            end
        end
    end
    cat(3, discretize(cur_liberties, maxval=3), discretize(op_liberties, maxval=3))
end

# Liberties after move
# Capture size of move
# self atari size of move
function after_move_features(board::Board; maxval=8)
    legal = legal_moves(board)
    liberties = zeros(Int, N, N)
    groupsize = zeros(Int, N, N)
    capture_size = zeros(Int, N, N)
    color = current_player(board)
    for move in legal
        group = board.groups[move]
        (friends, foes, empties) = get_friend_foe(board, move, color)
        liberties[move] = length(union(group.liberties, [f.liberties for f in friends]...))
        groupsize[move] = 1 + sum([length(f.members) for f in friends])
        for foe in foes
            if length(foe.liberties) == 1
                capture_size[move] = length(foe.members)
            end
        end
    end
    cat(3, discretize(liberties, maxval=maxval),
        discretize(groupsize, maxval=maxval),
        discretize(capture_size, maxval=maxval))
end

function ladder_capture(board::Board)
end

function ladder_escape(board::Board)
end

function sensibleness(board::Board)
end

function player_color(board::Board)
    is_black = current_player(board) == BLACK
    fill!(BitArray(N, N), is_black)
end

#TODO: Provide features memory to write into to avoid unnecessary copies
function get_features(board::Board; features::Vector{Function}=DEFAULT_FEATURES)
    processed = Vector{BitArray}()
    for feature in features
        push!(processed, feature(board))
    end
    return cat(3, processed...)
end

# Run feature extraction on an empty board to get the dimensionality
function get_input_size(features::Vector{Function})
    size(get_features(Board(), features=features))
end

# TODO: Make a FeatureSet type to hold this
LIBERTIES = [perplayer_liberties]
DEFAULT_FEATURES = [stone_color, liberties]
ALL_FEATURES = [stone_color, liberties, ones, zeros, turns_since, after_move_features, player_color]

ext: jl


text:
const EMPTYSTRING = ""

##################################################
#
##################################################
## Getters
MOI.get(model::Optimizer, ::MOI.NumberOfConstraints) = ctx_m(model.ctx)
MOI.get(model::Optimizer, ::MOI.NumberOfConstraints{MOI.SingleVariable, MOI.ZeroOne}) =
    rhp_get_nb_bin_var(model.ctx)
MOI.get(model::Optimizer, ::MOI.NumberOfConstraints{MOI.SingleVariable, MOI.Integer}) =
    rhp_get_nb_int_var(model.ctx)
MOI.get(model::Optimizer, ::MOI.NumberOfConstraints{MOI.SingleVariable, S}) where S <: SS =
    rhp_get_nb_var(model.ctx, S)
# TODO: a bit hacky, but that should work for MOI Test.
MOI.get(model::Optimizer, ::MOI.NumberOfConstraints{VAF, MOI.Nonnegatives}) =
    sum(typeof.(collect(keys(model.vaf_mapping))) .== MOI.ConstraintIndex{VAF, MOI.Nonnegatives})
MOI.get(model::Optimizer, ::MOI.NumberOfConstraints{VAF, MOI.Nonpositives}) =
    sum(typeof.(collect(keys(model.vaf_mapping))) .== MOI.ConstraintIndex{VAF, MOI.Nonpositives})
MOI.get(model::Optimizer, ::MOI.NumberOfConstraints{VAF, MOI.Zeros}) =
    sum(typeof.(collect(keys(model.vaf_mapping))) .== MOI.ConstraintIndex{VAF, MOI.Zeros})
MOI.get(model::Optimizer, ::MOI.NumberOfConstraints{VAF, MOI.SecondOrderCone}) =
    sum(typeof.(collect(keys(model.vov_mapping))) .== MOI.ConstraintIndex{VAF, MOI.SecondOrderCone})
MOI.get(model::Optimizer, ::MOI.NumberOfConstraints{VOV, T}) where T <: VLS =
    sum(typeof.(collect(keys(model.vov_mapping))) .== MOI.ConstraintIndex{VOV, T})
MOI.get(model::Optimizer, ::MOI.NumberOfConstraints{VOV, S}) where S <: Union{MOI.SOS1{Float64}, MOI.SOS2{Float64}} =
    sum(typeof.(collect(keys(model.vov_mapping))) .== MOI.ConstraintIndex{VOV, S})
MOI.get(model::Optimizer, ::MOI.NumberOfConstraints{MOI.ScalarAffineFunction{Float64}, S}) where S <: LS  =
    rhp_get_nb_lequ(model.ctx, S)
MOI.get(model::Optimizer, ::MOI.NumberOfConstraints{MOI.ScalarQuadraticFunction{Float64}, S}) where S <: LS  =
    sum(typeof.(collect(keys(model.quadfn_mapping))) .== MOI.ConstraintIndex{MOI.ScalarQuadraticFunction{Float64}, S})

function MOI.get(model::Optimizer, ::MOI.ConstraintSet, ci::MOI.ConstraintIndex{VOV, S}) where S <: Union{MOI.SOS1{Float64}, MOI.SOS2{Float64}}
    return model.sos_sets[ci]
end

function MOI.get(model::Optimizer, ::MOI.VariableName, ci::MOI.ConstraintIndex)
    check_inbounds(model, ci)
    return ctx_getequname(model.ctx, ci.value-1)
end

_chk_citype = Dict(
     MOI.LessThan{Float64} => (false, true),
     MOI.GreaterThan{Float64} => (true, false),
     MOI.EqualTo{Float64} => (true, true),
     MOI.Interval{Float64} => (true, true),
    )

_cone_moi_to_rhp = Dict(
                        MOI.LessThan{Float64} => CONE_R_MINUS,
                        MOI.GreaterThan{Float64} => CONE_R_PLUS,
                        MOI.EqualTo{Float64} => CONE_0
                       )

function chk_equ_citype(ctx, eidx, ::Type{MOI.ConstraintIndex{MOI.SingleVariable, S}}) where S<:SS
    (lb, ub) = ctx_getvarbounds(ctx, eidx)
    has_lb = isfinite(lb)
    has_ub = isfinite(ub)
    return (has_lb, has_ub) == _chk_citype[S]
end

function chk_equ_citype(ctx, eidx, ::Type{MOI.ConstraintIndex{F, S}}) where  {F<:SF, S<:LS}
    type, cone = reshop_get_equtype(ctx, eidx)
    res = (type == 2) && (cone == get(_cone_moi_to_rhp, S, CONE_NONE))
    return res
end

function MOI.get(model::Optimizer, ci_type::Type{MOI.ConstraintIndex{MOI.SingleVariable, S}}, cstr::String) where S<:SS
    eidx = ctx_getequbyname(model.ctx, cstr)
    if eidx >= 0
        if (!chk_equ_citype(model.ctx, eidx, ci_type))
            return nothing
        end
        return MOI.ConstraintIndex{MOI.SingleVariable, S}(eidx+1)
    elseif eidx == -2
        error("Multiple variables have the name $(name).")
    else
        return nothing
    end
end

function MOI.get(model::Optimizer, ci_type::Type{MOI.ConstraintIndex{F, S}}, cstr::String) where {F<:SF, S<:LS}
    eidx = ctx_getequbyname(model.ctx, cstr)
    if eidx >= 0
        chk_equ_citype(model.ctx, eidx, ci_type)
        return MOI.ConstraintIndex{F, S}(eidx+1)
    elseif eidx == -2
        error("Multiple variables have the name $(name).")
    else
        return nothing
    end
end


text:
#=
@static if isdefined(Base.Experimental, Symbol("@optlevel"))
    Base.Experimental.@optlevel 0
end
=#

struct ApplicationError <: Exception
    showerror::Any
    summary::String
end

Base.showerror(io::IO, e::ApplicationError) = e.showerror(io)

function Base.show(io::IO, e::ApplicationError)
    show(io, ApplicationError)
    print(io, "(")
    show(io, e.summary)
    print(io, ", …)")
end

function (cli::Types.CLI)(args...; verbose::Integer = 1, kwargs...)
    # TODO: set logger level using `verbose`
    ans = try
        cli.api(args...; kwargs...)
    catch err
        if err isa ApplicationError
            showerror(stderr, err)
            println(stderr)
            exit(1)
        else
            rethrow()
        end
    end
    if ans !== nothing
        if verbose > 0
            display(ans)
            println()
        end
    end
end

function check_toplevel(dir)
    cmd = setenv(`git rev-parse --show-toplevel`; dir = dir)
    toplevel = try
        read(cmd, String)
    catch err
        @debug("Failed to call `read(cmd, String)` in `check_toplevel($(repr(dir)))`", cmd)
        err isa Union{Base.IOError,ProcessFailedException} || rethrow()
        nothing
    end
    if toplevel !== nothing
        realpath(chomp(toplevel)) == realpath(dir) && return
    end
    git_dir = joinpath(dir, ".git")
    git_dir_exist = isfile(git_dir) || isdir(git_dir)
    err = ApplicationError("failed: check_toplevel") do io
        print(io, "Not a Git repository: ", dir)
        if !git_dir_exist
            println(io)
            print(io, "`.git` directory does not exist")
        end
    end
    throw(err)
end

default_message() = "Update @$(gethostname())"

function git_is_clean(dir)
    output = read(setenv(`git status --porcelain`; dir = dir), String)
    return all(isspace, output)
end

function commit_impl(dir; message::AbstractString)
    git(cmd) = setenv(`$git $cmd`; dir = dir)
    check_toplevel(dir)
    if git_is_clean(dir)
        @info "Git repository `$dir` is clean. Nothing to commit."
        return false
    end
    run(git(`add .`))
    run(git(`commit --allow-empty-message --message $message`))
    return true
end

function API.commit(dir; message::AbstractString = default_message())
    commit_impl(dir; message = message)
    return
end

function API.upload(dir; merge::Bool = true, message::AbstractString = default_message())
    git(cmd) = setenv(`$git $cmd`; dir = dir)
    commit_impl(dir; message = message) || return
    if merge
        run(git(`pull --strategy=ours --no-edit`))
    end
    run(git(`push`))
    return
end

function eachstring(f, buf::AbstractVector{UInt8})
    i = firstindex(buf)
    for j in eachindex(buf)
        if iszero(@inbounds buf[j])
            f(String(buf[i:j-1]))
            i = j + 1
        end
    end
    if i < lastindex(buf)
        f(String(buf[i:end]))
    end
    return
end

function API.empty(dir)
    check_toplevel(dir)
    git(cmd) = setenv(`$git $cmd`; dir = dir)
    buf = read(git(`ls-files -z`))
    eachstring(buf) do path
        fullpath = joinpath(dir, path)
        isfile(fullpath) || return
        println(stderr, "Remove: ", path)
        rm(fullpath)
    end
    return
end


text:
@testset "Monodromy" begin
    function toric_ed(A)
        d, n = size(A)
        @var t[1:d] y[1:n] u[1:n]

        φ = map(j -> prod(i -> t[i]^A[i, j], 1:d), 1:n)
        Dφ = [differentiate(φ[i], t[j]) for i = 1:n, j = 1:d]

        System([φ + y - u; Dφ' * y], parameters = u)
    end

    @testset "monodromy_solve" begin
        F = toric_ed([3 2 1 0; 0 1 2 3])

        result = monodromy_solve(
            F,
            target_solutions_count = 21,
            max_loops_no_progress = 20,
            threading = false,
        )
        @test is_success(result)
        @test length(solutions(result)) == 21
        @test isempty(multiplicities(solutions(result)))
        @test isempty(sprint(show, result)) == false

        # test seed
        result2 = monodromy_solve(
            F,
            target_solutions_count = 21,
            max_loops_no_progress = 20,
            threading = false,
            seed = result.seed,
        )
        @test result2.statistics.tracked_loops[] == result.statistics.tracked_loops[]

        result = monodromy_solve(
            F,
            target_solutions_count = 21,
            max_loops_no_progress = 20,
            threading = true,
        )
        @test is_success(result)
        @test length(solutions(result)) == 21

        # test that timeout works
        result_timeout = monodromy_solve(F, target_solutions_count = 21, timeout = 1e-12)
        @test length(solutions(result_timeout)) < 21

        # test input of length > 1
        x₀, p₀ = find_start_pair(F)
        result = monodromy_solve(
            F,
            [x₀ for _ = 1:30],
            p₀,
            target_solutions_count = 21,
            max_loops_no_progress = 20,
        )
        @test length(solutions(result)) == 21

        # test with false input
        @test_logs (:warn, "None of the provided solutions is a valid start solution.") (
            result = monodromy_solve(F, [rand(6) for _ = 1:10], p₀)
        )
        @test result.returncode == :invalid_startvalue
        result =
            monodromy_solve(F.expressions, [x₀, rand(6)], p₀, parameters = F.parameters)
        @test length(solutions(result)) == 21

        # different distance function
        result = monodromy_solve(F, x₀, p₀, distance = (x, y) -> 0.0)
        @test length(solutions(result)) == 1

        # Test stop heuristic with no target solutions count
        result = monodromy_solve(F, x₀, p₀)
        @test is_heuristic_stop(result)

        roots_of_unity(s) = begin
            t = cis(π * 2 / 3)
            t² = t * t
            (vcat(t * s[1], t * s[2], s[3:end]), vcat(t² * s[1], t² * s[2], s[3:end]))
        end
        result = monodromy_solve(
            F,
            x₀,
            p₀,
            target_solutions_count = 21,
            max_loops_no_progress = 100,
            equivalence_classes = false,
            group_action = roots_of_unity,
        )
        @test length(solutions(result)) == 21

        # equivalence classes
        result = monodromy_solve(
            F,
            x₀,
            p₀,
            equivalence_classes = true,
            target_solutions_count = 7,
            max_loops_no_progress = 20,
            group_actions = roots_of_unity,
        )
        @test nresults(result) == 7
        # Test that equivalence classes are on by default if we supply a group action
        result = monodromy_solve(
            F,
            x₀,
            p₀,
            group_action = roots_of_unity,
            max_loops_no_progress = 20,
        )
        @test nsolutions(result) == 7

        # AbstractSystem as input
        result = monodromy_solve(
            InterpretedSystem(F),
            group_action = roots_of_unity,
            target_solutions_count = 7,
            max_loops_no_progress = 20,
        )
        @test nsolutions(result) == 7

        # test reuse strategies
        result = monodromy_solve(
            F,
            x₀,
            p₀,
            group_action = roots_of_unity,
            target_solutions_count = 7,
            reuse_loops = :all,
            max_loops_no_progress = 200,
        )
        @test nresults(result) == 7

        result = monodromy_solve(
            F,
            x₀,
            p₀,
            group_action = roots_of_unity,
            target_solutions_count = 7,
            reuse_loops = :random,
            max_loops_no_progress = 200,
        )
        @test nresults(result) == 7

        result = monodromy_solve(
            F,
            x₀,
            p₀,
            group_action = roots_of_unity,
            target_solutions_count = 7,
            reuse_loops = :none,
            max_loops_no_progress = 200,
        )
        @test nresults(result) == 7
    end

    @testset "Method of moments" begin
        f = moments3()
        S₃ = SymmetricGroup(3)
        relabeling(v) = map(p -> [v[1:3][p]..., v[4:6][p]..., v[7:9][p]...], S₃)

        R = monodromy_solve(
            f;
            group_action = relabeling,
            show_progress = false,
            max_loops_no_progress = 1,
        )
        @test nsolutions(R) ≤ 225

        for threading in [true, false]
            R = monodromy_solve(
                f;
                group_action = relabeling,
                show_progress = false,
                max_loops_no_progress = 20,
                target_solutions_count = 225,
                threading = false,
            )
            @test nsolutions(R) == 225
        end
    end

    @testset "Projective + Group Actions" begin
        @var a[1:2] x[1:2] s[1:2] z m[0:5]

        f0 = a[1] + a[2]
        f1 = a[1] * x[1] + a[2] * x[2]
        f2 = a[1] * (x[1]^2 + s[1] * z) + a[2] * (x[2]^2 + s[2] * z)
        f3 = a[1] * (x[1]^3 + 3 * s[1] * x[1] * z) + a[2] * (x[2]^3 + 3 * s[2] * x[2] * z)
        f4 =
            a[1] * (x[1]^4 + 6 * s[1] * x[1]^2 * z + 3 * s[1]^2 * z^2) +
            a[2] * (x[2]^4 + 6 * s[2] * x[2]^2 * z + 3 * s[2]^2 * z^2)
        f5 =
            a[1] * (x[1]^5 + 10 * s[1] * x[1]^3 * z + 15 * x[1] * s[1]^2 * z^2) +
            a[2] * (x[2]^5 + 10 * s[2] * x[2]^3 * z + 15 * x[2] * s[2]^2 * z^2)

        M2 = System(
            [f0, f1, f2, f3, f4, f5] - m .* z .^ (1:6),
            variables = [a; x; s; z],
            parameters = m,
        )

        relabeling = let S₂ = SymmetricGroup(2)
            v -> map(p -> [v[1:2][p]..., v[3:4][p]..., v[5:6][p]..., v[7]], S₂)
        end

        R = monodromy_solve(
            M2;
            group_action = relabeling,
            show_progress = false,
            threading = false,
            max_loops_no_progress = 10,
        )

        @test nsolutions(R) == 9
    end

    @testset "permutations" begin
        @var x[1:2] a b c
        p₁ = x - [2; 0]
        p₂ = x - [-2; 0]
        c₁ = p₁[1]^2 + p₁[2]^2 - 1
        c₂ = p₂[1]^2 + p₂[2]^2 - 1

        F = System([c₁ * c₂; a * x[1] + b * x[2] - c]; parameters = [a, b, c])
        S = monodromy_solve(
            F,
            [[1.0, 0.0]],
            [1, 1, 1],
            permutations = true,
            max_loops_no_progress = 20,
        )
        A = permutations(S)
        B = permutations(S, reduced = false)

        @test size(A) == (2, 2)
        @test A == [1 2; 2 1] || A == [2 1; 1 2]
        @test size(B, 1) == 2
        @test size(B, 2) > 2

        start_sols = solve(F; target_parameters = [1, 1, 1], start_system = :total_degree)
        S = monodromy_solve(F, solutions(start_sols), [1, 1, 1]; permutations = true)
        C = permutations(S)
        @test size(C, 1) == 4
    end

    @testset "Linear subspaces" begin
        @var x[1:4]
        f1 = rand_poly(x, 6)
        F = System([f1], x)
        res = monodromy_solve(F; dim = 3, threading = false)
        @test nsolutions(res) == 6
        @test trace(res) < 1e-6
        @test is_success(res)

        res = monodromy_solve(F; dim = 3, trace_test = false, threading = false)
        @test nsolutions(res) == 6
        @test is_heuristic_stop(res)

        res = monodromy_solve(
            F;
            dim = 3,
            trace_test = true,
            trace_test_tol = 1e-50,
            threading = false,
        )
        @test nsolutions(res) == 6
        @test is_heuristic_stop(res)

        @var x[1:4]
        f1 = rand_poly(x, 6; homogeneous = true)
        F = System([f1], x)
        res = monodromy_solve(F; dim = 2, compile = false)
        @test nsolutions(res) == 6
        @test trace(res) < 1e-6
        @test is_success(res)

        F = System([f1, rand_poly(x, 3), rand_poly(x, 4)], x)
        res = monodromy_solve(F; codim = 3, compile = false)
        @test nsolutions(res) == 72
        @test trace(res) < 1e-6
        @test is_success(res)
    end

    @testset "Monodromy rational functions" begin
        @var x y
        @var u[1:4]
        f1 = u[1] / x^2 + u[2]
        f2 = u[3] / y^2 + u[4]
        F = System([f1, f2], variables = [x, y], parameters = u[1:4])
        r1 = monodromy_solve(
            F,
            compile = true,
            target_solutions_count = 4,
            max_loops_no_progress = 100,
        )
        r2 = monodromy_solve(
            F,
            compile = false,
            target_solutions_count = 4,
            max_loops_no_progress = 100,
        )
        @test nsolutions(r1) == nsolutions(r2) == 4

        @var A1[1:3, 1:4] A2[1:3, 1:4]
        @var x[1:3] u1[1:2] u2[1:2]
        y1 = A1 * [x; 1]
        y2 = A2 * [x; 1]
        f = sum((u1 - y1[1:2] ./ y1[3]) .^ 2) + sum((u2 - y2[1:2] ./ y2[3]) .^ 2)
        F = System(
            differentiate(f, x),
            variables = x,
            parameters = [u1; u2; vec(A1); vec(A2)],
        )
        r = monodromy_solve(
            F,
            compile = false,
            max_loops_no_progress = 100,
            target_solutions_count = 6,
        )
        @test nsolutions(r) == 6
    end

    @testset "parameter homotopy with monodromy result" begin
        F = toric_ed([3 2 1 0; 0 1 2 3])
        mres = monodromy_solve(
            F,
            target_solutions_count = 21,
            max_loops_no_progress = 20,
            threading = false,
        )
        @test nsolutions(solve(F, mres, target_parameters = randn(ComplexF64, 4))) == 21
    end

    @testset "Verify solution_completeness" begin
        @var x y a b c
        f = x^2 + y^2 - 1
        l = a * x + b * y + c
        sys = System([f, l]; parameters = [a, b, c])
        res = solve(sys, [-0.6 - 0.8im, -1.2 + 0.4im]; target_parameters = [1, 2, 3])
        @test verify_solution_completeness(sys, solutions(res), [1, 2, 3])
        @test verify_solution_completeness(
            sys,
            solutions(res),
            [1, 2, 3],
            show_progress = false,
        )
        @test verify_solution_completeness(
            sys,
            solutions(res),
            [1, 2, 3],
            monodromy_options = (compile = false,),
            parameter_homotopy_options = (compile = false,),
            show_progress = false,
        )
        @test false == verify_solution_completeness(
            sys,
            solutions(res)[1:1],
            [1, 2, 3],
            monodromy_options = (compile = false,),
            parameter_homotopy_options = (compile = false,),
            show_progress = true,
        )
        @test false == verify_solution_completeness(
            sys,
            solutions(res),
            [1, 2, 3],
            monodromy_options = (compile = false,),
            parameter_homotopy_options = (compile = false,),
            show_progress = false,
            trace_tol = 1e-60,
        )
    end
end


text:
# As you prepare to give the amphipods your solution, you notice that the
# diagram they handed you was actually folded up. As you unfold it, you
# discover an extra part of the diagram.
#
# Between the first and second lines of text that contain amphipod starting
# positions, insert the following lines:
#
#   #D#C#B#A#
#   #D#B#A#C#
#
# So, the above example now becomes:
#
# #############
# #...........#
# ###B#C#B#D###
#   #D#C#B#A#
#   #D#B#A#C#
#   #A#D#C#A#
#   #########
#
# The amphipods still want to be organized into rooms similar to before:
#
# #############
# #...........#
# ###A#B#C#D###
#   #A#B#C#D#
#   #A#B#C#D#
#   #A#B#C#D#
#   #########
#
# In this updated example, the least energy required to organize these
# amphipods is 44169:
#
# #############
# #...........#
# ###B#C#B#D###
#   #D#C#B#A#
#   #D#B#A#C#
#   #A#D#C#A#
#   #########
#
# #############
# #..........D#
# ###B#C#B#.###
#   #D#C#B#A#
#   #D#B#A#C#
#   #A#D#C#A#
#   #########
#
# #############
# #A.........D#
# ###B#C#B#.###
#   #D#C#B#.#
#   #D#B#A#C#
#   #A#D#C#A#
#   #########
#
# #############
# #A........BD#
# ###B#C#.#.###
#   #D#C#B#.#
#   #D#B#A#C#
#   #A#D#C#A#
#   #########
#
# #############
# #A......B.BD#
# ###B#C#.#.###
#   #D#C#.#.#
#   #D#B#A#C#
#   #A#D#C#A#
#   #########
#
# #############
# #AA.....B.BD#
# ###B#C#.#.###
#   #D#C#.#.#
#   #D#B#.#C#
#   #A#D#C#A#
#   #########
#
# #############
# #AA.....B.BD#
# ###B#.#.#.###
#   #D#C#.#.#
#   #D#B#C#C#
#   #A#D#C#A#
#   #########
#
# #############
# #AA.....B.BD#
# ###B#.#.#.###
#   #D#.#C#.#
#   #D#B#C#C#
#   #A#D#C#A#
#   #########
#
# #############
# #AA...B.B.BD#
# ###B#.#.#.###
#   #D#.#C#.#
#   #D#.#C#C#
#   #A#D#C#A#
#   #########
#
# #############
# #AA.D.B.B.BD#
# ###B#.#.#.###
#   #D#.#C#.#
#   #D#.#C#C#
#   #A#.#C#A#
#   #########
#
# #############
# #AA.D...B.BD#
# ###B#.#.#.###
#   #D#.#C#.#
#   #D#.#C#C#
#   #A#B#C#A#
#   #########
#
# #############
# #AA.D.....BD#
# ###B#.#.#.###
#   #D#.#C#.#
#   #D#B#C#C#
#   #A#B#C#A#
#   #########
#
# #############
# #AA.D......D#
# ###B#.#.#.###
#   #D#B#C#.#
#   #D#B#C#C#
#   #A#B#C#A#
#   #########
#
# #############
# #AA.D......D#
# ###B#.#C#.###
#   #D#B#C#.#
#   #D#B#C#.#
#   #A#B#C#A#
#   #########
#
# #############
# #AA.D.....AD#
# ###B#.#C#.###
#   #D#B#C#.#
#   #D#B#C#.#
#   #A#B#C#.#
#   #########
#
# #############
# #AA.......AD#
# ###B#.#C#.###
#   #D#B#C#.#
#   #D#B#C#.#
#   #A#B#C#D#
#   #########
#
# #############
# #AA.......AD#
# ###.#B#C#.###
#   #D#B#C#.#
#   #D#B#C#.#
#   #A#B#C#D#
#   #########
#
# #############
# #AA.......AD#
# ###.#B#C#.###
#   #.#B#C#.#
#   #D#B#C#D#
#   #A#B#C#D#
#   #########
#
# #############
# #AA.D.....AD#
# ###.#B#C#.###
#   #.#B#C#.#
#   #.#B#C#D#
#   #A#B#C#D#
#   #########
#
# #############
# #A..D.....AD#
# ###.#B#C#.###
#   #.#B#C#.#
#   #A#B#C#D#
#   #A#B#C#D#
#   #########
#
# #############
# #...D.....AD#
# ###.#B#C#.###
#   #A#B#C#.#
#   #A#B#C#D#
#   #A#B#C#D#
#   #########
#
# #############
# #.........AD#
# ###.#B#C#.###
#   #A#B#C#D#
#   #A#B#C#D#
#   #A#B#C#D#
#   #########
#
# #############
# #..........D#
# ###A#B#C#.###
#   #A#B#C#D#
#   #A#B#C#D#
#   #A#B#C#D#
#   #########
#
# #############
# #...........#
# ###A#B#C#D###
#   #A#B#C#D#
#   #A#B#C#D#
#   #A#B#C#D#
#   #########
#
# Using the initial configuration from the full diagram, what is the least
# energy required to organize the amphipods?

# using Pkg; Pkg.add("DataStructures")
using DataStructures

# cost of moving each amphipod
const amphipodcost = [1, 10, 100, 1000]

# position of the rooms from the hallway
const roomposl = 3:2:9

# win condition
const win = [1 1 1 1; 2 2 2 2; 3 3 3 3; 4 4 4 4]

# read input
rooms = zeros(Int8, 4, 4)
rooms[1:4, 2:3] = [4 4; 3 2; 2 1; 1 3]
roomposition = 1
foreach(eachline("input.txt")) do l
  m = match(r"([A-D])#([A-D])#([A-D])#([A-D])", l)
  isnothing(m) && return

  rooms[:, roomposition] = map(c -> c[1] - 'A' + 1, m)
  global roomposition += 3
end

# Attempt to move amphipods home - mutates rooms and hallway, returns cost
function movehome!(rooms, hallway)
  cost = 0

  # keep trying to move amphipods home until we can't
  while true
    progress = false

    # for each amphipod in the hallway, see if we can move them home
    for (i,c) in enumerate(hallway)
      # check for an amphipod, and see if their room is available
      c == 0 && continue
      any(r -> r != 0 && r != c, rooms[c,:]) && continue

      # assume amphipod is coming from the left of the room
      roompos = roomposl[c]
      hallwaypos = (i+1):roompos
      if i > roompos
        # amphipod is coming from the right
        hallwaypos = roompos:(i-1)
      end

      # check if the hallway is clear to move through
      any(!=(0), hallway[hallwaypos]) && continue

      # cost is the length of the hallway, plus 1 to represent moving into the
      # room, plus each space we need to move in the room, times the amphipod's
      # movement energy.
      roompos = findlast(isequal(0), rooms[c,:])
      cost += (length(hallwaypos) + roompos) * amphipodcost[c]

      # update hallway and room
      rooms[c,roompos] = c
      hallway[i] = 0

      # makin' progress
      progress = true
    end

    !progress && break
  end

  return cost
end

# adds a state to the queue, or, if such a state already exists, update the
# cost if the new cost is less
function enqueueorupdate!(queue, state, cost)
  if haskey(queue, state)
    if cost < queue[state]
      queue[state] = cost
    end
  else
    enqueue!(queue, state, cost)
  end
end

# Dijkstra's, sorta... if you consider each "neighbor" to be all possible moves
# from the current move.
State = Tuple{Matrix{Int8}, Vector{Int8}}
queue = PriorityQueue{State, Int}((rooms, zeros(Int8, 11)) => 0)
visited = Set{State}()
while length(queue) > 0
  total_cost = peek(queue).second
  state = dequeue!(queue)
  push!(visited, state)

  local (rooms,hallway) = copy.(state)

  # are we winning yet?
  if rooms == win
    println("Result: ", total_cost)
    break
  end

  # first try moving amphipods home
  total_cost += movehome!(rooms, hallway)

  # if we've reached a win condition, push it into the queue - there's nothing
  # more to do, but there may be other ways to win that are shorter
  if rooms == win
    enqueueorupdate!(queue, (rooms, hallway), total_cost)
    continue
  end

  # now find all possible moves from the current state
  for (i,room) in enumerate(eachrow(rooms))
    # if the room is empty or already has the right amphipods, skip it
    all(c -> c == 0 || c == i, room) && continue

    # find the amphipod that's moving
    newrooms = copy(rooms)
    roomidx = findfirst(!=(0), room)
    movingamphipod = room[roomidx]
    newrooms[i,roomidx] = 0
    movement = roomidx - 1

    # find where they can move to
    hallwayexit = roomposl[i]
    for pos in [1,2,4,6,8,10,11]
      # positions we move through in the hallway
      hallwaypos = pos <= hallwayexit ? (pos:hallwayexit) : (hallwayexit:pos)

      # check if hallway is empty
      any(!=(0), hallway[hallwaypos]) && continue

      # cost is number of moves to exit the room, plus number of hallway spaces
      cost = total_cost + (length(hallwaypos) + movement) * amphipodcost[movingamphipod]

      # update hallway
      newhallway = copy(hallway)
      newhallway[pos] = movingamphipod

      # update queue
      newstate = (newrooms, newhallway)
      !in(newstate, visited) && enqueueorupdate!(queue, newstate, cost)
    end
  end
end


text:
###############
# Branch Type #
###############

type Branch <: GitLabType
    name::Nullable{GitLabString}
    protected::Nullable{Bool}
    commit::Nullable{Commit}
#=
    label::Nullable{GitLabString}
    ref::Nullable{GitLabString}
    sha::Nullable{GitLabString}
    user::Nullable{Owner}
    repo::Nullable{Repo}
    _links::Nullable{Dict}
    protection::Nullable{Dict}
=#
end

Branch(data::Dict) = json2gitlab(Branch, data)
Branch(name::AbstractString) = Branch(Dict("name" => name))

## namefield(branch::Branch) = isnull(branch.name) ? branch.ref : branch.name
namefield(branch::Branch) = branch.name

###############
# API Methods #
###############

function branches(repo; options...)
    ## MDP results, page_data = gh_get_paged_json("/repos/$(name(repo))/branches"; options...)
    results, page_data = gh_get_paged_json("/api/v3/projects/$(get(repo.id))/repository/branches"; options...)
    return map(Branch, results), page_data
end

function branch(repo, branch; options...)
    ## result = gh_get_json("/repos/$(name(repo))/branches/$(name(branch))"; options...)
    result = gh_get_json("/api/v3/projects/$(get(repo.id))/repository/branches/$(name(branch))"; options...)
    return Branch(result)
end


text:
#using IntervalTrees
using LAR
#include("src/inters.jl")
#include("test/boundary.jl")

function boxes2lar(boxes)
	V = Array{Float64,1}[]
	EV = Array{Int64,1}[]
	FV = Array{Int64,1}[]
	for k=1:size(boxes,2)
		xm,ym,xM,yM = boxes[:,k]
		push!(V,[xm,ym])
		push!(V,[xM,yM])
		push!(V,[xm,yM])
		push!(V,[xM,ym])
		push!(EV,[ 4(k-1)+1, 4(k-1)+3 ])
		push!(EV,[ 4(k-1)+1, 4(k-1)+4 ])
		push!(EV,[ 4(k-1)+2, 4(k-1)+3 ])
		push!(EV,[ 4(k-1)+2, 4(k-1)+4 ])
		push!(FV,[ 4(k-1)+1, 4(k-1)+2, 4(k-1)+3, 4(k-1)+4])
	end
	return hcat(V...), hcat(EV...), hcat(FV...)
end

lineArray = randomLines(600,.3)

larview(lineArray...)
V,EV = lines2lar(lineArray)
larview(V,EV)
chains = boundary(V,EV)
operator = boundaryOp(EV,chains)
FV = [sort(collect(Set(vcat([EV[:,abs(e)] for e in face]...)))) for face in chains]
viewLarIndices(V,EV,FV,0.1)

boxes = lar2boxes(V,FV)
W,EW,FW = boxes2lar(boxes)
larview(W,EW)

buckets = boxBucketing(boxes)

cell = 1200
Z,EZ,_ = boxes2lar(lar2boxes(V,[FV[f] for f in buckets[cell]]))
pivot = p.COLOR(p.RED)(lar2hpc(boxes2lar(boxes[:,cell])[1:2]...))
bucket = lar2hpc(Z,EZ)
p.VIEW(p.STRUCT([bucket,pivot]))


text:
using DrWatson
@quickactivate
using ArgParse
using GenerativeAD
import StatsBase: fit!, predict
using StatsBase
using BSON
# because of vae and 2stage
using DataFrames
using CSV
using ValueHistories
using Flux
using ConditionalDists
using GenerativeModels
import GenerativeModels: VAE
using Distributions
using DistributionsAD

s = ArgParseSettings()
@add_arg_table! s begin
	"max_seed"
		default = 1
		arg_type = Int
		help = "max_seed"
	"category"
		default = "wood"
		arg_type = String
		help = "category"
	"tab_name"
		default = "vae_LOSS_images_mvtec"
		arg_type = String
		help = "name of tab -> example: vae_LOSS_images, wae-vamp_AUC_images"
	"mi_only"
		arg_type = Int
		default = -1
		help = "index of model in range 1 to 10 or -1 for all models"
    "contamination"
    	arg_type = Float64
    	help = "contamination rate of training data"
    	default = 0.0
end
parsed_args = parse_args(ARGS, s)
@unpack category, max_seed, tab_name, mi_only, contamination = parsed_args

#######################################################################################
################ THIS PART IS TO BE PROVIDED FOR EACH MODEL SEPARATELY ################

sp = split(tab_name, "_")
enc = sp[1]
criterion = lowercase(sp[2])

modelname = "$(enc)_knn"

function sample_params()
	par_vec = (1:2:101,)
	argnames = (:k,)
	return (;zip(argnames, map(x->sample(x, 1)[1], par_vec))...)
end

function fit(data, parameters, aux_info)
	# construct model - constructor should only accept kwargs
	model = GenerativeAD.Models.knn_constructor(;v=:kappa, parameters...)

	# fit train data
	try
		global info, fit_t, _, _, _ = @timed fit!(model, data[1][1])
	catch e
		# return an empty array if fit fails so nothing is computed
		return (fit_t = NaN,), []
	end

	# construct return information - put e.g. the model structure here for generative models
	training_info = (
		fit_t = fit_t,
		model = nothing
		)

	# now return the different scoring functions
	function knn_predict(model, x, v::Symbol)
		try
			return predict(model, x, v)
		catch e
			if isa(e, ArgumentError) # this happens in the case when k > number of points
				return NaN # or nothing?
			else
				rethrow(e)
			end
		end
	end
	parameters = merge(parameters, aux_info)
	training_info, [(x -> knn_predict(model, x, v), merge(parameters, (distance = v,))) for v in [:gamma, :kappa, :delta]]
end

####################################################################
################ THIS PART IS COMMON FOR ALL MODELS ################
# set a maximum for parameter sampling retries
if abspath(PROGRAM_FILE) == @__FILE__

	try_counter = 0
	max_tries = 10*max_seed
	cont_string = (contamination == 0.0) ? "" : "_contamination-$contamination"
	while try_counter < max_tries
		parameters = sample_params()

		for seed in 1:max_seed
			savepath = datadir("experiments/images_mvtec$cont_string/$(modelname)/$(category)/ac=1/seed=$(seed)")
			mkpath(savepath)

			mi_indexes = (mi_only == -1) ? [1:10...] : [mi_only]

			for mi = mi_indexes

				aux_info = (model_index=mi, criterion=criterion)
				data = GenerativeAD.load_data("MVTec-AD", seed=seed, category=category,
					contamination=contamination, img_size=128)
				data, encoding_name, encoder_params = GenerativeAD.Models.load_encoding(tab_name, data, 1, dataset=category, seed=seed, model_index=mi)

				@info "Trying to fit $modelname on $category with parameters $(parameters)..."
				@info "Train/validation/test splits: $(size(data[1][1])) | $(size(data[2][1])) | $(size(data[3][1]))"
				@info "Number of features: $(size(data[1][1]))"
				# here, check if a model with the same parameters was already tested
				@info "Trying to fit $modelname on $category with parameters $(parameters)..."
				if GenerativeAD.check_params(savepath, merge(parameters, aux_info))
					training_info, results = fit(data, parameters, aux_info)
					if training_info.model !== nothing
						tagsave(joinpath(savepath, savename("model", parameters, "bson", digits=5)),
							Dict("model"=>training_info.model),
							safe = true)
						training_info = merge(training_info, (model = nothing,))
					end
					save_entries = merge(training_info, (modelname = modelname, seed = seed,
						category = category,
						contamination=contamination))
					# now loop over all anomaly score funs
					for result in results
						GenerativeAD.experiment(result..., data, savepath; save_entries...)
					end
					global try_counter = max_tries + 1
				else
					@info "Model already present, sampling new hyperparameters..."
					global try_counter += 1
				end
			end
		end
	end
	(try_counter == max_tries) ? (@info "Reached $(max_tries) tries, giving up.") : nothing
end


text:
## This function computes the Generalized Laguerre polynomials using the Tridiagonal matrix
include("orthopoly_evaluate.jl")
function compute_laguerre(n, a, x, all=true)
    # Evaluate Orthonormalized Laguerre polynomials at x
    #  Laguerre:
    #     w(x)=x^a exp(-x) with integral gamma(a+1)  (a>-1)
    #     Equivalent Mathematica:LaguerreL[n,a,x]/Sqrt[Gamma[n+a+1]/Gamma[n+1]]
    ## Laguerre Bidiagonal matrix
    B = diagm(sqrt((a + 1 +(0:n)))) # diagonal
    B = B - diagm(sqrt((1:n)), 1)   # superdiagonal
    T = B'*B
    c = gamma(a+1) #Integral of weight function
    ## Default to 0 through n computation
    if all
        phi = orthopoly_evaluate_all(T, x)
    else
        phi = orthopoly_evaluate_n(T, x)
    end
    phi/sqrt(c)
end


text:
module TestBasic
using ZulipReminderBot
using ZulipReminderBot: ZulipRequest, process, OPTS, Message
using ZulipReminderBot: isfixedtz
using JSON3
using Test
using Dates
using Logging

global_logger(ConsoleLogger(stderr, Logging.Warn))

f = readlines(joinpath(@__DIR__, "data", "requests.json"))
const MENTION_MSG = f[1]
# const PRIVATE_MSG = f[1]
# const WRONG_TOKEN = f[3]
# const NO_TOKEN = f[4]
# const NO_SENDER = f[5]
# const NO_DATA = f[6]
# const PRIVATE_HELP = f[7]
# const MENTION_HELP = f[8]

content(x) = JSON3.read(x).content

const TOKEN = "abc123"
setupbot!(; token = TOKEN, email = "foo@bar", apikey = "123")

@testset "test time" begin
    # msg = JSON3.read(MENTION_MSG, ZulipRequest)
    # c = Channel(1)
    # reply = process(msg, c)
    # @test startswith(content(reply), "Message is scheduled on 20")
end

@testset "test timezone parsing" begin
    @test isfixedtz("UTC+1")
    @test isfixedtz("UTC+12")
    @test isfixedtz("UTC-1")
    @test isfixedtz("UTC-12")
    @test !isfixedtz("UTC+123")
    @test !isfixedtz("GMT-1")
end

end # module


text:
function MultivariateNormalDistribution(arg0::RandomGenerator, arg1::Vector{jdouble}, arg2::Vector{Vector{jdouble}})
    return MultivariateNormalDistribution((RandomGenerator, Vector{jdouble}, Vector{Vector{jdouble}}), arg0, arg1, arg2)
end

function MultivariateNormalDistribution(arg0::Vector{jdouble}, arg1::Vector{Vector{jdouble}})
    return MultivariateNormalDistribution((Vector{jdouble}, Vector{Vector{jdouble}}), arg0, arg1)
end

function density(obj::MultivariateNormalDistribution, arg0::Vector{jdouble})
    return jcall(obj, "density", jdouble, (Vector{jdouble},), arg0)
end

function get_covariances(obj::MultivariateNormalDistribution)
    return jcall(obj, "getCovariances", RealMatrix, ())
end

function get_means(obj::MultivariateNormalDistribution)
    return jcall(obj, "getMeans", Vector{jdouble}, ())
end

function get_standard_deviations(obj::MultivariateNormalDistribution)
    return jcall(obj, "getStandardDeviations", Vector{jdouble}, ())
end

function sample(obj::MultivariateNormalDistribution)
    return jcall(obj, "sample", Vector{jdouble}, ())
end


text:
# Nightmare inducing logic to combine/remove nodes, e.g layers and activation functions or reshape between recurrent layers and dense layers.

struct ActivationFunction end

struct ActivationLayer end

struct RecurrentLayer end

struct ElemwiseOp end

"""
   check_combine(gb::CompGraphBuilder, n::OnnxNode)

Return a potentially combined [`OnnxNode`](@ref) and its inputs.

Purpose is to handle e.g activation functions to layers which Flux has as member of the layer struct while ONNX has as a separate node.

Main motivation is an attempt to make serialization/deserialization "pure" in the sense that `g -> serialize -> deserialize === g`.

Another motivation is that mutation becomes a bit harder if things like global pooling and dropping of dimensions are separated into two vertices.
"""
check_combine(gb::CompGraphBuilder, n::OnnxNode) = check_combine(optrait(n), n, gb)

optrait(n::OnnxNode) = optrait(Val(optype(n)), n)
function optrait(vot::Val{ot}, n::OnnxNode) where ot
      ot in keys(actfuns) && return ActivationFunction()
      ot in keys(fluxrecurrentlayers) && return RecurrentLayer()
      ot in keys(fluxlayers) && return ActivationLayer()
      return vot
end
optrait(::Val{:Mul}, ::OnnxNode) = ElemwiseOp()
optrait(::Val{:Add}, ::OnnxNode) = ElemwiseOp()

retnode(n, gb) = n, innodes(n, gb)
check_combine(x, n::OnnxNode, gb::CompGraphBuilder) = retnode(n, gb)

# Activation function: If there is exactly one input and that input is a layer with an activation function then combine them into one node
check_combine(::ActivationFunction, n::OnnxNode, gb::CompGraphBuilder) = check_actfun(n, gb, innodes(n, gb)...)

check_actfun(nact::OnnxNode, gb::CompGraphBuilder, args...) = retnode(nact, gb)
check_actfun(nact::OnnxNode, gb::CompGraphBuilder, nlayer::OnnxNode) = check_actfun(nact, gb, nlayer, optrait(nlayer))
check_actfun(nact::OnnxNode, gb::CompGraphBuilder, nlayer::OnnxNode, ::ActivationLayer) = wrapfrom(nact, nlayer, gb, :activation, actfuns)

function wrapfrom(nwrap::OnnxNode, nwrapped::OnnxNode, gb::CompGraphBuilder, attrib, from::AbstractDict=invariantops)
   @debug "Merge node $nwrap and node $nwrapped"
   nwrapped.attribute[attrib] = from[optype(nwrap)](nwrap.attribute, params(nwrap)...)
   return retnode(nwrapped, gb)
end

# Reshape
# 1. Check if input node is a global pool and if so combine them into one op
# 2. Check if input is from recurrent layers and if so remove the reshape as it is not needed
check_combine(::Val{:Reshape}, n::OnnxNode, gb::CompGraphBuilder) = check_reshape(n, gb, innodes(n, gb)...)

check_reshape(nreshape::OnnxNode, gb::CompGraphBuilder, args...) = retnode(nreshape, gb)
check_reshape(nreshape::OnnxNode, gb::CompGraphBuilder, innode::OnnxNode) = check_reshape(nreshape, gb, innode, optrait(innode))

check_reshape(nreshape::OnnxNode, gb::CompGraphBuilder, innode::OnnxNode, ::Val{:GlobalAveragePool}) = wrapfrom(nreshape, innode, gb, :wrap, pseudotransparentops)
check_reshape(nreshape::OnnxNode, gb::CompGraphBuilder, innode::OnnxNode, ::Val{:GlobalMaxPool}) = wrapfrom(nreshape, innode, gb, :wrap, pseudotransparentops)

check_reshape(nreshape::OnnxNode, gb::CompGraphBuilder, nrec::OnnxNode, ::RecurrentLayer) = check_recurrent_reshape(nreshape, gb, nrec)
check_reshape(nreshape::OnnxNode, gb::CompGraphBuilder, nsqueeze::OnnxNode, ::Val{:Squeeze}) = check_recurrent_reshape(nreshape, gb, innodes(nsqueeze, gb)...)

check_recurrent_reshape(nreshape::OnnxNode, gb::CompGraphBuilder, args...) = retnode(nreshape, gb)
check_recurrent_reshape(nreshape::OnnxNode, gb::CompGraphBuilder, innode::OnnxNode) = check_recurrent_reshape(nreshape, gb, innode, optrait(innode))

function check_recurrent_reshape(nreshape::OnnxNode, gb::CompGraphBuilder, nrec::OnnxNode, ::RecurrentLayer)
   outs = outnodes(nreshape, gb)
   # Maybe check that the reshape really is doing 3D -> 2D, but what else could it do in this case?
   if length(outs) > 0 && all(check_recurrent_reshape, Val.(optype.(outs)))
      @debug "Remove reshape found between $(outs) and $nrec"
      return check_combine(gb, innodes(nreshape, gb)[1])
   end
end

check_recurrent_reshape(::Val) = false
check_recurrent_reshape(::Val{:Gemm}) = true

# Squeeze
# 1. Check if input node is a global pool and if so combine them into one op
# 2.  Check if input is from recurrent layers and if so remove the squeeze as it is not needed
check_combine(::Val{:Squeeze}, n::OnnxNode, gb::CompGraphBuilder) = check_squeeze(n, gb, innodes(n, gb)...)

check_squeeze(nsqueeze::OnnxNode, gb::CompGraphBuilder, args...) = retnode(nsqueeze, gb)
check_squeeze(nsqueeze::OnnxNode, gb::CompGraphBuilder, innode::OnnxNode) = check_squeeze(nsqueeze, gb, innode, optrait(innode))

check_squeeze(nsqueeze::OnnxNode, gb::CompGraphBuilder, innode::OnnxNode, ::Val{:GlobalAveragePool}) = wrapfrom(nsqueeze, innode, gb, :wrap)
check_squeeze(nsqueeze::OnnxNode, gb::CompGraphBuilder, innode::OnnxNode, ::Val{:GlobalMaxPool}) = wrapfrom(nsqueeze, innode, gb, :wrap)

function check_squeeze(nsqueeze::OnnxNode, gb::CompGraphBuilder, innode::OnnxNode, ::RecurrentLayer)
   @debug "Remove squeeze after $innode"
   return retnode(innode, gb)
end

# For element wise operations, check if one of the inputs is a constant
check_combine(::ElemwiseOp, n::OnnxNode, gb::CompGraphBuilder) = check_elemwise(n, gb, innodes(n, gb))

function check_elemwise(nelemwise::OnnxNode, gb::CompGraphBuilder, innodes::AbstractVector{OnnxNode})
   ins_to_process = filter(!isnothing, map(innode -> check_elemwise(nelemwise, gb, innode, optrait(innode)), innodes))
   return nelemwise, ins_to_process
end

check_elemwise(::OnnxNode, ::CompGraphBuilder, innode::OnnxNode, ot) = innode

function check_elemwise(nelemwise::OnnxNode, gb::CompGraphBuilder, innode::OnnxNode, ::Val{:Constant})
   consts = get!(nelemwise.attribute, :Constant, [])
   push!(consts, sources[optype(innode)](innode.attribute, params(innode)...))
   return nothing
end


text:
using JuMP
using GLPK

model = Model(GLPK.Optimizer)
@variable(model, x >= 0)
@variable(model, y >= 0)
@constraint(model, 6x + 8y >= 100)
@constraint(model, 7x + 12y >= 120)
@objective(model, Min, 12x + 20y)

optimize!(model)

@show value(x);
@show value(y);
@show objective_value(model);

using JuMP

using GLPK

model = Model(GLPK.Optimizer);

@variable(model, x >= 0)
@variable(model, y >= 0);

@constraint(model, 6x + 8y >= 100)
@constraint(model, 7x + 12y >= 120);

@objective(model, Min, 12x + 20y);

optimize!(model)

@show value(x);
@show value(y);
@show objective_value(model);