text:
using Test

using MathOptInterface
const MOI = MathOptInterface
const MOIT = MathOptInterface.Test
const MOIU = MathOptInterface.Utilities
const MOIB = MathOptInterface.Bridges

include("utilities.jl")

# Model not supporting RotatedSecondOrderCone
MOIU.@model(NoRSOCModel,
            (),
            (MOI.EqualTo, MOI.GreaterThan, MOI.LessThan),
            (MOI.Zeros, MOI.Nonnegatives, MOI.Nonpositives, MOI.SecondOrderCone,
             MOI.ExponentialCone, MOI.PositiveSemidefiniteConeTriangle),
            (),
            (),
            (MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction),
            (MOI.VectorOfVariables,),
            (MOI.VectorAffineFunction, MOI.VectorQuadraticFunction))

# Model not supporting VectorOfVariables and SingleVariable
MOIU.@model(NoVariableModel,
            (MOI.ZeroOne, MOI.Integer),
            (MOI.EqualTo, MOI.GreaterThan, MOI.LessThan),
            (MOI.Zeros, MOI.Nonnegatives, MOI.Nonpositives, MOI.SecondOrderCone),
            (),
            (),
            (MOI.ScalarAffineFunction,),
            (),
            (MOI.VectorAffineFunction,))
function MOI.supports_constraint(::NoVariableModel, ::Type{MOI.SingleVariable},
                                 ::Type{<:MOI.AbstractScalarSet})
    return false
end

# Only supports GreaterThan and Nonnegatives
MOIU.@model(GreaterNonnegModel,
            (),
            (MOI.GreaterThan,),
            (MOI.Nonnegatives,),
            (),
            (),
            (MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction),
            (MOI.VectorOfVariables,),
            (MOI.VectorAffineFunction, MOI.VectorQuadraticFunction))
function MOI.supports_constraint(
    ::GreaterNonnegModel{T}, ::Type{MOI.SingleVariable},
    ::Type{<:Union{MOI.EqualTo{T}, MOI.LessThan{T}, MOI.Interval{T}}}) where T
    return false
end

MOIU.@model(ModelNoVAFinSOC,
            (),
            (MOI.EqualTo, MOI.GreaterThan, MOI.LessThan, MOI.Interval),
            (MOI.Zeros, MOI.Nonnegatives, MOI.Nonpositives, MOI.SecondOrderCone,
             MOI.RotatedSecondOrderCone, MOI.GeometricMeanCone,
             MOI.PositiveSemidefiniteConeTriangle, MOI.ExponentialCone),
            (MOI.PowerCone, MOI.DualPowerCone),
            (),
            (MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction),
            (MOI.VectorOfVariables,),
            (MOI.VectorAffineFunction, MOI.VectorQuadraticFunction))

MOI.supports_constraint(::ModelNoVAFinSOC{Float64},
                        ::Type{MOI.VectorAffineFunction{Float64}},
                        ::Type{MOI.SecondOrderCone}) = false

# Model supporting nothing
MOIU.@model NothingModel () () () () () () () ()
function MOI.supports_constraint(
    ::NothingModel{T}, ::Type{MOI.SingleVariable},
    ::Type{<:Union{MOI.EqualTo{T}, MOI.GreaterThan{T}, MOI.LessThan{T},
                   MOI.Interval{T}, MOI.Integer, MOI.ZeroOne}}) where T
    return false
end

struct BridgeAddingNoConstraint{T} <: MOI.Bridges.AbstractBridge end
MOIB.added_constraint_types(::Type{BridgeAddingNoConstraint{T}}) where {T} = []
function MOI.supports_constraint(::Type{<:BridgeAddingNoConstraint},
                                 ::Type{MOI.SingleVariable},
                                 ::Type{MOI.Integer})
    return true
end
function MOIB.concrete_bridge_type(::Type{<:BridgeAddingNoConstraint{T}},
                                   ::Type{MOI.SingleVariable},
                                   ::Type{MOI.Integer}) where {T}
    return BridgeAddingNoConstraint{T}
end

@testset "Bridge adding no constraint" begin
    mock = MOIU.MockOptimizer(NothingModel{Int}())
    bridged = MOIB.LazyBridgeOptimizer(mock)
    MOI.Bridges.add_bridge(bridged, BridgeAddingNoConstraint{Float64})
    @test MOI.Bridges.supports_bridging_constraint(bridged,
                                                   MOI.SingleVariable,
                                                   MOI.Integer)
end

@testset "Unsupported constraint with cycles" begin
    # Test that `supports_constraint` works correctly when it is not
    # supported but the bridges forms a cycle
    mock = MOIU.MockOptimizer(NothingModel{Float64}())
    bridged = MOIB.full_bridge_optimizer(mock, Float64)
    @test !MOI.supports_constraint(
        bridged, MOI.SingleVariable, MOI.GreaterThan{Float64})
    @test !MOI.supports_constraint(
        bridged, MOI.VectorAffineFunction{Float64}, MOI.Nonpositives)
end

mock = MOIU.MockOptimizer(NoRSOCModel{Float64}())
bridged_mock = MOIB.LazyBridgeOptimizer(mock)

@testset "UnsupportedConstraint when it cannot be bridged" begin
    x = MOI.add_variables(bridged_mock, 4)
    err = MOI.UnsupportedConstraint{MOI.VectorOfVariables,
                                    MOI.RotatedSecondOrderCone}()
    @test_throws err begin
        MOI.add_constraint(bridged_mock, MOI.VectorOfVariables(x),
                           MOI.RotatedSecondOrderCone(4))
    end
end

MOIB.add_bridge(bridged_mock, MOIB.SplitIntervalBridge{Float64})
MOIB.add_bridge(bridged_mock, MOIB.RSOCtoPSDBridge{Float64})
MOIB.add_bridge(bridged_mock, MOIB.SOCtoPSDBridge{Float64})
MOIB.add_bridge(bridged_mock, MOIB.RSOCBridge{Float64})

@testset "Name test" begin
    MOIT.nametest(bridged_mock)
end

@testset "Copy test" begin
    MOIT.failcopytestc(bridged_mock)
    MOIT.failcopytestia(bridged_mock)
    MOIT.failcopytestva(bridged_mock)
    MOIT.failcopytestca(bridged_mock)
    MOIT.copytest(bridged_mock, NoRSOCModel{Float64}())
end

# Test that RSOCtoPSD is used instead of RSOC+SOCtoPSD as it is a shortest path.
@testset "Bridge selection" begin
    MOI.empty!(bridged_mock)
    @test !(MOI.supports_constraint(bridged_mock,
                                    MOI.VectorAffineFunction{Float64},
                                    MOI.LogDetConeTriangle))
    x = MOI.add_variables(bridged_mock, 3)
    err = MOI.UnsupportedConstraint{MOI.VectorAffineFunction{Float64},
                                    MOI.LogDetConeTriangle}()
    @test_throws err begin
        MOIB.bridge_type(bridged_mock, MOI.VectorAffineFunction{Float64},
                         MOI.LogDetConeTriangle)
    end
    c = MOI.add_constraint(bridged_mock, MOI.VectorOfVariables(x),
                           MOI.RotatedSecondOrderCone(3))
    @test MOIB.bridge_type(bridged_mock, MOI.VectorOfVariables,
                MOI.RotatedSecondOrderCone) == MOIB.RSOCtoPSDBridge{Float64}
    @test MOIB.bridge(bridged_mock, c) isa MOIB.RSOCtoPSDBridge
    @test bridged_mock.dist[(MOI.VectorOfVariables,
                            MOI.RotatedSecondOrderCone)] == 1
end

@testset "Supports" begin
    full_bridged_mock = MOIB.full_bridge_optimizer(mock, Float64)
    greater_nonneg_mock = MOIU.MockOptimizer(GreaterNonnegModel{Float64}())
    full_bridged_greater_nonneg = MOIB.full_bridge_optimizer(
        greater_nonneg_mock, Float64)
    for F in [MOI.SingleVariable, MOI.ScalarAffineFunction{Float64},
              MOI.ScalarQuadraticFunction{Float64}]
        @test MOI.supports_constraint(full_bridged_mock, F,
                                      MOI.Interval{Float64})
        @test !MOI.supports_constraint(
            greater_nonneg_mock, F, MOI.LessThan{Float64})
        @test MOI.supports_constraint(
            full_bridged_greater_nonneg, F, MOI.LessThan{Float64})
    end
    for F in [MOI.VectorOfVariables, MOI.VectorAffineFunction{Float64},
              MOI.VectorQuadraticFunction{Float64}]
        @test MOI.supports_constraint(full_bridged_mock, F,
                                      MOI.PositiveSemidefiniteConeSquare)
        @test MOI.supports_constraint(full_bridged_mock, F,
                                      MOI.GeometricMeanCone)
        @test !MOI.supports_constraint(
            greater_nonneg_mock, F, MOI.Nonpositives)
        @test MOI.supports_constraint(
            full_bridged_greater_nonneg, F, MOI.Nonnegatives)
    end
    for F in [MOI.VectorOfVariables, MOI.VectorAffineFunction{Float64}]
        # The bridges in this for loop do not support yet
        # VectorQuadraticFunction. See TODO's for the reason.
        # TODO: Missing vcat for quadratic for supporting quadratic.
        @test MOI.supports_constraint(full_bridged_mock, F,
                                      MOI.RotatedSecondOrderCone)
        # TODO: Det bridges need to use MOIU.operate to support quadratic.
        @test MOI.supports_constraint(full_bridged_mock, F,
                                      MOI.LogDetConeTriangle)
        @test MOI.supports_constraint(full_bridged_mock, F,
                                      MOI.RootDetConeTriangle)
    end
    mock2 = MOIU.MockOptimizer(ModelNoVAFinSOC{Float64}())
    @test !MOI.supports_constraint(mock2, MOI.VectorAffineFunction{Float64},
                                   MOI.SecondOrderCone)
    full_bridged_mock2 = MOIB.full_bridge_optimizer(mock2, Float64)
    @test MOI.supports_constraint(full_bridged_mock2, MOI.VectorAffineFunction{Float64},
                                  MOI.SecondOrderCone)
    @testset "Unslack" begin
        for T in [Float64, Int]
            no_variable_mock = MOIU.MockOptimizer(NoVariableModel{T}())
            full_bridged_no_variable = MOIB.full_bridge_optimizer(
                no_variable_mock, T)
            for S in [MOI.LessThan{T}, MOI.GreaterThan{T}, MOI.EqualTo{T},
                      MOI.ZeroOne, MOI.Integer]
                @test MOI.supports_constraint(
                    full_bridged_no_variable, MOI.SingleVariable, S)
            end
            for S in [MOI.Nonpositives, MOI.Nonnegatives,
                      MOI.Zeros, MOI.SecondOrderCone]
                @test MOI.supports_constraint(
                    full_bridged_no_variable, MOI.VectorOfVariables, S)
            end
        end
    end
end

@testset "Combining two briges" begin
    full_bridged_mock = MOIB.full_bridge_optimizer(mock, Float64)
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [1, 1, 0, 1, 1, 0, 1, √2])
    config = MOIT.TestConfig()
    MOIT.rootdett1vtest(full_bridged_mock, config)
    MOIT.rootdett1ftest(full_bridged_mock, config)
    # Dual is not yet implemented for RootDet and GeoMean bridges
    ci = first(MOI.get(full_bridged_mock, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{Float64}, MOI.RootDetConeTriangle}()))
    test_delete_bridge(full_bridged_mock, ci, 4, ((MOI.VectorAffineFunction{Float64}, MOI.RotatedSecondOrderCone, 0),
                                                (MOI.VectorAffineFunction{Float64}, MOI.GeometricMeanCone, 0),
                                                (MOI.VectorAffineFunction{Float64}, MOI.PositiveSemidefiniteConeTriangle, 0)),
                       used_bridges = 3)
end

@testset "Continuous Linear" begin
    exclude = ["partial_start"] # VariablePrimalStart not supported.
    MOIT.contlineartest(bridged_mock, MOIT.TestConfig(solve=false), exclude)
end

@testset "Continuous Conic" begin
    MOIT.contconictest(MOIB.full_bridge_optimizer(mock, Float64), MOIT.TestConfig(solve=false), ["logdets", "rootdets", "psds"])
end


text:
let
    filename = joinpath(dirname(pathof(LDrawParser)),"..","assets","ATTEWalker.mpd")
    # filename = joinpath(dirname(pathof(LDrawParser)),"..","assets","Millennium Falcon.mpd")
    model = parse_ldraw_file(filename)
    for (k,p) in model.parts
        @test p.populated.status == false
    end
    @test isempty(model.sub_parts)
    @test !isempty(model.models["20009 - AT-TE Walker.mpd"].steps[1].lines)
    # load geometry
    if isdir(get_part_library_dir())
        LDrawParser.populate_part_geometry!(model)
        for (k,p) in Base.Iterators.flatten((model.parts,model.sub_parts))
            @assert p.populated.status == true
        end
    else
        @warn "Skipping tests for populate_part_geometry!(...) because $(get_part_library_dir()) does not exist"
    end

end

ext: jl


text:
#=
Predicate functions to filter points in Clouds based on positions or attributes.
A point is defined by a PointCould and index.

onsistent use of pointfilters:
value in scripts can either be <nothing> or a function according to format below
--> true: the points are used in the function.
--> False: the points are skipped

Format should be according to:

    function pointfilter(cloud::Cloud, index::Int)
        true
    end

=#

