text:
"""
Single variable bound constraints
"""
# MOI.supports_constraint(::Optimizer, ::Type{SVF}, ::Type{AbstractVector{<:MOI.AbstractScalarSet}}) = true
MOI.supports_constraint(::Optimizer, ::Type{SVF}, ::Type{MOI.LessThan{T}}) where {T<:Real} =
    true

MOI.supports_constraint(
    ::Optimizer,
    ::Type{SVF},
    ::Type{MOI.GreaterThan{T}},
) where {T<:Real} = true

MOI.supports_constraint(::Optimizer, ::Type{SVF}, ::Type{MOI.EqualTo{T}}) where {T<:Real} =
    true

MOI.supports_constraint(::Optimizer, ::Type{SVF}, ::Type{MOI.Interval{T}}) where {T<:Real} =
    true

MOI.supports_constraint(::Optimizer, ::Type{SVF}, ::Type{Integers}) = true

"""
Binary/Integer variable support
"""
MOI.supports_constraint(::Optimizer, ::Type{SVF}, ::Type{<:VAR_TYPES}) = true

"""
    check_var_bounds(model)

Checks whether all variables are integral and bounded on both sides. If not an error is throws
"""
function check_var_bounds(model::Optimizer)
    for var in model.variable_info
        if !var.has_lower_bound || !var.has_upper_bound || !var.is_integer
            throw(ErrorException("Each variable must be an integer and bounded. Currently the variable index $(var.idx) doesn't fulfill this requirements."))
        end
    end
end

"""
    addupd_var_in_inner_model(model::Optimizer, vidx::Int)

Adds all variables to model.inner.search_space
"""
function addupd_var_in_inner_model(model::Optimizer, vidx::Int)
    if length(model.inner.search_space) < vidx
        push!(model.inner.search_space, model.variable_info[vidx])
    else
        model.inner.search_space[vidx] = model.variable_info[vidx]
    end
end

"""
    MOI.add_variable(model::Optimizer)
"""
function MOI.add_variable(model::Optimizer)
    vidx = length(model.variable_info) + 1
    push!(model.variable_info, Variable(vidx))
    changes = changes = Vector{Vector{Tuple{Symbol,Int,Int,Int}}}()
    push!(changes, Vector{Tuple{Symbol,Int,Int,Int}}())
    model.variable_info[vidx].changes = changes
    push!(model.inner.subscription, Int[])
    push!(model.inner.bt_infeasible, 0)
    push!(model.inner.var_in_obj, false)
    addupd_var_in_inner_model(model, vidx)
    return MOI.VariableIndex(vidx)
end

MOI.add_variables(model::Optimizer, n::Int) = [MOI.add_variable(model) for i in 1:n]

function check_inbounds(model::Optimizer, vi::VI)
    num_variables = length(model.variable_info)
    if !(1 <= vi.value <= num_variables)
        @error "Invalid variable index $vi. ($num_variables variables in the model.)"
    end
    return
end

check_inbounds(model::Optimizer, var::SVF) = check_inbounds(model, var.variable)

has_upper_bound(model::Optimizer, vi::VI) = model.variable_info[vi.value].has_upper_bound

has_lower_bound(model::Optimizer, vi::VI) = model.variable_info[vi.value].has_lower_bound

is_fixed(model::Optimizer, vi::VI) = model.variable_info[vi.value].is_fixed

function MOI.add_constraint(model::Optimizer, v::SVF, t::MOI.Integer)
    vi = v.variable
    model.variable_info[vi.value].is_integer = true

    cidx = length(model.var_constraints) + 1
    push!(model.var_constraints, (vi.value, :int, typemin(Int), typemax(Int)))

    addupd_var_in_inner_model(model, vi.value)
    return MOI.ConstraintIndex{SVF,MOI.Integer}(cidx)
end

function MOI.add_constraint(model::Optimizer, v::SVF, t::MOI.ZeroOne)
    vi = v.variable
    model.variable_info[vi.value].is_integer = true

    # this gets called after setting lower and upper bound
    # => make sure that it's not already set
    if !has_upper_bound(model, vi) || model.variable_info[vi.value].upper_bound > 1
        model.variable_info[vi.value].upper_bound = 1
        model.variable_info[vi.value].max = 1
        model.variable_info[vi.value].has_upper_bound = true
    end
    if !has_lower_bound(model, vi) || model.variable_info[vi.value].lower_bound < 0
        model.variable_info[vi.value].lower_bound = 0
        model.variable_info[vi.value].min = 0
        model.variable_info[vi.value].has_lower_bound = true
    end
    min_val, max_val = model.variable_info[vi.value].lower_bound,
                       model.variable_info[vi.value].upper_bound
    values = collect(min_val:max_val)
    model.variable_info[vi.value].values = values
    model.variable_info[vi.value].init_vals = copy(values)
    model.variable_info[vi.value].init_val_to_index = 1:length(values)
    model.variable_info[vi.value].offset = 1
    model.variable_info[vi.value].indices = 1:length(values)
    model.variable_info[vi.value].first_ptr = 1
    model.variable_info[vi.value].last_ptr = length(values)
    addupd_var_in_inner_model(model, vi.value)

    cidx = length(model.var_constraints) + 1
    push!(model.var_constraints, (vi.value, :bin, typemin(Int), typemax(Int)))
    return MOI.ConstraintIndex{SVF,MOI.ZeroOne}(cidx)
end

function MOI.add_constraint(model::Optimizer, v::SVF, t::Integers)
    vi = v.variable
    model.variable_info[vi.value].is_integer = true

    set_vals = t.values
    set_variable_from_integers!(model.variable_info[vi.value], set_vals)

    addupd_var_in_inner_model(model, vi.value)
    min_val = model.variable_info[vi.value].min
    max_val = model.variable_info[vi.value].max

    cidx = length(model.var_constraints) + 1
    # TODO: If we want to do something with var_constraints later we need to save the actual input values
    push!(model.var_constraints, (vi.value, :integers, min_val, max_val))
    return MOI.ConstraintIndex{SVF,Integers}(cidx)
end

#=
Populating Variable bounds
=#
function MOI.add_constraint(
    model::Optimizer,
    v::SVF,
    interval::MOI.Interval{T},
) where {T<:Real}
    vi = v.variable
    check_inbounds(model, vi)
    isnan(interval.upper) &&
        throw(ErrorException("The interval bounds can not contain NaN and must be an Integer. Currently it has an upper bound of $(interval.upper)"))
    has_upper_bound(model, vi) && @error "Upper bound on variable $vi exists already."
    isnan(interval.lower) &&
        throw(ErrorException("The interval bounds can not contain NaN and must be an Integer. Currently it has a lower bound of $(interval.lower)"))
    has_lower_bound(model, vi) && @error "Lower bound on variable $vi exists already."
    is_fixed(model, vi) && @error "Variable $vi is fixed. Cannot also set upper bound."

    model.variable_info[vi.value].upper_bound = interval.upper
    model.variable_info[vi.value].max = interval.upper
    model.variable_info[vi.value].has_upper_bound = true
    model.variable_info[vi.value].lower_bound = interval.lower
    model.variable_info[vi.value].offset = 1 - interval.lower
    model.variable_info[vi.value].min = interval.lower
    model.variable_info[vi.value].has_lower_bound = true

    model.variable_info[vi.value].values =
        (model.variable_info[vi.value].min):(model.variable_info[vi.value].max)
    num_vals = model.variable_info[vi.value].max - model.variable_info[vi.value].min + 1
    model.variable_info[vi.value].indices = 1:num_vals
    model.variable_info[vi.value].first_ptr = 1
    model.variable_info[vi.value].last_ptr = num_vals
    model.variable_info[vi.value].init_vals = copy(model.variable_info[vi.value].values)
    model.variable_info[vi.value].init_val_to_index =
        copy(model.variable_info[vi.value].indices)

    addupd_var_in_inner_model(model, vi.value)

    cidx = length(model.var_constraints) + 1
    push!(model.var_constraints, (vi.value, :interval, interval.lower, interval.upper))
    return MOI.ConstraintIndex{SVF,MOI.Interval{T}}(cidx)
end

function MOI.add_constraint(model::Optimizer, v::SVF, lt::MOI.LessThan{T}) where {T<:Real}
    vi = v.variable
    check_inbounds(model, vi)
    isnan(lt.upper) &&
        throw(ErrorException("The variable bounds can not contain NaN and must be an Integer. Currently it has an upper bound of $(lt.upper)"))
    has_upper_bound(model, vi) && @error "Upper bound on variable $vi already exists."
    is_fixed(model, vi) && @error "Variable $vi is fixed. Cannot also set upper bound."

    model.variable_info[vi.value].upper_bound = lt.upper
    model.variable_info[vi.value].max = lt.upper
    model.variable_info[vi.value].has_upper_bound = true

    if has_lower_bound(model, vi)
        model.variable_info[vi.value].values =
            (model.variable_info[vi.value].min):(model.variable_info[vi.value].max)
        num_vals = model.variable_info[vi.value].max - model.variable_info[vi.value].min + 1
        model.variable_info[vi.value].indices = 1:num_vals
        model.variable_info[vi.value].first_ptr = 1
        model.variable_info[vi.value].last_ptr = num_vals
        model.variable_info[vi.value].init_vals = copy(model.variable_info[vi.value].values)
        model.variable_info[vi.value].init_val_to_index =
            copy(model.variable_info[vi.value].indices)
    end

    addupd_var_in_inner_model(model, vi.value)

    cidx = length(model.var_constraints) + 1
    push!(model.var_constraints, (vi.value, :leq, typemin(Int), lt.upper))
    return MOI.ConstraintIndex{SVF,MOI.LessThan{T}}(cidx)
end

function MOI.add_constraint(
    model::Optimizer,
    v::SVF,
    gt::MOI.GreaterThan{T},
) where {T<:Real}
    vi = v.variable
    check_inbounds(model, vi)
    isnan(gt.lower) &&
        throw(ErrorException("The variable bounds can not contain NaN and must be an Integer. Currently it has an upper bound of $(gt.upper)"))
    has_lower_bound(model, vi) && @error "Lower bound on variable $vi already exists."
    is_fixed(model, vi) && @error "Variable $vi is fixed. Cannot also set lower bound."

    model.variable_info[vi.value].lower_bound = gt.lower
    model.variable_info[vi.value].min = gt.lower
    model.variable_info[vi.value].offset = 1 - gt.lower
    model.variable_info[vi.value].has_lower_bound = true

    if has_upper_bound(model, vi)
        model.variable_info[vi.value].values =
            (model.variable_info[vi.value].min):(model.variable_info[vi.value].max)
        num_vals = model.variable_info[vi.value].max - model.variable_info[vi.value].min + 1
        model.variable_info[vi.value].indices = 1:num_vals
        model.variable_info[vi.value].first_ptr = 1
        model.variable_info[vi.value].last_ptr = num_vals
        model.variable_info[vi.value].init_vals = copy(model.variable_info[vi.value].values)
        model.variable_info[vi.value].init_val_to_index =
            copy(model.variable_info[vi.value].indices)
    end
    addupd_var_in_inner_model(model, vi.value)

    cidx = length(model.var_constraints) + 1
    push!(model.var_constraints, (vi.value, :geq, gt.lower, typemax(Int)))
    return MOI.ConstraintIndex{SVF,MOI.GreaterThan{T}}(cidx)
end

function MOI.add_constraint(model::Optimizer, v::SVF, eq::MOI.EqualTo{T}) where {T<:Real}
    vi = v.variable
    check_inbounds(model, vi)
    isnan(eq.value) && @error "Invalid fixed value $(eq.value)."
    has_lower_bound(model, vi) && @error "Variable $vi has a lower bound. Cannot be fixed."
    has_upper_bound(model, vi) && @error "Variable $vi has an upper bound. Cannot be fixed."
    is_fixed(model, vi) && @error "Variable $vi is already fixed."
    model.variable_info[vi.value].lower_bound = eq.value
    model.variable_info[vi.value].upper_bound = eq.value
    model.variable_info[vi.value].is_fixed = true
    model.variable_info[vi.value].has_lower_bound = true
    model.variable_info[vi.value].has_upper_bound = true
    model.variable_info[vi.value].min = eq.value
    model.variable_info[vi.value].max = eq.value
    model.variable_info[vi.value].values = [eq.value]
    model.variable_info[vi.value].indices = [1]
    model.variable_info[vi.value].first_ptr = 1
    model.variable_info[vi.value].last_ptr = 1
    model.variable_info[vi.value].init_vals = copy(model.variable_info[vi.value].values)
    model.variable_info[vi.value].init_val_to_index =
        copy(model.variable_info[vi.value].indices)
    addupd_var_in_inner_model(model, vi.value)

    cidx = length(model.var_constraints) + 1
    push!(model.var_constraints, (vi.value, :eq, eq.value, eq.value))
    return MOI.ConstraintIndex{SVF,MOI.EqualTo{T}}(cidx)
end

function set_variable_from_integers!(var::Variable, set_vals)
    min_val, max_val = extrema(set_vals)
    range = max_val - min_val + 1
    vals = copy(set_vals)

    # fill the indices such that 1:length(set_vals) map to set_vals
    indices = zeros(Int, range)
    # .- needed for the offset
    # values that don't exist get an index out of the range of the values array
    indices[set_vals .- (min_val - 1)] = 1:length(set_vals)
    j = length(set_vals) + 1
    for i in 1:range
        if indices[i] == 0
            indices[i] = j
            j += 1
        end
    end

    var.upper_bound = max_val
    var.max = max_val
    var.has_upper_bound = true
    var.lower_bound = min_val
    var.min = min_val
    var.has_lower_bound = true
    var.values = vals
    var.init_vals = copy(vals)
    var.offset = 1 - min_val
    var.indices = indices
    # needs copy to be different
    var.init_val_to_index = copy(indices)
    var.first_ptr = 1
    var.last_ptr = length(set_vals)
end

