text:
mutable struct ParametersFP
    roundingMethod::Int
    projectionMethod::Int
    perturbationMethod::Int
    restartMethod::Int
    presolve::Int
    improve::Bool
    logLevel::Bool
    timelim::Float64
    recursiveRounding::Int
    TT::Bool
    TTmin::Int
    TTmax::Int
    feasibilityCheck::Int # 0: both check, 1: only isfeasible_xOverline, 2: only isfeasible_xTilde
    restart::Int
    checkCycle::Int

    function ParametersFP()
        param = new()
        param.roundingMethod = 0
        param.projectionMethod = 0
        param.perturbationMethod = 0
        param.restartMethod = 0
        param.presolve = 0
        param.improve = false
        param.logLevel = false
        param.timelim = 30.0
        param.recursiveRounding = 0
        param.TT = true
        param.TTmin = 10
        param.TTmax = 30
        param.feasibilityCheck = 1
        param.restart = 100
        param.checkCycle = 3
        return param
    end
end

#-------------------------------------------------------------------------------
#-----
#-------------------------------------------------------------------------------
function feasibilitypump(
        myModel::CPLEX.Model,
        param::ParametersFP
        ;
        compile::Bool = false
    )
    #=---------------------------------
     Allow to precompile this function
    ---------------------------------=#
    if compile
        return :COMPILE, [NaN], NaN
    end

    timeStart = time_ns()

    roundingmethod! = whichRounding[0]
    projectionmethod = whichProjection[0]
    perturbmethod! = whichPerturb[0]
    restartmethod! = whichRestart[0]
    try
        roundingmethod! = whichRounding[param.roundingMethod]
    catch
        @warn("Invalid parameter for roundingMethod, default used")
    end
    try
        projectionmethod = whichProjection[param.projectionMethod]
    catch
        @warn("Invalid parameter for projectionMethod, default used")
    end
    try
        perturbmethod! = whichPerturb[param.perturbationMethod]
    catch
        @warn("Invalid parameter for perturbationMethod, default used")
    end
    try
        restartmethod! = whichRestart[param.restartMethod]
    catch
        @warn("Invalid parameter for perturbationMethod, default used")
    end
    restart = param.restart
    pasRestart = 0
    checkCycle = param.checkCycle
    presolve = param.presolve
    improve = param.improve
    timeLim = param.timelim
    feasibilityCheck = param.feasibilityCheck

    getData = param.logLevel

    # For getData
    myData = Vector{Float64}([
        0.0, # Percentage of variables fixed with presolve      1
        0.0, # Time for presolve                                2
        0.0, # Time for first linear relaxation                 3
        0.0, # Average time for projection                      4
        0.0, # Number of perturbation                           5
        0.0, # Average time for perturbation                    6
        0.0, # Average time for rounding                        7
        0.0, # Average time for scatter search                  8
        0.0, # Average time for feasibility check               9
        0.0, # Number of iterations                            10
        0.0, # Time for improve                                11
        0.0, # Number of feasibility checks                    12
    ])

    LPMethod_init = CPLEX.get_param(myModel.env, "CPX_PARAM_LPMETHOD")
    CPLEX.set_param!(myModel.env, "CPX_PARAM_LPMETHOD", 1)
    nVars = CPLEX.num_var(myModel)
    varTypes_init = CPLEX.get_vartype(myModel)
    varLB_init = CPLEX.get_varLB(myModel)
    varLB = copy(varLB_init)
    varUB_init = CPLEX.get_varUB(myModel)
    varUB = copy(varUB_init)
    indices = Vector{Int}()
    for i in 1:nVars
        if varTypes_init[i] == 'B'
            if varLB[i] != varUB[i]
                push!(indices, i)
            else
                # Fixed variable
            end
        end
    end

    saveSol = Vector{Vector{Float64}}()
    for i in 1:(checkCycle-1)
        push!(saveSol, zeros(Float64, nVars))
    end
    saveVect = 0

    objectiveFunction_init = CPLEX.get_obj(myModel)
    objectiveFunction = copy(objectiveFunction_init)
    objSense_init = CPLEX.get_sense(myModel)
    if objSense_init == :Max
        objectiveFunction = -objectiveFunction
    end

    # This is only useful for recursiverounding!
    sortBinVariables = false
    sortBinVariablesOnce = false
    objectiveFunctionRR = zeros(Float64, nVars)
    if param.recursiveRounding == 1
        sortBinVariables = true
        sortBinVariablesOnce = true
    elseif param.recursiveRounding == 2
        sortBinVariables = true
        sortBinVariablesOnce = false
    elseif param.recursiveRounding == 3
        sortBinVariables = true
        sortBinVariablesOnce = true
        objectiveFunctionRR .= objectiveFunction_init
    elseif param.recursiveRounding == 4
        sortBinVariables = true
        sortBinVariablesOnce = false
        objectiveFunctionRR .= objectiveFunction_init
    elseif param.recursiveRounding == 0
        #Nothing
    else
        @warn "Parameter recursiveRounding has an unauthorized value, parameter ignored"
    end

    #=--------------------------------------------------------------------------
    Variables needed for presolve and to check feasibility
    --------------------------------------------------------------------------=#
    A = CPLEX.get_constr_matrix(myModel)
    senses = CPLEX.get_constr_senses(myModel)
    rhs = CPLEX.get_rhs(myModel)
    nConstr = length(rhs)

    if presolve != 0
        timeTmp = @timed begin
            nBin = length(indices)
            indicesTmp = falses(nBin)
            for ind in 1:indices
                i = indices[ind]
                fixed = true
                val = -1
                for j in A.rowval[nzrange(A, i)]
                    if senses[j] == 76
                        if val == -1
                            if A[j,i] < 0
                                val = 1
                            elseif A[j,i] > 0
                                val = 0
                            end
                        elseif val == 0
                            if A[j,i] < 0
                                fixed = false
                                break
                            end
                        elseif val == 1
                            if A[j,i] > 0
                                fixed = false
                                break
                            end
                        end
                    elseif senses[j] == 71
                        if val == -1
                            if A[j,i] > 0
                                val = 1
                            elseif A[j,i] < 0
                                val = 0
                            end
                        elseif val == 0
                            if A[j,i] > 0
                                fixed = false
                                break
                            end
                        elseif val == 1
                            if A[j,i] < 0
                                fixed = false
                                break
                            end
                        end
                    elseif senses[j] == 69
                        fixed = false
                        break
                    end
                end
                if fixed
                    if presolve == 1 || (val == 0 && objectiveFunction[i] >= 0) || (val == 1 && objectiveFunction[i] <= 0)
                        varLB[i] = val
                        varUB[i] = val
                        indicesTmp[ind] = true
                    end
                end
            end
            CPLEX.set_varLB!(myModel, varLB)
            CPLEX.set_varUB!(myModel, varUB)
            deleteat!(indices, indicesTmp)
            if getData
                myData[1] = (nBin-length(indices))*100/nBin
            end
        end
        if getData
            myData[2] = timeTmp[2]
        end
    end

    indicesRR = copy(indices)
    A = sparse(A') # This transpose makes feasibility check much faster

    if param.TT
        TTmin = param.TTmin
        TTmax = param.TTmax
    else
        TTmin = round(Int, length(indices)*param.TTmin/100, RoundDown)
        TTmax = round(Int, length(indices)*param.TTmax/100, RoundUp)
    end
    if TTmin > length(indices)
        TTmin = length(indices)
    end
    if TTmax > length(indices)
        TTmax = length(indices)
    end
    if TTmin > TTmax
        @warn "TTmin is greater than TTmax, parameters are switched"
        TTmin = TTmax
        TTmax = param.TTmin
    end

    problemType_init = CPLEX.get_prob_type(myModel)
    CPLEX.set_prob_type!(myModel, 0)
    has_int_init = myModel.has_int
    myModel.has_int = false

    #=--------------------------------------------------------------------------
    We initialize variables for variants:
      coef, alpha and alphaChange are parameters of objective
      feasibility pump. freqRound is useful for frequenciesrounding
    --------------------------------------------------------------------------=#
    coef = norm(objectiveFunction)
    if coef > 1e-8
        coef = sqrt(length(indices))/coef
    end
    alphaChange = 0.9
    alpha = 1.0

    freqRound = zeros(Int, nVars)

    #=--------------------------------------------------------------------------
    First linear relaxation, rounding and feasibility check
    --------------------------------------------------------------------------=#
    #CPLEX.set_obj!(myModel, zeros(Float64, nVars))
    timeLim_init = CPLEX.get_param(myModel.env, "CPX_PARAM_TILIM")
    timeTmp = @timed begin
        CPLEX.set_param!(myModel.env, "CPX_PARAM_TILIM", maximum([0, timeLim - (time_ns() - timeStart) / 1.0e9]))
        CPLEX.optimize!(myModel)

        xOverline = zeros(Float64, nVars)
        if CPLEX.get_status(myModel) == :CPX_STAT_OPTIMAL
            xOverline = CPLEX.get_solution(myModel)
        end
    end
    if getData
        myData[3] = timeTmp[2]
    end

    isFeasible = false
    if feasibilityCheck != 2
        isFeasible = isfeasible_xOverline(xOverline, indices, timeStart, timeLim)
    end
    if feasibilityCheck != 2
        if getData
            myData[9] += timeTmp[2]
            myData[12] += 1.0
        end
    end
    nIter = 0
    xTilde = zeros(Float64, nVars)

    if !isFeasible
        if sortBinVariables
            sortReducedCosts!(indicesRR, myModel)
            if sortBinVariablesOnce
                sortBinVariables = false
            end
        end

        CPLEX.set_sense!(myModel, :Min)

        timeTmp = @timed begin
            roundingmethod!(xTilde, xOverline, timeStart, timeLim, indices, myModel, indicesRR, objectiveFunctionRR)
        end
        xTildeBis = copy(xTilde)
        for i in indices
            if xTilde[i] > 0.5
                freqRound[i] += 1
            end
        end
        if getData
            myData[7] += timeTmp[2]
        end

        timeTmp = @timed if feasibilityCheck != 1
            isFeasible = isfeasible_xTilde(xTilde, A, senses, rhs, timeStart, timeLim)
        end
        if feasibilityCheck != 1
            if getData
                myData[9] += timeTmp[2]
                myData[12] += 1.0
            end
        end
    else
        xTilde .= xOverline
    end

    #=--------------------------------------------------------------------------
    The feasibility pump
    --------------------------------------------------------------------------=#
    while !isFeasible && (time_ns()-timeStart)/1.0e9 < timeLim
        nIter += 1
        pasRestart += 1
        alpha *= alphaChange
        timeTmp = @timed begin
            xOverline = projectionmethod(myModel, timeStart, timeLim, xTilde, indices, objectiveFunction, alpha, coef)
        end
        if getData
            myData[4] += timeTmp[2]
        end
        timeTmp = @timed if feasibilityCheck != 2
            isFeasible = isfeasible_xOverline(xOverline, indices, timeStart, timeLim)
        end
        if feasibilityCheck != 2
            if getData
                myData[9] += timeTmp[2]
                myData[12] += 1.0
            end
        end

        if !isFeasible
            if sortBinVariables
                sortReducedCosts!(indicesRR, myModel)
            end

            timeTmp = @timed begin
                roundingmethod!(xTildeBis, xOverline, timeStart, timeLim, indices, myModel, indicesRR, objectiveFunctionRR)
            end
            if getData
                myData[7] += timeTmp[2]
            end
            if pasRestart == restart
                pasRestart = 0
                restartmethod!(xTildeBis, xOverline, indices)
            end
            if xTildeBis == xTilde
                timeTmp = @timed begin
                    perturbmethod!(xTildeBis, xOverline, indices, TTmin, TTmax, freqRound, nIter)
                end
                if getData
                    myData[5] += 1
                    myData[6] += timeTmp[2]
                end
            else
                if xTildeBis in saveSol
                    pasRestart = 0
                    restartmethod!(xTildeBis, xOverline, indices)
                end
            end
            if checkCycle > 1
                saveSol[saveVect+1] .= xTildeBis
                saveVect = mod(saveVect+1, checkCycle-1)
            end
            xTilde = copy(xTildeBis)

            for i in indices
                if xTilde[i] > 0.5
                    freqRound[i] += 1
                end
            end
            timeTmp = @timed if feasibilityCheck != 1
                isFeasible = isfeasible_xTilde(xTilde, A, senses, rhs, timeStart, timeLim)
            end
            if feasibilityCheck != 1
                if getData
                    myData[9] += timeTmp[2]
                    myData[12] += 1.0
                end
            end
        else
            xTilde .= xOverline
        end
    end

    objVal = NaN
    status = :TIMEOUT
    if isFeasible
        status = :SOL_FOUND
        objVal = sum(objectiveFunction_init .* xTilde)
        timeTmp = @timed if improve
            CPLEX.set_sense!(myModel, objSense_init)
            CPLEX.set_obj!(myModel, objectiveFunction_init)
            varLB = CPLEX.get_varLB(myModel)
            varUB = CPLEX.get_varUB(myModel)
            for i in indices
                varLB[i] = round(xTilde[i])
                varUB[i] = varLB[i]
            end
            CPLEX.set_varLB!(myModel, varLB)
            CPLEX.set_varUB!(myModel, varUB)

            CPLEX.set_param!(myModel.env, "CPX_PARAM_TILIM", maximum([0, timeLim - (time_ns() - timeStart) / 1.0e9]))
            CPLEX.optimize!(myModel)
            if CPLEX.get_status(myModel) == :CPX_STAT_OPTIMAL
                xTilde = CPLEX.get_solution(myModel)
                objVal = CPLEX.get_objval(myModel)
            end
            if getData
                myData[11] = timeTmp[2]
            end
        end
    else
        xTilde = repeat([NaN], nVars)
    end

    if getData
        myData[4] /= nIter
        myData[6] /= myData[5]
        myData[7] /= nIter+1
        myData[8] /= nIter-1
        myData[9] /= myData[12]
        myData[10] = nIter

        println(myData)
    end

    #=--------------------------------------------------------------------------
    We changed the model to solve it, now we repair it
    --------------------------------------------------------------------------=#
    CPLEX.set_param!(myModel.env, "CPX_PARAM_LPMETHOD", LPMethod_init)
    CPLEX.set_vartype!(myModel, varTypes_init)
    CPLEX.set_vartype!(myModel, varTypes_init)
    CPLEX.set_varLB!(myModel, varLB_init)
    CPLEX.set_varUB!(myModel, varUB_init)
    CPLEX.set_obj!(myModel, objectiveFunction_init)
    CPLEX.set_prob_type!(myModel, problemType_init)
    CPLEX.set_sense!(myModel, objSense_init)
    myModel.has_int = has_int_init
    CPLEX.set_param!(myModel.env, "CPX_PARAM_TILIM", timeLim_init)

    return status, xTilde, objVal
end

feasibilitypump(myModel::CPLEX.Model; compile::Bool = false) =
    feasibilitypump(myModel, ParametersFP(), compile = compile)

function isfeasible_xTilde(
        xTilde::Vector{Float64},
        A::SparseArrays.SparseMatrixCSC{Float64,Int64},
        senses::Vector{Int8},
        rhs::Vector{Float64},
        timeStart::UInt,
        timeLim::Float64
        ;
        compile = false
    )
    if compile
        return false
    end

    nVars::Int64, nConstr::Int64 = size(A)
    isFeasible = true
    sum::Float64 = 0.0
    i = 1
    while isFeasible && i <= nConstr
    #for i in 1:nConstr
        sum = 0.0
        #for j in 1:nVars
        for j in A.rowval[nzrange(A, i)]
            sum += A[j,i] * xTilde[j]
        end
        if senses[i] == 76
            if sum > rhs[i]
                isFeasible = false
            end
        elseif senses[i] == 71
            if sum < rhs[i]
                isFeasible = false
            end
        elseif senses[i] == 69
            if !isapprox(sum, rhs[i], atol=1e-8)
                isFeasible = false
            end
        end
        i += 1
        if (time_ns()-timeStart)/1.0e9 > timeLim
            isFeasible = false
        end
    end

    return isFeasible
end

function isfeasible_xOverline(
        xOverline::Vector{Float64},
        indices::Vector{Int},
        timeStart::UInt,
        timeLim::Float64
        ;
        compile = false
    )
    if compile
        return false
    end

    nBin = length(indices)
    ind = 1
    isFeasible = true
    while ind <= nBin && isFeasible
        i = indices[ind]
        if abs(xOverline[i]-round(xOverline[i])) > 1e-8
            isFeasible = false
        end
        if (time_ns()-timeStart)/1.0e9 > timeLim
            isFeasible = false
        end
        ind += 1
    end

    return isFeasible
end

function initializeparametersfp()
    return ParametersFP()
end

function setparamfp!(paramList::ParametersFP, param::Int, value::Int)
    if param == 10
        paramList.roundingMethod = value
    elseif param == 20
        paramList.projectionMethod = value
    elseif param == 30
        paramList.perturbationMethod = value
    elseif param == 32
        paramList.TTmin = value
    elseif param == 33
        paramList.TTmax = value
    elseif param == 34
        paramList.checkCycle = value
    elseif param == 35
        paramList.restart = value
    elseif param == 2
        paramList.timelim = Float64(value)
    elseif param == 11
        paramList.recursiveRounding = value
    elseif param == 41
        paramList.presolve = value
    elseif param == 100
        paramList.feasibilityCheck = value
    elseif param == 1
        try
            paramList.logLevel = Bool(value)
        catch
            @warn("The parameter "*string(param)*" can take only 0 or 1")
        end
    elseif param == 31
        try
            paramList.TT = Bool(value)
        catch
            @warn("The parameter "*string(param)*" can take only 0 or 1")
        end
    elseif param == 42
        try
            paramList.improve = Bool(value)
        catch
            @warn("The parameter "*string(param)*" can take only 0 or 1")
        end
    else
        @warn("The parameter "*string(param)*" doesn't exist")
    end
    return paramList
end

function setparamfp!(paramList::ParametersFP, param::Int, value::Float64)
    if param == 10
        paramList.roundingMethod = round(Int, value)
    elseif param == 20
        paramList.projectionMethod = round(Int, value)
    elseif param == 30
        paramList.perturbationMethod = round(Int, value)
    elseif param == 32
        paramList.TTmin = round(Int, value)
    elseif param == 33
        paramList.TTmax = round(Int, value)
    elseif param == 34
        paramList.checkCycle = round(Int, value)
    elseif param == 35
        paramList.restart = round(Int, value)
    elseif param == 2
        paramList.timelim = value
    elseif param == 11
        paramList.recursiveRounding = round(Int, value)
    elseif param == 41
        paramList.presolve = round(Int, value)
    elseif param == 100
        paramList.feasibilityCheck = round(Int, value)
    elseif param == 1
        try
            paramList.logLevel = Bool(round(Int, value))
        catch
            @warn("The parameter "*string(param)*" can take only 0 or 1")
        end
    elseif param == 31
        try
            paramList.TT = Bool(round(Int, value))
        catch
            @warn("The parameter "*string(param)*" can take only 0 or 1")
        end
    elseif param == 42
        try
            paramList.improve = Bool(round(Int, value))
        catch
            @warn("The parameter "*string(param)*" can take only 0 or 1")
        end
    else
        @warn("The parameter "*string(param)*" doesn't exist")
    end
    return paramList
end

function setparamfp!(paramList::ParametersFP, param::String, value::Union{Float64, Int})
    if haskey(__paramToInt, param)
        setparamfp!(paramList, __paramToInt[param], value)
    else
        @warn("The parameter "*string(param)*" doesn't exist")
    end
    return paramList
end

function getparamfp(paramList::ParametersFP, param::Int)
    value = NaN
    if param == 10
        value = paramList.roundingMethod
    elseif param == 20
        value = paramList.projectionMethod
    elseif param == 30
        value = paramList.perturbationMethod
    elseif param == 31
        value = paramList.TT
    elseif param == 32
        value = paramList.TTmin
    elseif param == 33
        value = paramList.TTmax
    elseif param == 34
        value = paramList.checkCycle
    elseif param == 35
        value = paramList.restart
    elseif param == 1
        value = paramList.logLevel
    elseif param == 2
        value = paramList.timelim
    elseif param == 11
        value = paramList.recursiveRounding
    elseif param == 41
        value = paramList.presolve
    elseif param == 100
        value = paramList.feasibilityCheck
    elseif param == 42
        value = paramList.improve
    else
        @warn("The parameter "*string(param)*" doesn't exist")
    end
    return value
end

function getparamfp(paramList::ParametersFP, param::String)
    value = NaN
    if haskey(__paramToInt, param)
        value = getparamfp(paramList, __paramToInt[param])
    end
    return value
end

function getparamnamefp(param::Int)
    paramName = ""
    if haskey(__paramToString, param)
        paramName = __paramToString[param]
    else
        @warn("No parameter "*string(param))
    end
    return paramName
end

function __compilefp__()
    feasibilitypump(
        CPLEX.Model(CPLEX.Env()),
        compile = true
    )

    # --- hybridation.jl --- #

    # --- perturbation.jl --- #
    for perturbmethod! in values(whichPerturb)
        perturbmethod!(
            [NaN],
            [NaN],
            [0],
            0,
            0,
            [0],
            0,
            compile = true
        )
    end

    # --- projection.jl --- #
    for projectionmethod in values(whichProjection)
        projectionmethod(
            CPLEX.Model(CPLEX.Env()),
            UInt(0),
            NaN,
            [NaN],
            [0],
            [NaN],
            NaN,
            NaN,
            compile = true
        )
    end

    # --- rounding.jl --- #
    for roundingmethod! in values(whichRounding)
        roundingmethod!(
            [NaN],
            [NaN],
            UInt(0),
            NaN,
            [0],
            CPLEX.Model(CPLEX.Env()),
            [0],
            [NaN],
            compile = true
        )
    end

    sortReducedCosts!(
        [0],
        CPLEX.Model(CPLEX.Env()),
        compile = true
    )

    isfeasible_xTilde(
        [NaN],
        sparse([NaN NaN; NaN NaN]),
        Vector{Int8}([0]),
        [NaN],
        UInt(0),
        NaN,
        compile = true
    )
    isfeasible_xOverline(
        [NaN],
        [0],
        UInt(0),
        NaN,
        compile = true
    )

    return 0
end

max_stars_repo_name: remi-garcia/FeasibilityPump

max_issues_repo_name: remi-garcia/FeasibilityPump

max_forks_repo_name: remi-garcia/FeasibilityPump


text:
import GLPK

include("LPBenchmark.jl")

function run_glpk(fname::String)
    glpk = GLPK.Optimizer(method=GLPK.INTERIOR)
    bridged = MOIB.full_bridge_optimizer(glpk, Float64)

    load_problem!(bridged, fname)

    MOI.set(bridged, MOI.TimeLimitSec(), 10_000)
    MOI.set(bridged, MOI.RawParameter("msg_lev"), 3)

    # Solve the problem
    run(bridged)

    # Display timing info
    return nothing
end

if abspath(PROGRAM_FILE) == @__FILE__
    run_glpk(ARGS[1])
end

ext: jl

max_stars_repo_name: mtanneau/LPBenchmarks

max_issues_repo_name: mtanneau/LPBenchmarks

max_forks_repo_name: mtanneau/LPBenchmarks


text:
using DataFrames
using Gadfly

function normalize_dictionary(dictionary)
    dictionary ./ sqrt(sum(dictionary .^ 2, 1))
end

# Like findmax but works on the absolute value
function find_abs_max(itr)
    i = indmax(abs(itr))
    return (i, itr[i])
end

function matching_persuit(signal, dictionary, sparsity)
    residual = copy(signal)
    to_return = zeros(size(dictionary)[2])

    for i = 1:sparsity
        dots = dictionary' * residual
        (best_atom_index, best_dot_product) = find_abs_max(dots)
        to_return[best_atom_index] += best_dot_product
        residual = residual - (best_dot_product .* dictionary[:,best_atom_index])
    end

    return to_return
end

# Works sorta like matching_persuit, but updates the residual at each step to
# be perpendicular to the span of the dictionary items.
function orthogonal_matching_persuit(signal, dictionary, sparsity)
    residual = copy(signal)
    num_atoms = size(dictionary)[2]

    orthonormal_basis = zeros(length(signal), sparsity)
    atom_indexes = zeros(sparsity)

    for i = 1:sparsity
        dots = dictionary' * residual
        (best_atom_index, best_dot) = find_abs_max(dots)
        atom_indexes[i] = best_atom_index
        # We are going to add this new best atom to our growing subspace.
        # We want to update the residual though to just be the component of
        # residual orthogonal to the subspace. At the end of the algorithm
        # we solve the signal in terms of the chosen vectors.
        orthogonalized_atom = copy(dictionary[:,best_atom_index])
        for j = 1:(i-1)
            basis_vec = orthonormal_basis[:,j]
            orthogonalized_atom -= dot(basis_vec, orthogonalized_atom) .* basis_vec
        end
        orthogonalized_atom /= norm(orthogonalized_atom)
        orthonormal_basis[:,i] = orthogonalized_atom
        residual = residual - (dot(orthogonalized_atom, residual) .* orthogonalized_atom)
    end

    final_atoms = dictionary[:,atom_indexes]
    # Least squares solve with the chosen atoms
    final_coords = final_atoms \ signal
    # create the vector with the specified sparsity
    to_return = zeros(num_atoms)
    for (atom_index, atom_coord) in zip(atom_indexes, final_coords)
        to_return[atom_index] = atom_coord
    end
    to_return
end

function get_sparse_encoding(signals, dictionary, sparsity)
    num_atoms = size(dictionary)[2]
    num_signals = size(signals)[2]

    encodings = zeros((num_atoms, num_signals))

    for i = 1:num_signals
        encodings[:,i] = orthogonal_matching_persuit(signals[:,i], dictionary, sparsity)
    end
    encodings
end

# Given signals (each sample a column vector), this will
# iteratively refine a dictionary that can reproduce the signals
# with a given degree of sparsity (l_0 norm of the encoding vectors)
function k_svd(signals, num_iters, num_atoms, sparsity)
    (signal_dimension, num_signals) = size(signals)
    dictionary = normalize_dictionary(rand(signal_dimension, num_atoms) .- .5)
    k_svd(signals, num_iters, num_atoms, sparsity, dictionary)
end

# Work with a dictionary that aready exists (possibly to further refine it).
function k_svd(signals, num_iters, num_atoms, sparsity, dictionary)
    # We will randomly change the order that we visit atoms
    # in the ksvd step
    atom_indexes = [1:num_atoms]
    last_print = 0
    residual_norm = 0
    start = time()

    for i = 1:num_iters
        encodings = get_sparse_encoding(signals, dictionary, sparsity)
        total_residual = signals - dictionary * encodings
        residual_norm = norm(total_residual)
        if time() - last_print > 10
            last_print = time()
	    println("Total residual: $residual_norm step: $i time: $(round(last_print - start, 2))s")
        end
        # Can randomize here to not always optimize the signals in the same order
        shuffle!(atom_indexes)
        for k = atom_indexes
            signals_using_atom = find(encodings[k,:])
	    if length(signals_using_atom) == 0
	        continue
	    end
            # We find the residual as if we didn't even have that atom (so add it back in)
            restricted_residual = total_residual[:,signals_using_atom] + dictionary[:,k] * encodings[k, signals_using_atom]
            # Now we want to find a new atom that best captures that residual
            (U, s, V) = svd(restricted_residual)
            # First singular value is the one that most captures the residual. Set this as the new atom
            dictionary[:,k] = U[:,1]
            # The right hand side here is the length of the projection onto the first singualar value,
            # which is precisely the encoding.
            encodings[k, signals_using_atom] = s[1] * (V')[1,:]
            # We removed the original atom from the residual above. Now we update the total residual
            # by subtracting out the result of this new encoding from the residual
            total_residual[:,signals_using_atom] = restricted_residual - dictionary[:,k] * encodings[k, signals_using_atom]
        end
    end
    dictionary
end

# Generate some fake data from combinations of noised up sin waves
function generate_data(num_signals, signal_dimension, num_sin_waves, noise_level)
    x = linspace(0, 2*pi, signal_dimension)

    basis_waves = [sin(rand() .+ x .* (2*num_signals*rand())) for i in 1:num_sin_waves]
    basis_waves = hcat(basis_waves...)

    selector_matrix = rand(num_sin_waves, num_signals)

    return basis_waves * selector_matrix + (noise_level .* rand(signal_dimension, num_signals))
end

# Plot the signal, the sparse approximation, and the basis functions used
# in the sparse approximation
function draw_info(signal, dictionary, sparsity)
    x = linspace(0,1,length(signal))
    encoding = orthogonal_matching_persuit(signal, dictionary, sparsity)
    basis_frames = [DataFrame(x=x, y=D[:,i], label="basis_$(i)") for i in find(encoding)]
    encoded_vector = DataFrame(x=x, y=D*encoding, label="approximation")
    raw_signal = DataFrame(x=x, y=signal, label="raw_signal")
    all_frames = vcat(encoded_vector, raw_signal, basis_frames...)
    plot(all_frames, x="x", y="y", color="label", Geom.line)
end

# Generates some fake data, fits a dictionary, and then graphs
# how the sparse approximation of the signal looks compared to
# the signal.
function demo()
    sparsity = 5
    S = generate_data(300, 100, sparsity, .2);
    D = k_svd(S, 200, 80, sparsity);
    draw_info(S[:,41], D, sparsity)
end

max_stars_repo_name: justinvf-zz/algorithmshop

max_issues_repo_name: justinvf-zz/algorithmshop

max_forks_repo_name: justinvf-zz/algorithmshop


text:
@testset "distance" begin
    p1 = Position([-180, -90])
    p2 = Position([180, -90])
    @test distance(p1, p1) == 0

    p3 = Position([-75.343, 39.984])
    p4 = Position([-75.534, 39.123])

    @test distance(p3, p4, "radians") ≈ 0.015 atol=0.001
    @test distance(p3, p4, "degrees") ≈ 0.87 atol=0.01
    @test distance(p3, p4, "miles") ≈ 60.35 atol=0.01

    @test rhumb_distance(p3, p4, "miles") ≈ 60.35 atol=0.01
    @test rhumb_distance(p3, p4, "degrees") ≈ 0.87 atol=0.01
    @test rhumb_distance(p3, p4, "nauticalmiles") ≈ 52.44 atol=0.01

    @test rhumb_distance(Point([-75.343, 39.984]), Point([-75.534, 39.123]), "miles") == rhumb_distance(p3, p4, "miles")

    a =  Point([-75.833, 39.284])
    b =  Point([-75.6, 39.984])
    c =  Point([-75.221, 39.125])
    d =  Point([-75.9221, 39.27])
    e =  Point([-75.534, 39.123])
    f = Point([-75.33, 39.44])
    target = Point([-75.4, 39.4])

    @test nearestpoint(target, [a, b, c, d, e, f]).coordinates == [-75.33, 39.44]

    @test pnorm_distance(Point([2, 0]), Point([0, 0]), 2) == 2
    @test pnorm_distance(Point([1, 1]), Point([0, 0]), 1) == 2

    mid = midpoint(Point([0, 10]), Point([0, 0]))

    @test distance(Point([0, 10]), mid) ≈ distance(Point([0, 0]), mid)

    point = Point([0, 0])
    line = LineString([[1, 1],[-1, 1]])

    @test point_to_line_distance(point, line, "miles") ≈ 69.11854715938406 atol=0.05

end


text:
#!/usr/bin/env julia

#start_pkg = time()

include("distributed.jl")

include("common.jl")

#include("second-stage-fp.jl")
include("second-stage-soft-fp.jl")

using Memento
const LOGGER = Memento.getlogger(PowerModelsSecurityConstrained)

using JuMP
using PowerModels

using Ipopt
using Cbc

try
    using Gurobi
catch
    @info "Gurobi does not appear to be installed. The gurobi command line option will result in an error."
end

#@everywhere Memento.config("debug")
#@everywhere setlevel!(LOGGER, "debug")

#println("package load time: $(time() - start_pkg)")

#println("script startup time: $(time() - start_init)")

function compute_solution1(con_file::String, inl_file::String, raw_file::String, rop_file::String, time_limit::Int, scoring_method::Int, network_model::String; output_dir::String="", scenario_id::String="none", gurobi=false)
    time_start = time()
    info(LOGGER, "time remaining: $(time_limit)")

    goc_data = parse_goc_files(con_file, inl_file, raw_file, rop_file, scenario_id=scenario_id)
    network = build_pm_model(goc_data)
    network["gen_flow_cuts"] = []
    network["branch_flow_cuts"] = []

    correct_network_solution!(network)
    write_solution1(network, output_dir=output_dir)

    result = Dict(
        "termination_status" => LOCALLY_SOLVED,
        "solution" => extract_solution(network)
    )

    time_data = time() - time_start

    ###### Preparations for Solver ######

    time_solve_start = time()
    nlp_solver = with_optimizer(Ipopt.Optimizer, tol=1e-6, mu_init=1e1)
    nlp_solver_relaxed = with_optimizer(Ipopt.Optimizer, tol=1e-6, mu_init=1e1)

    if gurobi
        qp_solver = with_optimizer(Gurobi.Optimizer, OptimalityTol=1e-6, Method=2, Crossover=0)
        qp_solver_relaxed = with_optimizer(Gurobi.Optimizer, OptimalityTol=1e-6, Method=2, Crossover=0, BarConvTol=5e-3)
    else
        qp_solver = with_optimizer(Cbc.Optimizer)
        qp_solver_relaxed = qp_solver
    end
    lp_solver = qp_solver

    time_filter = 0.0

    ###### DC OPF Solve ######

    time_dc_opf_start = time()

    network_apo = deepcopy(network)

    result = run_opf_cheap(network_apo, DCPPowerModel, lp_solver)
    if !(result["termination_status"] == OPTIMAL || result["termination_status"] == LOCALLY_SOLVED || result["termination_status"] == ALMOST_LOCALLY_SOLVED)
        warn(LOGGER, "base case DC-OPF solve failed with status $(result["termination_status"]), try with relaxed convergence tolerance")

        result = run_opf_cheap(network_apo, DCPPowerModel, qp_solver_relaxed)
        if !(result["termination_status"] == OPTIMAL || result["termination_status"] == LOCALLY_SOLVED || result["termination_status"] == ALMOST_LOCALLY_SOLVED)
            warn(LOGGER, "relaxed base case DC-OPF solve failed with status $(result["termination_status"])")
            result = Dict(
                "termination_status" => LOCALLY_SOLVED,
                "solution" => extract_solution(network)
            )
        else
            warn(LOGGER, "relaxed base case DC-OPF solve status $(result["termination_status"])")
        end
    end

    update_active_power_data!(network_apo, result["solution"])
    update_active_power_data!(network, result["solution"])

    write_solution1(network_apo, output_dir=output_dir, solution_file="solution1_apo.txt")

    correct_network_solution!(network)
    write_solution1(network, output_dir=output_dir)

    for (i,bus) in network["bus"]
        bus["vm"] = (bus["vmax"] + bus["vmin"])/2.0 + 0.04
        bus["vm_start"] = bus["vm"]
        bus["va_start"] = bus["va"]

        bus["vr_start"] = bus["vm"]*cos(bus["va"])
        bus["vi_start"] = bus["vm"]*sin(bus["va"])
    end
    for (i,gen) in network["gen"]
        gen["qg_start"] = 0.0
        gen["pg_start"] = gen["pg"]
    end

    gen_cost = calc_gen_cost(network)
    info(LOGGER, "generation cost: $(gen_cost)")

    time_dc_opf = time() - time_dc_opf_start
    info(LOGGER, "dc opf phase time: $(time_dc_opf)")

    ###### AC OPF Solve ######

    time_ac_opf_start = time()

    #result = run_opf_cheap(network, ACPPowerModel, nlp_solver)
    # if !(result["termination_status"] == OPTIMAL || result["termination_status"] == LOCALLY_SOLVED || result["termination_status"] == ALMOST_LOCALLY_SOLVED || result["termination_status"] == :Suboptimal)
    #     error(LOGGER, "voltage profile solver failed")
    # end

    deactivate_rate_a!(network)
    activate_rate_a_violations!(network)

    line_flow_vio = true
    while line_flow_vio

        result = run_opf_cheap_lazy_acr(network, nlp_solver, solution_processors=[sol_data_model!])
        if !(result["termination_status"] == OPTIMAL || result["termination_status"] == LOCALLY_SOLVED)
            warn(LOGGER, "base case AC-OPF solve failed with status $(result["termination_status"]), try with relaxed convergence tolerance")
            break
            # result = run_opf_cheap_lazy_acr(network, nlp_solver_relaxed, solution_processors=[sol_data_model!])
            # if !(result["termination_status"] == OPTIMAL || result["termination_status"] == LOCALLY_SOLVED || result["termination_status"] == ALMOST_LOCALLY_SOLVED)
            #     warn(LOGGER, "relaxed base case AC-OPF solve failed with status $(result["termination_status"])")
            #     break
            # else
            #     warn(LOGGER, "relaxed base case AC-OPF solve status $(result["termination_status"])")
            # end
        end
        update_data!(network, result["solution"])

        for (i,bus) in network["bus"]
            bus["vm_start"] = bus["vm"]
            bus["va_start"] = bus["va"]

            bus["vr_start"] = bus["vm"]*cos(bus["va"])
            bus["vi_start"] = bus["vm"]*sin(bus["va"])
        end
        for (i,gen) in network["gen"]
            gen["qg_start"] = gen["qg"]
            gen["pg_start"] = gen["pg"]
        end
        line_flow_vio = activate_rate_a_violations!(network)
    end
    #activate_rate_a!(network)

    balance = compute_power_balance_deltas!(network)
    info(LOGGER, "power balance cost: $(balance)")

    #update_data!(network, result["solution"])
    correct_network_solution!(network)
    write_solution1(network, output_dir=output_dir)

    time_ac_opf = time() - time_ac_opf_start
    info(LOGGER, "ac opf phase time: $(time_ac_opf)")

    network["delta"] = 0.0
    for (i,bus) in network["bus"]
        bus["va_start"] = bus["va"]
        bus["vm_start"] = bus["vm"]
    end

    for (i,gen) in network["gen"]
        gen["pg_base"] = gen["pg"]
        gen["pg_start"] = gen["pg"]
        gen["qg_start"] = gen["qg"]
    end

    time_solve = time() - time_solve_start
    info(LOGGER, "total solve time: $(time_solve)")

    ###### Memory Check######
    scopf_comp = true
    inverse_size = compute_inverse_size(network)
    workers = Distributed.workers()

    if 2*length(workers)*inverse_size > Sys.total_memory()
        scopf_comp = false
        gb_scale = 2^30
        inverse_mem = trunc(inverse_size/gb_scale, digits=2)
        required_mem = trunc((2*length(workers)*inverse_size)/gb_scale, digits=2)
        system_mem = trunc(Sys.total_memory()/gb_scale, digits=2)
        info(LOGGER, "number of workers $(length(workers)), estimated memory per worker $(2*inverse_mem) Gb")
        info(LOGGER, "insufficient memory to run SCOPF cut generation, required $(required_mem) Gb, system $(system_mem) Gb")
    end

    if scopf_comp
        #return
        ###### Setup Workers ######

        time_worker_start = time()
        workers = Distributed.workers()

        info(LOGGER, "start warmup on $(length(workers)) workers")
        worker_futures = []
        for wid in workers
            future = remotecall(load_network_global, wid, con_file, inl_file, raw_file, rop_file, scenario_id)
            push!(worker_futures, future)
        end

        # setup for contigency solve
        gen_cont_total = length(network["gen_contingencies"])
        branch_cont_total = length(network["branch_contingencies"])
        cont_total = gen_cont_total + branch_cont_total
        cont_per_proc = cont_total/length(workers)

        cont_order = contingency_order(network)
        cont_range = []
        for p in 1:length(workers)
            cont_start = trunc(Int, ceil(1+(p-1)*cont_per_proc))
            cont_end = min(cont_total, trunc(Int,ceil(p*cont_per_proc)))
            push!(cont_range, cont_start:cont_end,)
        end

        for (i,rng) in enumerate(cont_range)
            info(LOGGER, "task $(i): $(length(rng)) / $(rng)")
        end
        #pmap(filter_network_global_contingencies, cont_range)
        output_dirs = [output_dir for i in 1:length(workers)]

        info(LOGGER, "waiting for worker warmup to complete: $(time() - time_start)")
        for future in worker_futures
            wait(future)
        end

        time_worker = time() - time_worker_start
        info(LOGGER, "total worker warmup time: $(time_worker)")

        #return
        ###### DC SCOPF Solve ######

        time_remaining = time_limit - (time() - time_start)
        info(LOGGER, "time remaining before cut enumeration: $(time_remaining)")

        time_filter = 0.0
        objective = result["objective"]
        sol = result["solution"]

        iteration = 0
        network_apo["gen_flow_cuts"] = []
        network_apo["branch_flow_cuts"] = []
        solution_file_apo = ["solution1_apo.txt" for p in 1:length(workers)]
        while true
            time_filter_start = time()
            time_start_iteration = time()
            iteration += 1
            info(LOGGER, "cut enumeration iteration: $(iteration)")

            write_active_flow_cuts(network_apo, output_dir=output_dir)
            #cuts = pmap(check_contingencies_branch_power_remote, cont_range, output_dirs, [iteration for p in 1:length(workers)], [true for p in 1:length(workers)], solution_file_apo)
            cuts = pmap(check_contingencies_branch_power_bpv_remote, cont_range, output_dirs, [iteration for p in 1:length(workers)], solution_file_apo)
            time_filter += time() - time_filter_start

            cuts_found = sum(length(c.gen_cuts)+length(c.branch_cuts) for c in cuts)
            if cuts_found <= 0
                info(LOGGER, "no violated cuts found scopf fixed-point reached")
                network["gen_flow_cuts"] = network_apo["gen_flow_cuts"]
                network["branch_flow_cuts"] = network_apo["branch_flow_cuts"]
                break
            else
                for cut in cuts
                    info(LOGGER, "found $(length(cut.gen_cuts) + length(cut.branch_cuts)) branch flow violations")
                end
            end

            for cut in cuts
                append!(network_apo["gen_flow_cuts"], cut.gen_cuts)
                append!(network_apo["branch_flow_cuts"], cut.branch_cuts)
            end

            info(LOGGER, "active cuts: gen $(length(network_apo["gen_flow_cuts"])), branch $(length(network_apo["branch_flow_cuts"]))")
            #for c in network_apo["gen_contingencies"]
            #    info(LOGGER, "  $(c.label)")
            #end
            #for c in network_apo["branch_contingencies"]
            #    info(LOGGER, "  $(c.label)")
            #end

            time_solve_start = time()
            #result = run_scopf_cuts_soft(network_apo, DCPPowerModel, qp_solver)
            result = run_scopf_cuts_soft_bpv(network_apo, DCPPowerModel, qp_solver)
            if !(result["termination_status"] == OPTIMAL || result["termination_status"] == LOCALLY_SOLVED || result["termination_status"] == ALMOST_LOCALLY_SOLVED)
                warn(LOGGER, "scopf solve failed with status $(result["termination_status"])")

                result = run_scopf_cuts_soft_bpv(network_apo, DCPPowerModel, qp_solver_relaxed)
                if !(result["termination_status"] == OPTIMAL || result["termination_status"] == LOCALLY_SOLVED || result["termination_status"] == ALMOST_LOCALLY_SOLVED)
                    warn(LOGGER, "relaxed scopf solve failed with status $(result["termination_status"])")
                    break
                else
                    warn(LOGGER, "relaxed scopf solve status $(result["termination_status"])")
                end
            end
            info(LOGGER, "objective: $(result["objective"])")
            time_solve += time() - time_solve_start
            update_active_power_data!(network_apo, result["solution"])
            update_active_power_data!(network, result["solution"])
            write_solution1(network_apo, output_dir=output_dir, solution_file="solution1_apo.txt")

            balance = compute_power_balance_deltas!(network)
            info(LOGGER, "power balance cost: $(balance)")

            gen_cost = calc_gen_cost(network)
            info(LOGGER, "generation cost: $(gen_cost)")

            time_iteration = time() - time_start_iteration
            info(LOGGER, "iteration time: $(time_iteration)")

            time_remaining = time_limit - (time() - time_start)
            info(LOGGER, "time remaining: $(time_remaining)")

            if time_remaining < time_iteration + time_ac_opf
                warn(LOGGER, "insufficent time for next iteration, time remaining $(time_remaining), iteration time $(time_iteration), ac solve time $(time_ac_opf)")
                break
            end
        end

        ###### AC OPF Recovery Solve ######

        # result = run_opf_cheap_target(network, ACPPowerModel, nlp_solver)
        # if !(result["termination_status"] == OPTIMAL || result["termination_status"] == LOCALLY_SOLVED || result["termination_status"] == ALMOST_LOCALLY_SOLVED || result["termination_status"] == :Suboptimal)
        #     error(LOGGER, "voltage profile solver failed")
        # end
        # sol = result["solution"]
        # update_data!(network, sol)

        line_flow_vio = true
        while line_flow_vio
            result = run_opf_cheap_lazy_acr(network, nlp_solver, solution_processors=[sol_data_model!])

            if !(result["termination_status"] == OPTIMAL || result["termination_status"] == LOCALLY_SOLVED)
                warn(LOGGER, "base case AC polish solve failed with status $(result["termination_status"])")
                break
                # result = run_opf_cheap_lazy_acr(network, nlp_solver_relaxed, solution_processors=[sol_data_model!])
                # if !(result["termination_status"] == OPTIMAL || result["termination_status"] == LOCALLY_SOLVED || result["termination_status"] == ALMOST_LOCALLY_SOLVED)
                #     warn(LOGGER, "relaxed base case AC-OPF solve failed with status $(result["termination_status"])")
                #     break
                # else
                #     warn(LOGGER, "relaxed base case AC-OPF solve status $(result["termination_status"])")
                # end
            end
            update_data!(network, result["solution"])

            for (i,bus) in network["bus"]
                bus["vm_start"] = bus["vm"]
                bus["va_start"] = bus["va"]

                bus["vr_start"] = bus["vm"]*cos(bus["va"])
                bus["vi_start"] = bus["vm"]*sin(bus["va"])
            end
            for (i,gen) in network["gen"]
                gen["qg_start"] = gen["qg"]
                gen["pg_start"] = gen["pg"]
            end
            line_flow_vio = activate_rate_a_violations!(network)
        end

        balance = compute_power_balance_deltas!(network)
        info(LOGGER, "power balance cost: $(balance)")

        gen_cost = calc_gen_cost(network)
        info(LOGGER, "generation cost: $(gen_cost)")

        correct_network_solution!(network)
        write_solution1(network, output_dir=output_dir)
    end

    ###### Results Summary ######

    active_contingencies = Set(cut.cont_label for cut in [network["gen_flow_cuts"]; network["branch_flow_cuts"]])
    network["gen_contingencies_active"] = [cont for cont in network["gen_contingencies"] if cont.label in active_contingencies]
    network["branch_contingencies_active"] = [cont for cont in network["branch_contingencies"] if cont.label in active_contingencies]

    total_cuts = length(network["gen_flow_cuts"]) + length(network["branch_flow_cuts"])

    write_contingencies(network, output_dir=output_dir)

    write_scopf_summary(goc_data.scenario, network, gen_cost, branch_flow_cuts=total_cuts, load_time=time_data, solve_time=time_solve, filter_time=time_filter, total_time = time() - time_start)
end


text:
# Autogenerated wrapper script for librttopo_jll for x86_64-apple-darwin
export librttopo

using GEOS_jll
JLLWrappers.@generate_wrapper_header("librttopo")
JLLWrappers.@declare_library_product(librttopo, "@rpath/librttopo.1.dylib")
function __init__()
    JLLWrappers.@generate_init_header(GEOS_jll)
    JLLWrappers.@init_library_product(
        librttopo,
        "lib/librttopo.1.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using MLSuite
using Statistics
using Pkg
using Test

cd(mktempdir())

F, N, T = 10, 100, 10
x = randn(Float32, F, N, T)
w = rand(Float32, N, T)
y = mean(x, dims = 1)

regressors = [
    GridRegressor(),
    VWRegressor(),
    ScikitRegressor(name = "ridge", alpha = 0.01),
    ScikitRegressor(name = "lasso", alpha = 0.01),
    ScikitRegressor(name = "ridgecv"),
    ScikitRegressor(name = "lassocv"),
    ScikitRegressor(name = "kernelridge"),
    ScikitRegressor(name = "knn"),
    ScikitRegressor(name = "linearsvr"),
    ScikitRegressor(name = "svr"),
    ScikitRegressor(name = "mlp"),
    ScikitRegressor(name = "lightgbm"),
    H2oRegressor(name = "glm"),
    H2oRegressor(name = "random_forest"),
    H2oRegressor(name = "gbm"),
    H2oRegressor(name = "xgboost"),
    H2oRegressor(name = "deeplearning"),
    H2oRegressor(name = "automl"),
    H2oRegressor(name = "autoglm"),
    GbmRegressor(name = "lightgbm"),
    GbmRegressor(name = "xgboost"),
    # GbmRegressor(name = "catboost"),
    LgbmRegressor(warm_start = false),
    LgbmRegressor(warm_start = true)
]

if "Elemental" ∈ keys(Pkg.installed())
    using Elemental
    push!(regressors, ElRegressor(name = "lstsq"))
    push!(regressors, ElRegressor(name = "ridge", alpha = 0.1))
end

for model in regressors
    print(model)
    MLSuite.fit!(model, x, y, w)
    ŷ = MLSuite.predict(model, x)
    r2 = r2_score(vec(ŷ), vec(y))
    @test r2 > 0.5
end

classifiers = [
    VWClassifier(),
    ScikitClassifier(name = "logistic"),
    ScikitClassifier(name = "linearsvc"),
    ScikitClassifier(name = "svc"),
    ScikitClassifier(name = "mlp"),
    H2oClassifier(name = "glm"),
    H2oClassifier(name = "psvm"),
    H2oClassifier(name = "random_forest"),
    H2oClassifier(name = "gbm"),
    H2oClassifier(name = "xgboost"),
    H2oClassifier(name = "deeplearning"),
    H2oClassifier(name = "autoglm"),
    H2oClassifier(name = "automl"),
    GbmClassifier(name = "lightgbm"),
    GbmClassifier(name = "xgboost"),
    # GbmClassifier(name = "catboost"),
    LgbmClassifier()
]

for binary in [true, false], model in classifiers
    print(model)
    !binary && !MLSuite.support_multiclass(model) && continue
    y′ = binary ? MLSuite.signone.(y) : @. ifelse(abs(y) > 0.1, sign(y) + 1.0, 1.0)
    MLSuite.fit!(model, x, y′, w)
    ŷ = MLSuite.predict(model, x)
    prob = MLSuite.predict_proba(model, x)
    acc = accuracy_score(vec(ŷ), vec(y′))
    @test acc > 0.6
end

rankers = [
    SvmRanker(),
    # TfRanker(loss = "pairwise_hinge_loss"),
    # TfRanker(loss = "pairwise_soft_zero_one_loss"),
    # TfRanker(loss = "softmax_loss"),
    # TfRanker(loss = "sigmoid_cross_entropy_loss"),
    # TfRanker(loss = "mean_squared_loss"),
    # TfRanker(loss = "list_mle_loss"),
    # TfRanker(loss = "approx_ndcg_loss"),
    RanklibRanker(ranker = 0),
    RanklibRanker(ranker = 1),
    RanklibRanker(ranker = 2),
    RanklibRanker(ranker = 3, metric2t = "NDCG@50"),
    RanklibRanker(ranker = 3, metric2t = "DCG@50"),
    RanklibRanker(ranker = 3, metric2t = "P@50"),
    RanklibRanker(ranker = 3, metric2t = "RR@50"),
    RanklibRanker(ranker = 3, metric2t = "ERR@50"),
    RanklibRanker(ranker = 4, metric2t = "NDCG@50"),
    RanklibRanker(ranker = 4, metric2t = "DCG@50"),
    RanklibRanker(ranker = 4, metric2t = "P@50"),
    RanklibRanker(ranker = 4, metric2t = "RR@50"),
    RanklibRanker(ranker = 4, metric2t = "ERR@50"),
    RanklibRanker(ranker = 6, metric2t = "NDCG@50"),
    RanklibRanker(ranker = 6, metric2t = "DCG@50"),
    RanklibRanker(ranker = 6, metric2t = "P@50"),
    RanklibRanker(ranker = 6, metric2t = "RR@50"),
    RanklibRanker(ranker = 6, metric2t = "ERR@50"),
    RanklibRanker(ranker = 7),
    RanklibRanker(ranker = 8),
    GbmRanker(name = "lightgbm", objective = "rank:ndcg"),
    GbmRanker(name = "xgboost", objective = "rank:pairwise"),
    # GbmRanker(name = "catboost", objective = "rank:pairlogit")
]

for model in rankers
    print(model)
    y′ = @. ifelse(y > 0, 1, 0)
    MLSuite.fit!(model, x, y′, w)
    ŷ = MLSuite.predict(model, x)
    ŷ .-= mean(ŷ)
    ŷ = @. ifelse(ŷ > 0, 1, 0)
    acc = accuracy_score(vec(ŷ), vec(y′))
end


text:
type Input <: Op; Input(;o...)=new(); end
Kenv.kdef(:input,Input)
ninputs(::Input)=0
canoverwrite(::Input)=false
back_reads_x(::Input)=false
back_reads_y(::Input)=false
infersize(::Input,y)=tuple(y)
back(::Input,y...;o...)=nothing


text:
#
# ComponentPath manipulation methods
#

Base.length(path::ComponentPath) = length(path.names)
Base.isempty(path::ComponentPath) = isempty(path.names)

head(path::ComponentPath) = (isempty(path) ? Symbol[] : path.names[1])
tail(path::ComponentPath) = ComponentPath(length(path) < 2 ? Symbol[] : path.names[2:end])

# The equivalent of ".." in the file system.
Base.parent(path::ComponentPath) = ComponentPath(path.names[1:end-1])

# Return a string like "/##ModelDef#367/top/Comp1"
function Base.string(path::ComponentPath)
    s = join(path.names, "/")
    return is_abspath(path) ? string("/", s) : s
end

Base.joinpath(p1::ComponentPath, p2::ComponentPath) = is_abspath(p2) ? p2 : ComponentPath(p1.names..., p2.names...)
Base.joinpath(p1::ComponentPath, other...) = joinpath(joinpath(p1, other[1]), other[2:end]...)

"""
    _fix_comp_path!(child::AbstractComponentDef, parent::AbstractCompositeComponentDef)

Set the ComponentPath of a child object to extend the path of its composite parent.
For composites, also update the component paths for all connections. For leaf components,
also update the ComponentPath for ParameterDefs and VariableDefs.
"""
function _fix_comp_path!(child::AbstractComponentDef, parent::AbstractCompositeComponentDef)
    parent_path = pathof(parent)
    child.comp_path = child_path = ComponentPath(parent_path, child.name)
    # @info "Setting path of child $(child.name) with parent $parent_path to $child_path"

    # First, fix up child's namespace objs. We later recurse down the hierarchy.
    root = get_root(parent)

    # recursively reset all comp_paths to their abspath equivalent
    if is_composite(child)
        # do same recursively
        for grandchild in compdefs(child)
            # @info "recursively fixing comp path: child: $(pathof(child)), grandchild: $(pathof(grandchild))"
            _fix_comp_path!(grandchild, child)
        end

        # Fix internal param conns
        conns = child.internal_param_conns
        for (i, conn) in enumerate(conns)
            src_path = ComponentPath(child_path, conn.src_comp_path.names[end])
            dst_path = ComponentPath(child_path, conn.dst_comp_path.names[end])

            # @info "Resetting IPC src in $child_path from $(conn.src_comp_path) to $src_path"
            # @info "Resetting IPC dst in $child_path from $(conn.dst_comp_path) to $dst_path"

            # InternalParameterConnections are immutable, but the vector holding them is not
            conns[i] = InternalParameterConnection(src_path, conn.src_var_name,
                                                   dst_path, conn.dst_par_name,
                                                   conn.ignoreunits, conn.backup;
                                                   backup_offset=conn.backup_offset)
        end

    else
        for datum in [variables(child)..., parameters(child)...]
            # @info "Resetting leaf IPC from $(datum.comp_path) to $child_path"
            datum.comp_path = child_path
        end
    end
end

"""
    fix_comp_paths!(md::AbstractModelDef)

Recursively set the ComponentPaths in a tree below a ModelDef to the absolute path equivalent.
This includes updating the component paths for all internal/external connections. For leaf components,
we also update the ComponentPath for ParameterDefs and VariableDefs.
"""
function fix_comp_paths!(md::AbstractModelDef)
    for child in compdefs(md)
        _fix_comp_path!(child, md)
    end
end

"""
    comp_path(node::AbstractCompositeComponentDef, path::AbstractString)

Convert a string describing a path from a node to a ComponentPath. The validity
of the path is not checked. If `path` starts with "/", the first element in the
returned component path is set to the root of the hierarchy containing `node`.
"""
function comp_path(node::AbstractCompositeComponentDef, path::AbstractString)
    # empty path means just select the node's path
    isempty(path) && return node.comp_path

    elts = split(path, "/")

    if elts[1] == ""    # path started with "/"
        root = get_root(node)
        elts[1] = string(nameof(root))
    end
    return ComponentPath([Symbol(elt) for elt in elts])
end

# For leaf components, we can only "find" the component itself
# when the path is empty.
function find_comp(obj::ComponentDef, path::ComponentPath)
    return isempty(path) ? obj : nothing
end

function find_comp(obj::AbstractComponentDef, name::Symbol)
    # N.B. test here since compdef doesn't check existence
    return has_comp(obj, name) ? compdef(obj, name) : nothing
end

function find_comp(obj::AbstractCompositeComponentDef, path::ComponentPath)
    # @info "find_comp($(obj.name), $path)"
    # @info "obj.parent = $(printable(obj.parent))"

    if isempty(path)
        return obj
    end

    # Convert "absolute" path from a root node to relative
    if is_abspath(path)
        path = rel_path(obj.comp_path, path)
        # @info "abspath converted to relpath is $path"

    elseif (child = find_comp(obj, head(path))) !== nothing
        # @info "path is unchanged: $path"

    elseif nameof(obj) == head(path)
        # @info "nameof(obj) == head(path); path: $(printable(path))"
        path = tail(path)
    else
        error("Cannot find path $(printable(path)) from component $(printable(obj.comp_id))")
    end

    names = path.names
    if has_comp(obj, names[1])
        return find_comp(compdef(obj, names[1]), ComponentPath(names[2:end]))
    end

    return nothing
end

function find_comp(obj::AbstractCompositeComponentDef, pathstr::AbstractString)
    path = comp_path(obj, pathstr)
    find_comp(obj, path)
end

find_comp(cr::AbstractComponentReference) = find_comp(parent(cr), pathof(cr))

@delegate find_comp(m::Model, path::ComponentPath) => md

"""
Return the relative path of `descendant` if is within the path of composite `ancestor` or
or nothing otherwise.
"""
function rel_path(ancestor_path::ComponentPath, descendant_path::ComponentPath)
    a_names = ancestor_path.names
    d_names = descendant_path.names

    if ((a_len = length(a_names)) >= (d_len = length(d_names)) || d_names[1:a_len] != a_names)
        # @info "rel_path($a_names, $d_names) returning nothing"
        return nothing
    end

    return ComponentPath(d_names[a_len+1:end])
end

rel_path(obj::AbstractComponentDef, descendant_path::ComponentPath) = rel_path(obj.comp_path, descendant_path)

"""
Return whether component `descendant` is within the composite structure of `ancestor` or
any of its descendants. If the comp_paths check out, the node is located within the
structure to ensure that the component is really where it says it is. (Trust but verify!)
"""
function is_descendant(ancestor::AbstractCompositeComponentDef, descendant::AbstractComponentDef)
    a_path = ancestor.comp_path
    d_path = descendant.comp_path
    if d_path === nothing || (relpath = rel_path(a_path, d_path)) === nothing
        return false
    end

    # @info "is_descendant calling find_comp($a_path, $relpath)"
    return find_comp(ancestor, relpath)
end

"""
    is_abspath(path::ComponentPath)

Return true if the path starts from a ModelDef, whose name is generated with
gensym("ModelDef") names look like Symbol("##ModelDef#123")
"""
function is_abspath(path::ComponentPath)
    return ! isempty(path) && match(r"^##ModelDef#\d+$", string(path.names[1])) !== nothing
end

# Returns a dictionary of the paths associated with all components, including composite components
function _get_all_paths(m::Model)
    all_paths = Dict{Symbol, ComponentPath}()
    for comp in components(m) # iterate over top level ComponentInstances
        _add_paths(m, comp, all_paths)
    end
    return all_paths
end

# a helper function to perform a preorder traversal of a given top-level component
# in model m and add that path, and all sub-component paths, to the paths array
function _add_paths(m::Model, comp::Union{CompositeComponentInstance, LeafComponentInstance}, paths::Dict{Symbol, ComponentPath})
    if isa(comp, CompositeComponentInstance)
        paths[comp.comp_name] = comp.comp_path
        for subcomp in values(comp.comps_dict)
            _add_paths(m, subcomp, paths)
        end
    else # LeafComponentInstance
        paths[comp.comp_name] = comp.comp_path
    end
    return paths
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Private functions.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

################################################################################
#                                    Macros
################################################################################

"""
    @parse_value(T, str, line_num)

Parse the string `str` using the type `T`. If it is not succeeded, then throw an
error indicating the line `line_num` with the problem.

"""
macro parse_value(T, str, line_num)
    quote
        local _T = $(esc(T))
        local _str = $(esc(str))
        local _line_num = $(esc(line_num))

        try
            parse(_T, _str)
        catch e
            throw(ErrorException("The TLE file is not valid (error in line $_line_num): $(e.msg)."))
        end
    end
end

################################################################################
#                                  Functions
################################################################################

function _get_mant_exp(n)
    if abs(n) < 1e-20
        return 0., 0
    end

    exp  = floor(Int, log10(n)) + 1
    mant = n/10.0^exp

    return mant,exp
end

function _parse_tle(io::IO, verify_checksum::Bool = true)
    # State machine to read the TLE. It has three possible states:
    #
    #   :name -> Satellite name.
    #   :l1   -> Line 1.
    #   :l2   -> Line 2.
    #
    # The transitions are:
    #
    #   :name --> :l1 --> :l2
    #     ^                |
    #     |                |
    #      ----------------
    #
    state = :name

    # Auxiliary variables that will store the TLE.
    name = ""
    sat_num = 0
    classification = Char(0)
    int_designator = ""
    epoch_year = 0
    epoch_day = 0.0
    epoch = 0.0
    dn_o2 = 0.0
    ddn_o6 = 0.0
    bstar = 0.0
    elem_set_number = 0
    checksum_l1 = 0
    i = 0.0
    Ω = 0.0
    e = 0.0
    ω = 0.0
    M = 0.0
    n = 0.0
    rev_num = 0
    checksum_l2 = 0

    # Output array with the TLEs found in the file.
    tles = Array{TLE}(undef,0)

    line_num = 0

    skip_line_read = false

    line = nothing

    while !eof(io)
        # Read the current line, strip white spaces, and skip if it is blank or
        # is a comment.
        if !skip_line_read
            line = strip(readline(io))
            line_num += 1
            (isempty(line))  && continue
            (line[1] == '#') && continue
        else
            skip_line_read = false
        end

        # Check the state of the reading.
        if state == :name
            # Check if the line seems to be the first line of the TLE. In this
            # case, maybe the user has not provided a name. Then, change the
            # state to `:l1`, and skip the line reading in the next loop.
            if (line[1:2] == "1 ") && (length(line) == 69)
                name = "UNDEFINED"
                state = :l1
                skip_line_read  = true
                continue
            end

            # Otherwise, if the line is not blank, then it must be the name of
            # the satellite.
            #
            # NOTE: The name should not be bigger than 24 characters. However,
            # we will not check this here.

            # Copy the name and change to state to wait for the 1st line.
            name  = line
            state = :l1

        # TLE Line 1
        # ==========
        elseif state == :l1
            # The next non-blank line must be the first line of the TLE.
            # Otherwise, the file is not valid.

            # The first line must start with "1 " and have 69 characters.
            if (line[1:2] != "1 ") || (length(line) != 69)
                throw(ErrorException("The TLE file is not valid (error in line $line_num): This is not a valid 1st line."))
            end

            # Verify the checksum.
            # ====================

            checksum_l1   = @parse_value(Int, line[69], line_num)
            checksum_line = compute_checksum(line[1:end-1])

            if (verify_checksum) && (checksum_l1 != checksum_line)
                throw(ErrorException("The TLE file is not valid (error in line $line_num): Expected checksum: $checksum_line, line checksum: $checksum_l1."))
            end

            # Satellite number and classification
            # ===================================
            sat_num        = @parse_value(Int, line[3:7], line_num)
            classification = Char(line[8])

            # International designator
            # ========================
            int_designator = line[10:17]

            # Epoch
            # =====
            epoch_year = @parse_value(Int,     line[19:20], line_num)
            epoch_day  = @parse_value(Float64, line[21:32], line_num)

            # Convert the TLE year and date to JD.
            # ------------------------------------

            # If the two-number year is higher than 75, then consider it in the
            # past (e.g. 50 = 2050, but 81 = 1981).
            if epoch_year > 75
                epoch = datetime2julian(DateTime(1900 + epoch_year, 1, 1, 0, 0, 0)) - 1 +
                        epoch_day
            else
                epoch = datetime2julian(DateTime(2000 + epoch_year, 1, 1, 0, 0, 0)) - 1 +
                        epoch_day
            end

            # Mean motion derivatives
            # =======================
            dn_o2  = @parse_value(Float64, line[34:43], line_num)

            aux = ( (line[45] == ' ') ? "+." : line[45:45] * "." )*line[46:50]

            ddn_o6_dec = @parse_value(Float64, aux, line_num)
            ddn_o6_exp = @parse_value(Float64, line[51:52], line_num)
            ddn_o6 = ddn_o6_dec*10^ddn_o6_exp

            # BSTAR
            # =====

            aux = ( (line[54] == ' ') ? "+." : line[54:54] * "." )*line[55:59]

            bstar_dec = @parse_value(Float64, aux, line_num)
            bstar_exp = @parse_value(Float64, line[60:61], line_num)
            bstar = bstar_dec*10^bstar_exp

            # Ephemeris type
            # ==============
            (line[63] != '0' && line[63] != ' ') &&
                warn("Warning in TLE file (line $line_num): Ephemeris type should be 0!")

            # Element number
            # ==============
            elem_set_number = @parse_value(Int, line[65:68], line_num)

            # Change the state to wait for the line 2.
            state = :l2

        # TLE Line 2
        # ==========
        elseif state == :l2
            # The next non-blank line must be the second line of the TLE.
            # Otherwise, the file is not valid.

            # The second line must start with "2 " and have 69 characters.
            if (line[1:2] != "2 ") || (length(line) != 69)
                throw(ErrorException("The TLE file is not valid (error in line $line_num): This is not a valid 2st line."))
            end

            # Verify the checksum.
            # ====================

            checksum_l2   = @parse_value(Int, line[69], line_num)
            checksum_line = compute_checksum(line[1:end-1])

            if (verify_checksum) && (checksum_l2 != checksum_line)
                throw(ErrorException("The TLE file is not valid (error in line $line_num): Expected checksum: $checksum_line, line checksum: $checksum_l2."))
            end

            # Check satellite number with the one in the first line.
            # ======================================================
            sat_num_l2 = @parse_value(Int, line[3:7], line_num)

            if sat_num_l2 != sat_num
                throw(ErrorException("The TLE file is not valid (error in line $line_num): Satellite number in line 2 is not equal to that in line 1."))
            end

            # Inclination
            # ===========
            i = @parse_value(Float64, line[9:16], line_num)

            # RAAN
            # ====
            Ω = @parse_value(Float64, line[18:25], line_num)

            # Eccentricity
            # ============
            e = @parse_value(Float64, "." * line[27:33], line_num)

            # Argument of perigee
            # ===================
            ω = @parse_value(Float64, line[35:42], line_num)

            # Mean anomaly
            # ============
            M = @parse_value(Float64, line[44:51], line_num)

            # Mean motion
            # ===========
            n = @parse_value(Float64, line[53:63], line_num)

            # Revolution number at epoch
            # ==========================
            rev_num = @parse_value(Int, line[64:68], line_num)

            # Now that we have all the information, we can store our TLE in the
            # output array.
            push!(tles,
                  TLE(name,
                      sat_num,
                      classification,
                      int_designator,
                      epoch_year,
                      epoch_day,
                      epoch,
                      dn_o2,
                      ddn_o6,
                      bstar,
                      elem_set_number,
                      checksum_l1,
                      i,
                      Ω,
                      e,
                      ω,
                      M,
                      n,
                      rev_num,
                      checksum_l2))

            # Change the state to wait for another TLE.
            state = :name
        end
    end

    # If the final state is not :name, then we have an incomplete TLE. Thus,
    # throw and exception because the file is not valid.
    if state != :name
        throw(ErrorException("The TLE file is not valid.\nThere is an incomplete TLE."))
    end

    tles
end

"""
    compute_checksum(str::AbstractString)

Compute the checksum of the line `str` modulo 10.

The algorithm is simple: add all the numbers in the line, ignoring letters,
spaces, periods, and plus signs, but assigning +1 to the minus signs. The
checksum is the remainder of the division by 10.

"""
function compute_checksum(str::AbstractString)
    checksum = 0

    for c in str
        # Check if `c` is a number.
        if isnumeric(c)
            checksum += parse(Int, c)

        # Check if `c` is a minus sign, which has value 1.
        elseif c == '-'
            checksum += 1
        end

        # Otherwise, the value of the character is 0.
    end

    # Return the checksum modulo 10.
    checksum % 10
end

"""
    _show_tle(io::IO, tle::TLE, color::Bool = true)

Show the TLE `tle` in the IO `io`.

If `color` is `true`, then the text will be printed using colors. If `color` is
omitted, then it defaults to `true`.

"""
function _show_tle(io::IO, tle::TLE, color::Bool = true)
    # Colors will be printed only for STDOUT.
    b = (color) ? _b : ""
    d = (color) ? _d : ""
    g = (color) ? _g : ""
    y = (color) ? _y : ""

    # Print the TLE information.
    print(io, "                             $(g)TLE$(d)\n")
    print(io, "$(y)    ==========================================================$(d)\n")
    print(io, "$(b)                            Name: $(d)"); @printf(io, "%s\n", tle.name)
    print(io, "$(b)                Satellite number: $(d)"); @printf(io, "%d\n", tle.sat_num)
    print(io, "$(b)        International designator: $(d)"); @printf(io, "%s\n", tle.int_designator)
    print(io, "$(b)                    Epoch (Year): $(d)"); @printf(io, "%d\n", tle.epoch_year)
    print(io, "$(b)                     Epoch (Day): $(d)"); @printf(io, "%12.8f\n", tle.epoch_day)
    print(io, "$(b)              Epoch (Julian Day): $(d)"); @printf(io, "%12.5f\n", tle.epoch)
    print(io, "$(b)              Element set number: $(d)"); @printf(io, "%d\n", tle.elem_set_number)
    print(io, "$(b)                    Eccentricity: $(d)"); @printf(io, "%12.8f deg\n", tle.e)
    print(io, "$(b)                     Inclination: $(d)"); @printf(io, "%12.8f deg\n", tle.i)
    print(io, "$(b)                            RAAN: $(d)"); @printf(io, "%12.8f deg\n", tle.Ω)
    print(io, "$(b)             Argument of perigee: $(d)"); @printf(io, "%12.8f deg\n", tle.ω)
    print(io, "$(b)                    Mean anomaly: $(d)"); @printf(io, "%12.8f deg\n", tle.M)
    print(io, "$(b)                 Mean motion (n): $(d)"); @printf(io, "%12.8f revs/day\n", tle.n)
    print(io, "$(b)               Revolution number: $(d)"); @printf(io, "%d\n", tle.rev_num)
    print(io, "\n")
    print(io, "$(b)                              B*: $(d)"); @printf(io, "%f 1/[er]\n", tle.bstar)
    print(io, "\n")
    print(io, "$(b)                        1   d\n$(d)")
    print(io, "$(b)                       ---.--- n: $(d)"); @printf(io, "%f rev/day²\n", tle.dn_o2)
    print(io, "$(b)                        2  dt\n$(d)")
    print(io, "\n")
    print(io, "$(b)                        1   d²\n$(d)")
    print(io, "$(b)                       ---.--- n: $(d)"); @printf(io, "%f rev/day³\n", tle.ddn_o6)
    print(io, "$(b)                        6  dt²\n$(d)")
    print(io, "$(y)    ==========================================================$(d)")

    nothing
end


text:
## Problem
struct Feasibility{T1<:ProximableFunction,T2<:ProximableFunction}
    S1::T1
    S2::T2
    n::Int64
end

mutable struct FeasibilitySolution <: AbstractSolution
    x::Array{Float64, 1}
    status::Symbol
end

# Model with all data
mutable struct FeasibilityModel{T1<:ProximableFunction,T2<:ProximableFunction} <: FOSModel
    S1::T1
    S2::T2
    n::Int64
    alg::FOSAlgorithm                 #This descides which algorithm to call
    data::FOSSolverData               # Solver specific data is stored here
    solve_stat::Symbol
    obj_val::Float64
    options::Dict{Symbol, Any}
    enditr::Int64
    status_generator::Function
    init_duration::UInt64           # In nano-seconds
    history::ValueHistories.MVHistory
end

function FeasibilityModel(problem::Feasibility, alg::FOSAlgorithm, kwargs...)
    t1 = time_ns()

    model = FeasibilityModel(problem.S1, problem.S2, problem.n, alg, FOSSolverDataPlaceholder(), :NotSolved,
                     0.0, Dict{Symbol,Any}(kwargs), -1,
                     x -> (@error "No status generator defined"),
                     UInt64(1), ValueHistories.MVHistory())

    data, status_generator = init_algorithm!(model.alg, model)
    model.status_generator = status_generator
    model.data = data
    t2 = time_ns()
    model.init_duration = t2-t1

    return model
end

function solve!(problem::Feasibility, alg::FOSAlgorithm; kwargs...)
    model = FeasibilityModel(problem, alg, kwargs...)
    solution = solve!(model)
end

getinitialvalue(model::FeasibilityModel) = zeros(model.n)
getinitialvalue(model::FeasibilityModel, alg, data) = zeros(model.n)

function populate_solution(model::FeasibilityModel, alg, data, x, status)
    endstatus = status.status
    if endstatus == :Continue
        endstatus = :Indeterminate
    end
    solution = FeasibilitySolution(x, endstatus)
end

"""
    `S1, S2, n, status_generator = get_sets_and_status(alg::FOSAlgorithm, model::FOSMathProgModel)`
 Get the sets S1, S2 and their size n
"""
function get_sets_and_status(alg::FOSAlgorithm, model::FeasibilityModel)
    direct = true # alg.direct # No support in this model for setting flag here
    #!alg.direct && @warn "Not possible to set direct/indirect solve in algorithm, use appropriate set from ProximalOperators"
    status_generator = (mo, checki, eps, verbose, debug) ->
        FeasibilityStatus(mo.n, 0, mo, fill(NaN, mo.n), :Continue, checki, eps, verbose, false, direct, time_ns(), mo.init_duration, debug)
    return model.S1, model.S2, model.n, status_generator
end


text:
using Test
using IncrementalInference

@testset "basic per clique stopping criteria" begin

fg = generateGraph_LineStep(1)
smtasks = Task[]
tree = solveTree!(fg, smtasks=smtasks, recordcliqs=[:x0;], limititercliqs=[(:x0=>2);])
hist = fetchCliqHistoryAll!(smtasks)

@test haskey(hist, 1)

@test hist[1] |> length == 2

#normal solve should have 11 states, update when more are added.
fg = generateGraph_LineStep(1)
smtasks = Task[]
tree = solveTree!(fg, smtasks=smtasks, recordcliqs=[:x0;]);
hist = fetchCliqHistoryAll!(smtasks)

@test haskey(hist, 1)

@test hist[1] |> length == 12

end

@testset "test endless cycle case, issue #754" begin

fg = generateGraph_LineStep(5;
                            poseEvery=1,
                            landmarkEvery=5,
                            posePriorsAt=[0,2],
                            sightDistance=4,
                            solverParams=SolverParams(algorithms=[:default, :parametric]))

getSolverParams(fg).graphinit = false
getSolverParams(fg).treeinit = true
getSolverParams(fg).limititers = 50
smtasks = Task[]
tree = solveTree!(fg; smtasks=smtasks, verbose=true, timeout=50, recordcliqs=ls(fg));

end

@testset "basic test for tree initialization functionality" begin

# small canonical factor graph, without graphinit
fg = generateGraph_CaesarRing1D(graphinit=false)
getSolverParams(fg).graphinit = false
getSolverParams(fg).treeinit = true
@show getLogPath(fg)

# # debug only
# getSolverParams(fg).drawtree = true
# getSolverParams(fg).showtree = true
# getSolverParams(fg).dbg = true

# mkpath(getLogPath(fg))
# verbosefid = open(joinLogPath(fg, "csmVerbose.log"),"w")

tree = solveTree!(fg, timeout=70) # , verbose=true, verbosefid=verbosefid)

# flush(verbosefid)
# close(verbosefid)
# open(joinLogPath(fg, "csmLogicalReconstructMax.log"),"w") do io
#   IIF.reconstructCSMHistoryLogical(getLogPath(fg), fid=io)
# end

end

@testset "basic tree initialization limittreeinit_iters" begin

# part of fg that can init
fg = generateGraph_LineStep(3; poseEvery=1)

good_vars = sortDFG(ls(fg))
# fg.solverParams.showtree = true
# fg.solverParams.drawtree = true

# part of fg that cannot init
addVariable!(fg, :s0, ContinuousScalar)
addVariable!(fg, :s1, ContinuousScalar)
addVariable!(fg, :s2, ContinuousScalar)

addFactor!(fg, [:s0;:s1], LinearRelative(Normal()))
addFactor!(fg, [:s1;:s2], LinearRelative(Normal()))

smtasks = Task[]

solveTree!(fg; smtasks=smtasks, verbose=true)

for var in good_vars
    sppe = getVariable(fg,var) |> getPPE |> IIF.getPPESuggested
    println("Testing ", var,": ", sppe)
    @test isapprox(sppe[1], parse(Int,string(var)[end]), atol=0.15)
end

end

#


text:
export coefficient
typealias Danacoefficient DanaRealParametric
type _coefficient
	function _coefficient(_::Dict{Symbol,Any})
		fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
		fields[:Brief]="General Coefficient"
		fields[:Default]=1
		fields[:Lower]=-50
		fields[:Upper]=50
		drive!(fields,_)
		new(fields)
	end
	value::Dict{Symbol,Any}
end
typealias coefficient DanaRealParametric{_coefficient}


text:
# module for drawing map
# can put an object

using Plots
pyplot()

include(joinpath(split(@__FILE__, "src")[1], "src/model/object/object.jl"))

mutable struct Map
    objects

    # init
    function Map()
        self = new()
        self.objects = []
        return self
    end
end

function add_object(self::Map, object)
    push!(self.objects, object)
end

function draw!(self::Map)
    for obj in self.objects
        draw!(obj)
    end
end

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy


text:
"""
    setCGPInfo(; nodesFile::String = "", maxTreeDepth::Integer = 3,
               genj::GeneticJuliaStructure = GenJ)

Sets the information about a Canonical Genetic Programming problem.
This function constructs a structure of type `CGPInfo` and adds it to the GP
experiment info field of the experiment info of the main structure.

# Keyword Arguments
- `nodesFile::String = ""`: file in which the information about terminal and
    non-terminal nodes is. If this is not specified, the system provides with a
    file that is read by default and contains the simples arithmetic operations and terminals.
- `maxTreeDepth::Integer = 3`: maximum tree depth permitted. If this argument is
    not specified, it will be 3 by default.
- `genj::GeneticJuliaStructure = GenJ`: the main structure. For `code user`: don't modify
    the default value unless multiple experiments are going to be run.

# Examples
```jldoctest
julia> setCGPInfo() #Arguments set by default
```

```jldoctest
julia> setCGPInfo(nodesFile="myNodesFile.txt", maxTreeDepth=4) #Setting a specific node information file and a maximum tree depth of 4
```

See also: [`setSTGPInfo`](@ref), [`setGEPInfo`](@ref), [`setGEInfo`](@ref),
[`setGPExperimentInfo`](@ref)
"""
function setCGPInfo(; nodesFile::String = "", maxTreeDepth::Integer = 3,
                    genj::GeneticJuliaStructure = GenJ)

    if nodesFile == ""
        functionSet, terminalSet = createNodes("GeneticJulia/src/utils/GeneticProgramming/Canonical/exampleNodesCGP.json")
    elseif isfile(nodesFile)
        functionSet, terminalSet = createNodes(nodesFile)
    else
        error("Node information file $nodesFile does not exist in the specified directory")
    end

    functionSet = unique(functionSet)

    if maxTreeDepth <= 0
        error("Maximum tree depth must be an integer number greater than 0")
    end

    # Variables' names obtention
    variables = getVariableName.(filter(x -> typeof(x) <: VariableNode, terminalSet))

    setGPExperimentInfo(CGPInfo(functionSet, terminalSet, maxTreeDepth, variables), genj=genj)
end # function

"""
    setSTGPInfo(; nodesFile::String = "", maxTreeDepth::Integer = 3,
                genj::GeneticJuliaStructure = GenJ)

Sets the information about a Strongly Typed Genetic Programming problem.
This function constructs a structure of type `STGPInfo` and adds it to the GP
experiment information field of the experiment information of the main structure.

# Keyword Arguments
- `nodesFile::String = ""`: file in which the information about terminal and
    non-terminal nodes is. If this is not specified, the system provides with a
    file that is read by default and contains the simples arithmetic operations and terminals.
- `maxTreeDepth::Integer = 3`: maximum tree depth permitted. If this argument is
    not specified, it will be 3 by default.
- `genj::GeneticJuliaStructure = GenJ`: the main structure. For `code user`: don't modify
    the default value unless multiple experiments are going to be run.

# Examples
```jldoctest
julia> #Arguments set by default
julia> setSTGPInfo()
```

```jldoctest
julia> #Setting a specific node information file and a maximum tree depth of 4
julia> setSTGPInfo(nodesFile="myNodesFile.txt", maxTreeDepth=4)
```

See also: [`setCGPInfo`](@ref), [`setGEPInfo`](@ref), [`setGEInfo`](@ref),
[`setGPExperimentInfo`](@ref)
"""
function setSTGPInfo(; nodesFile::String = "", maxTreeDepth::Integer = 3,
                     genj::GeneticJuliaStructure = GenJ)

    if nodesFile == ""
        functionSet, terminalSet = createNodes("ECJ/utils/GeneticProgramming/StronglyTyped/exampleNodesSTGP.json")
    elseif isfile(nodesFile)
        functionSet, terminalSet = createNodes(nodesFile)
    else
        error("Node information file $nodesFile does not exist in the specified directory")
    end

    functionSet = unique(functionSet)

    if maxTreeDepth <= 0
        error("Maximum tree depth must be an integer number greater than 0")
    end

    # Variables' names obtention
    variables = getVariableName.(filter(x -> typeof(x) <: VariableNode, terminalSet))

    setGPExperimentInfo(STGPInfo(functionSet, terminalSet, maxTreeDepth, variables), genj=genj)
end # function

"""
    setGEPInfo(; nodesFile::String = "", head::Integer = 7, genj::GeneticJuliaStructure = GenJ)

Sets the information about a Gene Expression Programming problem.
This function constructs a structure of type `GEPInfo` and adds it to the GP
experiment information field of the experiment information of the main structure.

# Keyword Arguments
- `nodesFile::String = ""`: file in which the information about terminal and
    non-terminal nodes is. If this is not specified, the system provides with a
    file that is read by default and contains the simples arithmetic and boolean
    operations and terminals.
- `head::Integer = 7`: length of the head of the genotype, which will determine
    the total length of the genotype, along with the maximum arity of the non-terminals.
    If this argument is not specified, it will be 7 by default.
- `genj::GeneticJuliaStructure = GenJ`: the main structure. For `code user`: don't modify
    the default value unless multiple experiments are going to be run.

# Examples
```jldoctest
julia> #Arguments set by default
julia> setGEPInfo()
```

```jldoctest
julia> #Setting a specific node information file and a head length of 5
julia> setSTGPInfo(nodesFile="myNodesFile.txt", head=5)
```

See also: [`setCGPInfo`](@ref), [`setSTGPInfo`](@ref), [`setGEInfo`](@ref),
[`setGPExperimentInfo`](@ref)
"""
function setGEPInfo(; nodesFile::String = "", head::Integer = 7, genj::GeneticJuliaStructure = GenJ)

    if nodesFile == ""
        functionSet, terminalSet = createNodes("ECJ/utils/GeneticProgramming/GeneExpressionProgramming/exampleNodesGEP.json")
    elseif isfile(nodesFile)
        functionSet, terminalSet = createNodes(nodesFile)
    else
        error("Node information file $nodesFile does not exist in the specified directory")
    end

    if head <= 0
        error("Head length must be an integer number greater than 0")
    end

    # Maximum arity obtention
    maxArity = 0
    for func in functionSet
        if maxArity < func._arity
            maxArity = func._arity
        end
    end

    # Tail length calculation
    tail = head * (maxArity - 1) + 1

    # Variables' names obtention
    variables = getVariableName.(filter(x -> typeof(x) <: VariableNode, terminalSet))

    setGPExperimentInfo(GEPInfo(functionSet, terminalSet, head, tail, variables), genj=genj)
end # function

"""
    setGEInfo(N::Array{String}, T::Array{String}, P::Array{String}, S::String,
              variables::Array{String}; genj::GeneticJuliaStructure = GenJ,
              maxProductions::Integer = 15, maxDepth::Integer = -1)

Sets the information about a Grammatical Evolution problem.
This function constructs a structure of type `GEInfo` and adds it to the GP
experiment information field of the experiment information of the main structure.

# Arguments
- `N::Array{String}`: array that contains the symbols of non-terminal grammatical
    expressions.
- `T::Array{String}`: array that contains the symbols of terminal grammatical
    expressions.
- `P::Array{String}`: array that contains the production rules.
- `S:String`: initial symbol from which the expressions are generated.
- `variables::Array{String}`: array that contains the name of the variables of the problem.

# Keyword Arguments
- `maxProductions::Integer = 15`: maximum number of productions permitted.
- `maxDepth::Integer = -1`: maximum tree depth permitted.
- `genj::GeneticJuliaStructure = GenJ`: the main structure. For `code user`: don't modify
    the default value unless multiple experiments are going to be run.

# Examples
```jldoctest
julia> #Declaring the Grammar

julia> N=["expr","op", "var"]
3-element Array{String,1}:
 "expr"
 "op"
 "var"

julia> T=["x", "y", "+", "-", "/", "*", "(", ")"]
8-element Array{String,1}:
 "x"
 "y"
 "+"
 "-"
 "/"
 "*"
 "("
 ")"

julia> P=[
    "<expr> -> (<expr> <op> <expr>)
             | <var>",

    "<op> -> +
           | -
           | /
           | *",

    "<var> -> x
            | y"
]
3-element Array{String,1}:
 "<expr> -> (<expr> <op> <expr>) | <var>"
 "<op> -> + | - | / | *"
 "<var> -> x | y"

julia> S="expr"
"expr"

#Arguments by default
julia> setGEInfo(N, T, P, S, ["x", "y"])
```

```jldoctest
julia> #Setting a maximum number of productions of 20
julia> setGEInfo(N, T, P, S, ["x", "y"], maxProductions=20)
```

```jldoctest
julia> #Setting a maximum tree depth of 5
julia> setGEInfo(N, T, P, S, ["x", "y"], maxProductions=0, maxDepth=5)
```

When setting a maximum number of production rules or a maximum tree depth,
take in consideration that those two cannot be set a the same time, otherwise
an error will be thrown.

```jldoctest
julia> #Setting a maximum tree depth of 5
julia> setGEInfo(N, T, P, S, ["x", "y"], maxProductions=20, maxDepth=5)
ERROR: Either maxDepth or maxProduction must be set, but not both. By default,
    maxProductions is 15, if you want to use maxDepth, set maxProductions to 0 or below.
    [...]
```

Also, either if you set `maxProductions` to a number lower than 6 or `maxDepth` to a
number lower than 3, a warning will be displayed saying that if the Grammar is
long, generated trees might be too similar with so few productions or small depth
or incorrect due to no termination.

```jldoctest
julia> #Setting a maximum tree depth of 5
julia> setGEInfo(N, T, P, S, ["x", "y"], maxProductions=4)
Warning: maxProductions have been set to 1, if your grammar is long, grammatical
    trees will be too similar or, in some cases, incorrect due to no termination
    [...]
```

See also: [`setCGPInfo`](@ref), [`setSTGPInfo`](@ref), [`setGEPInfo`](@ref),
[`setGPExperimentInfo`](@ref)
"""
function setGEInfo(N::Array{String}, T::Array{String}, P::Array{String}, S::String,
                   variables::Array{String}; genj::GeneticJuliaStructure = GenJ,
                   maxProductions::Integer = 15, maxDepth::Integer = -1)

    if maxDepth > 0 && maxProductions < 0
        error("If maxDepth is set, maxProduction must be set. By
               default, maxProductions is 15, if you want to use maxDepth, set
               maxProductions to a higher number.")
    end
    if maxDepth <= 0 && maxProductions <= 0
        error("Either maxDepth or maxProduction must be set. By
               default, maxProductions is 15, if you want to use maxDepth, set
               maxProductions to a higher number.")
    end

    if 0 < maxProductions < 6
        @warn "maxProductions have been set to $maxProductions, if your grammar
               is long, grammatical trees will be too similar or, in some cases,
               incorrect due to no termination"
    end

    if 0 < maxDepth < 3
        @warn "maxDepth have been set to $maxDepth, if your grammar is long,
               grammatical trees will be too similar or, in some cases, incorrect
               due to no termination"
    end

    # Existance of every variable as a terminal node is checked
    for i=1:length(variables)
        check = false
        for j=1:length(N)
            if variables[i] == T[j]
                check = true
                break
            end
        end
        if !check
            error("Variable $(variables[i]) isn't found in the terminal set")
        end
    end

    grammar = createGrammar(N, T, P, S)
    setGPExperimentInfo(GEInfo(grammar, maxProductions, maxDepth, variables), genj=genj)
end # function

"""
    setGPExperimentInfo(GPInfo::GPExperimentInfo; genj::GeneticJuliaStructure = GenJ)

Sets the information about a Genetic Programming problem to the given one.

# Arguments
- `GPInfo::GPExperimentInfo`: experiment information for a GP problem.

# Keyword Arguments
- `genj::GeneticJuliaStructure = GenJ`: the main structure. For `code user`: don't modify
    the default value unless multiple experiments are going to be run.

See also: [`setCGPInfo`](@ref), [`setSTGPInfo`](@ref), [`setGEPInfo`](@ref),
[`setGEInfo`](@ref)
"""
function setGPExperimentInfo(GPInfo::GPExperimentInfo; genj::GeneticJuliaStructure = GenJ)
    genj._experimentInfo._GPExperimentInfo = GPInfo
end # function

"""
    setGPExperimentInfo(; genj::GeneticJuliaStructure = GenJ)

Sets the default information about a Genetic Programming problem according to
the individual type specified.

# Keyword Arguments
- `genj::GeneticJuliaStructure = GenJ`: the main structure. For `code user`: don't modify
    the default value unless multiple experiments are going to be run.

See also: [`setCGPInfo`](@ref), [`setSTGPInfo`](@ref), [`setGEPInfo`](@ref),
[`setGEInfo`](@ref)
"""
function setGPExperimentInfo(; genj::GeneticJuliaStructure = GenJ)
    if genj._experimentInfo._individualType <: CGPGenotype
        setCGPInfo(genj=genj)
    elseif genj._experimentInfo._individualType <: STGPGenotype
        setSTGPInfo(genj=genj)
    elseif genj._experimentInfo._individualType <: GEPGenotype
        setGEPInfo(genj=genj)
    end
end # function

max_stars_repo_name: SergioGmezM/GeneticJulia

max_issues_repo_name: SergioGmezM/GeneticJulia

max_forks_repo_name: SergioGmezM/GeneticJulia


text:
load("distributions.jl")
using Distributions

function ind_max(x)
  findmax(x)[2]
end

abstract BanditAlgorithm
abstract BanditArm
abstract BanditEnvironment

load("arms/adversarial.jl")
load("arms/bernoulli.jl")
load("arms/normal.jl")

load("algorithms/epsilon_greedy/standard.jl")
load("algorithms/epsilon_greedy/annealing.jl")
load("algorithms/softmax/standard.jl")
load("algorithms/softmax/annealing.jl")
load("algorithms/ucb/ucb1.jl")
load("algorithms/exp3/exp3.jl")
load("algorithms/hedge/hedge.jl")

load("environments/2_seasonal_arms.jl")
load("environments/5_constant_arms.jl")
load("environments/5_drifting_arms.jl")
load("environments/5_near_uniform_arms.jl")
load("environments/5_step_function_arms.jl")

load("testing_framework/tests.jl")

max_stars_repo_name: rmminusrslash/BanditsBook

max_issues_repo_name: rmminusrslash/BanditsBook

max_forks_repo_name: rmminusrslash/BanditsBook