ground_water(cloud::Cloud, i::Integer) = (cloud[:classification][i] == 2) || (cloud[:classification][i] == 9)
ground(cloud::Cloud, i::Integer) = (cloud[:classification][i] == 2)
notground(cloud::Cloud, i::Integer) = (cloud[:classification][i] != 2)
water(cloud::Cloud, i::Integer) = (cloud[:classification][i] == 9)
notwater(cloud::Cloud, i::Integer) = (cloud[:classification][i] != 9)
notoutlier(cloud::Cloud, i::Integer) = cloud[:classification][i] != 7 && cloud[:classification][i] != 18
outlier(cloud::Cloud, i::Integer) = cloud[:classification][i] == 7 || cloud[:classification][i] == 18
unclassified(cloud::Cloud, i::Integer) = cloud[:classification][i] == 0
unclassified_lastreturn(cloud::Cloud, i::Integer) = lastreturn(cloud, i) && unclassified(cloud, i)
unclassified_firstreturn(cloud::Cloud, i::Integer) = firstreturn(cloud, i) && unclassified(cloud, i)

highscanangle(cloud::Cloud, i::Integer) = abs(cloud[:scan_angle][i]) > 20.0
lowscanangle(cloud::Cloud, i::Integer) = abs(cloud[:scan_angle][i]) <= 20.0

# dummy filters
filtertrue(::Cloud, ::Integer) = true
filterfalse(::Cloud, ::Integer) = false

# ASPRS, but only on normalized pointcloud(!)
low_veg(cloud::Cloud, i::Integer) = 0.5 < getz(cloud, i) <= 2.0
med_veg(cloud::Cloud, i::Integer) = 2.0 < getz(cloud, i) <= 5.0
high_veg(cloud::Cloud, i::Integer) =  5.0 < getz(cloud, i)

non_tree(cloud::Cloud, i::Integer) = getz(cloud, i) < 50.0  # higher trees don't exist ;)

# check if last return
lastreturn(cloud::Cloud, i::Integer) = cloud[:return_number][i] == cloud[:number_of_returns][i]
firstreturn(cloud::Cloud, i::Integer) = cloud[:return_number][i] == 1
singlereturn(cloud::Cloud, i::Integer) = cloud[:return_number][i] == cloud[:number_of_returns][i] == 1

# edge of flightline
flightedge(cloud::Cloud, i::Integer) = cloud[:edge_of_flight_line][i] == 1

# combined

"First returns, not ground and not on edge. Default for chm method"
function unclassifiedfirstreturnonedge(cloud::Cloud, i::Integer)
    unclassified_firstreturn(cloud, i) && !flightedge(cloud, i)
end

"First returns, not ground and not with high scan angle."
function unclassifiedfirstreturnlowscanangle(cloud::Cloud, i::Integer)
    unclassified_firstreturn(cloud, i) && lowscanangle(cloud, i)
end

"First returns, not ground, only trees and not with high scan angle."
function unclassifiedfirstreturnnontreelowscanangle(cloud::Cloud, i::Integer)
    unclassified_firstreturn(cloud, i) && lowscanangle(cloud, i) && non_tree(cloud, i)
end

"not outlier and is last return. default pointfilter for classify_ground! method"
function notoutlier_lastreturn(cloud::Cloud, i::Integer)
    lastreturn(cloud, i) && notoutlier(cloud, i)
end

"not outlier and is last return. default pointfilter for classify_ground! method"
function notoutlier_ground(cloud::Cloud, i::Integer)
    notoutlier(cloud, i) && ground(cloud, i)
end

"not outlier and is last return. default pointfilter for classify_ground! method"
function notoutlier_notwater_lastreturn(cloud::Cloud, i::Integer)
    lastreturn(cloud, i) && notoutlier(cloud, i) && notwater(cloud, i)
end


text:
include("BinDeps.jl")

using .BinDeps

@static if Base.VERSION < v"1.3-" || (Sys.iswindows() && Sys.WORD_SIZE == 32)
    @BinDeps.setup()

    using Libdl

    function validate_libz_version(name, handle)
        f = Libdl.dlsym_e(handle, "zlibVersion")
        f == C_NULL && return false
        ver = VersionNumber(unsafe_string(ccall(f, Cstring, ())))
        println("Version Libz is: $ver")
        # Version 1.2.8 or above
        return ver >= v"1.2.8"
    end

    function validate_openssl_version(name, handle)
        f = Libdl.dlsym_e(handle, "OpenSSL_version_num")
        f == C_NULL && return false
        v = ccall(f, Culong, ())
        println("Version OpenSSL is: $(string(v, base=16))")
        # Version 1.1.0 or above
        return v >= 0x1010000f
    end

    libz = library_dependency("libz", aliases=["libz", "libzlib", "zlib1"], validate=validate_libz_version)
    libcrypto = library_dependency("libcrypto", aliases=["libcrypto", "libcrypto-1_1-x64", "libcrypto-1_1"], validate=validate_openssl_version)

    prefix = joinpath(@__DIR__, "usr")

    if !Sys.iswindows()
        provides(Sources,
                 URI("https://github.com/madler/zlib/archive/v1.2.11.tar.gz"),
                 libz, unpacked_dir="zlib-1.2.11")
        provides(SimpleBuild,
                 (@build_steps begin
                  GetSources(libz)
                  @build_steps begin
                  ChangeDirectory(joinpath(BinDeps.depsdir(libz), "src", "zlib-1.2.11"))
                     `./configure --prefix=$prefix`
                     `make`
                     `make install`
                 end
              end), libz, os = :Unix)

        osslver = "1_1_0k"
        osslfn  = "OpenSSL_$(osslver)"
        ossldir = "openssl-$(osslfn)"
        osslpkg = "$(osslfn).tar.gz"

        provides(Sources,
                 URI("https://github.com/openssl/openssl/archive/$(osslpkg)"),
                 libcrypto, unpacked_dir="$(ossldir)")
        provides(SimpleBuild,
                 (@build_steps begin
                  GetSources(libcrypto)
                  @build_steps begin
                     ChangeDirectory(joinpath(BinDeps.depsdir(libcrypto), "src", "$(ossldir)"))
                     `./config --prefix=$prefix`
                     `make depend`
                     `make install`
                 end
              end), libcrypto, os = :Unix)

    else
        zlib_bn  = "zlib-1.2.11-win$(Sys.WORD_SIZE)-mingw"
        zlib_fn  = "$(zlib_bn).zip"
        zlib_uri = "https://bintray.com/vszakats/generic/download_file?file_path=$(zlib_fn)"
        provides(Binaries, URI(zlib_uri), libz, filename="$(zlib_fn)", unpacked_dir="$(zlib_bn)")

        openssl_bn   = "openssl-1.1.0i-win$(Sys.WORD_SIZE)-mingw"
        openssl_fn   = "$(openssl_bn).zip"
        openssl_uri  = "https://bintray.com/vszakats/generic/download_file?file_path=$(openssl_fn)"
        provides(Binaries, URI(openssl_uri), libcrypto, filename="$(openssl_fn)", unpacked_dir="$(openssl_bn)")
    end

    @BinDeps.install Dict([:libz => :libz, :libcrypto => :libcrypto])
end


text:
@testset "quadratic problems" begin
    @testset "qp1_test" begin
        #=
            min x^2 + xy + y^2 + yz + z^2
        s.t.
            x + 2y + 3z >= 4 (a)
            x +  y      >= 1 (b)
            x,y \in R

        dual
            max 4a +  b - w1^2 - w1 w2 - w2^2 - w2 w3 - w3^2
        s.t.
            - a -  b + 2w1 + 1w2        = 0
            -2a -  b + 1w1 + 2w2 + 1w3  = 0
            -3a            + 1w2 + 2w3  = 0
            a >= 0
            b >= 0
        =#
        primal_model = qp1_test()
        dual_model, primal_dual_map = dual_model_and_map(primal_model)

        @test MOI.get(dual_model, MOI.NumberOfVariables()) == 2 + 3
        list_of_cons =  MOI.get(dual_model, MOI.ListOfConstraints())
        @test Set(list_of_cons) == Set([
            (SVF, MOI.GreaterThan{Float64})
            (SAF{Float64}, MOI.EqualTo{Float64})
        ])
        @test MOI.get(dual_model, MOI.NumberOfConstraints{SVF, MOI.GreaterThan{Float64}}()) == 2
        @test MOI.get(dual_model, MOI.NumberOfConstraints{SAF{Float64}, MOI.EqualTo{Float64}}()) == 3
        obj_type = MOI.get(dual_model, MOI.ObjectiveFunctionType())
        @test obj_type == SQF{Float64}
        obj = MOI.get(dual_model, MOI.ObjectiveFunction{obj_type}())
        @test MOI.constant(obj) == 0.0
        @test MOI.coefficient.(obj.affine_terms) == [1.0; 4.0]
        @test MOI.coefficient.(obj.quadratic_terms) == [-2.0; -1.0; -2.0; -1.0; -2.0]

        eq_con1_fun = MOI.get(dual_model, MOI.ConstraintFunction(), CI{SAF{Float64}, MOI.EqualTo{Float64}}(1))
        eq_con1_set = MOI.get(dual_model, MOI.ConstraintSet(), CI{SAF{Float64}, MOI.EqualTo{Float64}}(1))
        @test MOI.coefficient.(eq_con1_fun.terms) == [1.0; 1.0; -2.0; -1.0]
        @test MOI.constant.(eq_con1_fun) == 0.0
        @test MOI.constant(eq_con1_set) == 0.0
        eq_con2_fun = MOI.get(dual_model, MOI.ConstraintFunction(), CI{SAF{Float64}, MOI.EqualTo{Float64}}(2))
        eq_con2_set = MOI.get(dual_model, MOI.ConstraintSet(), CI{SAF{Float64}, MOI.EqualTo{Float64}}(2))
        @test MOI.coefficient.(eq_con2_fun.terms) == [2.0; 1.0; -1.0; -2.0; -1.0]
        @test MOI.constant.(eq_con2_fun) == 0.0
        @test MOI.constant(eq_con2_set) == 0.0
        eq_con3_fun = MOI.get(dual_model, MOI.ConstraintFunction(), CI{SAF{Float64}, MOI.EqualTo{Float64}}(3))
        eq_con3_set = MOI.get(dual_model, MOI.ConstraintSet(), CI{SAF{Float64}, MOI.EqualTo{Float64}}(3))
        @test MOI.coefficient.(eq_con3_fun.terms) == [3.0; -1.0; -2.0]
        @test MOI.constant.(eq_con3_fun) == 0.0
        @test MOI.constant(eq_con3_set) == 0.0

        primal_con_dual_var = primal_dual_map.primal_con_dual_var
        @test primal_con_dual_var[CI{SAF{Float64}, MOI.GreaterThan{Float64}}(1)] == [VI(1)]
        @test primal_con_dual_var[CI{SAF{Float64}, MOI.GreaterThan{Float64}}(2)] == [VI(2)]

        primal_var_dual_con = primal_dual_map.primal_var_dual_con
        @test primal_var_dual_con[VI(1)] == CI{SAF{Float64}, MOI.EqualTo{Float64}}(1)
        @test primal_var_dual_con[VI(2)] == CI{SAF{Float64}, MOI.EqualTo{Float64}}(2)
        @test primal_var_dual_con[VI(3)] == CI{SAF{Float64}, MOI.EqualTo{Float64}}(3)

        primal_var_dual_quad_slack = primal_dual_map.primal_var_dual_quad_slack
        @test primal_var_dual_quad_slack[VI(1)] == VI(2+1)
        @test primal_var_dual_quad_slack[VI(2)] == VI(2+2)
        @test primal_var_dual_quad_slack[VI(3)] == VI(2+3)
    end
    @testset "qp2_test" begin
        #=
            min 2 x^2 + y^2 + xy + x + y + 1
        s.t.
            x + y = 1 (a)
            x    >= 0 (b)
               y >= 0 (c)

        dual
            max a + 1 - 2 w1^2 - w2^2 - w1 w2
        s.t.
            - a - b     + 4w1 + 1w2 + 1 = 0
            - a     - c + 1w1 + 2w2 + 1 = 0
            a \in R
            b >= 0
            c >= 0
        =#
        primal_model = qp2_test()
        dual_model, primal_dual_map = dual_model_and_map(primal_model)

        @test MOI.get(dual_model, MOI.NumberOfVariables()) == 3 + 2
        list_of_cons =  MOI.get(dual_model, MOI.ListOfConstraints())
        @test Set(list_of_cons) == Set([
            (SVF, MOI.GreaterThan{Float64})
            (SAF{Float64}, MOI.EqualTo{Float64})
        ])
        @test MOI.get(dual_model, MOI.NumberOfConstraints{SVF, MOI.GreaterThan{Float64}}()) == 2
        @test MOI.get(dual_model, MOI.NumberOfConstraints{SAF{Float64}, MOI.EqualTo{Float64}}()) == 2
        obj_type = MOI.get(dual_model, MOI.ObjectiveFunctionType())
        @test obj_type == SQF{Float64}
        obj = MOI.get(dual_model, MOI.ObjectiveFunction{obj_type}())
        @test MOI.constant(obj) == 1.0
        @test MOI.coefficient.(obj.affine_terms) == [1.0]
        @test MOI.coefficient.(obj.quadratic_terms) == [-4.0; -1.0; -2.0]

        eq_con1_fun = MOI.get(dual_model, MOI.ConstraintFunction(), CI{SAF{Float64}, MOI.EqualTo{Float64}}(1))
        eq_con1_set = MOI.get(dual_model, MOI.ConstraintSet(), CI{SAF{Float64}, MOI.EqualTo{Float64}}(1))
        @test MOI.coefficient.(eq_con1_fun.terms) == [1.0; 1.0; -4.0; -1.0]
        @test MOI.constant.(eq_con1_fun) == 0.0
        @test MOI.constant(eq_con1_set) == 1.0
        eq_con2_fun = MOI.get(dual_model, MOI.ConstraintFunction(), CI{SAF{Float64}, MOI.EqualTo{Float64}}(2))
        eq_con2_set = MOI.get(dual_model, MOI.ConstraintSet(), CI{SAF{Float64}, MOI.EqualTo{Float64}}(2))
        @test MOI.coefficient.(eq_con2_fun.terms) == [1.0; 1.0; -1.0; -2.0]
        @test MOI.constant.(eq_con2_fun) == 0.0
        @test MOI.constant(eq_con2_set) == 1.0

        primal_con_dual_var = primal_dual_map.primal_con_dual_var
        @test primal_con_dual_var[CI{SAF{Float64}, MOI.EqualTo{Float64}}(1)] == [VI(1)]
        @test primal_con_dual_var[CI{SVF, MOI.GreaterThan{Float64}}(1)] == [VI(2)]
        @test primal_con_dual_var[CI{SVF, MOI.GreaterThan{Float64}}(2)] == [VI(3)]

        primal_var_dual_con = primal_dual_map.primal_var_dual_con
        @test primal_var_dual_con[VI(1)] == CI{SAF{Float64}, MOI.EqualTo{Float64}}(1)
        @test primal_var_dual_con[VI(2)] == CI{SAF{Float64}, MOI.EqualTo{Float64}}(2)

        primal_var_dual_quad_slack = primal_dual_map.primal_var_dual_quad_slack
        @test primal_var_dual_quad_slack[VI(1)] == VI(3+1)
        @test primal_var_dual_quad_slack[VI(2)] == VI(3+2)
    end
