text:
function _bls_by_shape(bls::AbstractBL{SIDE},i::Interval) where SIDE
	if SIDE == 1
		stop = i.start
		start = stop + bls.depth
	else
		stop = i.stop
		start = stop - bls.depth
	end
	return bls(start=start,stop=stop)
end

################################################################################
# Plotting

@recipe function f(bnd::Boundary{1})
	bnd.shape
end


text:
# This file contains the functions and types related to the models used to select which case to run in the A-DA algorithm

abstract CaseModel

using DecisionTree
using GLM

# types for the models
type BiasedCoin <: CaseModel
  p::Vector # vector with the probabilites for selecting case 1,2,3, and 4
end

type LogisticRegression <: CaseModel
  β_for_model1or3::Vector # parameter vector for the logistic regression model for case 1 or 3
  β_for_model2or4::Vector # parameter vector for the logistic regression model  for case 1 or 3
  mean_posterior::Vector
end

type DT <: CaseModel
  decisiontree1or3::DecisionTree.Node # DT model  for case 1 or 3
  decisiontree2or4::DecisionTree.Node # DT model  for case 1 or 3
end

# functions for selecting case 1 or 3
"""
    selectcase1or3(model::BiasedCoin)

Returns 1 if case 1 is selected.
"""
function selectcase1or3(model::BiasedCoin, theta::Vector, prediction_sample_ml_star::Real, prediction_sample_ml_old::Real)

  if rand(Bernoulli(model.p[1])) == 1
    return 1
  else
    return 0
  end

end

"""
    selectcase1or3(model::LogisticRegression)

Returns 1 if case 1 is selected.
"""
function selectcase1or3(model::LogisticRegression, theta::Vector, prediction_sample_ml_star::Real, prediction_sample_ml_old::Real)

  # create new obs using theta

  theta_log_reg_mod = zeros(length(theta)+2)
  model.mean_posterior
  for i = 1:length(theta)
    theta_log_reg_mod[i] = sqrt((model.mean_posterior[i] - theta[i]).^2)
  end

  theta_log_reg_mod[length(theta)+1] = sqrt(sum((mean_posterior - theta).^2))
  theta_log_reg_mod[length(theta)+2] = prediction_sample_ml_star/prediction_sample_ml_old

  prob_case_1 = _predict(model.β_for_model1or3, [1;theta_log_reg_mod])
  if rand(Bernoulli(prob_case_1)) == 1
    return 1
  else
    return 0
  end
end

"""
    selectcase1or3(model::DT)

Returns 1 if case 1 is selected.
"""
function selectcase1or3(model::DT, theta::Vector, prediction_sample_ml_star::Real, prediction_sample_ml_old::Real)

  theta_tree_model = zeros(length(theta)+1)
  theta_tree_model = [theta; prediction_sample_ml_star/prediction_sample_ml_old]

  if apply_tree(model.decisiontree1or3, theta_tree_model) == "case 1"
    return 1
  else
    return 0
  end

end

# functions for selecting case 2 or 4
"""
    selectcase2or4(model::BiasedCoin)

Returns 1 if case 2 is selected.
"""
function selectcase2or4(model::BiasedCoin, theta::Vector, prediction_sample_ml_star::Real, prediction_sample_ml_old::Real)

  if rand(Bernoulli(model.p[2])) == 1
    return 1
  else
    return 0
  end

end

"""
    selectcase2or4(model::LogisticRegression)

Returns 1 if case 2 is selected.
"""
function selectcase2or4(model::LogisticRegression, theta::Vector, prediction_sample_ml_star::Real, prediction_sample_ml_old::Real)

  # create new obs using theta

  theta_log_reg_mod = zeros(length(theta)+2)
  model.mean_posterior
  for i = 1:length(theta)
    theta_log_reg_mod[i] = sqrt((model.mean_posterior[i] - theta[i]).^2)
  end

  theta_log_reg_mod[length(theta)+1] = sqrt(sum((mean_posterior - theta).^2))
  theta_log_reg_mod[length(theta)+2] = prediction_sample_ml_star/prediction_sample_ml_old

  prob_case_2 = _predict(model.β_for_model2or4, [1;theta_log_reg_mod])
  if rand(Bernoulli(prob_case_2)) == 1
    return 1
  else
    return 0
  end
end

"""
    selectcase2or4(model::DT)

Returns 1 if case 2 is selected.
"""
function selectcase2or4(model::DT, theta::Vector,prediction_sample_ml_star::Real, prediction_sample_ml_old::Real)

  theta_tree_model = zeros(length(theta)+1)
  theta_tree_model = [theta; prediction_sample_ml_star/prediction_sample_ml_old]

  if apply_tree(model.decisiontree2or4, theta_tree_model) == "case 2"
    return 1
  else
    return 0
  end

end

# help functions

doc"""
Prediction for the logistic regression model at theta.
"""
function _predict(beta::Vector, theta::Vector)

  p_hat =  exp(theta'*beta)./(1+exp(theta'*beta))
  return p_hat[1]

end

ext: jl

max_stars_repo_path: select case/selectcase.jl

max_stars_repo_name: SamuelWiqvist/adamcmcpaper

max_issues_repo_path: select case/selectcase.jl

max_issues_repo_name: SamuelWiqvist/adamcmcpaper

max_forks_repo_path: select case/selectcase.jl

max_forks_repo_name: SamuelWiqvist/adamcmcpaper


text:
using OrthogonalPolynomialsQuasi, ContinuumArrays, DomainSets, Test
import OrthogonalPolynomialsQuasi: Hilbert, StieltjesPoint