"""
    set_init_fixes!(com::CS.CoM)

Redefines variables if they are fixed to a specific value.
Return feasibility
"""
function set_init_fixes!(com::CS.CoM)
    for fixes in com.init_fixes
        vidx = fixes[1]
        val = fixes[2]
        var = com.search_space[vidx]
        !fix!(com, var, val) && return false
    end
    return true
end


text:
###############################################################################
# Functions for calculating draft rules and estimating abilities
# for simulated and actual sports seasons
# Author: Evan Munro
# Date: January, 2020
###############################################################################

using Random, Combinatorics, Distributions, StatsBase, JuMP, Ipopt, MosekTools, SpecialFunctions, Suppressor

############################################################################
# Functions for calculating abilities
#############################################################################

#calculate updated abilities depending on if you win or lose the match
function contingent_abilities(matches,winners,losers,N,lambda)
    m = length(winners)
    t1,t2 = matches[m,:]
    winners[m] = t1
    losers[m] = t2
    ahatw = estimate_alpha(matches,winners,losers,N,lambda)
    winners[m] = t2
    losers[m] = t1
    ahatl = estimate_alpha(matches,winners,losers,N,lambda)
    winners[m] = 0
    losers[m] = 0
    return ahatw, ahatl
end

#MLE estimation of abilities based on records
function estimate_alpha(matches,wins,loses,N,lambda=10)
    model=Model(Ipopt.Optimizer)
    @variable(model,a[1:N],start=1/N)
    ncdf(x) = (1+erf(x/sqrt(2)))/2
    register(model,:ncdf,1,ncdf,autodiff=true)
    @constraint(model, sum(a)==1)
    @constraint(model, a .>=0)
    @NLobjective(model,Max,sum(log(1-ncdf((a[j]-a[i])/sqrt(2))) for (i,j) in zip(wins,loses)) -lambda*sum((a[i] - 1/N)^2 for i in 1:N))
    @suppress_out begin
        status = optimize!(model);
    end
    return value.(a)
end

#to make parallel
#export JULIA_NUM_THREADS=4

############################################################################
# Functions for calculating expected values of league objective
#############################################################################

function expected_objective_sim(s,abilities,calc_ability=false,mult=10,sims=1000)
    exp_obj = 0.0
    exp_stop = 0.0
    p = game_prob(abilities)
    M = size(s.matches,1)
    for i in 1:sims
        winners = zeros(Int64,M)
        for m in 1:M
            match = s.matches[m,:]
            t1,t2 = match
            weights = ProbabilityWeights([p[t1,t2], 1-p[t1,t2]])
            gwin = sample(match,weights)
            winners[m]=gwin
        end
        record = calculate_draft_rule(s, winners, calc_ability, abilities,mult)
        exp_obj += objective(record)
        exp_stop += record.stop
    end
    return exp_obj/sims, exp_stop/sims
end
function expected_objective(s,mult=10)

    M = size(s.matches,1)
    match_list = fill(Int[], M)
    for m in 1:M
        match_list[m] = s.matches[m,:]
    end
    all_outcomes = collect(Iterators.product(match_list...))
    exp_obj = 0
    exp_stop = 0
    p_total = 0
    for out=1:length(all_outcomes)
        winners = [all_outcomes[out]...]
        p_outcome = prob_winners(s, winners, s.abilities)
        record = calculate_draft_rule(s, winners, mult)
        exp_obj += p_outcome* objective(record)
        exp_stop += p_outcome*record.stop
        p_total += p_outcome
    end
    return exp_obj, exp_stop#, p_total
end

############################################################################
# Functions for calculating win and loss probabilities in single games
# and in seasons
#############################################################################

#probability of a certain history
function prob_winners(s, winners, abilities)
    prob = 1
    pwins = game_prob(abilities)
    for m in 1:size(s.matches,1)
        match = s.matches[m,:]
        w = winners[m]
        l = match[match.!=w][1]
        prob = prob*pwins[w,l]
    end
    return prob
end

#exact calculation of win probability based on enumerating all possible
#results from the season
function win_prob_ex(existing_wins,teams,matches,win_prob,top_k=1)

    all_outcomes = collect(Iterators.product(matches...))
    victory_prob = zeros(length(teams))
    for out =1:length(all_outcomes)
        prob = 1
        outcome = all_outcomes[out]
        wins=[count(x->x==i,outcome) for i in 1:length(teams)] + existing_wins
        for m =1:length(outcome)
            match = matches[m]
            winner = outcome[m]
            loser = match[findall(x->x!=winner,match)[1]]
            prob = prob*win_prob[winner,loser]
        end
        win_record = maximum(wins)
        winners = findall(x->x==win_record,wins)
        for j=1:length(winners)
            victory_prob[winners[j]] += prob/length(winners)
        end
    end
    return victory_prob
end

#approximate calculation of win probability based on simulation
function win_prob_sim(existing_wins, teams, matches, p, top_k, sims=1000)
    win_record = zeros(length(teams))
    bottom_record = zeros(length(teams))
    M = size(matches,1)
    for s in 1:sims
        wins = zeros(length(teams))
        wins += existing_wins
        for m in 1:M
            match =matches[m,:]
            t1,t2 = match
            weights = ProbabilityWeights([p[t1,t2], 1-p[t1,t2]])
            gwin = sample(match,weights)
            wins[gwin]+=1
        end
        cutoff = sort(wins,rev=true)[top_k]
        bottom = sort(wins)[1]
        losers = (wins .<= bottom)
        winners = (wins .>= cutoff)
        win_record = win_record .+ winners/(sum(winners)/top_k)
        bottom_record = bottom_record .+ losers/(sum(losers))
    end
    win_record = win_record./sims
    loss_record = bottom_record ./sims
    #re-normalize in case of ties
    return win_record, loss_record
end

function game_prob(a,sigma=1)
    d = Normal(0,1)
    N = length(a)
    p = zeros(N,N)
    for i in 1:N
        for j in 1:N
            p[i,j] = 1-cdf(d,(a[j]-a[i])/(sqrt(2)*sigma))
        end
    end
    return p
end

#############################################################################
# Structures to hold details and results for a season
#############################################################################

struct Season
    teams::Array{Int}
    matches::Array{Int,2}
    abilities::Array{Float64}
    pwin::Array{Float64,2}
end

function Season(N::Int, G::Int, abilities::Array{Float64})
    teams = Array(1:N)
    team_combos = collect(Iterators.flatten(combinations(teams, 2)))
    team_combos = transpose(reshape(team_combos, (2, Integer(length(team_combos)/2))))
    matches = repeat(team_combos, G)
    pwin = game_prob(abilities)
    return Season(teams, matches, abilities, pwin)
end

struct Record
    draft_prob::Array{Float64}
    stop::Int
    cwins::Array
    abilities::Array{Float64}
end

function summarize(r::Record)
    println("Final Draft probabilities: ")
    println(r.draft_prob[:,size(r.draft_prob,2)])
    println("Stopped at match ", r.stop)
    println("Final record: ")
    println(r.cwins[:,size(r.cwins,2)])
    println("Ability estimates: ")
    println(r.abilities[:,size(r.abilities,2)])
end

function objective(r::Record)
    d = draft(r)
    #adjust = d[r.cwins[:,size(r.cwins,2)].==maximum(r.cwins[:,size(r.cwins,2)])][1]
    obj = d[r.cwins[:,size(r.cwins,2)].==minimum(r.cwins[:,size(r.cwins,2)])][1]
    return obj #+ adjust/2
end

function draft(r::Record)
    return r.draft_prob[:,size(r.draft_prob,2)]
end

#function ability(r::Record)

function add_win(wins,team::Int)::Array{Int}
    new_wins = copy(wins)
    new_wins[team] +=1
    return new_wins
end

function get_losers(winners::Array, matches::Array{Int,2})::Array{Int}
    losers = zeros(Int,length(winners))
    for i in 1:length(winners)
        match = matches[i,:]
        w = winners[i]
        l = match[match.!=w][1]
    end
    return losers
end

#############################################################################
# Functions for calculating evolution of rule
##############################################################################

function calculate_draft_rule(s::Season, winners=nothing,
                              multiple=10, top_k = 1, seed=0)
    #Random.seed!(seed)
    #default is to not calculate ability updates
    adjust = true
    matches = s.matches
    teams = s.teams
    abilities = s.abilities
    M = size(matches,1)
    N = length(teams)

    if winners == nothing
        winners = zeros(Int64,M)
        losers = zeros(Int64,M)
    else
        losers = get_losers(winners,matches)
    end

    draft_prob = zeros(N,M+1) .+ 1/N
    cwins = zeros(Int64,N,M)

    stop_time = M
    for m in 1:M
        if m>1
            lag= m-1
        else
            lag = 1
        end
        match = matches[m,:]
        t1, t2 = match

        #if simulating season, then sample result
        if winners[m] == 0
            win_prob = game_prob(abilities)
            weights = ProbabilityWeights([win_prob[match...], 1-win_prob[match...]])
            w = sample(match,weights)
            winners[m] = w
        end

        if adjust
            probw1,probl1 = win_prob_sim(add_win(cwins[:,lag],t1),
                                        teams, matches[(m+1):M,:],
                                        s.pwin,top_k)
            probw2,probl2 = win_prob_sim(add_win(cwins[:,lag],t2),
                                        teams, matches[(m+1):M,:],
                                        s.pwin,top_k)

            diff1 = max(multiple*(probw1[t1]-probw2[t1]),0)
            diff2 = max(multiple*(probw2[t2]-probw1[t2]),0)

            if  (diff1 <  (probl2[t1]-probl1[t1])) || (diff2 < (probl1[t2]-probl2[t2]))
                adjust = false
                stop_time = m
            else
                if winners[m]==t1
                    draft_prob[:,(m+1):(M+1)] .= probl1
                else
                    draft_prob[:,(m+1):(M+1)] .= probl2
                end
            end
        end

        cwins[:,m] = cwins[:,lag]
        cwins[winners[m],m] += 1
        losers[m] = match[match.!=winners[m]][1]
    end
    return Record(draft_prob, stop_time, cwins, abilities)
end

##OLDWITHH ABILITIES
#estimate evolution of draft rule based on real data
function calculate_draft_rule(s::Season, winners, multiple,
                                top_k, calc_ability, abilities,
                                lambda=10, seed=0)
    #Random.seed!(seed)
    #default is to not calculate ability updates
    println(multiple)
    println(top_k)
    adjust = true
    matches = s.matches
    teams = s.teams
    M = size(matches,1)
    N = length(teams)
    ahat = zeros(N,M) .+ 1/N

    if abilities==nothing && winners == nothing
        error("To simulate a season, must provide abilities to simulate from")
    end

    if abilities == nothing
        abilities = zeros(N) .+ 1/N
    end

    if calc_ability==false
        ahat =  zeros(N,M) .+ abilities
        println("not calculating")
    end
    if winners == nothing
        winners = zeros(Int64,M)
        losers = zeros(Int64,M)
    else
        losers = get_losers(winners,matches)
    end

    draft_prob = zeros(N,M+1) .+ 1/N
    cwins = zeros(Int64,N,M)

    stop_time = M
    for m in 1:M
        if m>1
            lag= m-1
        else
            lag = 1
        end
        match = matches[m,:]
        t1, t2 = match
        if calc_ability
            ahat1,ahat2 = contingent_abilities(matches[1:m,:],
                                              copy(winners[1:m]),
                                              copy(losers[1:m]), N,lambda)
        else
            ahat1 = abilities
            ahat2 = abilities
        end
        #if simulating season, then sample result
        if winners[m] == 0
            win_prob = game_prob(abilities)
            weights = ProbabilityWeights([win_prob[match...], 1-win_prob[match...]])
            w = sample(match,weights)
            winners[m] = w
        end

        if adjust
            probw1,probl1 = win_prob_sim(add_win(cwins[:,lag],t1),
                                        teams, matches[(m+1):M,:],
                                        game_prob(ahat1),top_k)
            probw2,probl2 = win_prob_sim(add_win(cwins[:,lag],t2),
                                        teams, matches[(m+1):M,:],
                                        game_prob(ahat2),top_k)

            diff1 = max(multiple*(probw1[t1]-probw2[t1]),0)
            diff2 = max(multiple*(probw2[t2]-probw1[t2]),0)
            if t1==1
                incentives[m] = diff1/multiple
            elseif t2==1
                incentives[m] = diff2/multiple
            else
                incentives[m] = incentives[lag]
            end
            if  (diff1 < (probl2[t1]-probl1[t1])) || (diff2 <(probl1[t2]-probl2[t2]))
                adjust = false
                stop_time = m
            else
                if winners[m]==t1
                    draft_prob[:,(m+1):(M+1)] .= probl1
                else
                    draft_prob[:,(m+1):(M+1)] .= probl2
                end
            end
        end

        if calc_ability
            if winners[m]==t1
                ahat[:,m] = ahat1
            else
                ahat[:,m] = ahat2
            end
        end
        cwins[:,m] = cwins[:,lag]
        cwins[winners[m],m] += 1
        losers[m] = match[match.!=winners[m]][1]
    end
    return Record(draft_prob, stop_time, cwins, ahat)
end

#############################################################################
# Functions for calculating globally optimal rule
##############################################################################