end


text:
import PyCall

# , Conda

# Conda.add("git")
# Conda.add("pip")
# Conda.add("cython")
# Conda.add("numpy")

# const pip = joinpath(Conda.BINDIR, "pip")

# proxy_arg = String[]
# if haskey(ENV, "http_proxy")
#     push!(proxy_arg, "--proxy")
#     push!(proxy_arg, ENV["http_proxy"])
# end

pip = `$(PyCall.python) -m pip`
run(`$pip install --user --upgrade pip setuptools`)
run(`$pip install --user cython numpy`)

# Install pyconcorde (The Concorde TSP Solver + QSOpt LP Library)
# https://github.com/jvkersch/pyconcorde
run(`$pip install git+git://github.com/jvkersch/pyconcorde.git`)

# Install elkai (The LKH Solver)
# https://github.com/filipArena/elkai
run(`$pip install elkai`)


text:
@polly function kernel_seidel_2d(tsteps, A)
    n = size(A,1)

    for t = 1:tsteps, i = 2:(n-1), j = 2:(n-1)
        A[i,j] = (A[i-1,j-1] + A[i-1,j] + A[i-1,j+1] + A[i,j-1] + A[i,j] + A[i,j+1] + A[i+1,j-1] + A[i+1,j] + A[i+1,j+1])/9.0
    end
end

let
    tsteps = 500
    n = 2000

    A = zeros(Float32, n, n)

    for i = 1:n, j = 1:n
        A[i,j] = (i*(j+2) + 2) / n
    end

    SUITE["seidel-2d"] = @benchmarkable kernel_seidel_2d($tsteps, A) setup = (A = copy($A))
end


text:
export width, height, maxwidth, minwidth, maxheight, minheight,
       boxsize, hbox, vbox, hskip, vskip, floating, wrap, pad,
       padinner, alignitems, aligncontent, justifycontent, clip,
       container, alignself, floating

for (fn, prop) in [:width => "width",
                   :height => "height",
                   :minwidth => "min-width",
                   :maxwidth => "max-width",
                   :minheight => "min-height",
                   :maxheight => "max-height",
                  ]
    @eval function $fn(val, elem=nothing)
        style(elem, $prop => val)
    end
end

function boxsize(w, h, elem=nothing)
    style(elem, "width"=>w, "height"=>h)
end

function flex(elem=nothing)
    style(elem, "display"=>"flex")
end

function container(xs::AbstractVector)
    dom"div"(xs...)
end

container(xs...) = container([xs...])

"""
    hbox(el...)

    Horizontally align mulitple web components.
"""
function hbox(elems::AbstractVector)
    container(elems)(style("display" => "flex", "flex-direction"=>"row"))
end
hbox(xs...) = hbox([xs...])

"""
    vbox(el...)

    Vertically align mulitple web components.
"""
function vbox(elems::AbstractVector)
    container(elems)(style("display" => "flex", "flex-direction"=>"column"))
end
vbox(xs...) = vbox([xs...])

hskip(x) = boxsize(x, 0px, empty())
vskip(y) = boxsize(0px, y, empty())

function wrap(elem=nothing; reverse=false)
    style(elem, "flex-wrap" => reverse ? "wrap-reverse" : "wrap")
end

function floating(direction, elem=nothing)
    style(elem, float=direction)
end

function justifycontent(justification, elem=nothing)
    style(elem, "justify-content" => justification)
end

function alignitems(alignment, elem=nothing)
    style(elem, "align-items" => alignment)
end

function alignself(alignment, elem=nothing)
    style(elem, "align-self" => alignment)
end

function aligncontent(alignment, elem=nothing)
    style(elem, "align-content" => alignment)
end

function _pad(side::Union{String,Symbol}, length)
    assertoneof(side, ["top", "bottom", "left", "right"], "side")
    "padding-$side" => length
end

function padinner(side::Union{String,Symbol}, length, elem=nothing)
    style(elem, _pad(side, length))
end

function padinner(sides::Union{Vector, Tuple}, length, elem=nothing)
    style(elem, (_pad(side, length) for side in sides)...)
end

function padinner(length, elem=nothing)
    style(elem, "padding" => length)
end

function pad(sides, length, elem)
    padinner(sides, length, container(elem))
end

function pad(length, elem)
    padinner(length, container(elem))
end

function clip(overflow, elem)
    assertoneof(overflow, ["hidden", "visible", "scroll", "auto"], "overflow")
    style(container(elem), "overflow" => overflow)
end


text:
import NumericExtensions

################################################################################
#
#   Type for keeping track of nearby pairs
#
################################################################################
type MolecularDynamicsIndex
    pairs::Array{(Int64,Int64),1}               # tuple vector of nearby pairs
    cells::Array{Array{Int64,1},3}              # grid containers for particles
    δ::Float64                                  # max distance worth including
    λ::Int64                                    # steps before index regen
    Λ::Int64                                    # current step
end

################################################################################
#
#   Convenience constructor for MolecularDynamicsIndex
#
################################################################################
MolecularDynamicsIndex( δ, λ ) =
    MolecularDynamicsIndex(
        Array((Int64,Int64),0),
        Array(Array{Int64,1},0,0,0),
        δ,
        λ,
        10 )

################################################################################
#
#   Function for finding the shortest vector twixt two points
#
#       Calculate minimum squared distance, accounting for periodicity of the
#       containing box.
#
################################################################################
function Δ( p2::Array{Float64,1}, p1::Array{Float64,1}, l::Float64, Δp::Array{Float64,1} )
    Δp=p2-p1
    lo2=l÷2
    for i in 1:3
        if Δp[i] > lo2
            Δp[i] -= l
        elseif Δp[i] < lo2
            Δp[i] += l
        end
    end
    return Δp
end

################################################################################
#
#   Simple function for finding nearby pairs
#
#       Iterates through all pairs and, if a pair is sufficiently close, adds
#       it to the list of interacting pairs as a tuple of the particle indices.
#
#       1) If the index counter Λ has reached the planned index lifetime λ,
#          rebuild the index.
#       2) Find the square of the max interaction distance δ; call it δsq
#       3) Clear the list of interacting pairs and reinitialize it as empty.
#       4) Loop through all ordered pairs of particles (j,k), skipping j==k
#       5) Find Δr⃗, the shortest distance between particles j and k
#       6) If the squared magnitude of Δr⃗ is less that δsq, the particles are
#          close enough to interact; push the tuple (j,k) into the pairs list.
#       7) Increment Λ so that it rolls back to 1 on passing λ.
#
################################################################################
function index!( y::Array{Float64, 2}, i::MolecularDynamicsIndex, l::Float64 )
    if i.Λ == i.λ
        print("Rebuilding the index... ")
        δsq = i.δ^2
        i.pairs = Array((Int64,Int64),0)
        for j in 1:size(y)[2]
            for k in 1:size(y)[2]
                if j == k
                    continue
                else
                    Δr⃗ = Δ(y[:,j],y[:,k],l)
                    if δsq >= dot(Δr⃗,Δr⃗)
                        push!( i.pairs, (j,k) )
        end; end; end; end
        print("done. ")
    end
    i.Λ = ( i.Λ % i.λ ) + 1
end

################################################################################
#
#   Function for finding f⃗, PE, PEtotal, and the Virial at step n
#
#       Calculate kinetic energy by iterating through the interaction list,
#       calculating rSquared, PE, and F for each interacting body. The
#       interaction list holds a list of (Int64,Int64) specifying the indices
#       of the current particle as well as a nearby neighbor.
#
#       Use Δ(r1,r2,l) to find the shortest distance between pairs, taking
#       account of periodic boundary conditions. Then calculate dot(Δr⃗,Δr⃗),
#       which gives the square of the distance between the pairs, which
#       conveniently is the only value we need for subsequent calculations.
#
#       Increment f⃗, PE, PEtot, and the Virial by appropriate amounts. Save a
#       bit of time by holding off on multiplying the virial by 1/2 until all
#       the contributions to it have been added in.
#
################################################################################
function potentialenergyandforce!( r::MolecularDynamicsTrial, n::Int64 )
    m = r.numBodies
    l = r.L
    x = r.y[:,:,n]

    ΔR⃗ = reshape(x,3,m) .- reshape(x,3,1,m) # All interatomic distances, aperiodic
    ΔR⃗ *= (1/l)                             # Scale out box length and shift to
    ΔR⃗ += 0.5                               # middle, so that 0 < ΔR⃗min < 1.
    ΔR⃗ = mod(ΔR⃗,1)                          # Pick values that match this criterion.
    ΔR⃗ -= 0.5                               # Shift values back into place, and
    ΔR⃗ *= l                                 # scale them back up.

    Rn2 = NumericExtensions.sumsq(ΔR⃗,1)     # All the minimum square distances;
    Rn2 = 1 ./ Rn2                          # We will use Rn2 = R^-2, hence the name
    for i in 1:m                            # Avoid infinite self-interactions
        Rn2[1,i,i] = 0.0
    end
    Rn6 = Rn2 .* Rn2 .* Rn2                 # Rn6 = Rn^-6, which will be handy

    ΔPE = 4(Rn6 .* Rn6 - Rn6)
    Δvir = 0.25Rn6 + 0.125ΔPE
    Δf⃗ = 2Rn2 .* (Δvir .* ΔR⃗)

    r.pe[:,n] = sum(ΔPE,2)
    r.pet[n] = sum(r.pe[:,n])
    r.vir[n] = sum(Δvir)
    r.f[:,:,n] = sum(Δf⃗,3)
end

################################################################################
#
#   Calculate just PE and the virial
#
################################################################################
function potentialenergy!( r::MolecularDynamicsTrial, n::Int64 )
    m = r.numBodies
    l = r.L
    x = r.y[:,:,n]

    ΔR⃗ = reshape(x,3,m) .- reshape(x,3,1,m) # All interatomic distances, aperiodic
    ΔR⃗ *= (1/l)                             # Scale out box length and shift to
    ΔR⃗ += 0.5                               # middle, so that 0 < ΔR⃗min < 1.
    ΔR⃗ = mod(ΔR⃗,1)                          # Pick values that match this criterion.
    ΔR⃗ -= 0.5                               # Shift values back into place, and
    ΔR⃗ *= l                                 # scale them back up.

    Rn2 = NumericExtensions.sumsq(ΔR⃗,1)     # All the minimum square distances;
    Rn2 = 1 ./ Rn2                          # We will use Rn2 = R^-2, hence the name
    for i in 1:m                            # Avoid infinite self-interactions
        Rn2[1,i,i] = 0.0
    end
    Rn6 = Rn2 .* Rn2 .* Rn2                 # Rn6 = Rn^-6, which will be handy

    ΔPE = 4(Rn6 .* Rn6 - Rn6)
    Δvir = 0.25Rn6 + 0.125ΔPE

    r.pe[:,n] = sum(ΔPE,2)
    r.pet[n] = sum(r.pe[:,n])
    r.vir[n] = sum(Δvir)
end

################################################################################
#
#   Function for finding KE, KEtotal, and T at step n
#
#       In Verlet's units, KE = 24v^2 and T = (2/3)KEavg = (2/3)KEtotal/N.
#
#       1) Calculate KE for each particle
#       2) Sum up the KE to get KEtotal
#       3) Use KEtotal to find the temperature, T
#
################################################################################
function kineticenergy!( r::MolecularDynamicsTrial, n::Int64 )
    print("Calculating kinetic energies... ")
    for i=1:r.numBodies
        r.ke[i,n] = 0.5dot( r.v[:,i,n], r.v[:,i,n] )
    end
    r.ket[n] = sum( r.ke[:,n] )
    r.T[n] = (2/3)r.ket[n] / r.numBodies
    print("done. ")