@testset "Stieltjes" begin
    T = Chebyshev()
    wT = ChebyshevWeight() .* T
    x = axes(wT,1)
    z = 0.1+0.2im
    S = inv.(z .- x')
    @test S isa StieltjesPoint{ComplexF64,Float64,ChebyshevInterval{Float64}}
    f = wT * [[1,2,3]; zeros(∞)];
    J = T \ (x .* T)
    @test π*((z*I-J) \ f.args[2])[1,1] ≈ (S*f)[1]

    x = Inclusion(0..1)
    y = 2x .- 1
    wT2 = wT[y,:]
    S = inv.(z .- x')
    f = wT2 * [[1,2,3]; zeros(∞)];

    @test (π/2*(((z-1/2)*I - J/2) \ f.args[2]))[1] ≈ (S*f)[1]
end

@testset "Hilbert" begin
    wT = ChebyshevWeight() .* Chebyshev()
    wU = UltrasphericalWeight(1) .*  Ultraspherical(1)
    x = axes(wT,1)
    H = inv.(x .- x')
    @test H isa Hilbert{Float64,ChebyshevInterval{Float64}}

    @test (Ultraspherical(1) \ (H*wT))[1:10,1:10] == diagm(1 => fill(-π,9))
    @test (Chebyshev() \ (H*wU))[1:10,1:10] == diagm(-1 => fill(1.0π,9))

    # check consistency
    @test (Ultraspherical(1) \ (H*wT) * (wT \ wU))[1:10,1:10] ==
                ((Ultraspherical(1) \ Chebyshev()) * (Chebyshev() \ (H*wU)))[1:10,1:10]

    # Other axes
    x = Inclusion(0..1)
    y = 2x .- 1
    H = inv.(x .- x')

    wT2 = wT[y,:]
    wU2 = wU[y,:]
    @test (Ultraspherical(1)[y,:]\(H*wT2))[1:10,1:10] == diagm(1 => fill(-π,9))
    @test (Chebyshev()[y,:]\(H*wU2))[1:10,1:10] == diagm(-1 => fill(1.0π,9))
end

@testset "Log kernel" begin
    T = Chebyshev()
    wT = ChebyshevWeight() .* Chebyshev()
    x = axes(wT,1)
    L = log.(abs.(x .- x'))
    D = T \ (L * wT)
    @test (L * (wT * (T \ exp.(x))))[0.] ≈ -2.3347795490945797  # Mathematica

    x = Inclusion(-1..1)
    T = Chebyshev()[1x, :]
    L = log.(abs.(x .- x'))
    wT = (ChebyshevWeight() .* Chebyshev())[1x, :]
    @test (T \ (L*wT))[1:10,1:10] ≈ D[1:10,1:10]

    x = Inclusion(0..1)
    T = Chebyshev()[2x.-1, :]
    wT = (ChebyshevWeight() .* Chebyshev())[2x .- 1, :]
    L = log.(abs.(x .- x'))
    u =  wT * (2 *(T \ exp.(x)))
    @test u[0.1] ≈ exp(0.1)/sqrt(0.1-0.1^2)
    @test (L * u)[0.5] ≈ -7.471469928754152 # Mathematica
end


text:
using FIX
using DataStructures

import FIX: onFIXMessage

export logout
export onFIXMessage

function fixconnect(user::GDAXUser, handler::H) where {H <: AbstractGDAXMessageHandler}
    client = FIXClient(MbedTLS.connect("127.0.01", 4198),
                        handler,
                        Dict(8=>"FIX.4.2",
                             49=>user.key,
                             56=>"Coinbase"),
                             RateLimit(50, Dates.Second(1)))
    start(client)
    send_message(client, login_message(user.key, user.secret, user.passphrase))
    return client
end

function logout(client::FIXClient{TCPSocket, H}) where {H <: AbstractGDAXMessageHandler}
    send_message(client, logout_message())
end

function login_message(api_key::String, api_secret::String, passphrase::String)::Dict{Int64, String}
    t = Dates.format(now(Dates.UTC), "yyyy-mm-ddTHH:MM:SS.sss")
    seq_num = "1"
    message   = join([t , "A" , seq_num, api_key, "Coinbase", passphrase], Char(1))
    hmac_key  = base64decode(api_secret)
    h = HMACState("sha256", hmac_key)
    Nettle.update!(h, message)
    signature = base64encode(Nettle.digest!(h))

    return Dict{Int64, String}(35  => "A",
                             98  => "0",
                            108  => "30",
                            52 => t,
                            554  => passphrase,
                            96   => signature,
                            8013 => "Y" )
end

function logout_message()
    return Dict{Int64, String}(35=>"5")
end

function placeOrder_fake(this::FIXClient{TCPSocket, H},
                        side::String,
                        instrument::String,
                        lots::Float64,
                        price::Float64) where {H <: AbstractGDAXMessageHandler}
    id = string(Base.Random.uuid4());
    ord = Dict{Int64, String}(35=>"D",
                        21 => "1",
                        11 => id,
                        55 => instrument,
                        54 => side,
                        44 => string(signif(price, 6)),
                        38 => string(signif(lots, 8)),
                        40 => "2",
                        59 => "1",
                        7928 => "B")
    FIX.send_message_fake(this, ord)
end

function placeOrder(this::FIXClient{TCPSocket, H},
                        side::String,
                        instrument::String,
                        lots::Float64,
                        price::Float64) where {H <: AbstractGDAXMessageHandler}
    id = string(Base.Random.uuid4())
    ord = Dict{Int64, String}(35=>"D",
                        21 => "1",
                        11 => id,
                        55 => instrument,
                        54 => side,
                        44 => string(signif(price, 6)),
                        38 => string(signif(lots, 8)),
                        40 => "2",
                        59 => "1",
                        7928 => "B")
    send_message(this, ord)
end

function cancelAll(this::FIXClient{TCPSocket, H}) where {H <: AbstractGDAXMessageHandler}
    for order in FIX.getOpenOrders(this)
        c = Dict{Int64, String}()
        c[11] = string(Base.Random.uuid4())
        c[37] = order[37]
        c[41] = order[11]
        c[55] = order[55]
        c[35] = "F"
        send_message(this, c)
    end
end

mutable struct FIXOrderChain
    cl::Pair{String, Dict{Int64, String}}
    cl_ex::Pair{String, Dict{Int64, String}}
    rm::String
    rm_ex::Bool
    rejected::Bool
    function FIXOrderChain()
        return new(Pair("", Dict{Int64, String}()),
                    Pair("", Dict{Int64, String}()),
                    "",
                    false,
                    false)
    end
end

mutable struct FIXQuote <: AbstractGDAXMessageHandler
    side::String
    instrument::String
    lots::Float64
    price::Float64
    lots_tolerance::Float64
    price_tolerance::Float64
    chains::OrderedDict{String, FIXOrderChain}
    storage::Vector{FIXOrderChain}
    client::FIXClient
    m_lock::ReentrantLock
    function FIXQuote(side::String,
                    instrument::String,
                    lots::Float64,
                    price::Float64;
                    lots_tolerance::Float64 = 1e-4+1e-6,
                    price_tolerance::Float64 = 1.0)
        this = new()
        this.side = side
        this.instrument = instrument
        this.lots = lots
        this.price = price
        this.lots_tolerance = lots_tolerance
        this.price_tolerance = price_tolerance
        this.chains = OrderedDict{String, FIXOrderChain}()
        this.storage = FIXOrderChain[]
        this.m_lock = ReentrantLock()
        return this
    end
end

function onWSMessage(this::FIXQuote, msg::Dict{String, Any})
    maintain(this)
end

function onFIXMessage(this::FIXQuote, msg::OrderedDict{Int64, String})
    msg_type = msg[35]
    if msg_type == "8"
        exec_type = msg[150]
        for chain in values(this.chains)
            if exec_type == "0" #new order single
                if msg[11] == chain.cl.first
                    chain.cl_ex = Pair(msg[37], msg)
                end
            elseif exec_type == "1" #fill
                if msg[37] == chain.cl_ex.first
                    filled_size = parse(Float64, msg[32])
                    got_size = parse(Float64, chain.cl_ex.second[38])
                    chain.cl_ex.second[38] = string(signif(got_size - filled_size, 8))
                end
            elseif exec_type == "3" #reject
                if msg[37] == chain.cl_ex.first #newOrderSingle reject
                    chain.rejected = true
                end
            elseif exec_type == "4" #cancel
                if msg[37] == chain.cl_ex.first
                    chain.rm_ex = true
                end
            elseif exec_type == "9" #order cancel reject
                cancel_request_client_id = msg[11]
                for (id, chain) in this.chains
                    if chain.rm == cancel_request_client_id
                        chain.rm = ""
                        chain.rm_ex = false
                    end
                end
            end
        end
    end
    maintain(this)
end

function maintain(this::FIXQuote)
    lock(this.m_lock)
    for (id, chain) in this.chains
        if chain.rm_ex || chain.rejected
            push!(this.storage, chain)
            delete!(this.chains, id)
        end
    end

    if isempty(this.chains)
        this.lots > 0.0 && !isnan(this.price) && addNewChain(this)
    else
        for (id, chain) in this.chains
            if isempty(chain.cl_ex.first)
                continue
            end
            if !isempty(chain.rm)
                continue
            end
            cur_price = parse(Float64, chain.cl_ex.second[44])
            cur_lots  = parse(Float64, chain.cl_ex.second[38])

            price_ok = abs(cur_price - this.price) < this.price_tolerance
            lots_ok = abs(cur_lots - this.lots) < this.lots_tolerance

            if !price_ok || !lots_ok
                closeChain(this, chain)
                if this.lots > 0.0 && !isnan(this.price)
                    addNewChain(this)
                    break #maintain(this)
                end
            end
        end
    end
    unlock(this.m_lock)
end

function getPrice(this::FIXOrderChain)
    if isempty(this.cl_ex.first)
        return NaN
    end
    return parse(Float64, this.cl_ex.second[44])
end

function getLots(this::FIXOrderChain)
    if isempty(this.cl_ex.first)
        return NaN
    end
    return parse(Float64, this.cl_ex.second[38])
end

function closeChain(this::FIXQuote, chain::FIXOrderChain)
    client = this.client
    c = Dict{Int64, String}()
    c[11] = string(Base.Random.uuid4())
    c[37] = chain.cl_ex.first
    c[41] = chain.cl.first
    c[55] = this.instrument
    c[35] = "F"
    (msg, msg_str) = send_message(client, c)
    chain.rm = msg[11]
    return nothing
end

function addNewChain(this::FIXQuote)
    numleft = FIX.numMsgsLeftToSend(this.client)
    if numleft < 10
        @printf("[%ls][FIXQuote][RateLimit]: blocked. %2.2f msg per second\n", now(), this.client.m_messages.outgoing.ratelimit())
        return nothing
    end
    (msg, msg_str) = placeOrder(this.client, this.side, this.instrument, this.lots, this.price)
    id = msg[11]
    this.chains[id] = FIXOrderChain()
    this.chains[id].cl = Pair(id, msg)
    return nothing
end

max_stars_repo_name: DatName/gdax-julia

max_issues_repo_name: DatName/gdax-julia

max_forks_repo_name: DatName/gdax-julia


text:
f_1p(x) = 5.0 + (x[1]-3.0)^2 # [100., missing, missing, missing]

f_2p_1im(x) = 5.0 + (x[1]-3.0)^2 + 0.0*x[2]  # [3., missing, missing, missing]

f_2p(x) = 5.0 + (x[1]-3.0)^2 + (x[2]-4.0)^2  # [3., 4., missing, missing]

f_3p_1im(x) = 5.0 + (x[1]-3.0)^2 + (x[2]/x[3]-4.0)^2 # [3., missing, missing, missing]

f_3p_1im_dep(x) = 5.0 + (x[1]-3.0)^2 + (x[1]-x[2]-1.0)^2 + 0*x[3]^2 # [3., 2., missing]

f_4p_2im(x) = 5.0 + (x[1]-3.0)^2 + (x[2]-4.0)^2 + 0.0*x[3] + 0.0*x[4] # [3., 4., missing, missing]

f_4p_3im(x) = 5.0 + (x[1]-3.0)^2 + (x[2]/x[3]-4.0)^2 + 0.0*x[4] # [3., missing, missing, missing]

f_1p_ex(x) = 5.0 + (x[1]-1e-8)^2 # [1e-8, missing, missing, missing]

f_5p_3im(x) = 5.0 + (x[1]-3.0)^2 + (exp(x[2])-1.0)^2 + (x[3]/x[4]-4.0)^2 + 0.0*x[5] # [3., 0., missing, missing, missing]

f_3p_im(x) = 5.0 + (x[1]-3.0)^2 + (exp(x[2])-1.0)^2 + 0.0*x[3] # [3.0, 0., missing]

# test each alg on all functions
function test_alg(
  alg::NamedTuple,
  func_dict::AbstractDict = test_funcs;
  bounds::Tuple{Float64,Float64} = (-Inf,Inf),
  tol::Float64 = 1e-2,
  kwargs...
)
  @testset "get_interval() for $(alg.algorithm)" begin
    for (f_name, f) in func_dict
      #println("Testing $f_name")
      @testset "Case $f_name" begin
        for i in eachindex(f.x0)
          ep = get_interval(
            f.x0,
            i,
            f.func,
            :CICO_ONE_PASS;
            theta_bounds=fill(bounds,length(f.x0)),
            scan_tol=1e-8,
            local_alg = alg.algorithm,
            loss_crit = f.loss_crit,
            silent = true,
            kwargs...
          )
          should_skip = f_name in alg.skip
          @test ep.result[1].status == f.status[i][1] skip = should_skip
          @test ep.result[2].status == f.status[i][2] skip = should_skip
          if isa(f.endpoints[i][1], Nothing)
            @test isa(ep.result[1].value, Nothing) skip = should_skip
          else
            @test isapprox(ep.result[1].value, f.endpoints[i][1], atol=tol) skip = should_skip
          end
          if isa(f.endpoints[i][2], Nothing)
            @test isa(ep.result[2].value, Nothing) skip = should_skip
          else
            @test isapprox(ep.result[2].value, f.endpoints[i][2], atol=tol) skip = should_skip
          end
        end
      end
    end
  end
end

# functions dict
test_funcs = Dict(
  :f_1p => (
    func = f_1p,
    x0 = [3.,1.],
    endpoints = [(1.,5.),(nothing,nothing)],
    status = [(:BORDER_FOUND_BY_SCAN_TOL,:BORDER_FOUND_BY_SCAN_TOL),(:SCAN_BOUND_REACHED,:SCAN_BOUND_REACHED)],
    loss_crit = 9.
  ),

  :f_2p_1im => (
    func = f_2p_1im,
    x0 = [3.,1.],
    endpoints = [(1.,5.),
                 (nothing,nothing)],
    status = [(:BORDER_FOUND_BY_SCAN_TOL,:BORDER_FOUND_BY_SCAN_TOL),
              (:SCAN_BOUND_REACHED,:SCAN_BOUND_REACHED)],
    loss_crit = 9.
  ),

  :f_2p => (
    func = f_2p,
    x0 = [3.,4.],
    endpoints = [(1.,5.),
                 (2.,6.)],
    status = [(:BORDER_FOUND_BY_SCAN_TOL,:BORDER_FOUND_BY_SCAN_TOL),
              (:BORDER_FOUND_BY_SCAN_TOL,:BORDER_FOUND_BY_SCAN_TOL)],
    loss_crit = 9.
  ),

  :f_3p_1im => (
    func = f_3p_1im,
    x0 = [3.,4.,1.1],
    endpoints = [(1.,5.),
                 (nothing,nothing),
                 (nothing,nothing)],
    status = [(:BORDER_FOUND_BY_SCAN_TOL,:BORDER_FOUND_BY_SCAN_TOL),
              (:SCAN_BOUND_REACHED,:SCAN_BOUND_REACHED),
              (:SCAN_BOUND_REACHED,:SCAN_BOUND_REACHED)],
    loss_crit = 9.
  ),

  :f_3p_1im_dep => (
    func = f_3p_1im_dep,
    x0 = [3., 2., 2.1],
    endpoints = [(1.,5.),
                 (2.0-2.0*sqrt(2.),2.0+2.0*sqrt(2.)),
                 (nothing,nothing)],
    status = [(:BORDER_FOUND_BY_SCAN_TOL,:BORDER_FOUND_BY_SCAN_TOL),
              (:BORDER_FOUND_BY_SCAN_TOL,:BORDER_FOUND_BY_SCAN_TOL),
              (:SCAN_BOUND_REACHED,:SCAN_BOUND_REACHED)],
    loss_crit = 9.
  ),

  :f_4p_2im => (
    func = f_4p_2im,
    x0 = [3.,4.,1.,1.],
    endpoints = [(1.,5.),
                 (2.,6.),
                 (nothing,nothing),
                 (nothing,nothing)],
    status = [(:BORDER_FOUND_BY_SCAN_TOL,:BORDER_FOUND_BY_SCAN_TOL),
              (:BORDER_FOUND_BY_SCAN_TOL,:BORDER_FOUND_BY_SCAN_TOL),
              (:SCAN_BOUND_REACHED,:SCAN_BOUND_REACHED),
              (:SCAN_BOUND_REACHED,:SCAN_BOUND_REACHED)],
    loss_crit = 9.
  ),

  :f_4p_3im => (
    func = f_4p_3im,
    x0 = [3.,4.,1.1,1.1],
    endpoints = [(1.,5.),
                 (nothing,nothing),
                 (nothing,nothing),
                 (nothing,nothing)],
    status = [(:BORDER_FOUND_BY_SCAN_TOL,:BORDER_FOUND_BY_SCAN_TOL),
              (:SCAN_BOUND_REACHED,:SCAN_BOUND_REACHED),
              (:SCAN_BOUND_REACHED,:SCAN_BOUND_REACHED),
              (:SCAN_BOUND_REACHED,:SCAN_BOUND_REACHED)],
    loss_crit = 9.
  ),

  :f_1p_ex => (
    func = f_1p_ex,
    x0 = [1.5, 2.],
    endpoints = [(-2+1e-8,2+1e-8), (nothing, nothing)],
    status = [(:BORDER_FOUND_BY_SCAN_TOL,:BORDER_FOUND_BY_SCAN_TOL),(:SCAN_BOUND_REACHED,:SCAN_BOUND_REACHED)],
    loss_crit = 9.
  ),

  :f_5p_3im => (
    func = f_5p_3im,
    x0 = [3., 0.5, 8., 2., 2.],
    endpoints = [(1.,5.),
                 (nothing,log(3)),
                 (nothing,nothing),
                 (nothing,nothing),
                 (nothing,nothing)],
    status = [(:BORDER_FOUND_BY_SCAN_TOL,:BORDER_FOUND_BY_SCAN_TOL),
              (:SCAN_BOUND_REACHED,:BORDER_FOUND_BY_SCAN_TOL),
              (:SCAN_BOUND_REACHED,:SCAN_BOUND_REACHED),
              (:SCAN_BOUND_REACHED,:SCAN_BOUND_REACHED),
              (:SCAN_BOUND_REACHED,:SCAN_BOUND_REACHED)],
    loss_crit = 9.
  ),

  :f_3p_im => (
    func = f_3p_im,
    x0 = [3.,1.,1,],
    endpoints = [(1.,5.),
                 (nothing,log(3)),
                 (nothing,nothing)],
    status = [(:BORDER_FOUND_BY_SCAN_TOL,:BORDER_FOUND_BY_SCAN_TOL),
              (:SCAN_BOUND_REACHED,:BORDER_FOUND_BY_SCAN_TOL),
              (:SCAN_BOUND_REACHED,:SCAN_BOUND_REACHED)],
    loss_crit = 9.
  )
)


text:
using Dash, DashHtmlComponents, DashCoreComponents

app = dash()

app.layout = html_div() do
    dcc_input(id = "input-4", value = "1", type = "text"),
    html_tr((html_td("x^2 ="), html_td(id = "square"))),
    html_tr((html_td("x^3 ="), html_td(id = "cube"))),
    html_tr((html_td("2^x ="), html_td(id = "twos"))),
    html_tr((html_td("3^x ="), html_td(id = "threes"))),
    html_tr((html_td("x^x ="), html_td(id = "xx")))
end

callback!(
    app,
    Output("square", "children"),
    Output("cube", "children"),
    Output("twos", "children"),
    Output("threes", "children"),
    Output("xx", "children"),
    Input("input-4", "value"),
) do x
    if x == "" || x == nothing
        return ("", "", "", "", "")
    end

    x = parse(Int64, x)
    return (x^2, x^3, 2^x, 3^x, x^x)
end

run_server(app, "0.0.0.0", debug=true)

max_stars_repo_name: joelostblom/dash-docs

max_issues_repo_name: joelostblom/dash-docs

max_forks_repo_name: joelostblom/dash-docs


text:
# Test functionality of PDDL axioms / derived predicates
@testset "axioms" begin

path = joinpath(dirname(pathof(PDDL)), "..", "test", "axioms")

domain = load_domain(joinpath(path, "domain.pddl"))
@test domain.name == Symbol("blocksworld-axioms")
@test convert(Term, domain.predicates[:above]) == pddl"(above ?x ?y)"
@test Clause(pddl"(handempty)", [pddl"(forall (?x) (not (holding ?x)))"]) in
      values(domain.axioms)

problem = load_problem(joinpath(path, "problem.pddl"))
@test problem.name == Symbol("blocksworld-problem")
@test problem.objects == @pddl("a", "b", "c")

state = initstate(domain, problem)
state = execute(domain, state, pddl"(pickup b)")
@test domain[state => pddl"(holding b)"] == true
state = execute(domain, state, pddl"(stack b c)")
@test domain[state => pddl"(on b c)"] == true
state = execute(domain, state, pddl"(pickup a)")
@test domain[state => pddl"(holding a)"] == true
state = execute(domain, state, pddl"(stack a b)")
@test domain[state => pddl"(above a c)"] == true

@test satisfy(domain, state, problem.goal) == true

end # axioms


text:
module BernoulliTest

using Test
using ForneyLab
using ForneyLab: outboundType, isApplicable, prod!, unsafeMean, unsafeVar, vague, dims, logPdf, naturalParams, standardDistribution
using ForneyLab: SPBernoulliOutNP, SPBernoulliIn1PN, SPBernoulliOutNB, VBBernoulliOut, VBBernoulliIn1

@testset "Bernoulli ProbabilityDistribution and Message construction" begin
    @test ProbabilityDistribution(Univariate, Bernoulli, p=0.2) == ProbabilityDistribution{Univariate, Bernoulli}(Dict(:p=>0.2))
    @test_throws Exception ProbabilityDistribution(Multivariate, Bernoulli)
    @test ProbabilityDistribution(Bernoulli, p=0.2) == ProbabilityDistribution{Univariate, Bernoulli}(Dict(:p=>0.2))
    @test ProbabilityDistribution(Bernoulli) == ProbabilityDistribution{Univariate, Bernoulli}(Dict(:p=>0.5))
    @test Message(Bernoulli) == Message{Bernoulli, Univariate}(ProbabilityDistribution{Univariate, Bernoulli}(Dict(:p=>0.5)))
    @test Message(Univariate, Bernoulli) == Message{Bernoulli, Univariate}(ProbabilityDistribution{Univariate, Bernoulli}(Dict(:p=>0.5)))
    @test_throws Exception Message(Multivariate, Bernoulli)
end

@testset "dims" begin
    @test dims(ProbabilityDistribution(Bernoulli, p=0.5)) == ()
end

@testset "vague" begin
    @test vague(Bernoulli) == ProbabilityDistribution(Bernoulli, p=0.5)
end

@testset "unsafe mean and variance" begin
    @test unsafeMean(ProbabilityDistribution(Bernoulli, p=0.2)) == 0.2
    @test unsafeVar(ProbabilityDistribution(Bernoulli, p=0.5)) == 0.25
end

@testset "log pmf" begin
    @test isapprox(logPdf(ProbabilityDistribution(Bernoulli, p=0.2),1), -1.6094379124341003)
end

@testset "prod!" begin
    @test ProbabilityDistribution(Bernoulli, p=0.2) * ProbabilityDistribution(Bernoulli, p=0.8) == ProbabilityDistribution(Bernoulli, p=0.5000000000000001)
    @test_throws Exception ProbabilityDistribution(Bernoulli, p=0.0) * ProbabilityDistribution(Bernoulli, p=1.0)
end

@testset "natural parameters" begin
    d = ProbabilityDistribution(Univariate, Bernoulli, p=0.2)
    η = naturalParams(d)
    s = standardDistribution(Univariate, Bernoulli, η=η)
    @test d.params[:p] == s.params[:p] # Test conversion consistency

    x = [0.0, 1.0]
    d_x = logPdf.([d], x)
    η_x = logPdf.(Univariate, Bernoulli, x; η=η)
    @test isapprox(d_x, η_x) # Test pdf consistency
end

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

@testset "SPBernoulliOutNP" begin
    @test SPBernoulliOutNP <: SumProductRule{Bernoulli}
    @test outboundType(SPBernoulliOutNP) == Message{Bernoulli}
    @test isApplicable(SPBernoulliOutNP, [Nothing, Message{PointMass}])

    @test ruleSPBernoulliOutNP(nothing, Message(Univariate, PointMass, m=1.0)) == Message(Univariate, Bernoulli, p=1.0)
end

@testset "SPBernoulliIn1PN" begin
    @test SPBernoulliIn1PN <: SumProductRule{Bernoulli}
    @test outboundType(SPBernoulliIn1PN) == Message{Beta}
    @test isApplicable(SPBernoulliIn1PN, [Message{PointMass}, Nothing])

    @test ruleSPBernoulliIn1PN(Message(Univariate, PointMass, m=1.0), nothing) == Message(Univariate, Beta, a=2.0, b=1.0)
end

@testset "SPBernoulliOutNB" begin
    @test SPBernoulliOutNB <: SumProductRule{Bernoulli}
    @test outboundType(SPBernoulliOutNB) == Message{Bernoulli}
    @test isApplicable(SPBernoulliOutNB, [Nothing, Message{Beta}])

    @test ruleSPBernoulliOutNB(nothing, Message(Univariate, Beta, a=1.0, b=1.0)) == Message(Univariate, Bernoulli, p=0.5)
end

@testset "VBBernoulliOut" begin
    @test VBBernoulliOut <: NaiveVariationalRule{Bernoulli}

    @test outboundType(VBBernoulliOut) == Message{Bernoulli}
    @test isApplicable(VBBernoulliOut, [Nothing, ProbabilityDistribution])
    @test !isApplicable(VBBernoulliOut, [ProbabilityDistribution, Nothing])

    @test ruleVBBernoulliOut(nothing, ProbabilityDistribution(Univariate, PointMass, m=0.2)) == Message(Univariate, Bernoulli, p=0.2)
    @test ruleVBBernoulliOut(nothing, ProbabilityDistribution(Univariate, Beta, a=1.0, b=1.0)) == Message(Univariate, Bernoulli, p=0.5)
end

@testset "VBBernoulliIn1" begin
    @test VBBernoulliIn1 <: NaiveVariationalRule{Bernoulli}
    @test outboundType(VBBernoulliIn1) == Message{Beta}
    @test isApplicable(VBBernoulliIn1, [ProbabilityDistribution, Nothing])

    @test ruleVBBernoulliIn1(ProbabilityDistribution(Univariate, Bernoulli, p=0.2), nothing) == Message(Univariate, Beta, a=1.2, b=1.8)
end

@testset "averageEnergy and differentialEntropy" begin
    @test differentialEntropy(ProbabilityDistribution(Univariate, Bernoulli, p=0.25)) == averageEnergy(Bernoulli, ProbabilityDistribution(Univariate, Bernoulli, p=0.25), ProbabilityDistribution(Univariate, PointMass, m=0.25))
end

end # module


text:
using SparseArrays
using LightGraphs

struct GraphIndex{T}
    n_nodes::T
    n_edges::T
    neighbors_per_node::Array{T}
    node_neighbors::Array{T, 2}
    population::Array{T}
    src_dst_to_edge::SparseMatrixCSC{T, T}
    edges::Array{T, 2}
end

function GraphIndex(graph::SimpleGraph, population::Array{Int64})::GraphIndex{Int64}
    n_nodes = nv(graph)
    n_edges = ne(graph)
    all_edges = zeros(Int64, 2, n_edges)
    cut_edges = BitSet(Int64[])

    src_dst_to_edge = spzeros(Int64, n_nodes, n_nodes)
    for (index, edge) in enumerate(edges(graph))
        all_edges[1, index] = src(edge)
        all_edges[2, index] = dst(edge)
        src_dst_to_edge[src(edge), dst(edge)] = index
        src_dst_to_edge[dst(edge), src(edge)] = index
    end

    node_neighbors = zeros(Int64, 5, n_nodes) # Every node in a planar graph has degree ≦5.
    neighbors_per_node = zeros(Int64, n_nodes)
    for index in 1:n_nodes
        for (neighbor_idx, neighbor) in enumerate(neighbors(graph, index))
            node_neighbors[neighbor_idx, index] = neighbor
            neighbors_per_node[index] += 1
        end
    end

    return GraphIndex{Int64}(n_nodes, n_edges, neighbors_per_node,
                      node_neighbors, population, src_dst_to_edge,
                      all_edges)
end

max_stars_repo_name: pjrule/gridwalk

max_issues_repo_name: pjrule/gridwalk

max_forks_repo_name: pjrule/gridwalk


text:
# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

using DataStructures, JSON, GZip

"""

    read_egret_solution(path::String)::OrderedDict

Read a JSON solution file produced by EGRET and transforms it into a
dictionary having the same structure as the one produced by
UnitCommitment.solution(model).
"""
function read_egret_solution(path::String)::OrderedDict
    egret = _read_json(path)
    T = length(egret["system"]["time_keys"])

    solution = OrderedDict()
    is_on = solution["Is on"] = OrderedDict()
    production = solution["Production (MW)"] = OrderedDict()
    reserve = solution["Reserve (MW)"] = OrderedDict()
    production_cost = solution["Production cost (\$)"] = OrderedDict()
    startup_cost = solution["Startup cost (\$)"] = OrderedDict()

    for (gen_name, gen_dict) in egret["elements"]["generator"]
        if endswith(gen_name, "_T") || endswith(gen_name, "_R")
            gen_name = gen_name[1:end-2]
        end
        if "commitment" in keys(gen_dict)
            is_on[gen_name] = gen_dict["commitment"]["values"]
        else
            is_on[gen_name] = ones(T)
        end
        production[gen_name] = gen_dict["pg"]["values"]
        if "rg" in keys(gen_dict)
            reserve[gen_name] = gen_dict["rg"]["values"]
        else
            reserve[gen_name] = zeros(T)
        end
        startup_cost[gen_name] = zeros(T)
        production_cost[gen_name] = zeros(T)
        if "commitment_cost" in keys(gen_dict)
            for t in 1:T
                x = gen_dict["commitment"]["values"][t]
                commitment_cost = gen_dict["commitment_cost"]["values"][t]
                prod_above_cost = gen_dict["production_cost"]["values"][t]
                prod_base_cost = gen_dict["p_cost"]["values"][1][2] * x
                startup_cost[gen_name][t] = commitment_cost - prod_base_cost
                production_cost[gen_name][t] = prod_above_cost + prod_base_cost
            end
        end
    end
    return solution
end


text:
using SententialDecisionDiagrams

# set up vtree and manager
vtree = SententialDecisionDiagrams.read_vtree("$(@__DIR__)/input/simple.vtree")
sdd = SententialDecisionDiagrams.sdd_manager(vtree)

println("Created an SententialDecisionDiagrams with $(SententialDecisionDiagrams.var_count(sdd)) variables")
root = SententialDecisionDiagrams.read_cnf("$(@__DIR__)/input/simple.cnf", sdd, compiler_options=SententialDecisionDiagrams.CompilerOptions(vtree_search_mode=-1))
# For DNF functions use `read_dnf_file`

# Model Counting
wmc = SententialDecisionDiagrams.wmc_manager(root, log_mode=true)
w = SententialDecisionDiagrams.propagate(wmc)
println("Model count: $(convert(Int32,exp(w)))")

# Weighted Model Counting
lits = [SententialDecisionDiagrams.literal(i,sdd) for i in 1:SententialDecisionDiagrams.var_count(sdd)]

# Positive literal weight
SententialDecisionDiagrams.set_literal_weight(lits[1], log(0.5), wmc)
# Negative literal weight
SententialDecisionDiagrams.set_literal_weight(~lits[1], log(0.5), wmc)

w = SententialDecisionDiagrams.propagate(wmc)
println("Weighted model count: $(exp(w))")

# Visualize SententialDecisionDiagrams and VTREE
println("saving sdd and vtree ... ")
SententialDecisionDiagrams.dot("$(@__DIR__)/output/simple-vtree.dot", vtree)
SententialDecisionDiagrams.dot("$(@__DIR__)/output/simple-sdd-cnf.dot", root)


text:
using Distributed
addprocs(5)

using MendelIHT
using SnpArrays
using DataFrames
using Distributions
using Random
using LinearAlgebra
using GLM
using CSV
using DelimitedFiles

# function for reading in result after cv
function read_mse_results()
	filename = "mse_model_"
	mses = zeros(50)
	for i in 1:50
		cur_file = filename * string(i)
		mses[i] = (isfile(cur_file) ? readdlm(cur_file)[1] : Inf)
	end
	return mses
end
mses = read_mse_results()
findmin(mses)

# using Plots
# plot(mses)

# function for converting nongenetic covariates
function make_nongenetic_covariates()
	z = CSV.read("Covariate_polrGWAS.csv")
	z = [z[:, 3:6] z[:, 8] z[:, 12:end]]
	rename!(z, :x1 => :bmi)
	CSV.write("Covariate_Final.csv", z)
end

# runs a small example so cross validation code compiles
function run_iht_once_to_force_compile(d, l)
	n = 1000
	p = 10000
	k = 10

	#set random seed
	Random.seed!(2019)

	#construct x matrix and non genetic covariate (intercept)
	T = Float64
	x = simulate_random_snparray(n, p, undef)
	xbm = SnpBitMatrix{T}(x, model=ADDITIVE_MODEL, center=true, scale=true);
	z = ones(T, n, 1)

	# simulate response, true model b, and the correct non-0 positions of b
	true_b = zeros(T, p)
	true_b[1:k] .= collect(0.1:0.1:1.0)
	true_c = [T.(4.0)]
	shuffle!(true_b)
	correct_position = findall(!iszero, true_b)

	#simulate phenotypes (e.g. vector y)
	if d == Normal || d == Bernoulli || d == Poisson
	    prob = linkinv.(l, xbm * true_b)
	    clamp!(prob, -20, 20)
	    y = [rand(d(i)) for i in prob]
	    # prob = linkinv.(l, x * true_b + z * true_c)
	    # clamp!(prob, -20, 20)
	    # y = [rand(d(i)) for i in prob]
	    # k = k + 1
	elseif d == NegativeBinomial
	    nn = 10
	    μ = linkinv.(l, xbm * true_b)
	    # μ = linkinv.(l, x * true_b + z * true_c)
	    # k = k + 1
	    clamp!(μ, -20, 20)
	    prob = 1 ./ (1 .+ μ ./ nn)
	    y = [rand(d(nn, Float64(i))) for i in prob] #number of failtures before nn success occurs
	elseif d == Gamma
	    μ = linkinv.(l, xbm * true_b)
	    β = 1 ./ μ # here β is the rate parameter for gamma distribution
	    y = [rand(d(α, i)) for i in β] # α is the shape parameter for gamma
	end
	y = T.(y)

	#specify path and folds
	path = collect(1:10)
	num_folds = 5
	folds = rand(1:num_folds, size(x, 1))

	# run IHT
	scratch_folder = "/u/scratch/b/biona001/"
	cv_iht_distribute_fold(d(), l, x, z, y, 1, path, num_folds, destin=scratch_folder, showinfo=false, folds=folds, debias=false, parallel=true)
end

#main function
function run(cur_path::Int64)
	# d = Normal
	d = Bernoulli
	l = canonicallink(d())
	run_iht_once_to_force_compile(d, l)

	#Import covariates (data is in $SCRATCH on hoffman2)
	filedir = "/u/scratch/b/biona001/ukbiobank/"
	x = SnpArray(filedir * "ukb.plink.filtered.imputed.bed")
	z = CSV.read(filedir * "Covariate_Final.csv")
	z = convert(Matrix{Float64}, z)
	standardize!(z)
	z = [ones(size(z, 1)) z]

	# quantitative phenotype
	# y = CSV.read(filedir * "Covariate_polrGWAS.csv")[:AveSBP] # average systolic blood pressure
	# y = convert(Vector{Float64}, y)

	# binary phenotype
	sbp = CSV.read(filedir * "Covariate_polrGWAS.csv")[:AveSBP] # average systolic blood pressure
	dbp = CSV.read(filedir * "Covariate_polrGWAS.csv")[:AveDBP] # average diastolic blood pressure
	y = zeros(Float64, length(sbp))
	y[sbp .>= 140] .= 1.0
	y[dbp .>= 90]  .= 1.0

	#set random seed
	Random.seed!(2019)

	#specify path and folds
	path = [cur_path]
	num_folds = 5
	folds = rand(1:num_folds, size(x, 1))

	# run cross validation
	scratch_folder = "/u/scratch/b/biona001/"
	mses = cv_iht_distribute_fold(d(), l, x, z, y, 1, path, num_folds, destin=scratch_folder, folds=folds, debias=false, parallel=true)
	writedlm("mse_model_$(path[1])", mses)
end

# arguments passed from command line on hoffman cluster
const cur_path = parse(Int64, ARGS[1])

run(cur_path)


text:
export newton, is_success, NewtonResult, NewtonCache

"""
    NewtonResult

Result returned by [`newton`](@ref).

## Fields

* `return_code::Symbol`: Can be `:success`, `:rejected` or `:max_iters`.
* `x::Vector{ComplexF64}`: The last value obtained.
* `accuracy::Float64`: Estimate of the distance of `x` to a true zero.
* `iters::Int` Number of iterations performed.
* `contraction_ratio::Float64`: The value `|xᵢ - xᵢ₋₁| / |xᵢ₋₁ - xᵢ₋₂| `.
"""
struct NewtonResult
    return_code::Symbol
    x::Vector{ComplexF64}
    accuracy::Float64
    residual::Float64
    iters::Int
    contraction_ratio::Float64
end

Base.show(io::IO, ::MIME"application/prs.juno.inline", r::NewtonResult) = r
Base.show(io::IO, result::NewtonResult) = print_fieldnames(io, result)

"""
    is_success(R::NewtonResult)

Returns `true` if the [`newton`](@ref) was successfull.
"""
is_success(R::NewtonResult) = R.return_code == :success

"""
    NewtonCache(F::AbstractSystem)

Pre-allocates the necessary memory for [`newton`](@ref).
"""
struct NewtonCache{M}
    x::Vector{ComplexF64}
    Δx::Vector{ComplexF64}
    x_ext::Vector{ComplexDF64}
    J::MatrixWorkspace{M}
    r::Vector{ComplexF64}
end
function NewtonCache(F::AbstractSystem)
    m, n = size(F)
    x = zeros(ComplexF64, n)
    Δx = copy(x)
    x_ext = zeros(ComplexDF64, n)
    J = MatrixWorkspace(m, n)
    r = zeros(ComplexF64, m)
    NewtonCache(x, Δx, x_ext, J, r)
end

"""
    newton(
        F::AbstractSystem,
        x₀::AbstractVector,
        norm::AbstractNorm = InfNorm(),
        cache::NewtonCache = NewtonCache(F, x₀);
        options...
    )

An implemenetation of a local Newton's method with various options to specify
convergence criteria.
Returns a [`NewtonResult`](@ref).
The computations are always performed in complex arithmetic with
double precision, i.e., using `Complex{Float64}`.
The optional `cache` argument pre-allocates the necessary memory.
This is useful if the method is called repeatedly.

### Options

* `abs_tol::Float64 = 1e-8`: The method is declared converged if `norm(xᵢ₊₁ - xᵢ) < abs_tol`.
* `rel_tol::Float64 = abs_tol`: The method is declared converged if
  `norm(xᵢ₊₁ - xᵢ) < abs_tol * norm(x₀)`.
* `max_iters::Int = 20`: The maximal number of iterations.
* `extended_precision::Bool = false`: An optional use of extended precision for the evaluation
  of `F(x)`. This can increase the achievable accuracy.
* `contraction_factor::Float64 = 1.0`: The Newton updates have to satisfy
  ``|xᵢ₊₁ - xᵢ| < a^{2^(i-1)}|x₁ - x₀|`` for ``i ≥ 1`` where ``a`` is `contraction_factor`.
* `min_contraction_iters::Int = typemax(Int)`:  The minimal number of iterations the
  `contraction_factor` has to be satisfied. If after `min_contraction_iters` many
  iterations the contraction factor is not satisfied the step is accepted anyway.
* `max_abs_norm_first_update::Float64 = Inf`: The initial guess `x₀` is rejected if
  `norm(x₁ - x₀) >  max_abs_norm_first_update`
* `max_rel_norm_first_update::Float64 = max_abs_norm_first_update`: The initial guess `x₀`
  is rejected if `norm(x₁ - x₀) >  max_rel_norm_first_update * norm(x₀)`
"""
newton(f::System, args...; kwargs...) = newton(ModelKitSystem(f), args...; kwargs...)
function newton(
    F::AbstractSystem,
    x₀::AbstractVector,
    norm::AbstractNorm = InfNorm(),
    cache::NewtonCache = NewtonCache(F);
    extended_precision::Bool = false,
    abs_tol::Float64 = 1e-8,
    rel_tol::Float64 = abs_tol,
    max_iters::Int = 20,
    contraction_factor::Float64 = 1.0,
    min_contraction_iters::Int = typemax(Int),
    max_abs_norm_first_update::Float64 = Inf,
    max_rel_norm_first_update::Float64 = max_abs_norm_first_update,
)
    @unpack x, Δx, x_ext, J, r = cache
    x .= x₀
    a = contraction_factor
    norm_x = norm(x)
    norm_Δxᵢ = norm_Δxᵢ₋₁ = NaN
    res = NaN
    for i = 1:max_iters
        evaluate_and_jacobian!(r, matrix(J), F, x)
        if extended_precision
            x_ext .= x
            evaluate!(r, F, x_ext)
        end

        updated!(J)
        LA.ldiv!(Δx, J, r)
        norm_Δxᵢ = norm(Δx)

        x .= x .- Δx

        if i > 1
            if norm_Δxᵢ < a * norm_Δxᵢ₋₁
                norm_Δxᵢ₋₁ = norm_Δxᵢ
                a *= a
            elseif i > min_contraction_iters ||
                   norm_Δxᵢ < rel_tol * norm_x ||
                   norm_Δxᵢ < abs_tol
                return NewtonResult(
                    :success,
                    x,
                    norm_Δxᵢ,
                    norm(r),
                    i,
                    norm_Δxᵢ / norm_Δxᵢ₋₁,
                )
            else
                return NewtonResult(
                    :rejected,
                    x,
                    norm_Δxᵢ,
                    norm(r),
                    i,
                    norm_Δxᵢ / norm_Δxᵢ₋₁,
                )
            end
        elseif i == 1 && (
            norm_Δxᵢ > max_rel_norm_first_update * norm_x ||
            norm_Δxᵢ > max_abs_norm_first_update
        )
            return NewtonResult(:rejected, x, norm(r), norm_Δxᵢ, i, NaN)
        else
            norm_Δxᵢ₋₁ = norm_Δxᵢ
        end
    end
    return NewtonResult(:max_iters, x, norm_Δxᵢ, norm(r), max_iters, norm_Δxᵢ / norm_Δxᵢ₋₁)
end


text:
__precompile__()
module Signals

export Signal

const _async_mode = Ref(true)
async_mode() = _async_mode.x
async_mode(b::Bool) = _async_mode.x = b

const _debug_mode = Ref(false)
debug_mode() = _debug_mode.x
debug_mode(b::Bool) = _debug_mode.x = b

include("PullAction.jl")
include("Signal.jl")
include("error.jl")
include("pushpull.jl")
include("eventloop.jl")
include("operators.jl")
include("bind.jl")
include("async_remote.jl")
include("time.jl")
include("TypedSignal.jl")

@doc """
    S = Signal(val; strict_push = false)

Create a source `Signal` with initial value `val`, setting `strict_push` to `true`
guarantees that every push to this `Signal` will be carried out independently.
Otherwise if updates occur faster than what the `eventloop` can process, then only
the last value before the `eventloop` kicks in will be used(*default*).

    S = Signal(f,args...; v0 = nothing)

Create a derived `Signal` whos value is `f(args...)`, args can be of any type,
`Signal` args get replaced by their value before calling `f(args...)`. reads best with
with `do` notation(see example).if `v0` is not `nothing` then `f(args...)` will not
be called directly after Signal creation instead the Signal will be initialized to have value v0.

# Syntax

    S[] = val
Set the value of `S` to `val` without propogating the change to the rest of the signal graph,
usefull in pull based paradigm.

    S()
`pull!` Signal, pulling any changes in the Signal graph that affects `S`.

    S(val)
Set the value of `S` to `val` and pushes the changes along the Signal graph.

    S[]
Get the current value stored in `S` without pulling changes from the graph.

# Examples

    julia> A = Signal(1) # source Signal
    Signal
    value: 1

    julia> B = 2 # non-Signal input
    2

    julia> C = Signal(A, B) do a, b # derived Signal
               a + b
           end

    Signal
    value: 3

    julia> A[] = 10 # set value without propogation
    10
    julia> C[] # reads current value
    3
    julia> C() # pull latest changes from the Signal graph
    12
    julia> A(100) # set value to a signal and propogate this change
    100
    julia> C[]
    102
""" Signal

end # module


text:
function hilbert2int(coords::Vector{Int64})
  nD = length(coords)
  coord_chunks = unpack_coords(coords)
  nChunks = length(coord_chunks)
  mask = 2^nD - 1
  start, finish = initial_start_end(nChunks, nD)
  #Int here?
  index_chunks = zeros(Int64,nChunks)
  for j = 0:(nChunks-1)
    i = gray_decode_travel(start, finish, mask, coord_chunks[j+1])
    index_chunks[j+1] = i
    start, finish = child_start_end(start, finish, mask, i)
  end
  return pack_index(index_chunks, nD)
end

function unpack_coords(coords::Vector{Int64})
  nD = length(coords)
  #biggest = reduce(max, coords)  # the max of all coords
  #nChunks = max(1, ceil(Int64,log(2, biggest + 1))) # max # of bits
  nChunks = max(1, ceil(Int64,log(2, maximum(coords) + 1))) # max # of bits
  return transpose_bits(coords, nChunks)
end

function pack_indexSlow(chunks::Vector{Int64}, nD::Int64)
  p = 2^nD  # Turn digits mod 2**nD back into a single number:
  #return reduce( lambda n, chunk: n * p + chunk, chunks )
  return reduce((n,chunk)->n*p+chunk, chunks)
end

function pack_index(chunks::Vector{Int64}, nD::Int64)
  p = 2^nD  # Turn digits mod 2**nD back into a single number:
  #return reduce( lambda n, chunk: n * p + chunk, chunks )
  sum = 0::Int64
  @simd for i = 1:length(chunks)
    sum *= p
    sum += chunks[i]
  end
  return sum
end

function initial_start_end(nChunks::Int64, nD::Int64)
  # This orients the largest cube so that
  # its start is the origin (0 corner), and
  # the first step is along the x axis, regardless of nD and nChunks:
  return 0, 2^mod(-nChunks-1,nD)  # in Python 0 <=  a % b  < b.
end

function gray_encode( bn::Int64)
    assert(bn >= 0)
    return bn $ fld(bn,2)
end

function gray_decode(n::Int64)
  sh = 1
  while true
    div = n >> sh
    n $= div
    if div <= 1
      return n
    end
    sh = sh << 1
  end
end

function gray_decode_travel(start::Int64, finish::Int64, mask::Int64, g::Int64)
  travel_bit = start $ finish
  modulus = mask + 1          # == 2**nBits
  rg = (g $ start)*fld(modulus ,travel_bit*2)
  return gray_decode((rg | fld(rg,modulus)) & mask)
end

## gray_encode_travel -- gray_encode given start and end using bit rotation.
#    Modified Gray code.  mask is 2**nbits - 1, the highest i value, so
#        gray_encode_travel( start, end, mask, 0 )    == start
#        gray_encode_travel( start, end, mask, mask ) == end
#        with a Gray-code-like walk in between.
#    This method takes the canonical Gray code, rotates the output word bits,
#    then xors ("^" in Python) with the start value.
#
function gray_encode_travel(start::Int64, finish::Int64, mask::Int64, i::Int64)
    travel_bit = start $ finish
    modulus = mask + 1          # == 2**nBits
    # travel_bit = 2**p, the bit we want to travel.
    # Canonical Gray code travels the top bit, 2**(nBits-1).
    # So we need to rotate by ( p - (nBits-1) ) == (p + 1) mod nBits.
    # We rotate by multiplying and dividing by powers of two:
    g = gray_encode( i ) * ( travel_bit * 2 )
    return ((g | fld(g, modulus)) & mask) $ start
end

## child_start_end( parent_start, parent_end, mask, i ) -- Get start & end for child.
#    i is the parent's step number, between 0 and mask.
#    Say that parent( i ) =
#           gray_encode_travel( parent_start, parent_end, mask, i ).
#    And child_start(i) and child_end(i) are what child_start_end()
#    should return -- the corners the child should travel between
#    while the parent is in this quadrant or child cube.
#      o  child_start( 0 ) == parent( 0 )       (start in a corner)
#      o  child_end( mask ) == parent( mask )   (end in a corner)
#      o  child_end(i) - child_start(i+1) == parent(i+1) - parent(i)
#         (when parent bit flips, same bit of child flips the opposite way)
#    Those constraints still leave choices when nD (# of bits in mask) > 2.
#    Here is how we resolve them when nD == 3 (mask == 111 binary),
#    for parent_start = 000 and parent_end = 100 (canonical Gray code):
#         i   parent(i)    child_
#         0     000        000   start(0)    = parent(0)
#                          001   end(0)                   = parent(1)
#                 ^ (flip)   v
#         1     001        000   start(1)    = parent(0)
#                          010   end(1)                   = parent(3)
#                ^          v
#         2     011        000   start(2)    = parent(0)
#                          010   end(2)                   = parent(3)
#                 v          ^
#         3     010        011   start(3)    = parent(2)
#                          111   end(3)                   = parent(5)
#               ^          v
#         4     110        011   start(4)    = parent(2)
#                          111   end(4)                   = parent(5)
#                 ^          v
#         5     111        110   start(5)    = parent(4)
#                          100   end(5)                   = parent(7)
#                v          ^
#         6     101        110   start(6)    = parent(4)
#                          100   end(6)                   = parent(7)
#                 v          ^
#         7     100        101   start(7)    = parent(6)
#                          100   end(7)                   = parent(7)
#    This pattern relies on the fact that gray_encode_travel()
#    always flips the same bit on the first, third, fifth, ... and last flip.
#    The pattern works for any nD >= 1.
#
function child_start_end(parent_start::Int64, parent_end::Int64, mask::Int64, i::Int64)
  start_i = max(0,    (i - 1) & ~1)  # next lower even number, or 0
  end_i =   min(mask, (i + 1) |  1)  # next higher odd number, or mask
  child_start = gray_encode_travel(parent_start, parent_end, mask, start_i)
  child_end   = gray_encode_travel(parent_start, parent_end, mask, end_i)
  return child_start, child_end
end

## transpose_bits --
#    Given nSrcs source ints each nDests bits long,
#    return nDests ints each nSrcs bits long.
#    Like a matrix transpose where ints are rows and bits are columns.
#    Earlier srcs become higher bits in dests;
#    earlier dests come from higher bits of srcs.

@inbounds function transpose_bits(srcs::Vector{Int64}, nDests::Int64)
  #srcs = copy(srcs)  # Make a copy we can modify safely.
  nSrcs = length(srcs)
  dests = Array{Int64}(nDests)#zeros(Int64,nDests)
  # Break srcs down least-significant bit first, shifting down:
  for j = (nDests - 1):-1:0
    # Put dests together most-significant first, shifting up:
    dest = 0
    for k = 0:(nSrcs-1)
      dest = dest*2 + mod(srcs[k+1],2)
      srcs[k+1] = fld(srcs[k+1],2)
    end
    dests[j+1] = dest
  end
  return dests
end

#     Overwrites srcs! and writes data to dests,
#apparently slower!!!!
@inbounds function transpose_bits!(srcs::Vector{Int64}, nDests::Int64, dests::Vector{Int64})
  nSrcs = length(srcs)
  # Break srcs down least-significant bit first, shifting down:
  for j = (nDests - 1):-1:0
    # Put dests together most-significant first, shifting up:
    dest = 0
    for k = 0:(nSrcs-1)
      dest = dest*2 + mod(srcs[k+1],2)
      srcs[k+1] = fld(srcs[k+1],2)
    end
    dests[j+1] = dest
  end
  return dests
end


text:
function NaturalRanking()
    return NaturalRanking(())
end

function NaturalRanking(arg0::NaNStrategy)
    return NaturalRanking((NaNStrategy,), arg0)
end

function NaturalRanking(arg0::NaNStrategy, arg1::RandomGenerator)
    return NaturalRanking((NaNStrategy, RandomGenerator), arg0, arg1)
end

function NaturalRanking(arg0::NaNStrategy, arg1::TiesStrategy)
    return NaturalRanking((NaNStrategy, TiesStrategy), arg0, arg1)
end

function NaturalRanking(arg0::RandomGenerator)
    return NaturalRanking((RandomGenerator,), arg0)
end

function NaturalRanking(arg0::TiesStrategy)
    return NaturalRanking((TiesStrategy,), arg0)
end

function equals(obj::Object, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_nan_strategy(obj::NaturalRanking)
    return jcall(obj, "getNanStrategy", NaNStrategy, ())
end

function get_ties_strategy(obj::NaturalRanking)
    return jcall(obj, "getTiesStrategy", TiesStrategy, ())
end

function hash_code(obj::Object)
    return jcall(obj, "hashCode", jint, ())
end

function rank(obj::NaturalRanking, arg0::Vector{jdouble})
    return jcall(obj, "rank", Vector{jdouble}, (Vector{jdouble},), arg0)
end

function to_string(obj::Object)
    return jcall(obj, "toString", JString, ())
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct VirtualMachineInstanceView <: SwaggerModel
    platformUpdateDomain::Any # spec type: Union{ Nothing, Int32 } # spec name: platformUpdateDomain
    platformFaultDomain::Any # spec type: Union{ Nothing, Int32 } # spec name: platformFaultDomain
    rdpThumbPrint::Any # spec type: Union{ Nothing, String } # spec name: rdpThumbPrint
    vmAgent::Any # spec type: Union{ Nothing, VirtualMachineAgentInstanceView } # spec name: vmAgent
    maintenanceRedeployStatus::Any # spec type: Union{ Nothing, MaintenanceRedeployStatus } # spec name: maintenanceRedeployStatus
    disks::Any # spec type: Union{ Nothing, Vector{DiskInstanceView} } # spec name: disks
    extensions::Any # spec type: Union{ Nothing, Vector{VirtualMachineExtensionInstanceView} } # spec name: extensions
    bootDiagnostics::Any # spec type: Union{ Nothing, BootDiagnosticsInstanceView } # spec name: bootDiagnostics
    statuses::Any # spec type: Union{ Nothing, Vector{InstanceViewStatus} } # spec name: statuses

    function VirtualMachineInstanceView(;platformUpdateDomain=nothing, platformFaultDomain=nothing, rdpThumbPrint=nothing, vmAgent=nothing, maintenanceRedeployStatus=nothing, disks=nothing, extensions=nothing, bootDiagnostics=nothing, statuses=nothing)
        o = new()
        validate_property(VirtualMachineInstanceView, Symbol("platformUpdateDomain"), platformUpdateDomain)
        setfield!(o, Symbol("platformUpdateDomain"), platformUpdateDomain)
        validate_property(VirtualMachineInstanceView, Symbol("platformFaultDomain"), platformFaultDomain)
        setfield!(o, Symbol("platformFaultDomain"), platformFaultDomain)
        validate_property(VirtualMachineInstanceView, Symbol("rdpThumbPrint"), rdpThumbPrint)
        setfield!(o, Symbol("rdpThumbPrint"), rdpThumbPrint)
        validate_property(VirtualMachineInstanceView, Symbol("vmAgent"), vmAgent)
        setfield!(o, Symbol("vmAgent"), vmAgent)
        validate_property(VirtualMachineInstanceView, Symbol("maintenanceRedeployStatus"), maintenanceRedeployStatus)
        setfield!(o, Symbol("maintenanceRedeployStatus"), maintenanceRedeployStatus)
        validate_property(VirtualMachineInstanceView, Symbol("disks"), disks)
        setfield!(o, Symbol("disks"), disks)
        validate_property(VirtualMachineInstanceView, Symbol("extensions"), extensions)
        setfield!(o, Symbol("extensions"), extensions)
        validate_property(VirtualMachineInstanceView, Symbol("bootDiagnostics"), bootDiagnostics)
        setfield!(o, Symbol("bootDiagnostics"), bootDiagnostics)
        validate_property(VirtualMachineInstanceView, Symbol("statuses"), statuses)
        setfield!(o, Symbol("statuses"), statuses)
        o
    end
end # type VirtualMachineInstanceView

const _property_map_VirtualMachineInstanceView = Dict{Symbol,Symbol}(Symbol("platformUpdateDomain")=>Symbol("platformUpdateDomain"), Symbol("platformFaultDomain")=>Symbol("platformFaultDomain"), Symbol("rdpThumbPrint")=>Symbol("rdpThumbPrint"), Symbol("vmAgent")=>Symbol("vmAgent"), Symbol("maintenanceRedeployStatus")=>Symbol("maintenanceRedeployStatus"), Symbol("disks")=>Symbol("disks"), Symbol("extensions")=>Symbol("extensions"), Symbol("bootDiagnostics")=>Symbol("bootDiagnostics"), Symbol("statuses")=>Symbol("statuses"))
const _property_types_VirtualMachineInstanceView = Dict{Symbol,String}(Symbol("platformUpdateDomain")=>"Int32", Symbol("platformFaultDomain")=>"Int32", Symbol("rdpThumbPrint")=>"String", Symbol("vmAgent")=>"VirtualMachineAgentInstanceView", Symbol("maintenanceRedeployStatus")=>"MaintenanceRedeployStatus", Symbol("disks")=>"Vector{DiskInstanceView}", Symbol("extensions")=>"Vector{VirtualMachineExtensionInstanceView}", Symbol("bootDiagnostics")=>"BootDiagnosticsInstanceView", Symbol("statuses")=>"Vector{InstanceViewStatus}")
Base.propertynames(::Type{ VirtualMachineInstanceView }) = collect(keys(_property_map_VirtualMachineInstanceView))
Swagger.property_type(::Type{ VirtualMachineInstanceView }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_VirtualMachineInstanceView[name]))}
Swagger.field_name(::Type{ VirtualMachineInstanceView }, property_name::Symbol) =  _property_map_VirtualMachineInstanceView[property_name]

function check_required(o::VirtualMachineInstanceView)
    true
end

function validate_property(::Type{ VirtualMachineInstanceView }, name::Symbol, val)
end


text:
### General Utility Functions ###

# Data Read: Reads the structured data from the given instance file
# Input:  - filename: A string containing the file name with its relative path.
# Output: - n: An Integer containing the number of tasks;
#         - tasks_times: An (n) Array of Integers with each task's time;
#         - prec_relations: A (n,n) 2-dimensional Binary Array representing
#                           each task's Precedence Relation.
#         - times_sum: An Integer representing the sum of all the tasks' times.
function data_read(filename::String)
    # Gets every line in the instance file for processing
    lines = []
    open(filename,"r") do f
        lines = readlines(filename)
    end

    # The number n of tasks is given in the first line
    n = parse(Int64, lines[1])

    # Gets each task's time in the next n lines and the sum of all times
    task_times = Int64[]
    times_sum = 0
    for i = 2:(n+1)
        curr_time = parse(Int64, lines[i])
        append!(task_times, curr_time)
        times_sum += curr_time
    end

    # Gets the i,j Direct Precedence Relations until the "-1,-1" end mark
    # and transforms it into a Binary Matrix of the Relations
    prec_relations = zeros(Int64, (n,n))
    curr_line = n+2
    # Julia has no "Do While" so a first verification is necessary to avoid a break
    rs = split(lines[curr_line], ",")
    r = [parse(Int64, rs[1]), parse(Int64, rs[2])]
    while (r[1] > 0)
        prec_relations[r[1], r[2]] = 1

        curr_line += 1
        rs = split(lines[curr_line], ",")
        r = [parse(Int64, rs[1]), parse(Int64, rs[2])]
    end

    println("Data successfully read from file \"", filename, "\" containing ", n, " tasks.")
    println("")

    # Returns the structured data
    return n, task_times, prec_relations, times_sum
end

# Solution Write: Writes the final solution and time in a file with stdout format
# Input:  - filename: A string with the name for the file that will be written;
#         - solution: The final best solution encountered during execution;
#         - time: The time it took to find the given solution.
# Output: None.
function solution_write(filename::String, solution, time)
    open(filename, "w") do f
        if (time == 0)
            println(f, "No feasible solution found.\n")

            println("- FINAL RESULT: No feasible solution found.")
        else
            println(f, "Best solution: ", solution[2], "\r\n")
            println(f, "With score: ", solution[1], "\r\n")
            println(f, "Found in: ", time, " seconds.\r\n")

            println("- FINAL RESULT: ")
            println("- Best solution: ", solution[2])
            println("- With score: ", solution[1])
            println("- Found in: ", time, " seconds.")
        end
    end

    println("Solution saved in file: ", filename)
    println("")
end

# Args Read: Reads any (or provides default) Argments passed through command line
# Input:  - args: An Array of Strings containing the argments from the command line.
# Output: - file: A string with the name of the file for the instance given;
#         - stations: An Integer with the number of stations available, defaults to 1;
#         - population: An Integer with the size of the population, defaults to 100.
function args_read(args = ["instances/HAHN.IN2", "1", "500"])
    if (length(args) == 0)
        println("WARNING: No argument provided, using default instance Hahn and stations and population sizes of 1 and 100, respectively")
        file = "instances/HAHN.IN2"
        stations = 1
        population = 500
    elseif (length(args) == 1)
        println("WARNING: Only one argument provided, using default stations and population sizes of 1 and 100, respectively")
        file = args[1]
        stations = 1
        population = 500
    elseif (length(args) == 2)
        println("WARNING: Population size not provided, using default of 500")
        file = args[1]
        stations = parse(Int64, args[2])
        population = 500
    else
        file = args[1]
        stations = parse(Int64, args[2])
        population = parse(Int64, args[3])
    end

    println("Arguments successfully read. Using these parameters:")
    println("File Path and Name: \"", file, "\".")
    println("Number of Stations: ", stations, ".")
    println("Size of Population: ", population, ".")
    println("")

    return file, stations, population
end

### Genetic Algorithm Functions ###

# Rand Gene: Generates a random gene containing a solution for the problem
# Input:  - s: An Integer representing the given number of stations available;
#         - n: An Integer representing the number of tasks to be assigned.
# Output: - gene: A (s,n) 2-dimensional Integer Array representing the sequence
#                 of tasks n assigned to each station s.
function rand_gene(s::Int64, n::Int64)
    gene = zeros(Int64, (s,n))
    task_counter = ones(Int64, s)

    for i = 1:n
        rnds = rand(1:s)
        gene[rnds, task_counter[rnds]] = i
        task_counter[rnds] += 1
    end

    # 0.1% chance of just generating a linear gene
    rnd_linear = rand(1:1000)
    if (rnd_linear == 500)
        gene = zeros(Int64, (s,n))
        task_counter = 1
        threshold = ceil(Int64, (n/s))

        for i = 1:s
            for j = 1:threshold
                if (task_counter <= n)
                    gene[i, j] = task_counter
                    task_counter += 1
                end
            end
        end
    end

    return gene
end

# Rand Population:
# Input:  - size: An Integer representing the size of the generated population;
#         - s, n, task_times, times_sum, prec_relations: Needed from gene and fitness functions.
# Output: - population: A random population containing "size"s [fitness, gene] pairs.
function rand_population(size::Int64, s::Int64, n::Int64, task_times, times_sum::Int64, prec_relations)
    population = []
    for i = 1:size
        gene = rand_gene(s, n)
        fit = fitness(gene, task_times, times_sum, prec_relations)
        curr_individual = [fit, gene]
        append!(population, [curr_individual])
    end

    return population
end

# Crossover: Creates a new Gene by doing a crossover between two parents
# Input:  - parent_1: A [fitness, gene] individual representing the first parent;
#         - parent_2: A [fitness, gene] individual representing the second parent.
#         - task_times, times_sum, prec_relations: Needed from fitness function.
# Output: - new_individual: A new [fitness, gene] individual.
function crossover(parent_1, parent_2, task_times, times_sum::Int64, prec_relations)
    gene1 = parent_1[2]
    gene2 = parent_2[2]

    s = size(gene1)[1]
    n = size(gene1)[2]

    genes = [gene1, gene2]
    task_counter = ones(Int64, s)
    new_gene = zeros(Int64, (s,n))

    for i = 1:n
        rnd_parent = rand(1:2)
        for j = 1:s
            discover_task = find(genes[rnd_parent][j, :] .== i)
            if (!isempty(discover_task))
                new_gene[j, task_counter[j]] = i
                task_counter[j] += 1
            end
        end
    end

    fit = fitness(new_gene, task_times, times_sum, prec_relations)
    new_individual = [fit, new_gene]

    return new_individual
end

# Fitness: Calculates the fitness of the given gene
#          In SALBP, it's the longest cycle among all stations
#          A cycle of a station is the sum the times of all the tasks given to it
# Input:  - gene: A (s,n) 2-dimensional Integer Array representing a gene;
#         - task_times: An (n) Array of Integers with each task's time.
#         - times_sum: An Integer with the sum of all tasks' times to use as
#                      "punishment" for each infeasibility in the gene.
#         - prec_relations: A (n,n) 2-dimensional Binary Array representing
#                           each task's Precedence Relation.
# Output: - max_cycle: An Integer representing the fitness of the given gene.
function fitness(gene, task_times, times_sum::Int64, prec_relations)
    max_cycle = 0
    s = size(gene)[1]
    n = size(gene)[2]
    tasks_not_exec = ones(Int64, n)

    # For each station
    for i = 1:s
        current_cycle = 0
        # And each task assigned to it
        for j = 1:n
            curr_task_id = gene[i, j]
            if (curr_task_id > 0)
                # Adds the current task time to the current cycle total
                current_cycle += task_times[curr_task_id]
                # This is for penalities addition on infeasibilities
                for t = 1:n
                    # If a task has not yet been executed
                    if (tasks_not_exec[t] == 1)
                        # And it's in the precedence list for the current task
                        if(prec_relations[t, curr_task_id] == 1)
                            # Adds the penalty
                            current_cycle += times_sum
                        end
                    end
                end
                # Removes the current task from the not executed yet list
                tasks_not_exec[curr_task_id] = 0
            end
        end
        # Keeps the maximum cycle found
        if (current_cycle > max_cycle)
            max_cycle = current_cycle
        end
    end

    return max_cycle
end

# Genetic Algorithm:
# Input:  - init_pop: An array containing the initial [fitness, gene] pairs population
#                     of the previously defined size;
#         - pop_size, s, n, task_times, times_sum, prec_relations: Needed from gene, population
#                                                                  and fitness functions;
#         - max_gen: An integer containing the stop criteria of the GA, max generations.
# Output: - best_solution: The [fitness, gene] pair of the best solution found;
#         - time_taken: The time it took to find the best solution found;
function genetic_algorithm(init_pop, pop_size::Int64, s::Int64, n::Int64, task_times, times_sum::Int64, max_gen::Int64, prec_relations)
    tic()
    population = init_pop
    best_found_fit = times_sum + 1
    best_solution = [0, 0]
    time_taken = 0

    println("Genetic Algorithm started.")

    # Until the end criteria is hit
    for i = 1:max_gen
        population = sort(population, lt=lexless)
        curr_best = population[1][1]

        if (curr_best <= times_sum)
            if (curr_best < best_found_fit)
                best_found_fit = curr_best
                best_solution = population[1]
                time_taken += toc()
                tic()
                println("IT [", i, "] - ", "Found a better solution with score: ", best_solution[1])
            end
        end

        # Creates a new population for the next generation
        class_A = ceil(Int64, pop_size * 0.2)
        class_C = ceil(Int64, pop_size * 0.5)

        # Class A stays the same, changes start in Class B
        for i = class_A:class_C-1
            # Class B crossover
            curr_parents = i - class_A + 1
            population[i] = crossover(population[curr_parents], population[curr_parents + 1], task_times, times_sum, prec_relations)
        end
        # Class C just replaces the worst half of the population for new random genes
        population[(class_C+1):end,:] = rand_population(class_C, s, n, task_times, times_sum, prec_relations)
    end

    return best_solution, time_taken
end

### Main Function ###
function main()
    max_gen = 500
    out_filename = "main_test.txt"

    println("- Reading arguments...")
    println("")
    file_name, num_stations, population_size = args_read(ARGS)

    # TEST #
    #file_name = "instances/TEST1.IN2"
    #file_name = "instances/HAHN.IN2"
    #num_stations = 10
    #population_size = 100
    # TEST #

    println("- Setting up data from given file...")
    println("")
    num_tasks, task_times, prec_relations, times_sum = data_read(file_name)

    println("- Starting Genetic Algorithm execution...")
    println("")
    initial_population = rand_population(population_size, num_stations, num_tasks, task_times, times_sum, prec_relations)
    best_solution, time_taken = genetic_algorithm(initial_population, population_size, num_stations, num_tasks, task_times, times_sum, max_gen, prec_relations)

    solution_write(out_filename, best_solution, time_taken)

    println("- Execution ended.")
    println("")
end

main()

########## TESTS ##########

#=
file1 = "instances/HAHN.IN2"
file2 = "instances/LUTZ3.IN2"
file3 = "instances/WEE-MAG.IN2"

ns , ts, rs, bigm = data_read(file1)
println(ns)
println(ts)
for i = 1:ns
    println(i, rs[i,:])
end
println(bigm)

fs, ms, ps = args_read(["3"])

genes0 = rand_gene(ms, ns)

fs, ms, ps = args_read(["1", "2"])

genes = rand_gene(ms, ns)
println(genes)

testfits = [0, 0]
for i = 1:size(genes)[2]
    if (genes[1, i] > 0)
        testfits[1] += ts[genes[1, i]]
    end
    if (genes[2, i] > 0)
        testfits[2] += ts[genes[2, i]]
    end
end
println(testfits)

fits = fitness(genes, ts, bigm)
println(fits)

new_genes = crossover(genes0, genes)
println(new_genes)
=#

max_stars_repo_name: AlexandreSSJr/SALBP-GA

max_issues_repo_name: AlexandreSSJr/SALBP-GA

max_forks_repo_name: AlexandreSSJr/SALBP-GA