function calc_objective(s, mult)
    println(mult)
    topk = 1
    N = length(s.teams)
    y = zeros(N) .+ 1/N
    M = size(s.matches, 1)
    ntd_val = zeros(2^M)
    ntd_wts = zeros(M)
    losers = ones(Int,2^M)
    adjust = true
    model = Model(Mosek.Optimizer)

    #draft rule for all partial histories
    @variable(model, 0<=x[i=0:M, j=1:2^i, k=1:N]<=1)

    #x = model[:x]
    #FAIR constraint
    for k in 1:N
        @constraint(model,x[0,1,k]==1/N)
    end

    #PROB constraint
    for i in 1:M
        for j in 1:2^i
            @constraint(model, sum(x[i,j,k] for k in 1:N)==1)
        end
    end

    #DC Constraint under equal ability
    for i in 0:(M-1)
        for j in 1:2^i
            for k in 1:N
                @constraint(model, x[i,j,k] == 1/2*x[i+1,2*j,k]+ 1/2*x[i+1,2*j-1,k])
            end
        end
    end

    history = zeros(N)
    step_forward(model, history, y, adjust, ntd_val, ntd_wts, losers,
                  s, 1, 1, mult, topk)

    @objective(model, Min, sum((1-x[M, j, losers[j]])^2 - (0-x[M, j, losers[j]])^2 +
                             sum((0-x[M, j, k])^2 for k in 1:N)  for j in 1:2^M)/2^M)
    #@NLobjective(model, Min, 1/2*sum(abs(1-x[M, j, losers[j]]) - abs(0-x[M, j, losers[j]]) +
    #                         sum(abs(0-x[M, j, k]) for k in 1:N)  for j in 1:2^M)/2^M)

    optimize!(model)
    println("here")
    opt_wts = zeros(M)
    for m in 1:M
        wt = 0
        for j in 1:2^(m-1)
            for n in 1:N
                wt += abs(value(x[m,2*j-1,n]) - value(x[m-1,j,n]))
                wt += abs(value(x[m,2*j,n]) -  value(x[m-1,j,n]))
            end
        end
        opt_wts[m] = wt/(N*2^m)
    end
    println("OBJECTIVE VALUE: ", objective_value(model))
    return objective_value(model),
            mean(ntd_val), opt_wts, ntd_wts

end

#adds NTD constraint and calculates draft rule as step forward through
#tree of all possible season outcomes
#only for equal abilities, for now
function step_forward(model, history, y, adjust, ntd_val, ntd_wts, losers, s, m, j, pi, topk)
    M = size(s.matches,1)
    t1,t2 = s.matches[m,:]
    probw1, probl1 = win_prob_sim(add_win(history,t1), s.teams,
                    s.matches[(m+1):M,:], s.pwin, topk)
    probw2, probl2 = win_prob_sim(add_win(history,t2), s.teams,
                    s.matches[(m+1):M,:], s.pwin, topk)

    diff1 = pi*(probw1[t1]-probw2[t1])
    diff2 = pi*(probw2[t2]-probw1[t2])

    x = model[:x]
    @constraint(model, x[m, 2*j, t1] - x[m, 2*j-1, t1] <= diff1)
    @constraint(model, x[m, 2*j-1, t2] - x[m, 2*j, t2] <= diff2)

    if !adjust || diff1 < (probl2[t1]-probl1[t1]) || diff2 < (probl1[t2]-probl2[t2])
        adjust = false
        y1 = copy(y)
        y2 = copy(y)
    else
        y1 = probl1
        y2 = probl2
    end

    ntd_wts[m] += (mean(abs.(y1.-y)) + mean(abs.(y2.-y)))/2^m

    if m < M
        step_forward(model, add_win(history, t1), y1, adjust,
                     ntd_val, ntd_wts, losers, s, m+1, 2*j-1, pi, topk)
        step_forward(model, add_win(history, t2), y2, adjust,
                     ntd_val, ntd_wts, losers, s, m+1, 2*j, pi, topk)
    else
        #which loser for the given history
        N = length(s.teams)
        l1 = argmin(add_win(history,t1))
        l2 = argmin(add_win(history,t2))
        losers[2*j-1] = l1
        losers[2*j] = l2
        ntd_val[2*j-1] = ((1-y1[l1])^2 - (0-y1[l1])^2 + sum( (0-y1[k])^2 for k in 1:N))
        ntd_val[2*j] = ((1-y2[l2])^2  - (0-y2[l2])^2 + sum((0-y2[k])^2 for k in 1:N))
    end
end

ext: jl

max_stars_repo_name: evanmunro/draft-policy

max_issues_repo_name: evanmunro/draft-policy

max_forks_repo_name: evanmunro/draft-policy


text:
using Test
using ImageCore
using IterTools
using ReferenceTests
using ImageDistances

# general distances should cover any combination of number_types and color_types unless it's special designed
include("testutils.jl")

include("hausdorff.jl")
include("metrics.jl")
include("ciede2000.jl")

nothing


text:
# unit test suite of JuliaAutonomy

module TestSuite
  using Test

  # test modules
  include(joinpath(split(@__FILE__, "runtests.jl")[1], "test_localization.jl"))
  include(joinpath(split(@__FILE__, "runtests.jl")[1], "test_common.jl"))
  include(joinpath(split(@__FILE__, "runtests.jl")[1], "test_model.jl"))
  include(joinpath(split(@__FILE__, "runtests.jl")[1], "test_prob_stats.jl"))
  include(joinpath(split(@__FILE__, "runtests.jl")[1], "test_slam.jl"))
  include(joinpath(split(@__FILE__, "runtests.jl")[1], "test_decision_making.jl"))
  include(joinpath(split(@__FILE__, "runtests.jl")[1], "test_inference.jl"))

  function main()
    TestCommon.main()
    TestLocalization.main()
    TestModel.main()
    TestProbStats.main()
    TestSlam.main()
    TestDecisionMaking.main()
    TestInference.main()
  end
end

if abspath(PROGRAM_FILE) == @__FILE__
  using .TestSuite
  TestSuite.main()
end

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy


text:
using Displaz

function example(N)
    t = linspace(0,1, N)'
    P = (10.+0.2*rand(size(t))) .* vcat(exp(5*t).*cos(100*pi*t),
                                        exp(5*t).*sin(100*pi*t),
                                        100*t);
    C = vcat(t, 1.-t, zeros(t));
    sz  = map(Float32, t)
    shape = mod(rand(UInt8, N), 6)
    plot3d!(P, color = C, label="Top series", markershape=shape)
    P[3,:] = -P[3,:]
    plot3d!(P, color = C, markersize = sz, markershape="-", label="Bottom series")
end

example(1_000_000)


text:
# R2_IRI_temporal_include.jl

include("../tcm_code/tcm_transactionsToCounts.jl")

include("../tcm_code/tcm_estimateMNL.jl")
include("../tcm_code/tcm_predictMNL.jl")

include("../tcm_code/tcm_estimateMMNL.jl")
include("../tcm_code/tcm_predictMMNL.jl")

include("../tcm_code/tcm_estimateHALOMNL.jl")
include("../tcm_code/tcm_predictHALOMNL.jl")

include("../tcm_code/tcm_estimateRanking.jl")
include("../tcm_code/tcm_predictRanking.jl")

include("../tcm_code/tcm_estimateForest.jl")

include("../tcm_code/tcm_predictForest.jl")

include("../tcm_code/tcm_evaluateKL.jl")

include("../tcm_code/tcm_leafStats.jl")
include("../tcm_code/tcm_rankingStats.jl")

include("R2_IRI_temporal_helperfn.jl")

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_licenses: MIT Unlicense


text:
### A Pluto.jl notebook ###
# v0.14.8

using Markdown
using InteractiveUtils

# ╔═╡ 0831131a-077d-11eb-2412-fbd0e7cc3436
using Pkg, DrWatson

# ╔═╡ 08314e98-077d-11eb-1d86-f9c9dd52e2c3
begin
	@quickactivate "StatisticalRethinkingTuring"
	using Turing
	using StatisticalRethinking
end

# ╔═╡ 041b4048-077d-11eb-20b8-6f8b8ee72626
md"## Clip-05-01-18t.jl"

# ╔═╡ bbbbd388-3ada-4108-bd5e-d9bc107f9d01
versioninfo()

# ╔═╡ 08320f9a-077d-11eb-2181-9bacc9603e02
begin
	df = CSV.read(sr_datadir("WaffleDivorce.csv"), DataFrame)
	df.D = zscore(df.Divorce)
	df.M = zscore(df.Marriage)
	df.A = zscore(df.MedianAgeMarriage)
end

# ╔═╡ 084594e8-077d-11eb-082b-6f19dc8a8e0f
std(df.MedianAgeMarriage)

# ╔═╡ 084663b4-077d-11eb-3586-73882ea2ceec
@model function m5_1_A(A, D)
    a ~ Normal(0, 0.2)
    bA ~ Normal(0, 0.5)
    σ ~ Exponential(1)
    μ = lin(a, A, bA)
    D ~ MvNormal(μ, σ)
end

# ╔═╡ 0857f156-077d-11eb-0829-799d64e773d2
begin
	m5_1_At = m5_1_A(df.A, df.D)
	prior5_1_At = sample(m5_1_At, Prior(), 50) |> DataFrame
	prior5_1_At = prior5_1_At[:, [:a, :bA, :σ]]
	Text(precis(prior5_1_At; io=String))
end