end


text:
function _unpack_kwargs(;
    dim_max=1,
    cutoff=0,
    modulus=2,
    field=Mod{modulus},
    verbose=false,
    alg=:cohomology,
    reps=alg == :cohomology ? false : 1:dim_max,
    implicit=alg != :homology,
    # deprecated
    progress=nothing,
    field_type=nothing,
    #
    filtration_kwargs...,
)
    if !isnothing(progress)
        @warn "`progress` is deprecated. Use `verbose` instead" maxlog = 1
        verbose = progress
    end
    if !isnothing(field_type)
        @warn "`field_type` is deprecated. Use `field` instead" maxlog = 1
        field = field_type
    end
    ripserer_kwargs = (
        dim_max=dim_max,
        cutoff=cutoff,
        field=field,
        verbose=verbose,
        alg=alg,
        reps=reps,
        implicit=implicit,
    )

    return ripserer_kwargs, filtration_kwargs
end

"""
    ripserer(Type{<:AbstractFiltration}, args...; kwargs...)
    ripserer(filtration::AbstractFiltration; kwargs...)

Compute the persistent homology of a filtration. The filtration can be given as an
[`AbstractFiltration`](@ref) type, followed by its arguments, or as an initialized object
(see examples below). If only data is given, `Rips` is used by default.

Returns a `Vector` of [`PersistenceDiagram`](@ref)s with (`dim_max` + 1) elements. The
diagrams are sorted by dimension; the first element of the result is the 0-dimensional
diagram, and the last is the (`dim_max`)-dimensional diagram.

# Keyword Arguments

* `dim_max`: compute persistent homology up to this dimension. Defaults to `1`.

* `modulus`: compute persistent homology with coefficients in the prime field of integers
  mod `modulus`. Defaults to `2`.

* `field`: use this type of field of coefficients. Defaults to
  [`Ripserer.Mod`](@ref)`{modulus}`.

* `threshold`: compute persistent homology up to diameter smaller than threshold. Note that
  this parameter is passed to the filtration constructor. When using low thresholds with
  Rips filtrations, consider setting `sparse=true` for optimal performance.

* `cutoff`: only keep intervals with `persistence(interval) > cutoff`. Defaults to `0`. When
  `cutoff < 0`, the result will also contain zero-length intervals.

* `reps`: if `true`, attach representative (co)cycles to persistence intervals. Can also be
  set to collection of integers to only find representatives in specified dimensions,
  e.g. `reps=1:2` will only find representatives in dimensions 1 and 2. This is useful for
  large filtrations (such as cubical) where calculating zero-dimensional representatives can
  be very slow. Defaults to `false` for cohomology and `1:dim_max` for homology.
  Representatives are wrapped in a [`Chain`](@ref).

* `verbose`: If `true`, show a verbose bar. Defaults to `false`.

* `alg`: select the algorithm used in computation. The options are:

  - `:cohomology`: Default and fastest algorithm. When `reps` is set, intervals are equipped
    with representative _co_cycles.

  - `:homology`: Significantly slower than `:cohomology`, but finds representative
    cycles. Does not find infinite intervals beyond dimension 0.

  - `:involuted`: Use cohomology result to compute representative cycles. Can be extremely
    efficient compared to `:homology`, especially with `Rips` filtrations. See [this
    paper](https://arxiv.org/abs/2105.03629) for more information.

* `implicit`: If `true`, an implicit reduction algorithm is used. Defaults to `true` for
  :cohomology and `:involuted`, and `false` for `:homology`. `implicit=false` is not
  recommended for `:cohomology` because it disables the emergent pairs optimization.

Other `kwargs...` are passed to the filtration.

# Examples

```jldoctest
julia> ts = range(0, 2π; length=20)[1:(end - 1)];

julia> X = [((2 + cos(θ)) * cos(φ), (2 + cos(θ)) * sin(φ), sin(θ)) for θ in ts for φ in ts];

julia> ripserer(X)
2-element Vector{PersistenceDiagramsBase.PersistenceDiagram}:
 361-element 0-dimensional PersistenceDiagram
 362-element 1-dimensional PersistenceDiagram

julia> ripserer(EdgeCollapsedRips, X; modulus=7, threshold=2)
2-element Vector{PersistenceDiagramsBase.PersistenceDiagram}:
 361-element 0-dimensional PersistenceDiagram
 362-element 1-dimensional PersistenceDiagram

julia> ripserer(Rips(X; threshold=1); alg=:involuted)
2-element Vector{PersistenceDiagramsBase.PersistenceDiagram}:
 361-element 0-dimensional PersistenceDiagram
 362-element 1-dimensional PersistenceDiagram

```

"""
function ripserer(F::Type, args...; kwargs...)
    ripserer_kwargs, filtration_kwargs = _unpack_kwargs(; kwargs...)
    start_time = time_ns()
    filtration = F(args...; verbose=ripserer_kwargs.verbose, filtration_kwargs...)
    return _ripserer(filtration, start_time; ripserer_kwargs...)
end

# Default: Rips
ripserer(dist; kwargs...) = ripserer(Rips, dist; kwargs...)

function ripserer(filtration::AbstractFiltration; kwargs...)
    ripserer_kwargs, other_kwargs = _unpack_kwargs(; kwargs...)
    start_time = time_ns()
    return _ripserer(filtration, start_time; ripserer_kwargs..., other_kwargs...)
end

_reps(reps::Bool, _) = reps
_reps(reps, dim) = dim in reps

function _ripserer(
    filtration::AbstractFiltration,
    start_time;
    dim_max,
    cutoff,
    field,
    verbose,
    alg,
    reps,
    implicit,
)
    if field <: Union{Signed,Unsigned,AbstractFloat}
        error("$field is not a field! Please try a differnet field type")
    end
    index_overflow_check(filtration, field, dim_max)
    result = _ripserer(
        Val(alg), filtration, cutoff, verbose, field, dim_max, reps, implicit
    )
    elapsed = round((time_ns() - start_time) / 1e9; digits=3)
    @prog_println verbose "Done. Time: " ProgressMeter.durationstring(elapsed)

    return result
end

function _ripserer(
    ::Val{:cohomology}, filtration, cutoff, verbose, field, dim_max, reps, implicit
)
    result = PersistenceDiagram[]
    zeroth, to_reduce, to_skip = zeroth_intervals(
        filtration, cutoff, verbose, field, _reps(reps, 0)
    )
    push!(result, zeroth)
    if dim_max > 0
        matrix = CoboundaryMatrix{implicit}(field, filtration, to_reduce, to_skip)
        for dim in 1:dim_max
            push!(result, compute_intervals!(matrix, cutoff, verbose, _reps(reps, dim)))
            if dim < dim_max
                matrix = next_matrix(matrix, verbose)
            end
        end
    end
    return result
end

function _ripserer(
    ::Val{:homology}, filtration, cutoff, verbose, field, dim_max, reps, implicit
)
    result = PersistenceDiagram[]
    zeroth, to_reduce, to_skip = zeroth_intervals(
        filtration, cutoff, verbose, field, _reps(reps, 0)
    )
    push!(result, zeroth)
    if dim_max > 0
        simplices = columns_to_reduce(filtration, Iterators.flatten((to_reduce, to_skip)))
        for dim in 1:dim_max
            matrix = BoundaryMatrix{implicit}(field, filtration, simplices)
            push!(result, compute_intervals!(matrix, cutoff, verbose, _reps(reps, dim)))
            if dim < dim_max
                simplices = columns_to_reduce(filtration, simplices)
            end
        end
    end
    return result
end

function _ripserer(
    ::Val{:involuted}, filtration, cutoff, verbose, field, dim_max, reps, implicit
)
    result = PersistenceDiagram[]
    zeroth, to_reduce, to_skip = zeroth_intervals(
        filtration, cutoff, verbose, field, _reps(reps, 0)
    )
    push!(result, zeroth)
    if dim_max > 0
        comatrix = CoboundaryMatrix{true}(field, filtration, to_reduce, to_skip)
        for dim in 1:dim_max
            columns, inf_births = compute_death_simplices!(comatrix, verbose, cutoff)
            matrix = BoundaryMatrix{implicit}(field, filtration, columns)
            diagram = compute_intervals!(matrix, cutoff, verbose, _reps(reps, dim))
            for birth_simplex in inf_births
                push!(
                    diagram.intervals,
                    interval(comatrix, birth_simplex, nothing, 0, _reps(reps, dim)),
                )
            end
            push!(result, diagram)
            if dim < dim_max
                comatrix = next_matrix(comatrix, verbose)
            end
        end
    end
    return result
end

function _ripserer(::Val{A}, args...) where {A}
    return throw(ArgumentError("unsupported alg=$A"))
end


text:
import RRSTExperiments: InputEdge

function solve_rec_st_hurwicz(n::Int,
    E::Vector{InputEdge},
    k::Int,
    λ::Float64)

    V = collect(1:n) # set of nodes
    Vminus1 = setdiff(V, [1]) # commodity nodes
    A = [(e.i, e.j) for e ∈ E] ∪ [(e.j, e.i) for e ∈ E]
    L = n-1-k

    # Model
    model = Model(CPLEX.Optimizer)
    set_silent(model)
    set_optimizer_attribute(model, "CPX_PARAM_EPINT", 1e-8)

    # Variables
    @variable(model, fx[A, Vminus1] ≥ 0)
    @variable(model, wx[A] ≥ 0)
    @variable(model, x[E] ≥ 0, Bin)

    @variable(model, fy¹[A, Vminus1] ≥ 0)
    @variable(model, wy¹[A] ≥ 0)
    @variable(model, y¹[E] ≥ 0, Bin)
    @variable(model, z¹[E] ≥ 0)

    @variable(model, fy²[A, Vminus1] ≥ 0)
    @variable(model, wy²[A] ≥ 0)
    @variable(model, y²[E] ≥ 0, Bin)
    @variable(model, z²[E] ≥ 0)

    # Objective
    @objective(model, Min, sum(e.C * x[e] for e ∈ E) + λ * sum(e.c * y¹[e] for e ∈ E) + (1-λ) * sum((e.c + e.d) * y²[e] for e ∈ E))

    #Constraints
    for k ∈ Vminus1 # sources
        @constraint(model, sum(fx[(j,1),k] for j ∈ filter(j -> (j,1) ∈ A, V)) -
        sum(fx[(1,j),k] for j ∈ filter(j -> (1,j) ∈ A,V)) == -1)
        @constraint(model, sum(fy¹[(j,1),k] for j ∈ filter(j -> (j,1) ∈ A, V)) -
        sum(fy¹[(1,j),k] for j ∈ filter(j -> (1,j) ∈ A,V)) == -1)
        @constraint(model, sum(fy²[(j,1),k] for j ∈ filter(j -> (j,1) ∈ A, V)) -
        sum(fy²[(1,j),k] for j ∈ filter(j -> (1,j) ∈ A,V)) == -1)
    end

    for k ∈ Vminus1, i ∈ Vminus1 # balances
        if i ≠ k
            @constraint(model, sum(fx[(j,i),k] for j ∈ filter(j -> (j,i) ∈ A, V)) -
            sum(fx[(i,j),k] for j ∈ filter(j -> (i,j) ∈ A, V)) == 0)
            @constraint(model, sum(fy¹[(j,i),k] for j ∈ filter(j -> (j,i) ∈ A, V)) -
            sum(fy¹[(i,j),k] for j ∈ filter(j -> (i,j) ∈ A, V)) == 0)
            @constraint(model, sum(fy²[(j,i),k] for j ∈ filter(j -> (j,i) ∈ A, V)) -
            sum(fy²[(i,j),k] for j ∈ filter(j -> (i,j) ∈ A, V)) == 0)
        end
    end

    for k ∈ Vminus1 # sinks
        @constraint(model, sum(fx[(j,k),k] for j ∈ filter(j -> (j,k) ∈ A,V)) -
        sum(fx[(k,j),k] for j ∈ filter(j -> (k,j) ∈ A,V)) == 1)
        @constraint(model, sum(fy¹[(j,k),k] for j ∈ filter(j -> (j,k) ∈ A,V)) -
        sum(fy¹[(k,j),k] for j ∈ filter(j -> (k,j) ∈ A,V)) == 1)
        @constraint(model, sum(fy²[(j,k),k] for j ∈ filter(j -> (j,k) ∈ A,V)) -
        sum(fy²[(k,j),k] for j ∈ filter(j -> (k,j) ∈ A,V)) == 1)
    end

    for k ∈ Vminus1, a ∈ A # capacity
        @constraint(model, fx[a,k] ≤ wx[a])
        @constraint(model, fy¹[a,k] ≤ wy¹[a])
        @constraint(model, fy²[a,k] ≤ wy²[a])
    end

    # tree
    @constraint(model, sum(wx[a] for a ∈ A) == n - 1)
    @constraint(model, sum(wy¹[a] for a ∈ A) == n - 1)
    @constraint(model, sum(wy²[a] for a ∈ A) == n - 1)

    for e ∈ E
        @constraint(model, x[e] == wx[(e.i, e.j)] + wx[(e.j, e.i)])
        @constraint(model, y¹[e] == wy¹[(e.i, e.j)] + wy¹[(e.j, e.i)])
        @constraint(model, y²[e] == wy²[(e.i, e.j)] + wy²[(e.j, e.i)])
    end

    for e ∈ E
        @constraint(model, x[e] ≥ z¹[e])
        @constraint(model, x[e] ≥ z²[e])

        @constraint(model, y¹[e] ≥ z¹[e])
        @constraint(model, y²[e] ≥ z²[e])
    end

    @constraint(model, sum(z¹[e] for e ∈ E) ≥ L)
    @constraint(model, sum(z²[e] for e ∈ E) ≥ L)

    # Solve
    optimize!(model)

    status=termination_status(model)
    obj_value = objective_value(model)
    if status == MOI.OPTIMAL
        return status, obj_value, get_first_stage_solution(Array(value.(x)), E)
    else
        return status, missing, missing
    end
end

if false
    n = 6
    E = [
        InputEdge(1, 2, 0.0, 9.0, 1.0) # 1
        InputEdge(2, 3, 2.0, 5.0, 1.0) # 2
        InputEdge(3, 4, 1.0, 1.0, 1.0) # 3
        InputEdge(1, 4, 1.0, 2.0, 1.0) # 4
        InputEdge(2, 4, 5.0, 4.0, 1.0) # 5
        InputEdge(1, 5, 8.0, 3.0, 1.0) # 6
        InputEdge(4, 6, 0.0, 6.0, 1.0) # 7
        InputEdge(5, 6, 3.0, 2.0, 1.0) # 8
        InputEdge(5, 7, 5.0, 9.0, 1.0) # 9
        InputEdge(6, 7, 7.0, 2.0, 1.0) # 10
        InputEdge(1, 6, 2.0, 3.0, 1.0) # 11
    ]
    k = 1
    λ = 1.0
    status, obj_value, x, y¹, y² = solve_rec_st_hurwicz(n, E, k, λ)
end


text:
# type definition lint code

function linttype(ex::Expr, ctx::LintContext)
    if ctx.quoteLvl > 0
        return  # do not lint types in quotes, see issue 166
    end

    @checktoplevel(ctx, "type")
    typeparams = Symbol[]

    # TODO: this duplicates the code in functions.jl
    processCurly = (sube)->begin
        for i in 2:length(sube.args)
            adt = sube.args[i]
            if isa(adt, Symbol)
                foundobj = lookup(ctx, adt)
                if foundobj !== nothing && foundobj.typeactual <: Type
                    msg(ctx, :I393, adt, "using an existing type as type parameter name is probably a typo")
                end
                # TODO: review all uses of this function
                push!(typeparams, adt)
            elseif isexpr(adt, :(<:))
                temptype = adt.args[1]
                typeconstraint = adt.args[2]

                if temptype != :T
                    foundobj = lookup(ctx, temptype)
                    if foundobj !== nothing && foundobj.typeactual <: Type
                        msg(ctx, :E538, temptype, "known type in parametric data type, " *
                            "use {T<:...}")
                    end
                end
                if istype(ctx, typeconstraint)
                    dt = parsetype(ctx, typeconstraint)
                    if isa(dt, Type) && isconcretetype(dt)
                        msg(ctx, :E513, adt, "leaf type as a type constraint makes no sense")
                    end
                end
                push!(typeparams, temptype)
            end
        end
    end

    tname = Symbol("")
    if isa(ex.args[2], Symbol)
        tname = ex.args[2]
    elseif isexpr(ex.args[2], :curly)
        tname = ex.args[2].args[1]
        processCurly(ex.args[2])
    elseif isexpr(ex.args[2], :(<:))
        if isa(ex.args[2].args[1], Symbol)
            tname = ex.args[2].args[1]
        elseif isexpr(ex.args[2].args[1], :curly)
            tname = ex.args[2].args[1].args[1]
            processCurly(ex.args[2].args[1])
        end
    end
    if tname != Symbol("")
        if islower(string(tname)[1])
            msg(ctx, :I771, tname, "type names should start with an upper case")
        end
        set!(ctx.current, tname, VarInfo(location(ctx), Type))
    end

    fields = Any[]
    funcs = Any[]

    for def in ex.args[3].args
        if isa(def, LineNumberNode)
            ctx.line = ctx.lineabs + def.line-1
        elseif isa(def, Symbol)
            # it means Any, probably not a very efficient choice
            if !pragmaexists("Ignore untyped field $(def)", ctx.current, deep=false)
                msg(ctx, :I691, def, "a type is not given to the field which can be slow")
            end
            push!(fields, (def, Any))
        elseif isexpr(def, :macrocall) && def.args[1] == Symbol("@lintpragma")
            lintlintpragma(def, ctx)
        elseif isexpr(def, :call) && def.args[1] == Symbol("lintpragma")
            lintlintpragma(def, ctx)
            msg(ctx, :E425, "use lintpragma macro inside type declaration")
        elseif isexpr(def, :(::))
            if isexpr(def.args[2], :curly) && def.args[2].args[1] == :Array && length(def.args[2].args) <= 2 &&
                !pragmaexists("Ignore dimensionless array field $(def.args[1])", ctx.current, deep=false)
                msg(ctx, :I692, def.args[1], "array field has no dimension which can be slow")
            end
            push!(fields, (def.args[1], def.args[2]))
        elseif isexpr(def, :(=)) && isexpr(def.args[1], :call) || isexpr(def, :function)
            # curly bracket doesn't belong here. catch it first before linting the rest of the function body
            if def.args[1].head == :tuple
                # if julia supports anonymous constructor syntactic sugar, remove this, and make sure ctx.scope is type name
                msg(ctx, :E417, "anonymous function inside type definition")
            end
            push!(funcs, (def, ctx.line))
        end
    end

    for f in funcs
        ctx.line = ctx.lineabs + f[2]
        lintfunction(f[1], ctx; ctorType = tname)
    end
end

function linttypealias(ex::Expr, ctx::LintContext)
    # TODO: make this just part of lintassignment
    sym = withincurly(ex.args[1])
    @checkisa(ctx, sym, Symbol)
    set!(ctx.current, sym, VarInfo(location(ctx), Type))
end

function lintabstract(ex::Expr, ctx::LintContext)
    sym = withincurly(ex.args[1])
    if isa(sym, Symbol)
        set!(ctx.current, sym, VarInfo(location(ctx), Type))
    elseif isexpr(ex.args[1], :(<:))
        sym = withincurly(ex.args[1].args[1])
        @checkisa(ctx, sym, Symbol)
        set!(ctx.current, sym, VarInfo(location(ctx), Type))
    end
end

function lintbitstype(ex::Expr, ctx::LintContext)
    @checkisa(ctx, ex.args[2], Symbol)
    set!(ctx.current, ex.args[2], VarInfo(location(ctx), Type))
end


text:
# ------------------------------------------------------------------
# Copyright (c) 2015, Júlio Hoffimann Mendes <juliohm@stanford.edu>
# Licensed under the ISC License. See LICENCE in the project root.
# ------------------------------------------------------------------

module ImageQuilting

using ImageFiltering
using ImageMorphology
using LightGraphs
using LightGraphsFlows
using Base: @nexprs, @nloops, @nref
using Primes: factor
using Combinatorics: nthperm!
using StatsBase: sample, weights
using ProgressMeter: Progress, next!
using FFTW: set_num_threads
using CpuId: cpucores
using RecipesBase
using SparseArrays: spzeros
using Random: shuffle!
using Statistics: mean, std

# GeoStats.jl interface
using GeoStatsBase
using GeoStatsDevTools

import GeoStatsBase: preprocess, solve_single

include("utils.jl")
# include("utils_gpu.jl")
include("datatypes.jl")
include("plot_recipes.jl")
include("imfilter_cpu.jl")
# include("imfilter_gpu.jl")
include("relaxation.jl")
include("tau_model.jl")
include("dijkstra_cut.jl")
include("boykov_kolmogorov_cut.jl")
include("voxel_reuse.jl")
include("iqsim.jl")

include("geostats_api.jl")

export
  # functions
  iqsim,
  voxelreuse,

  # data types
  HardData,

  # geostats solver
  ImgQuilt

end


text:
# # Factorizations and other fun
# Based on work by Andreas Noack
#
# ## Outline
#  - Factorizations
#  - Special matrix structures
#  - Generic linear algebra

#-

# Before we get started, let's set up a linear system and use `LinearAlgebra` to bring in the factorizations and special matrix structures.

using LinearAlgebra
A = rand(3, 3)
x = fill(1, (3,))
b = A * x

# ## Factorizations
#
# #### LU factorizations
# In Julia we can perform an LU factorization
# ```julia
# PA = LU
# ```
# where `P` is a permutation matrix, `L` is lower triangular unit diagonal and `U` is upper triangular, using `lufact`.
#
# Julia allows computing the LU factorization and defines a composite factorization type for storing it.

Alu = lu(A)

#-

typeof(Alu)

# The different parts of the factorization can be extracted by accessing their special properties

Alu.P

#-

Alu.L

#-

Alu.U

# Julia can dispatch methods on factorization objects.
#
# For example, we can solve the linear system using either the original matrix or the factorization object.

A\b

#-

Alu\b

# Similarly, we can calculate the determinant of `A` using either `A` or the factorization object

det(A) ≈ det(Alu)

# #### QR factorizations
#
# In Julia we can perform a QR factorization
# ```
# A=QR
# ```
#
# where `Q` is unitary/orthogonal and `R` is upper triangular, using `qrfact`.

Aqr = qr(A)

# Similarly to the LU factorization, the matrices `Q` and `R` can be extracted from the QR factorization object via

Aqr.Q

#-

Aqr.R

# #### Eigendecompositions

#-

# The results from eigendecompositions, singular value decompositions, Hessenberg factorizations, and Schur decompositions are all stored in `Factorization` types.
#
# The eigendecomposition can be computed

Asym = A + A'
AsymEig = eigen(Asym)

# The values and the vectors can be extracted from the Eigen type by special indexing

AsymEig.values

#-

AsymEig.vectors

# Once again, when the factorization is stored in a type, we can dispatch on it and write specialized methods that exploit the properties of the factorization, e.g. that $A^{-1}=(V\Lambda V^{-1})^{-1}=V\Lambda^{-1}V^{-1}$.

inv(AsymEig)*Asym

# ## Special matrix structures
# Matrix structure is very important in linear algebra. To see *how* important it is, let's work with a larger linear system

n = 1000
A = randn(n,n);

# Julia can often infer special matrix structure

Asym = A + A'
issymmetric(Asym)

# but sometimes floating point error might get in the way.

Asym_noisy = copy(Asym)
Asym_noisy[1,2] += 5eps()

#-

issymmetric(Asym_noisy)

# Luckily we can declare structure explicitly with, for example, `Diagonal`, `Triangular`, `Symmetric`, `Hermitian`, `Tridiagonal` and `SymTridiagonal`.

Asym_explicit = Symmetric(Asym_noisy);

# Let's compare how long it takes Julia to compute the eigenvalues of `Asym`, `Asym_noisy`, and `Asym_explicit`

@time eigvals(Asym);

#-

@time eigvals(Asym_noisy);

#-

@time eigvals(Asym_explicit);

# In this example, using `Symmetric()` on `Asym_noisy` made our calculations about `5x` more efficient :)

#-

# #### A big problem
# Using the `Tridiagonal` and `SymTridiagonal` types to store tridiagonal matrices makes it possible to work with potentially very large tridiagonal problems. The following problem would not be possible to solve on a laptop if the matrix had to be stored as a (dense) `Matrix` type.

n = 1_000_000;
A = SymTridiagonal(randn(n), randn(n-1));
@time eigmax(A)

# ## Generic linear algebra
# The usual way of adding support for numerical linear algebra is by wrapping BLAS and LAPACK subroutines. For matrices with elements of `Float32`, `Float64`, `Complex{Float32}` or `Complex{Float64}` this is also what Julia does.
#
# However, Julia also supports generic linear algebra, allowing you to, for example, work with matrices and vectors of rational numbers.

#-

# #### Rational numbers
# Julia has rational numbers built in. To construct a rational number, use double forward slashes:

1//2

# #### Example: Rational linear system of equations
# The following example shows how linear system of equations with rational elements can be solved without promoting to floating point element types. Overflow can easily become a problem when working with rational numbers so we use `BigInt`s.

Arational = Matrix{Rational{BigInt}}(rand(1:10, 3, 3))/10

#-

x = fill(1, 3)
b = Arational*x

#-

Arational\b

#-

lu(Arational)

# ### Exercises
#
# #### 11.1
# What are the eigenvalues of matrix A?
#
# ```
# A =
# [
#  140   97   74  168  131
#   97  106   89  131   36
#   74   89  152  144   71
#  168  131  144   54  142
#  131   36   71  142   36
# ]
# ```
# and assign it a variable `A_eigv`

using LinearAlgebra

#-

@assert A_eigv ==  [-128.49322764802145, -55.887784553056875, 42.7521672793189, 87.16111477514521, 542.4677301466143]

# #### 11.2
# Create a `Diagonal` matrix from the eigenvalues of `A`.

@assert A_diag ==  [-128.493    0.0      0.0      0.0       0.0;
    0.0    -55.8878   0.0      0.0       0.0;
    0.0      0.0     42.7522   0.0       0.0;
    0.0      0.0      0.0     87.1611    0.0;
    0.0 0.0      0.0      0.0     542.468]
## #### 11.3
## Create a `LowerTriangular` matrix from `A` and store it in `A_lowertri`
@assert A_lowertri ==  [140    0    0    0   0;
  97  106    0    0   0;
  74   89  152    0   0;
 168  131  144   54   0;
 131   36   71  142  36]