# ╔═╡ 35f34180-07d0-11eb-2c53-d1a2454235d5
begin
	x = -2:0.1:2
	plot(;dpi = 460)
	for r in eachrow(prior5_1_At)
		p = lin(r.a, x, r.bA)
		plot!(x, p, color = :black, alpha = 0.4)
	end
	plot!(legend = false)

	quap5_1_At = quap(m5_1_At)
	dfa5_1_At = DataFrame(rand(quap5_1_At.distr, 1000)', quap5_1_At.params)

	A_seq = range(-3, 3.2, length = 30)
	mu5_1_At = lin(dfa5_1_At.a', A_seq, dfa5_1_At.bA') |> meanlowerupper

	scatter(df.A, df.D, alpha = 0.4, legend = false)
	plot!(A_seq, mu5_1_At.mean, ribbon =
		(mu5_1_At.mean .- mu5_1_At.lower, mu5_1_At.upper .- mu5_1_At.mean);dpi = 460)
	vline!([0])
end

# ╔═╡ ffe4b816-07d5-11eb-351d-2fa77b983973
@model function m5_2_M(M, D)
    a ~ Normal(0, 0.2)
    bM ~ Normal(0, 0.5)
    σ ~ Exponential(1)
    μ = lin(a, M, bM)
    D ~ MvNormal(μ, σ)
end

# ╔═╡ 0d7bab3c-07d7-11eb-3406-f3eea1352ba4
begin
	m5_2_Mt = m5_2_M(df.M, df.D)
	quap5_2_Mt = quap(m5_2_Mt)
end

# ╔═╡ 0d7c05c8-07d7-11eb-2253-437844d8bb9b
begin
	dfa5_2_Mt = DataFrame(rand(quap5_2_Mt.distr, 1000)', quap5_2_Mt.params)
	M_seq = range(-3, 3.2, length = 30)
	mu5_2_Mt = lin(dfa5_2_Mt.a', M_seq, dfa5_2_Mt.bM') |> meanlowerupper
	scatter(df.M, df.D, alpha = 0.4, legend = false)
	plot!(M_seq, mu5_2_Mt.mean, ribbon =
		(mu5_2_Mt.mean .- mu5_2_Mt.lower, mu5_2_Mt.upper .- mu5_2_Mt.mean))
	vline!([0];dpi = 460)
end

# ╔═╡ 0db57a6a-07d7-11eb-2536-914facad3260
@model function m5_3_A_M(A, M, D)
    a ~ Normal(0, 0.2)
    bM ~ Normal(0, 0.5)
    bA ~ Normal(0, 0.5)
    σ ~ Exponential(1)
    μ = lin(a, A, bA, M, bM)
    D ~ MvNormal(μ, σ)
end

# ╔═╡ 0dc1c694-07d7-11eb-1b23-fff550d68b74
begin
	m5_3_A_Mt = m5_3_A_M(df.A, df.M, df.D)
	quap5_3_A_Mt = quap(m5_3_A_Mt)
end

# ╔═╡ fa977d88-08d8-11eb-1ba3-3dbba5377f58
begin
	dfa5_3_A_Mt = DataFrame(rand(quap5_3_A_Mt.distr, 1000)', quap5_3_A_Mt.params)
	mu5_3_A_Mt = lin(dfa5_3_A_Mt.a', M_seq, dfa5_3_A_Mt.bM') |> meanlowerupper
	scatter(df.M, df.D, alpha = 0.4, legend = false)
	plot!(M_seq, mu5_3_A_Mt.mean, ribbon =
		(mu5_3_A_Mt.mean .- mu5_3_A_Mt.lower, mu5_3_A_Mt.upper .- mu5_3_A_Mt.mean))
	vline!([0];dpi = 460)
end

# ╔═╡ 86a824ae-0991-11eb-3ee5-8ffb37cce3ec
begin
	quap5_1_3, fig = plotcoef(
		[m5_1_At, m5_2_Mt, m5_3_A_Mt],
		[:a, :bM, :bA, :σ];
		func=quap
	)
	plot(fig)
end

# ╔═╡ 15658b5c-0993-11eb-21a9-0fc78a4a6fdf
quap5_1_3

# ╔═╡ 11f58eb0-08da-11eb-0f7f-239cec1bfdd5
begin
	age = randn(50)
	mar = rand.(Normal.(-age))
	div = rand.(Normal.(age))
end;

# ╔═╡ 156073d8-08da-11eb-18d2-ef100b4f33b0
@model function m5_4_AM(A, M)
    a ~ Normal(0, 0.2)
    bAM ~ Normal(0, 0.5)
    σ ~ Exponential(1)
    μ = lin(a, A, bAM)
    M ~ MvNormal(μ, σ)
end

# ╔═╡ 1560bce4-08da-11eb-10b7-cf52ec6e7f33
begin
	sort!(df, :A)  # for nice looking graphs
	quap5_4_AMt = quap(m5_4_AM(df.A, df.M))
	post5_4_AMt = DataFrame(rand(quap5_4_AMt.distr, 1000)', quap5_4_AMt.params)
	mu5_4_AMt = lin(post5_4_AMt.a', df.A, post5_4_AMt.bAM') |> meanlowerupper
	resid = df.M .- mu5_4_AMt.mean

	scatter(df.A, df.M, legend = false)
	plot!(df.A, mu5_4_AMt.mean, ribbon =
		(mu5_4_AMt.mean .- mu5_4_AMt.lower, mu5_4_AMt.upper .- mu5_4_AMt.mean))

	post5_3_A_Mt = DataFrame(rand(quap5_3_A_Mt.distr, 1000)', quap5_3_A_Mt.params)
	mu5_3_A_M2t = lin(post5_3_A_Mt.a', df.A, post5_3_A_Mt.bA', df.M, post5_3_A_Mt.bM') |> meanlowerupper

	scatter(df.D, mu5_3_A_M2t.mean,
		yerror = (mu5_3_A_M2t.mean .- mu5_3_A_M2t.lower, mu5_3_A_M2t.upper .- mu5_3_A_M2t.mean),
		legend = false)
	plot!(identity, range(extrema(df.D)..., length = 10))
	plot!(xlabel = "Observed divorce", ylabel = "Predicted divorce")
end;

# ╔═╡ 701051c2-08da-11eb-33cb-41036a9f2724
begin
	xy = [df.D mu5_3_A_M2t.mean]
	for loc in ("ID", "UT", "ME", "RI")
		coord = xy[df.Loc .== loc, :]
		annotate!(coord..., loc, :bottom)
	end
	plot!(;dpi = 460)
end

# ╔═╡ c13c8b1e-08db-11eb-2725-b521e559429c
begin
	x_real = randn(100)
	x_spur = rand.(Normal.(x_real))
	y = rand.(Normal.(x_real))
	d = DataFrame((; y, x_real, x_spur))  # the semicolon turns it into a NamedTuple which then
										  # gives the DataFrame the names of the columns

	@df d corrplot([:y :x_real :x_spur])  # either
	corrplot(Matrix(d), label = names(d);dpi = 460) # or
end

# ╔═╡ 11d3317c-08df-11eb-1843-13d7b9b390ac
md"## End of clip-05-01-18t.jl"

# ╔═╡ Cell order:
# ╟─041b4048-077d-11eb-20b8-6f8b8ee72626
# ╠═bbbbd388-3ada-4108-bd5e-d9bc107f9d01
# ╠═0831131a-077d-11eb-2412-fbd0e7cc3436
# ╠═08314e98-077d-11eb-1d86-f9c9dd52e2c3
# ╠═08320f9a-077d-11eb-2181-9bacc9603e02
# ╠═084594e8-077d-11eb-082b-6f19dc8a8e0f
# ╠═084663b4-077d-11eb-3586-73882ea2ceec
# ╠═0857f156-077d-11eb-0829-799d64e773d2
# ╠═35f34180-07d0-11eb-2c53-d1a2454235d5
# ╠═ffe4b816-07d5-11eb-351d-2fa77b983973
# ╠═0d7bab3c-07d7-11eb-3406-f3eea1352ba4
# ╠═0d7c05c8-07d7-11eb-2253-437844d8bb9b
# ╠═0db57a6a-07d7-11eb-2536-914facad3260
# ╠═0dc1c694-07d7-11eb-1b23-fff550d68b74
# ╠═fa977d88-08d8-11eb-1ba3-3dbba5377f58
# ╠═86a824ae-0991-11eb-3ee5-8ffb37cce3ec
# ╠═15658b5c-0993-11eb-21a9-0fc78a4a6fdf
# ╠═11f58eb0-08da-11eb-0f7f-239cec1bfdd5
# ╠═156073d8-08da-11eb-18d2-ef100b4f33b0
# ╠═1560bce4-08da-11eb-10b7-cf52ec6e7f33
# ╠═701051c2-08da-11eb-33cb-41036a9f2724
# ╠═c13c8b1e-08db-11eb-2725-b521e559429c
# ╟─11d3317c-08df-11eb-1843-13d7b9b390ac


text:
using Test
using DandelionWebSockets: CloseRequest, SocketClosed, protocolstate
using DandelionWebSockets: PongMissed, STATE_CLOSING

@testset "Connection management  " begin
    @testset "the server initiates a closing handshake" begin
        @testset "state is CLOSING" begin
            logic, handler, writer = makeclientlogic()

            close_frame = closeframe_from_server()
            handle(logic, FrameFromServer(close_frame))

            @test protocolstate(logic) == STATE_CLOSING
        end

        @testset "handler is notified of state change" begin
            logic, handler, writer = makeclientlogic()

            close_frame = closeframe_from_server()
            handle(logic, FrameFromServer(close_frame))

            @test handler.state == STATE_CLOSING
        end

        @testset "a close frame is sent in reply" begin
            logic, handler, writer = makeclientlogic()

            close_frame = closeframe_from_server()
            handle(logic, FrameFromServer(close_frame))

            frame = getframe(writer, 1)
            @test frame.opcode == OPCODE_CLOSE
        end
    end

    @testset "the client initiates a closing handshake" begin
        @testset "state is CLOSING" begin
            # Requirement
            # @7_1_2-1 Start the closing handshake

            logic, handler, writer = makeclientlogic()

            handle(logic, CloseRequest())

            @test protocolstate(logic) == STATE_CLOSING
        end

        @testset "handler is notified of state change" begin
            logic, handler, writer = makeclientlogic()

            handle(logic, CloseRequest())

            @test handler.state == STATE_CLOSING
        end

        @testset "a close frame is sent" begin
            logic, handler, writer = makeclientlogic()

            handle(logic, CloseRequest())

            frame = getframe(writer, 1)
            @test frame.opcode == OPCODE_CLOSE
        end
    end

    @testset "the server replies to a client initiated closing handshake" begin
        @testset "state is CLOSING" begin
            logic, handler, writer = makeclientlogic()

            handle(logic, CloseRequest())
            close_frame = closeframe_from_server()
            handle(logic, FrameFromServer(close_frame))

            @test protocolstate(logic) == STATE_CLOSING
        end
    end

    @testset "the socket is closed cleanly" begin
        @testset "the state is CLOSED" begin
            logic, handler, writer = makeclientlogic()

            handle(logic, CloseRequest())
            close_frame = closeframe_from_server()
            handle(logic, FrameFromServer(close_frame))
            handle(logic, SocketClosed())

            @test protocolstate(logic) == STATE_CLOSED
        end

        @testset "the handler is notified of the state change" begin
            logic, handler, writer = makeclientlogic()

            handle(logic, CloseRequest())
            close_frame = closeframe_from_server()
            handle(logic, FrameFromServer(close_frame))
            handle(logic, SocketClosed())

            @test handler.state == STATE_CLOSED
        end

        @testset "the client cleanup function is called" begin
            was_client_cleanup_called = false
            client_cleanup = () -> was_client_cleanup_called = true
            # TODO Rewrite to put in correct closing state
            logic, handler, writer = makeclientlogic(client_cleanup=client_cleanup)

            handle(logic, CloseRequest())
            close_frame = closeframe_from_server()
            handle(logic, FrameFromServer(close_frame))
            handle(logic, SocketClosed())

            @test was_client_cleanup_called == true
        end
    end

    @testset "close connection if enough pongs have been missed" begin
        @testset "state is closed" begin
            logic, handler, writer = makeclientlogic()

            handle(logic, PongMissed())

            @test protocolstate(logic) == STATE_CLOSED
        end

        @testset "handler is notified of the state change" begin
            logic, handler, writer = makeclientlogic()

            handle(logic, PongMissed())

            @test handler.state == STATE_CLOSED
        end
    end

    @testset "User request to close the connection during a close; Another Close frame is not sent" begin
        # Arrange
        logic, handler, writer = makeclientlogic()
        close_frame = closeframe_from_server()
        handle(logic, FrameFromServer(close_frame))
        # A Close frame is sent here.

        # Act
        handle(logic, CloseRequest())
        # A Close frame is _not_ sent here.

        # Assert
        # Only one close frame should have been sent.
        @test get_no_of_frames_written(writer) == 1
    end

    @testset "Socket unexpectedly closed; Connection is failed" begin
        logic, handler, writer = makeclientlogic()

        handle(logic, SocketClosed())

        @test protocolstate(logic) == STATE_CLOSED
    end

    @testset "Socket closed after client close request; handler is not notified twice" begin
        logic, handler, writer = makeclientlogic()

        handle(logic, CloseRequest())
        close_frame = closeframe_from_server()
        handle(logic, FrameFromServer(close_frame))
        handle(logic, SocketClosed())

        @test handler.statesequence == [STATE_CLOSING, STATE_CLOSED]
    end
end


text:
# Used "Conjugate Bayesian analysis of the Gaussian distribution" by Murphy as
# a reference.  Note that there were some typos in that document so the code
# here may not correspond exactly.

struct NormalWishart{T<:Real,V<:AbstractVector{T},M<:AbstractMatrix{T}} <: ContinuousMultivariateDistribution
    dim::Int
    zeromean::Bool
    mu::V
    kappa::T
    Tchol::Cholesky{T,M}  # Precision matrix (well, sqrt of one)
    nu::T

    function NormalWishart{T}(mu::AbstractVector{T}, kappa::T,
                              Tchol::Cholesky{T,M}, nu::T) where {T<:Real, M<:AbstractMatrix{T}}
        # Probably should put some error checking in here
        d = length(mu)
        zmean::Bool = true
        for i = 1:d
            if !iszero(mu[i])
                zmean = false
                break
            end
        end
        new{T,typeof(mu),M}(d, zmean, mu, T(kappa), Tchol, T(nu))
    end
end

function NormalWishart(mu::AbstractVector{U}, kappa::Real,
                       Tchol::Cholesky{S}, nu::Real) where {S<:Real, U<:Real}
    T = promote_type(U,S,typeof(kappa), typeof(nu))
    return NormalWishart{T}(convert(AbstractVector{T}, mu),T(kappa),Cholesky{T}(Tchol), T(nu))
end

function NormalWishart(mu::AbstractVector{T}, kappa::T,
                       Tmat::AbstractMatrix{T}, nu::T) where T<:Real
    NormalWishart{T}(mu, kappa, cholesky(Tmat), nu)
end

function insupport(::Type{NormalWishart}, x::Vector{T}, Lam::Matrix{T}) where T<:Real
    return (all(isfinite, x) &&
           size(Lam, 1) == size(Lam, 2) &&
           isApproxSymmmetric(Lam) &&
           size(Lam, 1) == length(x) &&
           hasCholesky(Lam))
end

pdf(nw::NormalWishart, x::Vector{T}, Lam::Matrix{S}) where T<:Real where S<:Real =
        exp(logpdf(nw, x, Lam))

function logpdf(nw::NormalWishart, x::Vector{T}, Lam::Matrix{T}) where T<:Real
    if !insupport(NormalWishart, x, Lam)
        return -Inf
    else
        p = length(x)

        nu = nw.nu
        kappa = nw.kappa
        mu = nw.mu
        Tchol = nw.Tchol
        hnu = 0.5 * nu
        hp = 0.5 * p

        # Normalization
        logp = hp*(log(kappa) - Float64(log2π))
        logp -= hnu * logdet(Tchol)
        logp -= hnu * p * log(2.)
        logp -= logmvgamma(p, hnu)

        # Wishart (MvNormal contributes 0.5 as well)
        logp += (hnu - hp) * logdet(Lam)
        logp -= 0.5 * tr(Tchol \ Lam)

        # Normal
        z = nw.zeromean ? x : x - mu
        logp -= 0.5 * kappa * dot(z, Lam * z)

        return logp

    end
end

function rand(nw::NormalWishart)
    Lam = rand(Wishart(nw.nu, nw.Tchol))
    Lsym = PDMat(Symmetric(inv(Lam) ./ nw.kappa))
    mu = rand(MvNormal(nw.mu, Lsym))
    return (mu, Lam)
end


text:
# # Data Structures
#
#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__notebooks/02_DataStructures.ipynb)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__notebooks/02_DataStructures.ipynb)
#
# ## Arrays
#
# Julia has a nice and flexible array interface. Arrays can have an arbitrary
# number of dimensions. Let's define a one-dimetional array (i.e. a vector):

vector = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0]

# The first index of an array in Julia is `1`:

vector[1]

# You can use `end` to access the last element of an array:

vector[end]

# Use ranges (`start:end`) to get a slice of the array:

vector[2:4]

# Ranges in Julia are iterable objects:

indexes = 2:4

#-

for i in indexes
	@show i
end

# Julia arrays, like the strings and ranges, are also iterables:

for element in vector
	println(element)
end

# #### Exercise 1
#
# Write a function to return the distance between two three dimensional points,
# i.e. two vector of three elements. You should use a `for` loop over a range
# and index the vectors.

## function distance(a, b...

#-

using Test
A = [1.25, 2.0, 3.6]
B = [-3.5, 4.7, 5.0]
@test distance(A, B) ≈ hypot((A - B)...)

# `...` "splats" the values contained in an iterable collection into a
# function call as individual arguments, e.g:

vector = [1, 2, 3]
hypot(vector...)  # hypot(1, 2, 3)

# You can use `push!` to add one element to the end of an array

vector = [1,2,3]

#-
push!(vector, 4)

# There are other useful
# [dequeues functions](https://docs.julialang.org/en/v1/base/collections/#Dequeues-1)
# defined in Julia, e.g. `pop!`, `append!`.
#
# In Julia, by convention, all the functions that modify their arguments
# should end with a bang or exclamation mark, `!`, see the
# [style guide](https://docs.julialang.org/en/v1/manual/style-guide/index.html#Append-!-to-names-of-functions-that-modify-their-arguments-1).

# ### Vectorized operations
#
# You can use a dot, `.`,  to indicate that a function, e.g. `log.(x)`, or
# operator, e.g. `x .^ y`,  should be applied element by element, see
# [dot syntax](https://docs.julialang.org/en/v1/manual/functions/#man-vectorized-1):

a = [1, -2, -3]
b = [-2, -4, 0]
a .* b

# This notation allows vectorizing any function, even element-wise functions
# defined by the user:

f(x) = 3.45x + 4.76

f.(sin.(a))

# Multiple vectorized operations get
# [fused in a single loop](https://julialang.org/blog/2017/01/moredots)
# without temporal arrays.

# ### Comprehensions
#
# You can use [comprehensions](https://docs.julialang.org/en/v1/manual/arrays/#Comprehensions-1)
# to create arrays and perform some operation

[ 2x for x in 1:10 ]

#-

result = [ 2x for x in 1:10 if x % 2 == 0 ]

# #### Exercise 3
#
# Write the equivalent of the previous expression using a `for` loop and `push!`.

## result = []
## for ...

# ### Matrices
#
# Matrices, bidimentional arrays, can be defined with the following notation:

matrix = [ 1.0 4.0 7.0
           2.0 5.0 8.0
		   3.0 6.0 9.0 ]

# You can use linear indexing (Julia arrays are stored in column major order)
# to access an element

matrix[2]

# Or using one index by dimension, i.e. `matrix[row_index, col_index]` :

matrix[2, 1]

# You can also use ranges and `end`. The colon, `:`, means that all the indices
# from that dimension should be used:

matrix[2:end, :]

# #### Comprehensions
#
# You can also use comprehensions to create matrices. In fact, you can create
# array of any desired dimension:

[ x + y for x in 1:5, y in 1:10 ]

# ## Dictionaries and pairs
#
# Dictionaries (hash tables) stores key => values pairs:

dictionary = Dict('A' => 'T', 'C' => 'G', 'T' => 'A', 'G' => 'C')

# You can get a value by indexing with the key:

dictionary['A'] # get(dictionary, 'A')

# If the key is not present in the dictionary, an error is raised:

dictionary['N'] # get(dictionary, 'N')

# The function `get` allows to specify a default value that is returned if
# the key is absent in the dictionary:

get(dictionary, 'N', '-')

# A nice thing about hash tables (dictionary keys, sets) is that test
# membership is $O(1)$ while it is $O(N)$ in lists/vectors/arrays:

'N' in keys(dictionary)

# A dictionary gives pairs when it is iterated:

for pair in dictionary
	println("pair: ", pair)  # each pair is key => value
	println("key: ", pair.first)  # pair.first == pair[1]
	println("value: ", pair.second)  # pair.second == pair[2]
end

# ## Tuples
#
# Tuples are immutable collections, while arrays are mutable:

point = [1.0, 2.0, 3.0]  # vector
point[1] = 10.0
point

#-

point = (1.0, 2.0, 3.0)  # tuple

#-

point[1] = 10.0

# You can index a tuple, like a vector, to get the stored element(s):

point[1:2]

# Tuples, vectors, pairs and other iterables can be easily unpacked using an
# assignation:

x, y, z = point
y

# You can use this unpacking when iterating a dictionary:

for (key, value) in dictionary
	println("key: ", key, " value: ", value)
end

# #### Exercise 3
#
# Write a function to return the reverse complement of a DNA sequence (string)
# using a dictionary, the `join` function and the `Base.Iterators.reverse`
# iterator. It should use a 'N' as complementary of any base different from
# 'A', 'C', 'T' or 'G':

## function reverse_complement(...

#-

using Test
@test reverse_complement("ACTGGTCCCNT") == "ANGGGACCAGT"

# ### Named tuples
#
# They can be an easy and fast way to store data:

point = (x=1.0, y=2.0, z=3.0)  # named tuple

# You can use `namedtuple.name` to access a particular element:

point.y

# ## Sets
#
# You can use `Set` to represent a set of unique elements:

set = Set([1, 2, 2, 3, 3, 3, 4, 4, 4, 4])

# Test membership is $O(1)$

4 in set

# You can get the intersection of two sets using `intersect` or
# `∩` (`\cap<TAB>`)

set_a = Set([1, 2, 3])
set_b = Set([2, 3, 4])
set_a ∩ set_b  # intersect(set, set_b)

# And the unioin of to sets using `union` or `∪` (`\cup<TAB>`)

set_a ∪ set_b  # union(set, set_b)

# The symmetric difference, i.e. disjunctive union, of two sets

symdiff(set_a, set_b)

max_stars_repo_name: diegozea/JuliaForBioinformatics

max_issues_repo_name: diegozea/JuliaForBioinformatics

max_forks_repo_name: diegozea/JuliaForBioinformatics


text:
module Dimensionless

using Unitful, LinearAlgebra

export DimBasis, dim_matrix
export number_of_dimensions, number_of_dimensionless, dimensionless, dimensionful, change_basis
export print_dimensionless

include("dim_basis.jl")
include("dim_basis_ops.jl")
include("utils.jl")

end # module


text:
using AnalogDiscoveryNeuromore

const fs = 512.0
const runForMin = 180
const twoChannels = true
const lowVolt = false

go(fs,runForMin,twoChannels,lowVolt)


text:
abstract type ChromoType end
abstract type Chromosome end
abstract type FeedForward <: ChromoType end

mutable struct MLPChromosome <: Chromosome
    id::Int64
    input_size::Int64
    output_size::Int64
    n_hidden_layers::Int64
    layers::Vector{Layer}
    fitness::Float64
    species_id::Int64
    parent1_id::Int64
    parent2_id::Int64
    genes_mapping::Dict{Int64, Tuple{Int64, Int64}}
    function MLPChromosome(g::Global)
        @assert g.cf.initial_hidden_nodes >= 1
        ntype = :MLP
        layers = Layer[MLP_layer(g, :HIDDEN, 1, g.cf.nn_activation, g.cf.input_nodes, g.cf.initial_hidden_nodes),
                       MLP_layer(g, :OUTPUT, 2, :identity, g.cf.initial_hidden_nodes, g.cf.output_nodes)]
        genes_mapping = merge(Dict([]), get_layer_dictionary(layers[1]))
        genes_mapping = merge!(genes_mapping, get_layer_dictionary(layers[2]))
        new(incChromeId!(g),
            g.cf.input_nodes, g.cf.output_nodes,
            1,
            layers, # array of layers
            0.0, # stub for fitness function
            0, # species_id
            # parents ids: they help in tracking chromosome's genealogy
            0, 0,
            genes_mapping
            )
    end
    function MLPChromosome(g::Global, parent1::MLPChromosome, parent2::MLPChromosome)
        @assert parent1.species_id == parent2.species_id
        new(incChromeId!(g),
            g.cf.input_nodes, g.cf.output_nodes,
            1,
            Layer[], # array of layers
            0.0, # stub for fitness function
            parent1.species_id, # species_id
            # parents ids: they help in tracking chromosome's genealogy
            parent1.id, parent2.id,
            Dict([])
            )
    end
end

function Base.show(io::IO, ch::MLPChromosome)
    println(io, "id: $(ch.id)")
    println(io, "nº hidden layers: $(ch.n_hidden_layers)")
    nlayers = length(ch.layers)
    println(io, "Layer 0 (INPUT) -> nodes: $(ch.input_size)")
    topology_str = "$(ch.input_size)"
    for l = 1:nlayers
        layer = ch.layers[l]
        ninputs = layer.nodes[1].ninputs
        nnodes = count((n) -> n.enable_bit, layer.nodes)
        l == nlayers ? print(io, "Layer $l (OUTPUT) -> ") : print(io, "Layer $l -> ")
        print(io, "id: $(layer.layer_id), position: $(layer.position), ")
        print(io, "nodes: $(layer.nodes_number)")
        #print(", weights: $(nnodes)x$(ninputs)")
        topology_str = string(topology_str, " -> ", nnodes)
        #if l != nlayers  println() end
        print(io, "\n")
    end
    print(io, "Topology: $topology_str")
end

function get_parameters(ch::MLPChromosome)
    number_layers = length(ch.layers)
    layers_params = Vector{Tuple{Matrix{Float64}, Vector{Float64}, Symbol, Int64}}(undef, number_layers)
    previous_enabled_nodes = ones(Bool, ch.layers[1].nodes[1].ninputs)
    for l = 1:number_layers
        #=if l == 1
            enabled_nodes = ones(Bool, ch.layers[l].nodes[1].ninputs)
            layers_params[l] = get_layer_params(ch.layers[l], enabled_nodes)
        else

        end=#
        layers_params[l] = get_layer_params(ch.layers[l], previous_enabled_nodes)
        previous_enabled_nodes = map((n) -> n.enable_bit, ch.layers[l].nodes)
    end
    #[get_layer_params(l) for l in ch.layers]
    layers_params
end

get_layers(ch::MLPChromosome) = ch.layers

function incChromeId!(g::Global)
    g.chromosomeCnt += 1
    return g.chromosomeCnt
end

create_minimal_chromossome(g::Global)::MLPChromosome = MLPChromosome(g)

function number_nodes(ch::MLPChromosome)
    if length(ch.layers) > 0
        return sum(map((l) -> length(l.nodes), ch.layers))
    end
    return 0
end

function maxInnov(ch::MLPChromosome)
    if length(ch.layers) > 0
        return map((l) -> maxInnov_layer(l), ch.layers) |> maximum
    end
    error("There is no layers in the chromossome.")
end

function mutate_add_node!(ch::MLPChromosome, g::Global)
    layer_pos = rand(1:ch.n_hidden_layers)
    number_inputs = ch.layers[layer_pos].nodes[1].ninputs
    new_node = NodeGene(get_new_innov_number!(g), :HIDDEN, number_inputs, randn(), g.cf.nn_activation)
    push!(ch.layers[layer_pos].nodes, new_node)
    add_node_weight!(ch.layers[layer_pos+1])
    ch.layers[layer_pos].nodes_number += 1
    merge!(ch.genes_mapping, Dict(new_node.id => (layer_pos, ch.layers[layer_pos].nodes_number)))
    return
end

function mutate_enable_bit!(ch::MLPChromosome)
    layer_pos = rand(1:ch.n_hidden_layers)
    l = ch.layers[layer_pos]
    nodes = l.nodes
    len = length(nodes)
    if len > 1
        enabled_nodes = get_enabled_nodes(l)
        len_enabled = length(enabled_nodes)
        if len_enabled <= 1 && len > len_enabled
            not_enabled = get_not_enabled_nodes(l)
            node_index = rand(1:length(not_enabled))
            mutate_enable_bit!(not_enabled[node_index])
        elseif len_enabled > 1
            node_index = rand(1:len)
            mutate_enable_bit!(nodes[node_index])
        end
    end
end

function mutate_remove_node!(ch::MLPChromosome)
    layer_pos = rand(1:ch.n_hidden_layers)
    nodes = ch.layers[layer_pos].nodes
    len = length(nodes)
    if  len > 1
        node_pos = rand(1:len)
        delete!(ch.genes_mapping, nodes[node_pos].id)
        deleteat!(nodes, node_pos)
        remove_node_weight!(ch.layers[layer_pos+1], node_pos)
        ch.layers[layer_pos].nodes_number -= 1
        merge!(ch.genes_mapping, get_layer_dictionary(ch.layers[layer_pos]))
    end
    return
end

function mutant_layer(g::Global, position::Int, ninputs::Int, nodes_number::Int, activation::Symbol)
    nodes = Vector{Node}(undef, nodes_number)
    for i = 1:nodes_number
        #=weights = zeros(Float64, ninputs)
        if i <= ninputs
            weights[i] = 1.0
        end=#
        weights = randn(Float64, ninputs)
        nodes[i] = NodeGene(get_new_innov_number!(g), :HIDDEN, weights, 0.0, activation)
    end
    MLP_layer(g, :HIDDEN, position, activation, nodes)
end

function mutate_add_layer!(ch::MLPChromosome, g::Global)
    #number_layers = length(get_layers(ch))
    number_layers = ch.n_hidden_layers
    new_layer_pos = rand(1:number_layers)
    number_nodes_nl = (Int ∘ floor)((ch.layers[new_layer_pos].nodes_number + ch.layers[new_layer_pos+1].nodes_number)/2)
    #new_layer = MLP_layer(g, :MLP, new_layer_pos+1, :tanh, ch.layers[new_layer_pos].nodes_number, number_nodes_nl)
    new_layer = mutant_layer(g, new_layer_pos+1, ch.layers[new_layer_pos].nodes_number, number_nodes_nl, g.cf.nn_activation)

    for n in ch.layers[new_layer_pos+1].nodes
        old_ninputs = n.ninputs
        n.ninputs = number_nodes_nl
        if number_nodes_nl <= old_ninputs
            n.weights = n.weights[1:number_nodes_nl]
        elseif number_nodes_nl > old_ninputs
            old_weights = n.weights
            new_weights = zeros(Float64, number_nodes_nl)
            new_weights[1:old_ninputs] = old_weights
            n.weights = new_weights
        end
    end
    insert!(ch.layers, new_layer_pos+1, new_layer)
    for i = (new_layer_pos+1):length(ch.layers)
        ch.layers[i].position = i
        merge!(ch.genes_mapping, get_layer_dictionary(ch.layers[i]))
    end
    ch.n_hidden_layers += 1
    return
end

function mutate_weights!(ch::MLPChromosome, g::Global)
    for l in ch.layers
        mutate_node_weight!(g, l)
    end
    return
end

function mutate!(ch::MLPChromosome, g::Global)
    # Mutates the chromosome
    if rand() < g.cf.prob_structural_mutation
        if rand() < g.cf.prob_addnode
            mutate_add_node!(ch, g)
            #println("Add node mutation")
        #=elseif rand() < g.cf.prob_removenode
            mutate_remove_node!(ch)
            #println("Remove node mutation")
        =#
        elseif rand() < g.cf.prob_addlayer
            mutate_add_layer!(ch, g)
            #println("Add layer mutation")
        end

    elseif rand() < g.cf.prob_enable_mutation
        mutate_enable_bit!(ch) # mutate node's enable bit
    else
        mutate_weights!(ch, g) # mutate node's weights and bias
    end
    return ch
end

function convert_to_FluxNet(ch::MLPChromosome, labels::Vector{Float64})::MLPClassifier
    params = get_parameters(ch)
    weights = getindex.(params, 1)
    biases = getindex.(params, 2)
    activations = getindex.(params, 3)
    number_nodes_layers = getindex.(params, 4)
    MLPClassifier(ch.input_size, labels, collect(zip(weights,biases)), activations, number_nodes_layers[1:end-1])
end

function get_node_from_dict(ch::MLPChromosome, key::Int)
    layer, node_pos = ch.genes_mapping[key]
    ch.layers[layer].nodes[node_pos]
end

function nodes_innov_numbers(ch::MLPChromosome)
    innov_numbers = Int[]
    for l in ch.layers
        append!(innov_numbers, map((n) -> n.id, l.nodes))
    end
    innov_numbers
end

# compatibility function
function distance(self::MLPChromosome, other::MLPChromosome, cf::Config)

    # Returns the distance between this chromosome and the other.
    chromo1, chromo2 = number_nodes(self) > number_nodes(other) ? (self,other) : (other,self)

    weight_diff = 0
    matching = 0
    disjoint = 0
    excess = 0

    #biggest_genome = max(length(chromo1.genes_mapping), length(chromo2.genes_mapping))

    max_innov_chromo2 = maxInnov(chromo2)
    max_innov_chromo1 = maxInnov(chromo1)

    max_id = -1

    for (k, g1) in chromo1.genes_mapping
        node1 = get_node_from_dict(chromo1, k)
        if haskey(chromo2.genes_mapping, k)
            # Homologous genes
            node2 = get_node_from_dict(chromo2, k)

            max_id = max(max_id, k)

            node1, node2 = length(node1.weights) > length(node2.weights) ? (node1, node2) : (node2, node1)
            max_len, min_len = length(node1.weights), length(node2.weights)
            for i = 1:max_len
                if i <= min_len
                    weight_diff += abs(node1.weights[i] - node2.weights[i])
                else
                    weight_diff += abs(node1.weights[i])
                end
            end
            matching += 1
        end
    end

    all_innov_numbers = union(nodes_innov_numbers(chromo1), nodes_innov_numbers(chromo2))

    disjoint = sum(map((n) -> n <= max_id, all_innov_numbers)) - matching
    excess = sum(map((n) -> n > max_id, all_innov_numbers))

    d = (cf.excess_coeficient * excess + cf.disjoint_coeficient * disjoint)#/biggest_genome
    #println("match: $matching, exc: $excess, disj: $disjoint e weight diff: $weight_diff.")

    return matching > 0 ? d + cf.weight_coeficient * weight_diff / matching : d
end

#=
function crossover(g::Global, self::MLPChromosome, other::MLPChromosome)
    # Crosses over parents' chromosomes and returns a child

    # This can't happen! Parents must belong to the same species.
    @assert self.species_id == other.species_id

    # If they're of equal fitnesses, choose the shortest
    parent1, parent2 = self.fitness >= other.fitness ? (self,other) : (other,self)
    if parent1.fitness == parent2.fitness
        parent1, parent2 = number_nodes(parent1) >= number_nodes(parent2) ?
                           (parent1, parent2) : (parent2, parent1)
    end

    # create a new child
    #child = MLPChromosome(incChromeId!(g), other.id, self.node_gene_type, self.conn_gene_type)
    child = deepcopy(parent1)
    inherit_genes!(g, child, parent1, parent2)
    child.species_id = parent1.species_id
    child.fitness = 0 #parent1.fitness
    #child._input_nodes = parent1._input_nodes

    return child
end
=#

function crossover(g::Global, self::MLPChromosome, other::MLPChromosome)
    # Crosses over parents' chromosomes and returns a child

    # This can't happen! Parents must belong to the same species.
    @assert self.species_id == other.species_id

    local child
    if self.fitness == other.fitness
        child = generate_and_inherit_genes(g, self, other)
    else
        parent1, parent2 = self.fitness > other.fitness ? (self,other) : (other,self)
        child = MLPChromosome(g, parent1, parent2)
        child.layers = deepcopy(parent1.layers)
        child.n_hidden_layers = length(parent1.layers)-1
        child.genes_mapping = deepcopy(parent1.genes_mapping)
        inherit_genes!(child, g, parent1, parent2)
    end
    child.fitness = 0

    return child
end

function inherit_genes!(child::MLPChromosome, g::Global, parent1::MLPChromosome, parent2::MLPChromosome)
    # Applies the crossover operator.
    @assert parent1.fitness >= parent2.fitness

    # Crossover node genes
    for (key, add) in parent1.genes_mapping
        node1 = get_node_from_dict(parent1, key)
        layer, node_pos = add
        if haskey(parent2.genes_mapping, key)
            node2 = get_node_from_dict(parent2, key)
            gene, parent_number = is_same_innov(node1, node2) ? get_child(g, node1, node2) : (deepcopy(node1), 1)
            child.layers[layer].nodes[node_pos] = gene
        else # Copy excess or disjoint genes from the fittest parent
            child.layers[layer].nodes[node_pos] = deepcopy(node1)
        end
    end
end

function generate_and_inherit_genes(g::Global, parent1::MLPChromosome, parent2::MLPChromosome)
    child = MLPChromosome(g, parent1, parent2)
    layer_nodes, parent_mapping = get_nodes_from_parents(g, parent1, parent2)
    layers = build_layers(g, layer_nodes, parent_mapping, parent1, parent2)
    child.layers = layers
    child.n_hidden_layers = length(layers)-1

    for l in child.layers
        merge!(child.genes_mapping, get_layer_dictionary(l))
    end

    child
end

function get_nodes_from_parents(g::Global, parent1::MLPChromosome, parent2::MLPChromosome)
    parent_mapping = Dict{Int, Int}([])
    p1_number_layers, p2_number_layers = length(parent1.layers), length(parent2.layers)
    number_layers = max(p1_number_layers, p2_number_layers)
    layer_nodes = [NodeGene[] for i=1:number_layers]
    # Crossover node genes
    for (key, add) in parent1.genes_mapping
        is_output_layer = false
        node1 = get_node_from_dict(parent1, key)
        layer_position, node_pos = add
        is_output_layer = layer_position == p1_number_layers
        node_gene = node1
        parent_number = 1
        if haskey(parent2.genes_mapping, key)
            node2 = get_node_from_dict(parent2, key)
            if is_same_innov(node1, node2)
                node_gene, parent_number = get_child(g, node1, node2)
            end
            layer_position_node2, _ = parent2.genes_mapping[key]
            is_output_layer = layer_position_node2 == p2_number_layers
        end
        if is_output_layer
            push!(layer_nodes[end], deepcopy(node_gene))
        else
            push!(layer_nodes[layer_position], deepcopy(node_gene))
        end
        parent_mapping[key] = parent_number
    end
    parent_number = 2
    for (key, add) in parent2.genes_mapping
        if !haskey(parent1.genes_mapping, key)
            layer_position, node_pos = add
            node_gene = get_node_from_dict(parent2, key)
            push!(layer_nodes[layer_position], deepcopy(node_gene))
            parent_mapping[key] = parent_number
        end
    end

    #=for l = 2:length(layer_nodes)
        length_previous_layer = length(layer_nodes[l-1])
        for n in layer_nodes[l] n.ninputs = length_previous_layer  end
    end=#

    sort!(layer_nodes[end], by = x -> x.id)
    layer_nodes, parent_mapping
end

function build_layers(g::Global, layers_nodes::Vector{Vector{NodeGene}}, parent_mapping::Dict, parent1::MLPChromosome, parent2::MLPChromosome)
    number_layers = length(layers_nodes)
    for current_layer = 2:number_layers
        number_nodes_last_layer = length(layers_nodes[current_layer-1])
        for node in layers_nodes[current_layer]
            node.ninputs = number_nodes_last_layer
            node.weights = randn(number_nodes_last_layer)
            parent_number = parent_mapping[node.id]
            node_mapping = parent_number == 1 ? parent1.genes_mapping : parent2.genes_mapping
            parent = parent_number == 1 ? parent1 : parent2
            layer_position, node_pos = node_mapping[node.id]
            last_layer_nodes = filter((n) -> parent_mapping[n[2].id] == parent_number, collect(enumerate(layers_nodes[current_layer-1])))
            for (weight_pos, last_layer_node) in last_layer_nodes
                _ , last_layer_node_pos = node_mapping[last_layer_node.id]
                node.weights[weight_pos] = parent.layers[layer_position].nodes[node_pos].weights[last_layer_node_pos]
            end
        end
    end
    layers = [MLP_layer(g, :HIDDEN, i, g.cf.nn_activation, layers_nodes[i]) for i = 1:(number_layers-1)]
    push!(layers, MLP_layer(g, :OUTPUT, number_layers, :identity, layers_nodes[number_layers]))
    layers
end

function get_child(g::Global, n1::NodeGene, n2::NodeGene)
    weights = deepcopy(n1.weights)
    inputs = min(n1.ninputs, n2.ninputs)
    parent = rand((1:2))
    bias = 0
    act = :sigm
    if parent == 1
        weights[1:inputs] = n1.weights[1:inputs]
        bias = n1.bias
        act = n1.activation
    else
        weights[1:inputs] = n2.weights[1:inputs]
        bias = n2.bias
        act = n2.activation
    end
    NodeGene(n1.id, :HIDDEN, weights, bias, act), parent
end

function chromosome_complexity(ch::MLPChromosome)
    number_weights = 0
    number_nodes = ch.input_size
    for l in ch.layers
        layer_nodes = count((n) -> n.enable_bit, l.nodes)
        number_weights += number_nodes*layer_nodes # + layer_nodes # bias weights
        number_nodes = layer_nodes
    end
    number_weights
end

max_stars_repo_name: Lucas-Muniz/MNEAT

max_issues_repo_name: Lucas-Muniz/MNEAT

max_forks_repo_name: Lucas-Muniz/MNEAT


text:
module StartMacros
using  ..MacroDatabase
using  ..MacroExpander
using  ..ShellCompiler
using  ..ShellInterop
export add_macro

function macro_maker(ARGS...; input="")
  name = ARGS[1]

  i  = 2
  ev = false
  mt = false
  fn = nothing
  jl = false

  only_one_error = "a macro can only have one transformation method"
  function check_multiple_methods()
    @assert input == "" && fn == nothing && !jl only_one_error
  end

  while i <= length(ARGS)
    if ARGS[i] == "-f"
      check_multiple_methods()
      fn = shellfunc(ARGS[i+1])
      i += 1
    elseif ARGS[i] == "-i"
      check_multiple_methods()
      fn = parse_shellfunc(name, ARGS[i+1])
      i += 1
    elseif ARGS[i] == "-j"
      jl = true
    elseif ARGS[i] == "-e"
      ev = true
    elseif ARGS[i] == "-m"
      mt = true
    else
      error("invalid argument $(ARGS[i])")
    end
    i += 1
  end

  if fn == nothing
    if input == ""
      error("no transformation methods specified")
    else
      fn = jl ? parse_juliafunc(name,input) : parse_shellfunc(name,input)
    end
  end

  Macros["@$name"] = Macro("@$name", fn, ev, mt)
  return ""
end
Macros["@macro"] = Macro("@macro", macro_maker, false, true)

function heredoc_dlm(name)
  replace(String(gensym("fxcrd_sh_julia_$name")), "#", "%")
end

function echo_code(ARGS...; input="")
  dlm = heredoc_dlm("code")
  if isempty(ARGS)
    return "cat - <<$dlm\n$input\n$dlm"
  elseif ARGS == ("-s",)
    return "cat - <<'$dlm'\n$input\n$dlm"
  else
    error("Invalid arguments to @code macro")
  end
end
Macros["@code"] = Macro("@code", echo_code, false, false)

function add_macro(str,input="")
  args = (x->x.token).(tokenize(str))
  macro_maker(args...,input=input)
end

module ShellEnv end

function parse_juliafunc(name,input)
  eval(ShellEnv, parse("""
    function $name(ARGS...; input="")
      $input
    end
  """))
end

function parse_shellfunc(name, f_input)
  expanded_input = compile_str(f_input)
  out = ShellCompiler.convert_to_julia(expanded_input)
  if out == nothing
    shell() do sh
      eval_sh(sh, "__fxcrd_MACROFUNC_$name(){\n$expanded_input\n}")
    end
    shellfunc("__fxcrd_MACROFUNC_$name")
  else
    return ShellEnv.eval(out)
  end
end

function shellfunc(name)
  function f(ARGS...; input="")::String
    dlm=heredoc_dlm("function")
    heredoc="<<'$dlm'\n$input\n$dlm"

    shell() do sh
      read_sh(sh, "$name $(ARGS...) $heredoc")
    end
  end
end

end

max_stars_repo_name: fcard/Fxcrd

max_issues_repo_name: fcard/Fxcrd

max_forks_repo_name: fcard/Fxcrd


text:
"""
A test for gage(ProjectedRotationAngle,...).

Arguments
----------
:Type{ProjectedRotationAngle}
    must be type ``ProjectedRotationAngle``

input:AbstractArray
    the vector to be rotated

target:AbstractArray
    the target vector to be aligned to

ref_axis:AbstractArray
    the reference axis around which input vector will be rotated

solution:AbstractFloat
    answer key

msg="":AbstractString
    optional test message

tol_near_zero=1e-15:AbstractFloat
    (optional) tolerance for the length of the input vectors to be close to zero
"""
function test_gage(::Type{ProjectedRotationAngle},
    input::AbstractArray, target::AbstractArray, ref_axis::AbstractArray,
    solution::AbstractFloat, msg::AbstractString="", tol_near_zero::AbstractFloat=1e-15)
    print_dashed_line(80)
    print_with_color(:blue, "Test gage(ProjectedRotationAngle,...)\n\n")
    if msg != ""
        print_with_color(:blue, "$(msg)\n\n")
    end

    @time answer = gage(ProjectedRotationAngle, input, target, ref_axis, tol_near_zero)
    @test_approx_eq answer solution

    print_with_color(:green, "VERIFIED!\n")
    print_dashed_line(80)
    return true
end

test_gage(ProjectedRotationAngle, [1.,0.,0.], [0.,1.,0.], [0., 0., 2.], pi/2,"test 1")
test_gage(ProjectedRotationAngle, [1.,0.,0.], [0.,0.,1.], [0., 1., 0.], pi/2,"test 2")
test_gage(ProjectedRotationAngle, [1./sqrt(2),1./sqrt(2),1.], [0.,0.,1.], [1., -1., 0.], pi/4,"test 3")


text:
include("header.jl")
include("simulator.jl")
include("visual.jl")

# load mech
gas = CreateSolution(mech);
const ns = gas.n_species;
const nr = gas.n_reactions;
const nu = ns + 2;
const np = nr;
const ind_diag = diagind(ones(nu - 1, nu - 1));
const ones_nu = ones(nu - 1);

include("sensitivity.jl")

# set condition
phi = 1.0;          # equivalence ratio
P = 40.0 * one_atm; # pressure, atm
T0 = 1200.0;        # initial temperature, K
get_Tcurve(phi, P, T0, zeros(nr); dT=dT, doplot=true, dTabort=dTabort);

# sampling for sensitivity
rng = Random.MersenneTwister(0x7777777);
n_sample = Int64(ceil(10 * log(np)));
p_sample = rand(rng, n_sample, np) .* 1.0 .- 0.5;
τ_sample = zeros(n_sample); # IDT sample
∇f_sample = similar(p_sample); # sensitivities sample

epochs = ProgressBar(1:n_sample);
for i in epochs
    p = p_sample[i, :];

    ts, pred = get_Tcurve(phi, P, T0, p; dT=dT, dTabort=dTabort);
    idt = interpx(ts, pred[end,:], pred[end,1] + dT);
    τ_sample[i] = deepcopy(idt);

    set_description(
        epochs,
        @sprintf("sample %d with %s", i, method)
    );

    if method == "sensBVP_mthread" # under testing
        # ts, pred = downsampling(ts, pred; dT=0.5, verbose=false);
        ∇f_sample[i, :] = sensBVP_mthread(ts, pred, p);
        # ∇f_sample[i, :] = sensBVP(ts, pred, p);

    elseif method == "sensBF_mthread" # unsuitable for large mechanism
        ∇f_sample[i, :] = sensBF_mthread(phi, P, T0, p; dT=dT, dTabort=dTabort, pdiff=5e-3);

    elseif method == "sensBFSA" # by ForwardDiff
        ∇f_sample[i, :] = sensBFSA(phi, P, T0, p; dT=dT, dTabort=dTabort);

    else
        @show "Wrong sensitivity method.";
    end
end

# eigen decomposition
C = ∇f_sample' * ∇f_sample;
eigs, eigvec = eigen(C);
eigs = reverse(eigs);
eigvec = reverse(eigvec, dims=2);

@save string(exp_path, "/eigs_$method.bson") eigs eigvec ∇f_sample p_sample τ_sample;
# @load string(exp_path, "/eigs_$method.bson") eigs eigvec ∇f_sample p_sample τ_sample;

# show resutls
#pyplot()
l_plt = []
plt = plot(xlabel="Index", ylabel="Eigenvalues", legend=false);
plot!(1:12, eigs[1:12], lw=3, marker=:circle, yscale=:log10);
xlims!((1, 12));
xticks!(1:1:12);
push!(l_plt, plt)

plt = plot(xlabel="Active variable", ylabel="IDT [s]", legend=false);
scatter!(p_sample * eigvec[:,1], τ_sample, yscale=:log10, msize=8);
push!(l_plt, plt)

plt = plot(xlabel="Active variable 1", ylabel="Active variable 2", legend=true);
z = log10.(τ_sample);
scatter!(p_sample * eigvec[:,1], p_sample * eigvec[:,2], marker_z=z, msize=8);
push!(l_plt, plt)

# For paper
h = plot(l_plt..., legend=false, framestyle = :box,
                      size = (1000, 900));
png(h, string(exp_path, "/eigs_$method.png"));

max_stars_repo_name: DENG-MIT/ArrheniusActiveSubspace

max_issues_repo_name: DENG-MIT/ArrheniusActiveSubspace

max_forks_repo_name: DENG-MIT/ArrheniusActiveSubspace


text:
export produce_or_load, tagsave

"""
    produce_or_load([prefix="",] c, f; kwargs...) -> file
Let `s = savename(prefix, c, suffix)`.
If a file named `s` exists then load it and return it.

If the file does not exist then call `file = f(c)`, save `file` as
`s` and then return the `file`.
The function `f` must return a dictionary.
The macros [`@dict`](@ref) and [`strdict`](@ref) can help with that.

## Keywords
* `tag = false` : Add the Git commit of the project in the saved file.
* `projectpath = projectdir()` : Path to search for a Git repo.
* `suffix = "bson"` : Used in `savename`.
* `force = false` : If `true` then don't check if file `s` exists and produce
  it and save it anyway.
* `kwargs...` : All other keywords are propagated to `savename`.

See also [`savename`](@ref) and [`tag!`](@ref).
"""
produce_or_load(c, f; kwargs...) = produce_or_load("", c, f; kwargs...)
function produce_or_load(prefix::String, c, f;
    tag::Bool = false, projectpath = projectdir(),
    suffix = "bson", force = false, kwargs...)

    s = savename(prefix, c, suffix; kwargs...)

    if !force && isfile(s)
        file = wload(s)
        return file
    else
        if force
            @info "Producing file $s now..."
        else
            @info "File $s does not exist. Producing it now..."
        end
        file = f(c)
        try
            mkpath(dirname(s))
            if tag
                tagsave(s, file, projectpath)
            else
            wsave(s, copy(file))
        end
        catch er
            @warn "Could not save file, got error $er. "*
            "\nReturning the file nontheless."
        end
        return file
    end
end

"""
    tagsave(file::String, d::Dict [, path = projectdir()])
First [`tag!`](@ref) dictionary `d` using the project in `path`,
and then save `d` in `file`.
"""
function tagsave(file, d, path = projectdir())
    d2 = tag!(d, path)
    wsave(copy(d2))
    return d2
end


text:
using Oceananigans.Architectures: AbstractCPUArchitecture
using Oceananigans.Fields: AbstractField, ReducedField

import Oceananigans.Fields: set!

const CubedSphereCPUField = CubedSphereField{X, Y, Z, <:AbstractCPUArchitecture} where {X, Y, Z}
const CubedSphereGPUField = CubedSphereField{X, Y, Z, <:AbstractGPUArchitecture} where {X, Y, Z}

function set!(u::CubedSphereCPUField , v::CubedSphereCPUField)
    for (u_face, v_face) in zip(faces(u), faces(v))
        @. u_face.data.parent = v_face.data.parent
    end
    return nothing
end

set!(field::CubedSphereCPUField, f::Function) = [set_face_field!(field_face, f) for field_face in faces(field)]
set!(field::CubedSphereCPUField, f::Number) = [set_face_field!(field_face, f) for field_face in faces(field)]
set!(field::CubedSphereGPUField, f::Function) = [set_face_field!(field_face, f) for field_face in faces(field)]
set!(field::CubedSphereGPUField, f::Number) = [set_face_field!(field_face, f) for field_face in faces(field)]

set_face_field!(field, a) = set!(field, a)


text:
using Gadfly
using Cairo
using TSne
using Statistics: mean, std
using CSV
using DataFrames: DataFrame
using DelimitedFiles: writedlm

"""
Normalize `A` columns, so that the mean and standard deviation
of each column are 0 and 1, resp.
"""
function rescale(A, dim::Integer=1)
    res = A .- mean(A, dims=dim)
    res ./= map(x -> x > 0.0 ? x : 1.0, std(A, dims=dim))
    res
end

if length(ARGS)==0
    println("usage:\n\tjulia demo.jl iris\n\tjulia demo.jl mnist")
    exit(0)
end

if ARGS[1] == "iris"
    using RDatasets
    println("Using Iris dataset.")
    iris = dataset("datasets","iris")
    X = convert(Matrix{Float64}, iris[:, 1:4])
    labels = iris[:, 5]
    plotname = "iris"
    initial_dims = -1
    iterations = 1500
    perplexity = 15
elseif ARGS[1] == "mnist"
    using MLDatasets
    println("Using MNIST dataset.")
    X, labels = MNIST.traindata(Float64);
    npts = min(2500, size(X, 3), length(labels))
    labels = string.(labels[1:npts])
    X = rescale(transpose(reshape(X[:, :, 1:npts], (28*28, :))))
    plotname = "mnist"
    initial_dims = 50
    iterations = 1000
    perplexity = 20
else
    error("Unknown dataset \"", ARGS[1], "\"")
end

println("X dimensions are: ", size(X))
Y = tsne(X, 2, initial_dims, iterations, perplexity)
println("Y dimensions are: ", size(Y))

CSV.write(plotname*"_tsne_out.csv", DataFrame(Y), writeheader=false)
open("labels.txt", "w") do io
    writedlm(io, (labels,), '\n')
end

theplot = plot(x=Y[:,1], y=Y[:,2], color=labels)
draw(PDF(plotname*".pdf", 4inch, 3inch), theplot)
#draw(SVG(plotname*".svg", 4inch, 3inch), theplot)


text:
export relative_genebodycoverage

"""
    relative_genebodycoverage(path_bam::String, path_bed12::String;
							output_prefix::String="", transcript_length_cut::Int=100, max_depth::Int=0, N_bin::Int = 100)

Calculates the relative gene body coverage.

# Arguments
---------
- `path_bam::String`: the path to a BAM file.
- `path_bed12::String`: the path to a reference gene model file with BED12 format.
- `output_prefix::String`: the prefix for output files. If this keyword is not specified or set to "" (defalut), no output files are saved.
- `transcript_length_cut`: the minimum length of transcripts to be used for calculation (Defalut: 100)
- `max_depth::Int`: the maximum depth for a position. Depth more than `max_depth` is cut to `max_depth`.
					If this keyword is set to 0 (defalut), no cut is occurred.
- `N_bin::Int`: the number of bins for a transcript (Defalut: 100).
"""
function relative_genebodycoverage(path_bam::String, path_bed12::String;
									output_prefix::String="", transcript_length_cut::Int=100, max_depth::Int=0, N_bin::Int = 100)
	if output_prefix != ""
		path_out = output_prefix * ".geneBodyCoverage.txt"
		path_out_plot = output_prefix * ".geneBodyCoverage.pdf"
		println(@sprintf "Start calculating relative gene body coverage...\n- bam: %s\n- bed12: %s\n- output: %s\n          %s\n- transcript_length_cut: %d\n- max_depth: %d\n- N_bin: %d\n" path_bam path_bed12 path_out path_out_plot transcript_length_cut max_depth N_bin)
	end

	# File check
	if !isfile(path_bam)
		error(@sprintf "No such file: %s\n" path_bam)
	end
	if !isfile(path_bed12)
		error(@sprintf "No such file: %s\n" path_bed12)
	end
	if output_prefix != ""
		if !uv_access_writable(dirname(output_prefix))
			error(@sprintf "Output files are not writable to: %s\n" (dirname(output_prefix) != "" ? dirname(output_prefix) : "."))
		end
	end

	# Loads transcripts information from BED12-format file
	transcripts = load_transcript(path_bed12)
	println(@sprintf "transcripts: %d" length(transcripts))

	# Define relative gene body coverage output
	relcov = zeros(Float64, N_bin)

	# Open BAM file
	open(BAM.Reader, path_bam, index=path_bam*".bai") do bam_reader

		# For each transcript
		for t in transcripts

			# Skip
			if sum(BED.blocksizes(t)) < transcript_length_cut
				continue
			end

			if BED.strand(t) != STRAND_POS && BED.strand(t) != STRAND_NEG  # STRAND_NA, STRAND_BOTH
				continue
			end

			# Skip the transcript whose region does not exist in BAM
			if ! exists_region_bam(bam_reader, BED.chrom(t), BED.chromstart(t), BED.chromend(t))
				continue
			end

			# Get read coverage of the percentile positions on a transcript
			if max_depth == 0
				relcov += coverage_transcript_percentile(bam_reader, t)
			else
				# The maximum number for (imcomplete) comaptibility with `--max-depth` parametr of `pileup()` used in RSeQC
				relcov += min.(max_depth, coverage_transcript_percentile(bam_reader, t))
			end
		end
    end

	if output_prefix != ""
		# Write result
		sample_id = replace(basename(path_bam), r"\.[^\.]+$" => s"")
		write_relative_genebodycoverage(relcov, path_out, sample_id)

		# Save plot
		plot_relative_coverage(relcov, out_path=path_out_plot)

		# Message
		println(@sprintf "Finished! Check output files:\n- %s\n- %s" path_out path_out_plot)
	end

	return relcov
end

"""
    coverage_transcript_percentile(bam_reader::BAM.Reader, t::BED.Record)

Gets read coverage on the percentile positions in a transcript
1 read for 1 transcript
"""
function coverage_transcript_percentile(bam_reader::BAM.Reader, t::BED.Record)

	# Get block_sizes  (= exon lengths) and block_starts (=start position of exons relative to chromStart of t)
	block_sizes  = BED.blocksizes(t)
	block_starts = BED.blockstarts(t) # Returns 1-based coordinate (though 0-based in BED12 format)
	transcript_start = BED.chromstart(t)  # Returns 1-based coordinate(though 0-based in BED12 format)
	transcript_end = BED.chromend(t)

	# Collect within-exon positions relative to the transcript start
	L = sum(block_sizes)
	exon_coordinates = zeros(Int, L)
	offset = 0
	for i in 1:length(block_sizes) # For each exon
		for j in 0:(block_sizes[i]-1)
			exon_coordinates[j+offset+1] = block_starts[i] + j
		end
		offset += block_sizes[i]
	end

	# Select percentile positions (resulting in 100 positions)
	percentile_positions = percentile_list(exon_coordinates)

	# Note: Read coverage is saved in `cov` in a 'Left justified' manner,
	#       and left is 5'-end.
	# If the strand of a transcript is '-' STRAND_NEG, read coverage was reversed.
	if BED.strand(t) == STRAND_NEG
		percentile_positions = reverse(percentile_positions)
	end

	return readcoverage_bam_base(bam_reader, BED.chrom(t), transcript_start, transcript_end)[percentile_positions]
end

"""
    write_relative_genebodycoverage(relative_coverage::Array{Float64,1}, path_out::String, sample_id::String)

Write result of `relative_genebodycoverage()`
"""
function  write_relative_genebodycoverage(relative_coverage::Array{Float64,1}, path_out::String, sample_id::String)
	open(path_out, "w") do fw
		# Header
		write(fw, "Percentile\t" * join(string.(collect(1:100)), "\t") * "\n")
		# Coverage
		write(fw, sample_id * "\t" * join(string.(relative_coverage), "\t") * "\n")
	end
end

"""
    convert2relativecoverage(converage::Array{Float64,1}, bin_number::Int)
Converts read coverage to relative gene body coverage

*future work to do (incomplete)
"""
function convert2relativecoverage(converage::Array{Float64,1}, bin_number::Int)
	relative_coverage = zeros(Float64, bin_number)
	L = length(coverage)

	for b in 1:bin_number

	end
end

"""
    percentile_list(N::Array{Int,1})

Finds the percentile of a list of values and returns the list of percentile of the values.

# Arguments
- `N::Array{Int,1}`: the list of values. Note N MUST BE already sorted.

Reference: https://github.com/MonashBioinformaticsPlatform/RSeQC/blob/cb42bd90afa8d131875b18d54f70821142b8ea79/rseqc/qcmodule/mystat.py#L156
"""
function percentile_list(N::Array{Int,1})
	if length(N) < 100
		return N
	end

	percentile = zeros(Int, 100)

	for i in 1:100
		k = (length(N)-1) * i/100.0
		f = floor(k)
		c = ceil(k)
		if f == c
			percentile[i] = N[Int(k)+1]
		else
			d0 = N[Int(f)+1] * (c-k)
			d1 = N[Int(c)+1] * (k-f)
			percentile[i] = Int(round(d0+d1))
		end
	end

	return percentile
end


text:
# paper figure for product of circular densities

using LinearAlgebra
# using DocStringExtensions
using KernelDensityEstimate

using Cairo
using Gadfly, Colors
using Random

using ApproxManifoldProducts

using TransformUtils
const TU = TransformUtils

setForceEvalDirect!(true)

## TOP LEFT two equally spaced

pts1 = TU.wrapRad.(0.1*randn(100).-pi/2.0)
pts2 = TU.wrapRad.(0.1*randn(100).+pi/2.0)

# pc1 = kde!_CircularNaiveCV(pts1)
# pc2 = kde!_CircularNaiveCV(pts2)
pc1 = kde!(pts1, [0.1;], (addtheta,), (difftheta,))
pc2 = kde!(pts2, [0.1;], (addtheta,), (difftheta,))

dummy = kde!(rand(100),[1.0;], (addtheta,), (difftheta,));

pGM, = prodAppxMSGibbsS(dummy, [pc1; pc2], nothing, nothing, Niter=2, addop=(addtheta,), diffop=(difftheta,), getMu=(getCircMu,));

pc12 = kde!_CircularNaiveCV(pGM[:])

pltl = plotKDECircular([pc1;pc2; pc12], scale=0.07, c=["green";"blue";"orange"])

pltl.coord = Coord.Cartesian(xmin = -1.6,
xmax = 1.6,
ymin = -1.6,
ymax = 1.6,
aspect_ratio=1.0)

pltl

##

pltl |> PDF("topleft.pdf",8cm,8cm)
pltl |> SVG("topleft.svg",8cm,8cm)
# @async run(`evince topright.pdf`)

## TOP RIGHT

pts1 = TU.wrapRad.(0.1*randn(100).-pi.+0.6)
pts2 = TU.wrapRad.(0.1*randn(100).+pi.-0.6)

# pc1 = kde!_CircularNaiveCV(pts1)
# pc2 = kde!_CircularNaiveCV(pts2)
pc1 = kde!(pts1, [0.1;], (addtheta,) , (difftheta,))
pc2 = kde!(pts2, [0.1;], (addtheta,) , (difftheta,))

# TODO: make circular KDE
dummy = kde!(rand(100),[1.0;], (addtheta,), (difftheta,));

pGM, = prodAppxMSGibbsS(dummy, [pc1; pc2], nothing, nothing, Niter=2,
                          addop=(addtheta,), diffop=(difftheta,), getMu=(getCircMu,));

pc12 = kde!_CircularNaiveCV(pGM[:])
getBW(pc12)[1,1]
# pc12 = kde!(pGM, [0.05;], (addtheta,), (difftheta,))

pltr = plotKDECircular([pc1;pc2; pc12], scale=0.07, c=["green";"blue";"orange"])

pltr.coord = Coord.Cartesian(xmin = -1.6,
xmax = 1.6,
ymin = -1.6,
ymax = 1.6,
aspect_ratio=1.0)

pltr

##

pltr |> PDF("topright.pdf",8cm,8cm)
pltr |> SVG("topright.svg",8cm,8cm)
# @async run(`evince topright.pdf`)

## BOTTOM LEFT  handrolically trying product of 3

pts1 = TU.wrapRad.(0.1*randn(100))
pts2 = TU.wrapRad.(0.1*randn(100).-2pi/3.0.+1.0)
pts3 = TU.wrapRad.(0.1*randn(100).+2pi/3.0)

# pc1 = kde!_CircularNaiveCV(pts1)
# pc2 = kde!_CircularNaiveCV(pts2)
pc1 = kde!(pts1, [0.1;], (addtheta,), (difftheta,))
pc2 = kde!(pts2, [0.1;], (addtheta,), (difftheta,))
pc3 = kde!(pts3, [0.1;], (addtheta,), (difftheta,))

# pl = plotKDECircular([pc1;pc2;pc3])

dummy = kde!(rand(100),[1.0;], (addtheta,), (difftheta,));

pGM, = prodAppxMSGibbsS(dummy, [pc1; pc2; pc3], nothing, nothing, Niter=1, addop=(addtheta,), diffop=(difftheta,), getMu=(getCircMu,));

pc123 = kde!_CircularNaiveCV(pGM[:])

plbl = plotKDECircular([pc1;pc2;pc3; pc123], scale=0.07, c=["green";"blue";"magenta";"orange"])

plbl.coord = Coord.Cartesian(xmin = -1.6,
xmax = 1.6,
ymin = -1.6,
ymax = 1.6,
aspect_ratio=1.0)

plbl

##

plbl |> PDF("bottomleft.pdf",8cm,8cm)
plbl |> SVG("bottomleft.svg",8cm,8cm)
# @async run(`evince topright.pdf`)

## BOTTOM RIGHT product of 3 equally spaced

pts1 = TU.wrapRad.(0.1*randn(100))
pts2 = TU.wrapRad.(0.1*randn(100).-2pi/3.0)
pts3 = TU.wrapRad.(0.1*randn(100).+2pi/3.0)

# pc1 = kde!_CircularNaiveCV(pts1)
# pc2 = kde!_CircularNaiveCV(pts2)
pc1 = kde!(pts1, [0.1;], (addtheta,), (difftheta,))
pc2 = kde!(pts2, [0.1;], (addtheta,), (difftheta,))
pc3 = kde!(pts3, [0.1;], (addtheta,), (difftheta,))

# pl = plotKDECircular([pc1;pc2;pc3])

dummy = kde!(rand(100),[1.0;], (addtheta,), (difftheta,));

pGM, = prodAppxMSGibbsS(dummy, [pc1; pc2; pc3], nothing, nothing, Niter=1, addop=(addtheta,), diffop=(difftheta,), getMu=(getCircMu,));

pc123 = kde!_CircularNaiveCV(pGM[:])

plbr = plotKDECircular([pc1;pc2;pc3; pc123], scale=0.07, c=["green";"blue";"magenta";"orange"])

plbr.coord = Coord.Cartesian(xmin = -1.6,
xmax = 1.6,
ymin = -1.6,
ymax = 1.6,
aspect_ratio=1.0)

plbr

##

plbr |> PDF("bottomright.pdf",8cm,8cm)
plbr |> SVG("bottomright.svg",8cm,8cm)
# @async run(`evince topright.pdf`)

##  Stack

pls = vstack(hstack(pltl,pltr), hstack(plbl,plbr))

##

pls |> PDF("plstack.pdf",12cm,12cm)
pls |> SVG("plstack.svg",12cm,12cm)


text:
"""Based on the desired aspect ratio L_ratio of the domain returns the number of
grid cells in y-direction and the length of that side of the domain."""
function domain_ratio(nx::Int,Lx::Real,L_ratio::Real)
    Δ = Lx / nx                        # grid spacing
    ny = Int(round(Lx / L_ratio / Δ))  # number of grid cells in y-direction
    Ly = ny * Δ                        # length of domain in y-direction
    return Δ,ny,Ly
end

"""Similar to the numpy meshgrid function:
repeats x length(y)-times and vice versa. Returns two matrices xx,yy of same shape so that
each row of xx is x and each column of yy is y."""
function meshgrid(x,y)
    # preallocate preserving the data type of x,y
    xx = zeros(typeof(x[1]),length(x),length(y))
    yy = zeros(typeof(y[1]),length(x),length(y))

    for i in 1:length(x)
        xx[i,:] .= x[i]
    end

    for i in 1:length(y)
        yy[:,i] .= y[i]
    end

    return xx,yy
end

""" Returns the time step dt,Δt,dtint and the total number of steps to integrate nt
based on the CFL criterion to resolve shallow water gravity waves. """
function timestep()
    # shallow water gravity wave phase speed
    c_phase = sqrt(gravity*water_depth)

    # the model timestep dt based on cfl stability criterion to resolve gravity waves
    # converting to integer, i.e. rounding up (ceil)
    dtint = Int(floor(cfl*Δ/c_phase))      # make the timestep slightly shorter due to floor
    nt = Int(ceil(Ndays*3600*24/dtint))    # number of time steps to integrate
    dt = Numtype(dtint)                    # convert to Numtype for multiplication of the RHS
    Δt = Numtype(dtint/Δ)                  # timestep combined with grid spacing Δ
    nΔt_diff = Numtype(nstep_diff)*Δt      # timestep for diffusive terms
    return dt,Δt,nΔt_diff,dtint,nt
end

""" Returns the tracer advection time step dtadv and the number of timesteps nadvstep after which
one evaluation of the tracer advection is computed. """
function adv_timestep()
    # round down to make sure nadvstep is integer, at least 1
    nadvstep = max(1,Int(floor(Δ/Uadv/dtint)))
    nadvstep_half = nadvstep ÷ 2    #
    # recompute the tracer advection time step to fit the rounding
    dtadvint = nadvstep*dtint
    dtadvu = Numtype(dtadvint*nx/Lx)    # [s/m] for dimensionless advection grid
    dtadvv = Numtype(dtadvint*ny/Ly)
    return dtadvu,dtadvv,dtadvint,nadvstep,nadvstep_half
end

const Δ,ny,Ly = domain_ratio(nx,Lx,L_ratio)

# number of grid points for u,v,q-grid in either x or y-direction
const nux = if (bc_x == "periodic") nx else nx-1 end
const nuy = ny
const nvx,nvy = nx,ny-1
const nqx = if (bc_x == "periodic") nx else nx+1 end
const nqy = ny+1

# total number of T,u,v,q-points
const nT = nx*ny
const nu = nux*nuy
const nv = nvx*nvy
const nq = nqx*nqy

# grid vectors
const x_T = Δ*Array(1:nx) .- Δ/2
const y_T = Δ*Array(1:ny) .- Δ/2

const x_u = if (bc_x == "periodic") Δ*Array(0:nx-1) else Δ*Array(1:nx-1) end
const y_u = y_T

const x_v = x_T
const y_v = Δ*Array(1:ny-1)

const x_q = if bc_x == "periodic" x_u else Δ*Array(1:nx+1) .- Δ end
const y_q = Δ*Array(1:ny+1) .- Δ

# halo of ghost points - don't change.
const halo = 2          # 2 because of biharmonic diffusion operator (wider stencil)
const haloη = 1
const halosstx = 1      #TODO change for domain decomposition
const halossty = 0

# is there a point on the left edge? ep - egde points
# used in some functions of rhs.jl to avoid an if
const ep = if bc_x == "periodic" 1 else 0 end
#TODO possibly different edge point variables needed on subdomains

# halo versions with additional ghost points
const x_T_halo = Δ*Array(0:nx+1) .- Δ/2
const y_T_halo = Δ*Array(0:ny+1) .- Δ/2

const x_u_halo = if (bc_x == "periodic") Δ*Array(-2:nx+1) else Δ*Array(-1:nx+1) end
const y_u_halo = Δ*Array(-1:ny+2) .- Δ/2

const x_v_halo = Δ*Array(-1:nx+2) .- Δ/2
const y_v_halo = Δ*Array(-1:ny+1)

# also two halo for q
const x_q_halo = if bc_x == "periodic" x_u_halo else Δ*Array(-1:nx+3) .- Δ end
const y_q_halo = Δ*Array(-1:ny+3) .- Δ

# time and output
const dt,Δt,nΔt_diff,dtint,nt = timestep()
const nout = Int(floor(output_dt*3600/dtint))   # output every nout time steps
const nout_total = (nt ÷ nout)+1                # total number of time steps for output
const t_vec = Array(0:nout_total-1)*dtint       # time vector for output

# advection time step
const dtadvu,dtadvv,dtadvint,nadvstep,nadvstep_half = adv_timestep()

println((dtint,dtint*nstep_diff,dtint*nstep_advcor,dtadvint,output_dt*3600))

max_stars_repo_name: eschnett/juls

max_issues_repo_name: eschnett/juls

max_forks_repo_name: eschnett/juls