max_stars_repo_name: fercarozzi/JuliaAcademyMaterials

max_issues_repo_name: fercarozzi/JuliaAcademyMaterials

max_forks_repo_name: fercarozzi/JuliaAcademyMaterials


text:
function loglikelihood(tree::Tree,
                       mod::T,
                       node_data::N;
                       output_calculations::Bool=false) where {T <: NASM, N <: Union{NodeDNA, NodeRNA}}
  # Error checking
  if !all(map(x -> x in keys(node_data), findleaves(tree)))
    error("Some leaves are missing sequence data")
  elseif length(findroots(tree)) > 1
    error("More than one root detected")
  end

  # Create a Dict to store likelihood calculations
  calculations = Dict{Int64, Array{Float64, 2}}()

  # Find node visit order for postorder traversal
  visit_order = postorder(tree)
  for i in visit_order
    if isleaf(tree, i)
      calculations[i] = node_data[i].data
    else
      branches = tree.nodes[i].out
      for j in branches
        branch_length = tree.branches[j].length
        child_node = tree.branches[j].target
        p = P(mod, branch_length)
        if !haskey(calculations, i)
          calculations[i] = p * calculations[child_node]
        else
          calculations[i] .*= p * calculations[child_node]
        end
      end
    end
  end
  if output_calculations
    return sum(log.(_π(mod)' * calculations[visit_order[end]])), calculations, visit_order
  else
    return sum(log.(_π(mod)' * calculations[visit_order[end]]))
  end
end


text:
using NicePipes
using Documenter

makedocs(;
    modules=[NicePipes],
    authors="Simeon Schaub <simeondavidschaub99@gmail.com> and contributors",
    repo="https://github.com/simeonschaub/NicePipes.jl/blob/{commit}{path}#L{line}",
    sitename="NicePipes.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://simeonschaub.github.io/NicePipes.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/simeonschaub/NicePipes.jl",
)


text:
export OnePlusEvo, evaluate, populate

"""

OnePlusEvo implements a 1+λ evolutionary algorithm. To do this, a new subtype of
AbstractEvolution is created, OnePlusEvo, and the following functions are
defined for this type:

evaluate
populate

the function generation can also be defined. Without specification, the default
generation function (which is void) is applied.
"""
mutable struct OnePlusEvo{T} <: AbstractEvolution
    config::NamedTuple
    logger::CambrianLogger
    population::Array{T}
    fitness::Function
    gen::Int
end

function OnePlusEvo{T}(cfg::NamedTuple, fitness::Function;
                      logfile=string("logs/", cfg.id, ".csv")) where T
    logger = CambrianLogger(logfile)
    population = initialize(T, cfg)
    OnePlusEvo(cfg, logger, population, fitness, 0)
end

function oneplus_populate(e::AbstractEvolution)
    p1 = max_selection(e.population)
    e.population[1] = p1
    for i in 2:e.config.n_population
        e.population[i] = mutate(p1)
    end
end

evaluate(e::OnePlusEvo) = fitness_evaluate(e, e.fitness)
populate(e::OnePlusEvo) = oneplus_populate(e)


text:
# ---
# title: 627. Swap Salary
# id: problem627
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories:
# link: <https://leetcode.com/problems/swap-salary/description/>
# hidden: true
# ---
#
# Given a table `salary`, such as the one below, that has m=male and f=female
# values. Swap all f and m values (i.e., change all f values to m and vice
# versa) with a **single update statement** and no intermediate temp table.
#
# Note that you must write a single update statement, **DO NOT** write any
# select statement for this problem.
#
#
#
# **Example:**
#
#
#
#     | id | name | sex | salary |
#     |----|------|-----|--------|
#     | 1  | A    | m   | 2500   |
#     | 2  | B    | f   | 1500   |
#     | 3  | C    | m   | 5500   |
#     | 4  | D    | f   | 500    |
#
#
# After running your **update** statement, the above salary table should have
# the following rows:
#
#
#
#     | id | name | sex | salary |
#     |----|------|-----|--------|
#     | 1  | A    | f   | 2500   |
#     | 2  | B    | m   | 1500   |
#     | 3  | C    | f   | 5500   |
#     | 4  | D    | m   | 500    |
#
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
__precompile__()

module ImagineFormat

using Images, FileIO, Compat
using SIUnits, SIUnits.ShortUnits
# using FileIO: skipmagic, stream, @format_str, Stream

export imagine2nrrd, Micron

Micron = SIUnits.NonSIUnit{typeof(Meter),:µm}()
convert(::Type{SIUnits.SIQuantity},::typeof(Micron)) = Micro*Meter

function load(f::File{format"Imagine"}; mode="r")
    fabs = File(format"Imagine", abspath(f.filename))
    open(fabs) do s
        skipmagic(s)
        load(s, mode=mode)
    end
end

function load(io::Stream{format"Imagine"}; mode="r")
    s = stream(io)
    h = parse_header(s)
    filename = s.name[7:end-1]
    basename, ext = splitext(filename)
    camfilename = basename*".cam"
    T = h["pixel data type"]
    sz = [h["image width"], h["image height"], h["frames per stack"], h["nStacks"]]
    if sz[4] == 1
        sz = sz[1:3]
        if sz[3] == 1
            sz = sz[1:2]
        end
    end
    havez = h["frames per stack"] > 1
    havet = h["nStacks"] > 1
    # Check that the file size is consistent with the expected size
    if !isfile(camfilename)
        warn("Cannot open ", camfilename)
        data = Array(T, sz[1], sz[2], sz[3], 0)
    else
        fsz = filesize(camfilename)
        n_stacks = sz[end]
        if fsz != sizeof(T)*prod(map(Int64,sz))  # guard against overflow on 32bit
            warn("Size of image file is different from expected value")
            n_stacks = ifloor(fsz / sizeof(T) / prod(sz[1:end-1]))
        end
        if sizeof(T)*prod(map(Int64,sz[1:end-1]))*n_stacks > typemax(UInt)
            warn("File size is too big to mmap on 32bit")
            n_stacks = ifloor(fsz / sizeof(T) / typemax(UInt))
        end
        if n_stacks < sz[end]
            println("Truncating to ", n_stacks, length(sz) == 4 ? " stacks" : " frames")
            sz[end] = n_stacks
        end
        data = SharedArray(camfilename, T, tuple(sz...), mode=mode)
    end
    um_per_pixel = h["um per pixel"]*µm
    pstart = h["piezo"]["stop position"]
    pstop = h["piezo"]["start position"]
    if length(sz)>2
        dz = abs(pstart - pstop)/sz[3]
    else dz = 0.0 end

    props = Dict(
        "spatialorder" => (havez ? ["x", "l", "z"] : ["x", "l"]),
        "colorspace" => "Gray",
        "pixelspacing" => (havez ? [um_per_pixel, um_per_pixel, dz] : [um_per_pixel, um_per_pixel]),
        "limits" => (UInt16(0), UInt16(2^h["original image depth"]-1)),
        "imagineheader" => h,
        "suppress" => Set(Any["imagineheader"]))
    if havet
        props["timedim"] = havez ? 4 : 3
    end
    Image(data, props)
end

abstract Endian
type LittleEndian <: Endian; end
type BigEndian <: Endian; end
const endian_dict = Dict("l"=>LittleEndian, "b"=>BigEndian)
const nrrd_endian_dict = Dict(LittleEndian=>"little",BigEndian=>"big")
parse_endian(s::AbstractString) = endian_dict[lowercase(s)]

function parse_vector_int(s::AbstractString)
    ss = split(s, r"[ ,;]", keep=false)
    v = Array(Int, length(ss))
    for i = 1:length(ss)
        v[i] = parse(Int,ss[i])
    end
    return v
end

const bitname_dict = Dict(
  "int8"      => Int8,
  "uint8"     => UInt8,
  "int16"     => Int16,
  "uint16"    => UInt16,
  "int32"     => Int32,
  "uint32"    => UInt32,
  "int64"     => Int64,
  "uint64"    => UInt64,
  "float16"   => Float16,
  "float32"   => Float32,
  "single"    => Float32,
  "float64"   => Float64,
  "double"    => Float64)

parse_bittypename(s::AbstractString) = bitname_dict[lowercase(s)]

function float64_or_empty(s::AbstractString)
    if isempty(s)
        return NaN
    else
        return parse(Float64,s)
    end
end

function parse_quantity_or_empty(s::AbstractString)
    if isempty(s)
        return NaN
    else
        return parse_quantity(s)
    end
end

_unit_string_dict = Dict("um" => Micro*Meter, "s" => Second, "us" => Micro*Second, "MHz" => Mega*Hertz)
function parse_quantity(s::AbstractString, strict::Bool = true)
    # Find the last character of the numeric component
    m = match(r"[0-9\.\+-](?![0-9\.\+-])", s)
    if m == nothing
        error("AbstractString does not have a 'value unit' structure")
    end
    val = parse(Float64, s[1:m.offset])
    ustr = strip(s[m.offset+1:end])
    if isempty(ustr)
        if strict
            error("AbstractString does not have a 'value unit' structure")
        else
            return val
        end
    end
    val * _unit_string_dict[ustr]
end

# Read and parse a *.imagine file (an Imagine header file)
const compound_fields = Any["piezo", "binning"]
const field_key_dict = Dict{AbstractString,Function}(
    "header version"               => x->parse(Float64,x),
    "app version"                  => identity,
    "date and time"                => identity,
    "rig"                          => identity,
    "byte order"                   => parse_endian,
    "stimulus file content"        => identity,  # stimulus info parsed separately
    "comment"                      => identity,
    "ai data file"                 => identity,
    "image data file"              => identity,
    "start position"               => parse_quantity,
    "stop position"                => parse_quantity,
    "bidirection"                  => x->parse(Int,x) != 0,
    "output scan rate"             => x->parse_quantity(x, false),
    "nscans"                       => x->parse(Int,x),
    "channel list"                 => parse_vector_int,
    "label list"                   => identity,
    "scan rate"                    => x->parse_quantity(x, false),
    "min sample"                   => x->parse(Int,x),
    "max sample"                   => x->parse(Int,x),
    "min input"                    => x->parse(Float64,x),
    "max input"                    => x->parse(Float64,x),
    "original image depth"         => x->parse(Int,x),
    "saved image depth"            => x->parse(Int,x),
    "image width"                  => x->parse(Int,x),
    "image height"                 => x->parse(Int,x),
    "number of frames requested"   => x->parse(Int,x),
    "nStacks"                      => x->parse(Int,x),
    "idle time between stacks"     => parse_quantity,
    "pre amp gain"                 => float64_or_empty,
    "EM gain"                      => float64_or_empty,
    "gain"                         => float64_or_empty,
    "exposure time"                => parse_quantity,
    "vertical shift speed"         => parse_quantity_or_empty,
    "vertical clock vol amp"       => x->parse(Float64,x),
    "readout rate"                 => parse_quantity_or_empty,
    "pixel order"                  => identity,
    "frame index offset"           => x->parse(Int,x),
    "frames per stack"             => x->parse(Int,x),
    "pixel data type"              => parse_bittypename,
    "camera"                       => identity,
    "um per pixel"                 => x->parse(Float64,x),
    "hbin"                         => x->parse(Int,x),
    "vbin"                         => x->parse(Int,x),
    "hstart"                       => x->parse(Int,x),
    "hend"                         => x->parse(Int,x),
    "vstart"                       => x->parse(Int,x),
    "vend"                         => x->parse(Int,x),
    "angle from horizontal (deg)"  => float64_or_empty)

function parse_header(s::IOStream)
    headerdict = Dict{Compat.ASCIIString, Any}()
    for this_line = eachline(s)
        this_line = strip(this_line)
        if !isempty(this_line) && !ismatch(r"\[.*\]", this_line)
            # Split on =
            m = match(r"=", this_line)
            if m.offset < 2
                error("Line does not contain =")
            end
            k = this_line[1:m.offset-1]
            v = this_line[m.offset+1:end]
            if in(k, compound_fields)
                thisdict = Dict{Compat.ASCIIString, Any}()
                # Split on semicolon
                strs = split(v, r";")
                for i = 1:length(strs)
                    substrs = split(strs[i], r":")
                    @assert length(substrs) == 2
                    k2 = strip(substrs[1])
                    func = field_key_dict[k2]
                    v2 = strip(substrs[2])
                    try
                        thisdict[k2] = func(v2)
                    catch err
                        println("Error processing key '", k2, "' with value ", v2)
                        rethrow(err)
                    end
                end
                headerdict[k] = thisdict
            else
                func = field_key_dict[k]
                try
                    headerdict[k] = func(v)
                catch err
                    println("Error processing key ", k, " with value ", v)
                    rethrow(err)
                end
            end
        end
    end
    return headerdict
end
function parse_header(io::Stream{format"Imagine"})
    skipmagic(io)
    parse_header(stream(io))
end
function parse_header(f::File{format"Imagine"})
    open(f) do io
        parse_header(io)
    end
end
parse_header(filename::AbstractString) = parse_header(query(filename))

function imagine2nrrd(sheader::IO, h::Dict{Compat.ASCIIString, Any}, datafilename = nothing)
    println(sheader, "NRRD0001")
    T = h["pixel data type"]
    if T<:AbstractFloat
        println(sheader, "type: ", (T == Float32) ? "float" : "double")
    else
        println(sheader, "type: ", lowercase(string(T)))
    end
    sz = [h["image width"], h["image height"], h["frames per stack"], h["nStacks"]]
    kinds = ["space", "space", "space", "time"]
    if sz[end] == 1
        sz = sz[1:3]
        kinds = kinds[[1,2,4]]
    end
    println(sheader, "dimension: ", length(sz))
    print(sheader, "sizes:")
    for z in sz
        print(sheader, " ", z)
    end
    print(sheader, "\nkinds:")
    for k in kinds
        print(sheader, " ", k)
    end
    print(sheader, "\n")
    println(sheader, "encoding: raw")
    println(sheader, "endian: ", nrrd_endian_dict[h["byte order"]])
    if isa(datafilename, AbstractString)
        println(sheader, "data file: ", datafilename)
    end
    sheader
end

function imagine2nrrd(nrrdname::AbstractString, h::Dict{Compat.ASCIIString, Any}, datafilename = nothing)
    sheader = open(nrrdname, "w")
    imagine2nrrd(sheader, h, datafilename)
    close(sheader)
end

"""
`save_header(filename, header)` writes a header dictionary in Imagine format.

`save_header(destname, srcname, img::AbstractArray, [T::Type =
eltype(img)])` writes a `.imagine` file with name `destname`, using
the `.imagine` file `srcname` as a template. Size and element type
fields are updated from `img` and `T`, respectively.
"""
function save_header(filename::AbstractString, h::Dict{Compat.ASCIIString, Any})
    open(filename, "w") do io
        write(io, magic(format"Imagine"))
        println(io, "\n[general]")
        writekv(io, h, ("header version", "app version", "date and time", "byte order", "rig"))
        println(io, "\n[misc params]")
        writekv(io, h, ("stimulus file content", "comment", "ai data file","image data file", "piezo"))
        println(io, "\n[ai]")
        writekv(io, h, ("nscans", "channel list", "label list", "scan rate", "min sample", "max sample", "min input", "max input"))
        println(io, "\n[camera]")
        writekv(io, h, ("original image depth", "saved image depth", "image width", "image height", "number of frames requested", "nStacks", "idle time between stacks", "pre amp gain", "gain", "exposure time", "vertical shift speed", "vertical clock vol amp", "readout rate", "pixel order", "frame index offset", "frames per stack", "pixel data type", "camera", "um per pixel", "binning", "angle from horizontal (deg)"))
    end
    nothing
end

function save_header(dest::AbstractString, src::AbstractString, img::AbstractArray, T::Type = eltype(img))
    h = parse_header(src)
    fillsize!(h, img)
    h["pixel data type"] = lowercase(string(T))
    h["byte order"] = ENDIAN_BOM == 0x04030201 ? "l" : "b"
    save_header(dest, h)
end

function fillsize!(h, img::AbstractImage)
    h["image width"] = size(img,"x")
    h["image height"] = size(img,"l")
    h["frames per stack"] = size(img,"z")
    h["nStacks"] = size(img,"t")
    h
end

function fillsize!{T}(h, img::AbstractArray{T,4})
    h["image width"] = size(img,1)
    h["image height"] = size(img,2)
    h["frames per stack"] = size(img,3)
    h["nStacks"] = size(img,4)
    h
end

function writekv(io, h, fieldnames)
    for fn in fieldnames
        if haskey(h, fn)
            writefield(io, fn, h[fn])
        end
    end
end
writefield(io, fn, v) = println(io, fn, "=", v)
function writefield(io, fn, v::Number)
    print(io, fn, "=")
    if haskey(write_dict, fn)
        write_dict[fn](io, v)
    else
        println(io, isnan(v) ? "" : v)
    end
end
writefield(io, fn, v::AbstractVector) = println(io, fn, "=", join(map(string,v)," "))

function writefield(io, fn, dct::Dict)
    print(io, fn, "=")
    ks = collect(keys(dct))
    vs = collect(values(dct))
    for i = 1:length(ks)
        k, v = ks[i], vs[i]
        print(io, k, ':')
        if haskey(write_dict, k)
            write_dict[k](io, v)
        else
            print(io, v)
        end
        print(io, i == length(ks) ? '\n' : ';')
    end
end

writeum(io,x) = print(io, round(Int, 10^6*(x/Meter)), " um")
writeus(io,x::SIUnits.SIQuantity) = print(io, round(Int, 10^6*(x/Second)), " us")
writeus(io,x) = isnan(x) || print(io, round(Int, 10^6*(x/Second)), " us")
writeMHz(io,x::SIUnits.SIQuantity) = print(io, round(Int, 1e-6*(x*Second)), " MHz")
writeMHz(io,x) = isnan(x) || print(io, round(Int, 1e-6*(x*Second)), " MHz")
const write_dict = Dict{Compat.ASCIIString,Function}(
    "bidirection"                  => (io,x)->x ? print(io, 1) : print(io, 0),
    "start position"               => writeum,
    "stop position"                => writeum,
    "vertical shift speed"         => (io,x)->(writeus(io,x); print(io,'\n')),
    "readout rate"                 => (io,x)->(writeMHz(io,x); print(io,'\n')),
)

function __init__()
    Base.rehash!(nrrd_endian_dict)
end

end


text:
struct Basis{N,T<:AbstractMatrix}
    dims::NTuple{N}
    B::T
end

function construct_basis(bf, xs::Tuple, params)
    dims = length.(xs)
    nr = prod(dims)
    nc = length(params)
    B = zeros(nr,nc)
    w = trues(nc)
    for (i,p) in enumerate(params)
        Bc = vec(bf(xs,p))
        if sum(Bc) != 0.0
            B[:,i] .= vec(bf(xs,p))
            w[i] = true
        end
    end
    return Basis(dims,B[:,w])
end

function construct_basis(bf, xs, params)
    construct_basis(bf,(xs,),params)
end

Base.:(*)(B::Basis,A::AbstractMatrix) = B.B*A
Base.:(*)(A::AbstractMatrix,B::Basis) = A*B.B

function evaluate(B::Basis, coeffs)
    if length(coeffs) != size(B.B,2)
        @error "Wrong Number of Coefficients for basis"
    end
    reshape(max.(B.B*coeffs,0.0),B.dims)
end


text:
include("utils.jl")

@testset "Minkowski Metric" begin
    N = Euclidean(3)
    M = MetricManifold(N, MinkowskiMetric())
    @test local_metric(M, zeros(3)) == Diagonal([1.0, 1.0, -1.0])
end


text:
function tau(n::Integer)::Integer
    ans = 0
    for i = 1:n
        if n % i == 0
            ans = 1 + ans
        end
    end
    return ans
end

τf = tau

function sigma(n::Integer)::Integer
    ans = 0
    for i = 1:n
        if n % i == 0
            ans = i + ans
        end
    end
    return ans
end

σf = sigma

function s(n::Integer)::Integer
    return sigma(n) - n
end

# Adding cache to this function can make thing easier
function getMappingFromSFunction(n::Integer)::Array{Bool}
    sn = s(n)
    mappingSofNProperty = [sn < n, sn == n, sn > n]
    return mappingSofNProperty
end

function isDeficient(n::Integer)::Bool
    return getMappingFromSFunction(n)[1]
end

function isPerfect(n::Integer)::Bool
    return getMappingFromSFunction(n)[2]
end

function isAbundant(n::Integer)::Bool
    return getMappingFromSFunction(n)[3]
end

# Including N
function pic(n::Integer)::Integer
    numberOfPrimesTillN = 0
    for i = 2:n
        isPrime = true
        for j = 1:i
            if (i % j == 0 && i != j) || j == 1
                isPrime = false
                break
            end
        end
        if isPrime
            numberOfPrimesTillN += 1
        end
    end
    return numberOfPrimesTillN
end

# Not Including N
function pi(n::Integer)::Integer
    numberOfPrimesTillN = 0
    for i = 2:n-1
        isPrime = true
        for j = 2:i
            if (i % j == 0 && i != j)
                isPrime = false
                break
            end
        end
        if isPrime && (i != n - 1)
            numberOfPrimesTillN += 1
        end
    end
    return numberOfPrimesTillN
end

πf = pi
πfc = pic

function 𝞋(n::Integer)::Float64
    value = 0.0
    for i = 2:n
        isPrime = true
        for j = 2:i
            if (i % j == 0 && i != j)
                isPrime = false
                break
            end
        end
        if isPrime && (i != n - 1)
            value += log(ℯ, p)
        end
    end
    return value
end

chebyshev_theta = 𝞋

function ψ(n::Integer)::Float64
    value = 0.00
    for p = 2:n
        isPrime = true
        for j = 2:p
            if (p % j == 0 && p != j)
                isPrime = false
                break
            end
        end
        if isPrime
            pk = p^1
            sum = 0.0
            while pk <= n
                sum += log(ℯ, p)
                pk *= p
                println("ln $p")
            end
            value += sum
        end
    end
    return value
end

chebyshev_psi = ψ

max_stars_repo_name: EllianCarlos/number-theory-but-in-julia

max_issues_repo_name: EllianCarlos/number-theory-but-in-julia

max_forks_repo_name: EllianCarlos/number-theory-but-in-julia


text:
# Synthesizer

abstract Envelope
type ADSREnvelope <: Envelope
    A::TIME_TYPE      #time
    D::TIME_TYPE      #time
    S::AMP_TYPE       #amp
    R::TIME_TYPE      #time
    amp::AMP_TYPE     #attack
    ADSREnvelope(A,D,S,R,amp) = new(A,D,S,R,amp)
end
const ENVELOPE_BOX = ADSREnvelope(0,0,1,0,1)
const ENVELOPE_STANDARD = ADSREnvelope(3f-2,1f-1,7f-1,5f-1,9f-1)
releasetime(env::ADSREnvelope) = env.R

function apply!(w::Array, env::ADSREnvelope, hold::Real, fs::Integer)
    hold = totimetype(hold)
    t1 = env.A
    t2 = env.A + env.D
    t3 = hold
    t4 = hold + env.R
    if t1 < hold < t2
        t2 = hold
    elseif hold < env.A
        t1 = hold
        t2 = hold
    end
    applylinspace!(w, fs, nsamples(w),  0, t1, toamptype(0), env.amp)
    applylinspace!(w, fs, nsamples(w), t1, t2, env.amp, env.S)
    applylinspace!(w, fs, nsamples(w), t2, t3, env.S, env.S)
    applylinspace!(w, fs, nsamples(w), t3, t4, env.S, toamptype(0))
    return w
end
function applylinspace!(w, fs, n, t1, t2, A1, A2)
    r = seconds2saferange(fs, n, t1, t2)
    m = ifelse( length(r)==1, A2, (A2 - A1)/(length(r) - 1) )
    @inbounds for c in 1:nchannels(w)
        for i = 1:length(r)
            w[r[i],c] *= A1 + (i-1)*m
        end
    end
end

# waves are periodic functions with period 1, and : [0, inf) -> [-1,1]
const twopi = totimetype(2*pi)
const twodpi = totimetype(2/pi)
const halft = totimetype(0.5)
nozerosign(x) = ifelse(x >= 0, one(x), oftype(x,-1))

sinewave(x) = sin(twopi*x)
squarewave(x) = nozerosign(sin(twopi*x))
sawtooth(x) = 2*(x + halft - floor(x + halft)) - 1
trianglewave(x) = twodpi*asin(sin(twopi*x))
whitenoise(x) = 2*rand(AMP_TYPE) - 1
binarynoise(x) = 2*round(rand(AMP_TYPE)) - 1
chirpf(x, k::Real=1, f0::Real = 0) = x.*(f0 + k*x)
ichirpf(x) = chirpf(x, -1, 20000)
expchirpf(x, k::Real=2, f0::Real = 1) = x.*(1f-2.*k.^(x))

type WaveForm
    wave::Function              # wave function
    env::Envelope               # envelope
    nchannels::Int              # number of channels
    phase::Vector{TIME_TYPE}    # phase shift of each channel
    function WaveForm(wave ,env, nchannels, phase)
        isa(phase, Real) && ( phase = repmat([phase], nchannels) )
        @assert nchannels == length(phase)
        return new(wave, env, nchannels, phase)
    end
end
WaveForm() = WaveForm(sinewave, ENVELOPE_BOX, 1, [0])
function WaveForm(; wave::Function=sinewave,
                env::Envelope=ENVELOPE_BOX,
                nchannels::Integer=1,
                phase::Union(Real,Vector)=0)
    WaveForm(wave, env, nchannels, phase)
end

# if wave is a periodic function with period 1
# returns enveloped waveform with frequency freq
function waveform(freq::Real, time::Real, fs::Integer, form::WaveForm=WaveForm())
    @assert time >= 0
    r = (0:(fs-1)*(totimetype(time)+releasetime(form.env)))*totimetype(freq/fs)
    w = time2wave(r, form)
    apply!(w, form.env, time, fs)
    return MCAudio(w, fs)
end
waveform(tone::ASCIIString, args...) = waveform(notefreq(tone), args...)

# return an AMP_TYPE Array
function time2wave(r::Range, form::WaveForm)
    n = length(r)
    if form.nchannels == 1
        w = Array(AMP_TYPE, n)
    else
        w = Array(AMP_TYPE, n, form.nchannels)
    end
    ww = Array(AMP_TYPE, n)
    return setwave!(w, ww, r, form.wave, form.phase)
end
function setwave!(w, ww, r, wave, fp)
    n = length(r)
    for c in 1:nchannels(w)
        rp = r + fp[c]
        broadcast!(wave, ww, rp)        # ww = wave(rp)
        copy!(w, 1 + (c-1)*n, ww, 1, n) # w[:,c] = ww
    end
    return w
end

type Note
	freq::TIME_TYPE     # frequency
	start::TIME_TYPE	# starting time
	S::TIME_TYPE		# sustain time
	amp::AMP_TYPE		# amplitude multiplier
	Note(freq, start, S, amp=1) = new(freq, start, S, amp)
end
Note(tone::ASCIIString, args...) = Note(notefreq(tone), args...)

type Notes
	notes::Vector{Note}
	wf::WaveForm
	fs::Int
    Notes(notes, wf, fs) = new(notes, wf, fs)
end

function Base.scale!(n::Note, from::Real, to::Real)
    n.freq *= to/from
    return n
end
function Base.scale!(ns::Notes, from::Real, to::Real)
    for i in 1:length(ns.notes)
        scale!(ns.notes[i], from, to)
    end
    return ns
end

# convert Notes to MCAudio
function render(ns::Notes)
    maxt = findmaxtime(ns.notes)
    r = unsafe_seconds2range(ns.fs, 0, 0.0, maxt)
    fs = ns.fs
    wf = ns.wf
    s = MCAudio(zeros(AMP_TYPE, r[end], ns.wf.nchannels), fs)

    for i in 1:length(ns.notes)
        w = waveform(ns.notes[i].freq, ns.notes[i].S, fs, wf)
        scale!(ns.notes[i].amp, w.v)
        addto!(s, ns.notes[i].start, w)
    end
    return s
end
function findmaxtime(ns::Vector{Note})
    maxt = 0
    @inbounds for i in 1:length(ns)
        @assert ns[i].start >= 0
        ti = ns[i].start + ns[i].S
        ti > maxt && (maxt = ti)
    end
    return maxt
end

immutable Tuning
    value2freq::Function    # convert frequency value to Hz (Real -> TIME_TYPE)
    note2value::Function    # convert note name to frequency value (ASCIIString -> Real)
    Tuning(value2freq, note2value) = new(value2freq, note2value)
end

function midinotenumber2freq(d::Real)
    return totimetype( 2^((d-69)/12)*440 )
end
function standardnote2value(n::ASCIIString)
    v = 60
    baseoct = 4
    note, oct = noteoct(n, baseoct)
    v = v + (oct - baseoct)*12
    if note == "C"
        # v += 0
    elseif note == "C#" || note == "Db"
        v += 1
    elseif note == "D"
        v += 2
    elseif note == "D#" || note == "Eb"
        v += 3
    elseif note == "E"
        v += 4
    elseif note == "F"
        v += 5
    elseif note == "F#" || note == "Gb"
        v += 6
    elseif note == "G"
        v += 7
    elseif note == "G#" || note == "Ab"
        v += 8
    elseif note == "A"
        v += 9
    elseif note == "A#" || note == "Bb"
        v += 10
    elseif note == "B"
        v += 11
    else
        error("unknown note ", string(note))
    end
    return v
end
function noteoct(n::ASCIIString, baseoct::Integer)
    ind = 0
    for i in 1:length(n)
        if isdigit(n[i])
            ind = i
            break
        end
    end
    if ind == 0
        oct = baseoct
    elseif ind == 1
        error("no note supplied: ", string(n))
    else
        oct = int(n[ind:end])
    end
    note = n[1:ind-1]
    return note, oct
end

const TUNING_STANDARD = Tuning(midinotenumber2freq, standardnote2value)

notefreq(x::ASCIIString, t::Tuning=TUNING_STANDARD) = notefreq(t.note2value(x), t)
notefreq(x::Real, t::Tuning=TUNING_STANDARD) = t.value2freq(x)

# convert bar to seconds
function bar(bt::Tuple, bpm::Real=120, sig::NTuple{2,Integer}=(4,4))
    n = length(bt)
    if n > 1
        for i in 2:n
            @assert 1 <= bt[i] < sig[1] + 1
        end
    end
    @assert sig == (4,4)  # other not implemented
    spb = 60/totimetype(bpm)    # seconds per beat
    t = totimetype(0)
    if n >= 1
        t += (totimetype(bt[1])-1)*totimetype(sig[1])*spb
    end
    if n >= 2
        t += (totimetype(bt[2])-1)*spb
    end
    if n >= 3
        t += (totimetype(bt[3])-1)*spb/totimetype(sig[1])
    end
    n >= 4 && error("tuple to long")

    return t::TIME_TYPE
end


text:
using Pkg

const PKGS = String.(keys(Pkg.installed()))

function precompile_one(pkg::String)
    eval(Meta.parse("using $pkg"))
end

function precompile_all(useprogress::Bool=true)
    if useprogress
        progress = Progress(length(PKGS), 1, "Precompiling Packages")
    end
    @inbounds for (i, pkg) in enumerate(PKGS)
        if useprogress
            next!(progress)
        else
            if i == 1
                println("precompiling packages:")
            end
            status = "$pkg ($i/$(length(PKGS)))"
            println("\t$status")
        end
        precompile_one(pkg)
    end
end

# first update everything
Pkg.update()

# determine if a progress bar is possible
useprogress = false
if "ProgressMeter" in PKGS
    using ProgressMeter
    useprogress = true
end

# precompile all packages iteratively
precompile_all(useprogress)

max_stars_repo_name: dysonance/Preferences

max_issues_repo_name: dysonance/Preferences

max_forks_repo_name: dysonance/Preferences


text:
using Justitia
using MathOptInterface
const MOI = MathOptInterface
using QuadraticRelaxations

mutable struct BoxQPInstance <: Justitia.AbstractInstance
    model::MOI.ModelLike
    Q::Matrix{Float64}
    c::Vector{Float64}
    primal_factory::Any
    cutoff::Float64

    # ASSUMPTION: BoxQP instance is minimization
    BoxQPInstance(model::MOI.ModelLike, Q::Matrix, c::Vector, primal_factory) =
        new(model, Q, c, primal_factory, Inf)
end

function Justitia.prep_instance!(instance::BoxQPInstance)
    primal_model = MOI.instantiate(instance.primal_factory)
    MOI.copy_to(primal_model, instance.model)
    MOI.optimize!(primal_model)
    @assert MOI.get(primal_model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
    instance.cutoff = MOI.get(primal_model, MOI.ObjectiveValue())
    return nothing
end

function _create_instance(family::String, instance::String, primal_factory, epigraph::Bool)
    Q, c = parse_box_qp(joinpath(INSTANCE_DIR, family, instance))
    return BoxQPInstance(formulate_box_qp(Q, c, epigraph = epigraph), Q, c, primal_factory)
end

struct BoxQPApproach <: Justitia.AbstractApproach
    factory::Any
    transformation::Function
end

function reformulate_with(relaxation, shift, soc_lower_bound)
    return moi_model -> QuadraticRelaxations.reformulate_quadratics(
        moi_model,
        QuadraticRelaxations.Reformulation(relaxation, shift, soc_lower_bound),
    )
end

function Justitia.build_model(
    approach::BoxQPApproach,
    instance::BoxQPInstance,
    config::Dict{String},
)
    @assert collect(keys(config)) == ["factory_config"]
    src_model = approach.transformation(instance.model)
    dest_model =
        MOI.instantiate(approach.factory(config["factory_config"], instance.cutoff))
    MOI.copy_to(dest_model, src_model, copy_names = false)
    return Justitia.MOIModel(dest_model)
end

Base.@kwdef mutable struct BoxQPResult <: Justitia.AbstractResult
    termination_status::MOI.TerminationStatusCode
    primal_status::MOI.ResultStatusCode
    solve_time_sec::Float64
    x::Union{Nothing,Vector{Float64}} = nothing
    primal_bound::Float64 = NaN
    dual_bound::Float64 = NaN
    true_objective::Float64 = NaN
    cutoff::Float64
    node_count::Int = -1
end

const _BOXQP_RESULT_FIELDS = [
    :termination_status,
    :primal_status,
    :primal_bound,
    :dual_bound,
    :true_objective,
    :cutoff,
    :solve_time_sec,
    :node_count,
]

function Justitia.CSVRecord(filename::String, ::Type{BoxQPResult})
    # Don't write feasible solution to CSV
    return Justitia.CSVRecord(
        filename,
        vcat(["instance", "approach"], [string(field) for field in _BOXQP_RESULT_FIELDS]),
    )
end

function Justitia.record_result!(
    table::Justitia.CSVRecord,
    result::BoxQPResult,
    instance_name::String,
    approach_name::String,
)
    println(
        table.fp,
        join(
            vcat(
                [instance_name, approach_name],
                [getfield(result, field) for field in _BOXQP_RESULT_FIELDS],
            ),
            ",",
        ),
    )
    return flush(table.fp)
end

function Justitia.tear_down(
    model::Justitia.MOIModel,
    instance::BoxQPInstance,
    ::Type{BoxQPResult},
)
    opt = model.opt
    result = BoxQPResult(
        termination_status = MOI.get(opt, MOI.TerminationStatus()),
        primal_status = MOI.get(opt, MOI.PrimalStatus()),
        solve_time_sec = MOI.get(opt, MOI.SolveTime()),
        cutoff = instance.cutoff,
    )
    if result.primal_status == MOI.FEASIBLE_POINT
        num_vars = MOI.get(opt, MOI.NumberOfVariables())
        all_var_vals =
            [MOI.get(opt, MOI.VariablePrimal(), MOI.VariableIndex(i)) for i = 1:num_vars]
        result.x = all_var_vals
        n = length(instance.c)
        @assert n == size(instance.Q, 1) == size(instance.Q, 2)
        x_val = all_var_vals[1:n]
        result.true_objective = -0.5dot(x_val, instance.Q * x_val) - dot(instance.c, x_val)
    end
    try
        primal_bound = MOI.get(opt, MOI.ObjectiveValue())
        result.primal_bound = primal_bound
    catch ArgumentError
        @warn "Primal bound not available; proceeding anyway."
    end
    try
        dual_bound = MOI.get(opt, MOI.ObjectiveBound())
        result.dual_bound = dual_bound
    catch ArgumentError
        @warn "Dual bound not available; proceeding anyway."
    end
    try
        node_count = MOI.get(opt, MOI.NodeCount())
        result.node_count = node_count
    catch ArgumentError
        @warn "Node count not available; proceeding anyway."
    end
    return result
end


text:
module WhipRest

using HttpServer

include("beans.jl")

  export ServerConfig, ResponseConfig, WhipRestServer, RestMessage, GET, POST, PUT, PATCH, DELETE, createServer, startServer

  function createServer(config::ServerConfig)
    restServer = WhipRestServer(config.port, Dict{String,Any}(),  Dict{String,Any}(), Dict{String,Any}(), Dict{String,Any}(), Dict{String,Any}())
    #Iterate over all the request response mappings
    for i = 1:length(config.responseConfigurations)
        responseConfiguration = config.responseConfigurations[i]
        #get the correct mappings set based on HTTP method
        mappings = getMappings(restServer, responseConfiguration.method)
        #Put the response handler into map for the given path
        mappings[responseConfiguration.path] = responseConfiguration.responseHandler
    end
    return restServer
  end

  function startServer(whipServer::WhipRestServer)
    #Create a new HTTP handler
    http = HttpHandler() do req::Request, res::Response
      #Get the mappings of path to action based on HTTP method
      mappings = getMappings(whipServer, getHTTP_METHOD(req.method))
      try
        #Extract the path from the request resource
        path = getPath(req.resource)
        #Try to find the function for the mapping
        responseHandler = mappings[String(path)]
        #Call the handler for this path with the content after the path
        params = getParams(req.resource)
        Response(responseHandler(RestMessage(params, req.data)))
      catch error
        if isa(error, KeyError)
           #No mapping was found for this path
           Response("404")
        end
      end

    end
    server = Server(http)
    run(server, whipServer.port)
  end

  function getPath(resource::String)
    #Get the entire resource from the first / to the ? or the end of string
    paramIndex = search(resource, '?')
    if paramIndex==0
      #Trim the / off the end if it exists
      if resource[end] == '/'
        path = resource[1:end-1]
      else
        path = resource[1:end]
      end
    else
      path = resource[1:paramIndex-1]
    end
    return path
  end

  function getParams(resource::String)
    paramIndex = search(resource, '?')
    paramsMap = Dict{String,String}()
    if paramIndex==0
      return paramsMap
    else
      params = split(resource[paramIndex+1:end], '&')
      for i = 1:length(params)
        pair = split(params[i], '=')
        paramsMap[pair[1]] = String(pair[2])
      end

      return paramsMap
    end
  end

  function getHTTP_METHOD(method::String)
    if method == "GET"
      return GET
    elseif method== "POST"
      return POST
    elseif method == "PUT"
      return PUT
    elseif method == "PATCH"
      return PATCH
    elseif method == "DELETE"
      return DELETE
    end
  end

  function getMappings(whipServer::WhipRestServer, method::HTTP_METHOD)
    if GET == method
      return whipServer.getMappings
    elseif POST == method
      return whipServer.postMappings
    elseif PUT == method
      return whipServer.putMappings
    elseif PACTH == method
      return whipServer.patchMappings
    elseif DELETE == method
      return whipServer.deleteMappings
    end
  end

end #end module WhipRest

max_stars_repo_name: WhippyCleric/WhipRest

max_issues_repo_name: WhippyCleric/WhipRest

max_forks_repo_name: WhippyCleric/WhipRest