text:
include("../src/MSA.jl")

t1 = "AACGTGATTGAC"
t2 = "TCGAGTGCTTTACAGT"
t3 = "GCCGTGCTAGTCG"
t4 = "TTCAGTGGACGTGGTA"
t5 = "GGTGCAGACC"

input_online = [t1, t2, t3, t4, t5]

#-------- Test input
input1 = "MMVHLTPMMKSAVTALWGKVNVDMVGGMALGRLLVVYPWTQRFFMSFGDLSTPDAVM"
input2 = "MMGLSDGMWQLVLNVWGKVMADIPGHGQMVLIRLFKGHPMTLMKFDKFKHLKSMDMMKAS"
input3 = "ALVMDNNAVAVSFSMMQMALVLKSWAILKKDSANIALRFFLKIFMVAPS"
input4 = "MMRPMPMLIRQSWRAVSRSPLMHGTVLFARLFALMPDLLPLFQYNCRQFSSPMD"

input_paperTest_sequence = [ input1, input2, input3, input4]

input10 = "AABB"
input11 = "ACA"
input12 = "CC"

input_shortTest_sequence = [ input10, input11, input12 ]

input13 = "AATGGTCATAGCGAGATGAAGCCACGTGATGGATAATATTGTGCAAACGACCTTATTAGCTATTGACCGTCGATGTCCAACGAGACAATT"
input14 = "GAATCTGTATTCTTCAAGCTTCAACTCCATGCACTACGAACGGTAGTGGTTCACATTGACCGTG"
input15 = "TTGGGCGCATTGACCGTCCTTCCTAGCGTATCATCAAACTTGTGATTCTCTATCTAGAGCAAAATGCGGTGTCCGCTATATGGAGATCTATTTCAAAATA"
input16 = "GTTACATATCAGAATCATTAGAAACGCTCTTAATGGGGTTAAGCAGAGACTTAGTAAGGATTAACTCCCAAGATGATTGACCGTGCTC"

input_project_example = [ input13, input14, input15, input16 ]

input_rand_seq = []
for i = 1:rand(4:4)
    push!(input_rand_seq, randstring("ACGT", rand(10:25)))
end
#--------

# start with input strings
# we want to end up with a useable chromosome representation for the string, which is the positions of the gaps
# we want an initialization for the population first, which gives us this initial set of chromosomes
# we can always use this to rebuild from the original string, so the original string should be carried over
# a fitness function is off of the original string with the gaps in, so we need to reconstruct from the chromosome and the string, the gapped string
# then use the fitness function

# hyper_parameters
score_list = [1, 0, 0]        # match, mismatch, gap penalty
gap_ext_penalty = -2          # if extending a gap, what is the penalty
initial_population_size = 100     # currently always the maintained population
gap_growth = 1.13                # after max input length is found, this is used to provide gaps for the max as well, ie [3, 4, 7] would be 7 * 1.2 = 9 (rounded up)
                                # so 2 gaps put in the largest one, 5 in the next, and the 3 sized input would get 6 gaps
elitism_prop = 0.3             # proportion of best to select to keep
num_crossover = 3               # number of crossovers to take as parents for generating against rest
children_cap = 100          # how many children can be the result of a crossover, "nothing" means there is no cap
generation_count = 1000          # how many times to run through children
mutation_strength = 0.03        # how much do the gaps shuffle, 1 would be full shuffle
crossover_criteria = 1          # save children from crossover when they are greater than 1 = average, 2 = worst, 3 = best of the 2 parents
crossover_version = 2           # which type of crossover to use, 1 = random cuts and always left/right, 2 = random cuts and random left/right mix
return_early = false            # if true, returns when children cap is met, otherwise when trying to add a new child, sort and pop off the worst
mutation_chance = 0.5            # the % chance each child has to undergo mutation after crossover
score_type = 2                  # 1 = basic match/mismatch/gap, 2 = pam250 and input gap penalty
printout = false                 # whether to print out step-wise operation of the MSA to see how it works

#input_sequence = input_rand_seq
#input_sequence = input_paperTest_sequence
#input_sequence = input_online
#input_sequence = input_shortTest_sequence
#input_sequence = input_project_example
input_sequence = []

score_list, input = msa_from_file("lidy.txt")

for i = 1:length(input)
    push!(input_sequence, input[i][2])
end

for i = 1:length(input_sequence)
    println(input_sequence[i])
end
println()

results = MSA(input_sequence, score_list, initial_population_size, gap_growth, elitism_prop, num_crossover, children_cap, generation_count, mutation_strength,
              crossover_criteria, crossover_version, return_early, mutation_chance, printout, score_type, gap_ext_penalty)

score = results[1]
results = build_strings(input_sequence, results[2])

open("out.txt", "w") do io
    write(io, string("Score: ", score, '\n'))
    for i = 1:length(input)
        write(io, input[i][1], '\n')
        write(io, results[i], '\n')
    end
end;

# code for checking an output score with the gaps inserted - as in the project description or a built string
# MSA outputs to a file the best result of the run at the end
#=
function score_sequence(input_sequence_set, sequence_length, score_list)
    score = 0
    score_buffer = 0
    item_one = nothing
    item_two = nothing
    num_items = length(input_sequence_set)

    # go through sequence length, to index by column, all input sequences + gaps should = sequence length
    for i = 1:sequence_length
        score_buffer = 0

        for j = 1:num_items
            for k = j + 1:num_items
                # get either a blank space or the character
                if input_sequence_set[j][i] == '_' || input_sequence_set[k][i] == '_' || input_sequence_set[j][i] == '-' || input_sequence_set[k][i] == '-'
                    if input_sequence_set[j][i] != input_sequence_set[k][i]
                        score_buffer += score_list[3]
                    end
                elseif input_sequence_set[j][i] == input_sequence_set[k][i]
                    score_buffer += score_list[1]
                else
                    score_buffer += score_list[2]
                end
            end
        end

        score += score_buffer
    end

    return score
end

i1 = "AATGGTCATAGC------------G-----AGATGA--AG-CCACGTGATGGATAATATTGTGCAAACGACCTTATTAGCT---ATTGACCGTCGATGTCCAACGAGACAAT---T"
i2 = "GAATC-------------------T-----GTATTCTTC--AA-GCTTCAACTC-----CAT--------------GCACT---ACGAACGGTAGTGGTTCACATTGACCGT---G"
i3 = "TTGGG-CGCATTGACCGTCCTTCCTAGCGTATCATC--AAACT-TGTGATTCTCTATCTAGAGCAAAATGCGGTGTCCGCT---AT-----AT-GGAGATCTATTTCAAAAT---A"
i4 = "GTTA--CATATC------------A---GAATCATT--AG-AAACGCTCTTAATGGGGTTAAGCAGAG-AC-TTAGTAAGGATTAACTCCCAA-GATGATT-----GACCGTGCTC"

i5 = "_AATG_GTCA_TAGCGA_G_A_T__G__AAGC_CACG_TGATGGATA__AT_ATTGTGCAAAC__G_ACCTTATTA_GCTATT_GA_CCG_TCGAT_G__TCCAACGAGA_CAATT"
i6 = "G_A_A__TC__T_G_TAT_T_CTT__C_AAG__C__T_TCA__ACTCC_ATG___CA_C__T_____ACG__A__A__CGG_TAG_T__G___G_T_T_CA_CA_TT_GACC_GTG"
i7 = "TTGGGCG_CATT_GAC_CGTCCTT_CCT_AG__C_GTATCATCAAACTTGTGATTCT_CTATCTAG_A_GC_AAAATGCGG_T_G_TCCGCTATATGGAGATCTATTTCAA_AATA"
i8 = "GT_TACAT_ATCAGA_ATC_A_TTAG__AAACGCTCT_TAATGG_G___GTTA___AGC_A__GAG_AC_TTAGTAAG_GATTAACTCC_CAAGAT_GA_TTGA_CC_GTGC__TC"

sequence = [i1, i2, i3, i4]
sequence2 = [i5, i6, i7, i8]

println(score_sequence(sequence, length(sequence[1]), score_list))
println(score_sequence(sequence2, length(sequence2[1]), score_list))
=#


text:
include("./logging.jl")

# BEGIN TYPE DescentSolver

"""Descent Solver

Args:
    inst (Instance)
    nb_test (Int): Nombre total de voisins testé
    nb_move (Int): nombre de voisins acceptés (améliorant ou non)
    nb_reject (Int): nombre de voisins refusés
    nb_cons_reject (Int): Nombre de refus consécutifs
    nb_cons_reject_max (Int): Nombre maxi de refus consécutifs

    duration (Float64): Durée réelle (mesurée) de l'exécution
    durationmax (Float64): Durée max de l'exécution (--duration)
    starttime (Float64): Heure de début d'une résolution

    cursol (Solution): Solution courante
    bestsol (Solution): Meilleure Solution rencontrée
    testsol (Solution): Nouvelle solution potentielle

    bestiter (Int)
    do_save_bestsol (Bool)

Notes:
    Amélioriations possibles:

        - revoir la gestion des options du solver (utiliser les params par clé-valeur)
        - plus besoin de gérer cursol (car ici, bestsol suffit contrainrement au ExploreSolver)

"""
mutable struct DescentSolver
    inst::Instance
    nb_test::Int          # Nombre total de voisins testé
    nb_move::Int          # nombre de voisins acceptés (améliorant ou non)
    nb_reject::Int        # nombre de voisins refusés
    nb_cons_reject::Int   # Nombre de refus consécutifs
    nb_cons_reject_max::Int # Nombre maxi de refus consécutifs

    duration::Float64     # durée réelle (mesurée) de l'exécution
    durationmax::Float64  # durée max de l'exécution (--duration)
    starttime::Float64    # heure de début d'une résolution

    cursol::Solution      # Solution courante
    bestsol::Solution     # meilleure Solution rencontrée
    testsol::Solution     # nouvelle solution potentielle

    bestiter::Int
    do_save_bestsol::Bool
end

function DescentSolver(inst::Instance;
                       startsol::Union{Nothing,Solution} = nothing)
    nb_test = 0
    nb_move = 0
    nb_reject = 0
    nb_cons_reject = 0
    nb_cons_reject_max = inst.nb_planes^2

    bestiter = 0

    durationmax = 366*24*3600 # 1 année par défaut !
    duration = 0.0 # juste pour initialisation
    starttime = 0.0 # juste pour initialisation

    if startsol === nothing
        # Pas de solution initiale => on en crée une
        cursol = Solution(inst)
    else
        cursol = startsol
        if lg2()
            println("Dans DescentSolver : cursol = opts[:startsol] ")
            println("cursol", to_s(cursol))
        end
    end

    bestsol = Solution(cursol)
    testsol = Solution(cursol)
    do_save_bestsol = true
    descent_solver = DescentSolver(inst, nb_test, nb_move, nb_reject,
                                   nb_cons_reject, nb_cons_reject_max, duration,
                                   durationmax, starttime, cursol, bestsol,
                                   testsol, bestiter, do_save_bestsol)
    return descent_solver
end

function solve(sv::DescentSolver, neighbour_operator!;
               nb_cons_reject_max::Int = 0,
               startsol::Union{Nothing, Solution} = nothing,
               durationmax::Int = 0,
               record_cost=false
               )
    ln2("BEGIN solve(DescentSolver)")
    if durationmax != 0
        sv.durationmax = durationmax
    end
    if nb_cons_reject_max != 0
        sv.nb_cons_reject_max = nb_cons_reject_max
    end

    if startsol !== nothing
        sv.cursol = startsol
        if sv.cursol.cost < sv.bestsol.cost
            copy!(sv.bestsol, sv.cursol) # cursol peut avoir été meilleure !
        end
        copy!(sv.testsol, sv.cursol)
        if lg2()
            println("Dans DescentSolver : sv.cursol = sv.opts[:startsol] ")
            println("sv.cursol", to_s(sv.cursol))
        end
    # else on garde la dernière solution sv.cursol
    end

    sv.starttime = time_ns()/1_000_000_000

    if lg3()
        println("Début de solve : get_stats(sv)=\n", get_stats(sv))
    end

    # improvement_costs = [sv.bestsol.cost]
    # improvement_steps = [0]
    while !finished(sv)
        copy!(sv.testsol, sv.cursol)
        neighbour_operator!(sv.testsol)
        sv.nb_test += 1

        # current_costs[sv.nb_test] = sv.testsol.cost

        degrad = sv.testsol.cost - sv.cursol.cost
        ln4("degrad=$(degrad)")
        if degrad < 0
            sv.nb_cons_reject = 0
            sv.nb_move += 1
            copy!(sv.cursol, sv.testsol)
            if sv.cursol.cost < sv.bestsol.cost
                copy!(sv.bestsol, sv.cursol)
                # push!(improvement_costs, sv.bestsol.cost)
                # push!(improvement_steps, sv.nb_test)
                sv.bestiter = sv.nb_test
                print_log(sv)
            end
        else
            sv.nb_reject += 1
            sv.nb_cons_reject += 1
        end
    end # fin while !finished
    ln2("END solve(DescentSolver)")
    # return improvement_costs, improvement_steps
end

"""Retourne true ssi l'état justifie l'arrêt de l'algorithme"""
function finished(sv::DescentSolver)
    sv.duration = time_ns()/1_000_000_000 - sv.starttime
    too_long = sv.duration >= sv.durationmax
    too_many_reject = (sv.nb_cons_reject >= sv.nb_cons_reject_max)
    stop = too_long || too_many_reject
    if stop && lg1()
        println("\nSTOP car :")
        println("     sv.nb_cons_reject=$(sv.nb_cons_reject)")
        println("     sv.nb_cons_reject_max=$(sv.nb_cons_reject_max)")
        println("     sv.duration=$(sv.duration)")
        println("     sv.durationmax=$(sv.durationmax)")
        println(get_stats(sv))
    end
    return stop
end

"""
Retourne un quadruplet d'indices destiné à affectuer deux shifts relativement
proches

Args:
    ecartmaxin est l'écart maxi au sein d'une paire d'indices
    ecartmaxout est l'écart maxi entre deux paires d'indices (cumulables)
    ecartmaxin et ecartmaxout sont imposés dans les bornes de l'instance
    abs(i2-i1) et abs(i4-i3) sont limités par ecartmaxin
    i3 est distant au maximum de ecartmaxout du couple (i1,i2)
"""
function sample_two_shifts(sol::Solution; ecartmaxin::Int=10, ecartmaxout::Int=-1)
    # Version stupide car trop large !
    i1, i2, i3, i4 = rand(1:sol.inst.nb_planes, 4)
    return (i1, i2, i3, i4)
end

function record_bestsol(sv::DescentSolver; movemsg="")
    copy!(sv.bestsol, sv.cursol)
    sv.bestiter = sv.nb_test
    if sv.do_save_bestsol
        write(sv.bestsol)
    end
    if lg3()
        print("\niter $(rpad(sv.nb_test, 4)):$(sv.nb_move)+$(sv.nb_reject) ")
        print("$movemsg ")
        print("bestsol=$(to_s(sv.bestsol))")
    elseif lg1()
        print("\niter $(rpad(sv.nb_test, 4)):$(sv.nb_move)+$(sv.nb_reject) ")
        print("$movemsg => bestcost=", sv.cursol.cost)
    end
end

function get_stats(sv::DescentSolver)
    # txt = <<-EOT.gsub /^ {4}/,''
    txt = "
    ==Etat de l'objet DescentSolver==
    sv.nb_test=$(sv.nb_test)
    sv.nb_move=$(sv.nb_move)
    sv.nb_cons_reject=$(sv.nb_cons_reject)
    sv.nb_cons_reject_max=$(sv.nb_cons_reject_max)

    sv.duration=$(sv.duration)
    sv.durationmax=$(sv.durationmax)

    sv.testsol.cost=$(sv.testsol.cost)
    sv.cursol.cost=$(sv.cursol.cost)
    sv.bestsol.cost=$(sv.bestsol.cost)
    sv.bestiter=$(sv.bestiter)
    sv.testsol.solver.nb_infeasable=$(sv.testsol.solver.nb_infeasable)
    "
    txt = replace(txt, r"^ {4}" => "")
end

# END TYPE DescentSolver

ext: jl


text:
import StatsBase: ProbabilityWeights
using Random, Distributions
using Test

rng = MersenneTwister(123)

@testset "Testing matrix-variates with $key" for (key, func) in
    Dict("rand(...)" => [rand, rand],
         "rand(rng, ...)" => [dist -> rand(rng, dist), (dist, n) -> rand(rng, dist, n)])

d = DiscreteNonParametric([40., 80., 120., -60.], [.4, .3, .1,  .2])

@test !(d ≈ DiscreteNonParametric([40., 80, 120, -60], [.4, .3, .1, .2], Distributions.NoArgCheck()))
@test d ≈ DiscreteNonParametric([-60., 40., 80, 120], [.2, .4, .3, .1], Distributions.NoArgCheck())

# Invalid probability
@test_throws ArgumentError DiscreteNonParametric([40., 80, 120, -60], [.5, .3, .1, .2])

# Invalid probability, but no arg check
DiscreteNonParametric([40., 80, 120, -60], [.5, .3, .1, .2], Distributions.NoArgCheck())

test_range(d)
vs = Distributions.get_evalsamples(d, 0.00001)
test_evaluation(d, vs, true)
test_stats(d, vs)
test_params(d)

@test func[1](d) ∈ [40., 80., 120., -60.]
@test func[1](sampler(d)) ∈ [40., 80., 120., -60.]

@test pdf(d, -100.) == 0.
@test pdf(d, -100) == 0.
@test pdf(d, -60.) == .2
@test pdf(d, -60) == .2
@test pdf(d, 100.) == 0.
@test pdf(d, 120.) == .1

@test cdf(d, -100.) == 0.
@test cdf(d, -100) == 0.
@test cdf(d, 0.) ≈ .2
@test cdf(d, 100.) ≈ .9
@test cdf(d, 100) ≈ .9
@test cdf(d, 150.) == 1.

@test ccdf(d, -100.) == 1.
@test ccdf(d, -100) == 1.
@test ccdf(d, 0.) ≈ .8
@test ccdf(d, 100.) ≈ .1
@test ccdf(d, 100) ≈ .1
@test ccdf(d, 150.) == 0

@test quantile(d, 0) == -60
@test quantile(d, .1) == -60
@test quantile(d, .3) == 40
@test quantile(d, .7) == 80
@test quantile(d, 1) == 120
@test minimum(d) == -60
@test maximum(d) == 120

@test insupport(d, -60)
@test insupport(d, -60.)
@test !insupport(d, 20)
@test !insupport(d, 20.)
@test insupport(d, 80)
@test !insupport(d, 150)

xs = support(d)
ws = ProbabilityWeights(probs(d))
@test mean(d) ≈ mean(xs, ws)
@test var(d) ≈ var(xs, ws, corrected=false)
@test skewness(d) ≈ skewness(xs, ws)
@test kurtosis(d) ≈ kurtosis(xs, ws)
@test entropy(d) ≈ 1.2798542258336676
@test mode(d) == 40
@test modes(d) == [40]
@test mgf(d, 0) ≈ 1.0
@test mgf(d, 0.17) ≈ 7.262034e7
@test cf(d, 0) ≈ 1.0
@test cf(d, 0.17) ≈ 0.3604521478 + 0.6953481124im

# Fitting
xs = [1, 2,3,4,5,4,3,2,3,2,1]
ws = [1.,2,1,0,1,1,2,3,2,6,1]

ss = suffstats(DiscreteNonParametric, xs)
@test ss isa Distributions.DiscreteNonParametricStats{Int,Float64,Vector{Int}}
@test ss.support == [1,2,3,4,5]
@test ss.freq == [2., 3., 3., 2., 1.]

ss2 = suffstats(DiscreteNonParametric, xs, ws)
@test ss2 isa Distributions.DiscreteNonParametricStats{Int,Float64,Vector{Int}}
@test ss2.support == [1,2,3,4,5]
@test ss2.freq == [2., 11., 5., 1., 1.]

d1 = fit_mle(DiscreteNonParametric, ss)
@test d1 isa DiscreteNonParametric{Int,Float64,Vector{Int}}
@test support(d1) == ss.support
@test probs(d1) ≈ ss.freq ./ 11

d2 = fit_mle(DiscreteNonParametric, xs)
@test typeof(d2) == typeof(d1)
@test support(d2) == support(d1)
@test probs(d2) ≈ probs(d1)

d3 = fit(DiscreteNonParametric, xs)
@test typeof(d2) == typeof(d1)
@test support(d3) == support(d1)
@test probs(d3) ≈ probs(d1)

# Numerical stability; see issue #872 and PR #926
p = [1 - eps(Float32), eps(Float32)]
d = Categorical(p)
@test ([rand(d) for _ = 1:100_000]; true)

end


text:
"""
Package: DataWrangler

    function d(x::{AbstractVector, AbstractArray},
               or::Int=1,
               la::Int=1;
               center::Bool=false)

Return Lagged differences of a given Vector or Array.

# Arguments
- `x`: Vector or Array of data.
- `or`: Order of the differences; number of recursive iterations on the same vector/array.
- `la`: Lag for the difference.
- `center`: Center the result in the response using Missing values.

# Returns
Laged differences Vector or Array of a given order.

# Examples
```julia-repl
julia> x = [1,2,3,4,5];
julia> d(x)
4-element Vector{Int64}:
 1
 1
 1
 1

julia> d(x,2)
3-element Vector{Int64}:
 0
 0
 0

julia> d(x,1,2)
3-element Vector{Int64}:
 2
 2
 2

julia> x = reshape(collect(1:20),10,2);

julia> d(x,2,2)
6×2 Matrix{Int64}:
 0  0
 0  0
 0  0
 0  0
 0  0
 0  0

julia> d(d(x,1,2),1,2) == d(x,2,2)
true
```
"""
function d(x::AbstractArray{<:Union{Missing,T}},
           order::Int=1,
           lag::Int=1;
           center::Bool=false) where T<:Real

    (lag == 0) | (order == 0) && (return x)

    n = size(x,1)
    nv = size(x,2)

    @assert 0 <= lag <= (n-1) "Lag must be larger or equal to size(x,1)"

    dx = x
    for i in order:-1:1
        dx = (dx .- circshift(dx,lag))[1+lag:end,:]
    end

    if center
        nl = n-size(dx,1)
        dx = vcat(dx, Array{Any}(missing,nl,size(x,2)))
        dx = circshift(dx,(div(nl,2),0))
    end

    return size(dx,2) == 1 ? dx[:,1] : dx

end

# Fractional Difference
function d(x::AbstractVector{<:T}, d::T;
           scale::T = T(size(x,1))) where T<:AbstractFloat

    N = size(x,1)
    np2 = nextpow(2,2*N-1)
    k = 1:N-1
    b = [T(1); cumprod((k.-d.-T(1))./k)]
    pad0 = zeros(T,np2-N)
    dx1 = fft(vcat(b,pad0))
    dx2 = fft(vcat(x,pad0))
    dx = ifft(dx2 .* dx1)
    vec(real(dx[1:N,:])) .* (N/scale)^d

end


text:
##################################################################################
#
#  LIFTING TRANSFORMS
#  Periodic boundaries
#
##################################################################################

# for split and merge
reqtmplength(x::AbstractArray) = (size(x,1)>>2) + (size(x,1)>>1)%2

# return scheme parameters adjusted for direction and type
function makescheme(::Type{T}, scheme::GLS, fw::Bool) where T<:Number
    n = length(scheme.step)
    stepseq = Vector{WT.LSStep{T}}(undef, n)
    for i = 1:n
        j = fw ? i : n+1-i
        stepseq[i] = WT.LSStep(scheme.step[j].steptype,
                            convert(Vector{T}, scheme.step[j].param.coef * (fw ? -1 : 1)),
                            scheme.step[j].param.shift)
    end
    norm1, norm2 =  convert(T, fw ? scheme.norm1 : 1/scheme.norm1),
                    convert(T, fw ? scheme.norm2 : 1/scheme.norm2)
    return stepseq, norm1, norm2
end

# 1-D
# inplace transform of y, no vector allocation
function _dwt!(y::AbstractVector{T}, scheme::GLS, L::Integer, fw::Bool,
        tmp::Vector{T} = Vector{T}(undef, reqtmplength(y))) where T<:Number

    n = length(y)
    0 <= L ||
        throw(ArgumentError("L must be positive"))
    sufficientpoweroftwo(y, L) ||
        throw(ArgumentError("size must have a sufficient power of 2 factor"))
    length(tmp) >= n>>2 ||
        throw(ArgumentError("length of tmp incorrect"))

    if L == 0
        return y
    end

    if fw
        lrange = 1:L
        ns = n
    else
        lrange = L:-1:1
        ns = detailn(n, L-1)
    end
    half = ns>>1
    s = y
    stepseq, norm1, norm2 = makescheme(T, scheme, fw)

    for l in lrange
        if fw
            Util.split!(s, ns, tmp)
            for step in stepseq
                lift!(s, half, step.param, step.steptype)
            end
            normalize!(s, half, ns, norm1, norm2)
            ns = ns>>1
            half = half>>1
        else
            normalize!(s, half, ns, norm1, norm2)
            for step in stepseq
                lift!(s, half, step.param, step.steptype)
            end
            Util.merge!(s, ns, tmp)        # inverse split
            ns = ns<<1
            half = half<<1
        end
    end
    return y
end
# 1-D, unsafe simple 1 level transform
# inplace transform of y, no vector allocation
# tmp: size at least n>>2
# oopc: use oop computation, if false iy and incy are assumed to be 1
# oopv: the out of place location
function unsafe_dwt1level!(y::AbstractArray{T}, iy::Integer, incy::Integer, oopc::Bool,
    oopv::FVector{T}, scheme::GLS, fw::Bool, stepseq::FVector, norm1::T, norm2::T,
    tmp::FVector{T}) where T<:Number
    if !oopc
        oopv = y
    end
    ns = length(oopv)
    half = ns>>1

    if fw
        if oopc
            Util.split!(oopv, y, iy, incy, ns)
        else
            Util.split!(oopv, ns, tmp)
        end
        for step in stepseq
            lift!(oopv, half, step.param, step.steptype)
        end
        if oopc
            normalize!(y, iy, incy, oopv, half, ns, norm1, norm2)
        else
            normalize!(oopv, half, ns, norm1, norm2)
        end
    else
        if oopc
            normalize!(oopv, y, iy, incy, half, ns, norm1, norm2)
        else
            normalize!(oopv, half, ns, norm1, norm2)
        end
        for step in stepseq
            lift!(oopv, half, step.param, step.steptype)
        end
        if oopc
            Util.merge!(y, iy, incy, oopv, ns)
        else
            Util.merge!(oopv, ns, tmp)
        end
    end

    return y
end

# 2-D
# inplace transform of y, no vector allocation
# tmp: size at least n>>2
# tmpvec: size at least n
function _dwt!(y::Matrix{T}, scheme::GLS, L::Integer, fw::Bool, tmp::Vector{T} = Vector{T}(undef, reqtmplength(y)), tmpvec::Vector{T} = Vector{T}(undef, size(y,1))) where T<:Number

    n = size(y,1)
    iscube(y) ||
        throw(ArgumentError("array must be square/cube"))
    0 <= L ||
        throw(ArgumentError("L must be positive"))
    sufficientpoweroftwo(y, L) ||
        throw(ArgumentError("size must have a sufficient power of 2 factor"))
    (length(tmp) >= n>>2) ||
        throw(ArgumentError("length of tmp incorrect"))
    (length(tmpvec) >= n) ||
        throw(ArgumentError("length of tmpvec incorrect"))

    if L == 0
        return y
    end
    row_stride = n
    plane_stride = n*n

    if fw
        lrange = 1:L
        nsub = n
    else
        lrange = L:-1:1
        nsub = div(n,2^(L-1))
    end
    stepseq, norm1, norm2 = makescheme(T, scheme, fw)

    # transforms with stride are out of place in a dense array for speed
    for l in lrange
        tmpsub = unsafe_vectorslice(tmpvec, 1, nsub)
        if fw
            # rows
            for i in 1:nsub
                xi = row_idx(i, n)
                unsafe_dwt1level!(y, xi, row_stride, true, tmpsub, scheme, fw,
                                    stepseq, norm1, norm2, tmp)
            end
            # columns
            for i in 1:nsub
                xi = col_idx(i, n)
                ya = unsafe_vectorslice(y, xi, nsub)
                unsafe_dwt1level!(ya, 1, 1, false, tmpsub, scheme, fw,
                                    stepseq, norm1, norm2, tmp)
            end
        else
            # columns
            for i in 1:nsub
                xi = col_idx(i, n)
                ya = unsafe_vectorslice(y, xi, nsub)
                unsafe_dwt1level!(ya, 1, 1, false, tmpsub, scheme, fw,
                                    stepseq, norm1, norm2, tmp)
            end
            # rows
            for i in 1:nsub
                xi = row_idx(i, n)
                unsafe_dwt1level!(y, xi, row_stride, true, tmpsub, scheme, fw,
                                    stepseq, norm1, norm2, tmp)
            end
        end

        nsub = (fw ? nsub>>1 : nsub<<1)
    end

    return y
end

# 3-D
# inplace transform of y, no vector allocation
# tmp: size at least n>>2
# tmpvec: size at least n
function _dwt!(y::Array{T,3}, scheme::GLS, L::Integer, fw::Bool, tmp::Vector{T} = Vector{T}(undef, reqtmplength(y)), tmpvec::Vector{T} = Vector{T}(undef, size(y,1))) where T<:Number

    n = size(y,1)
    iscube(y) ||
        throw(ArgumentError("array must be square/cube"))
    0 <= L ||
        throw(ArgumentError("L must be positive"))
    sufficientpoweroftwo(y, L) ||
        throw(ArgumentError("size must have a sufficient power of 2 factor"))
    (length(tmp) >= n>>2) ||
        throw(ArgumentError("length of tmp incorrect"))
    (length(tmpvec) >= n) ||
        throw(ArgumentError("length of tmpvec incorrect"))

    if L == 0
        return y
    end
    row_stride = n
    plane_stride = n*n

    if fw
        lrange = 1:L
        nsub = n
    else
        lrange = L:-1:1
        nsub = div(n,2^(L-1))
    end
    stepseq, norm1, norm2 = makescheme(T, scheme, fw)

    # transforms with stride are out of place in a dense array for speed
    for l in lrange
        tmpsub = unsafe_vectorslice(tmpvec, 1, nsub)
        if fw
            # planes
            for i in 1:nsub, j in 1:nsub
                xi = plane_idx(i, j, n)
                unsafe_dwt1level!(y, xi, plane_stride, true, tmpsub, scheme, fw,
                                    stepseq, norm1, norm2, tmp)
            end
            # rows
            for i in 1:nsub, j in 1:nsub
                xi = row_idx(i, j, n)
                unsafe_dwt1level!(y, xi, row_stride, true, tmpsub, scheme, fw,
                                    stepseq, norm1, norm2, tmp)
            end
            # columns
            for i in 1:nsub, j in 1:nsub
                xi = col_idx(i, j, n)
                ya = unsafe_vectorslice(y, xi, nsub)
                unsafe_dwt1level!(ya, 1, 1, false, tmpsub, scheme, fw,
                                    stepseq, norm1, norm2, tmp)
            end
        else
            # columns
            for i in 1:nsub, j in 1:nsub
                xi = col_idx(i, j, n)
                ya = unsafe_vectorslice(y, xi, nsub)
                unsafe_dwt1level!(ya, 1, 1, false, tmpsub, scheme, fw,
                                    stepseq, norm1, norm2, tmp)
            end
            # rows
            for i in 1:nsub, j in 1:nsub
                xi = row_idx(i, j, n)
                unsafe_dwt1level!(y, xi, row_stride, true, tmpsub, scheme, fw,
                                    stepseq, norm1, norm2, tmp)
            end
            # planes
            for i in 1:nsub, j in 1:nsub
                xi = plane_idx(i, j, n)
                unsafe_dwt1level!(y, xi, plane_stride, true, tmpsub, scheme, fw,
                                    stepseq, norm1, norm2, tmp)
            end
        end

        nsub = (fw ? nsub>>1 : nsub<<1)
    end

    return y
end

# WPT
# 1-D
function _wpt!(y::AbstractVector{T}, scheme::GLS, tree::BitVector, fw::Bool, tmp::Vector{T} = Vector{T}(undef, reqtmplength(y))) where T<:Number

    n = length(y)
    isvalidtree(y, tree) ||
        throw(ArgumentError("invalid tree"))
    length(tmp) >= n>>2 ||
        throw(ArgumentError("length of tmp incorrect"))

    if !tree[1]
        return y
    end

    stepseq, norm1, norm2 = makescheme(T, scheme, fw)

    Lmax = maxtransformlevels(n)
    L = Lmax
    while L > 0
        ix = 1
        k = 1
        fw  && (Lfw = Lmax-L)
        !fw && (Lfw = L-1)
        nj = detailn(n, Lfw)
        treeind = 2^(Lfw)-1

        while ix <= n
            if tree[treeind+k]
                dy = unsafe_vectorslice(y, ix, nj)
                unsafe_dwt1level!(dy, 1, 1, false, tmp, scheme, fw, stepseq, norm1, norm2, tmp)
            end
            ix += nj
            k += 1
        end
        L -= 1
    end

    return y
end

function normalize!(x::AbstractVector{T}, half::Int, ns::Int, n1::T, n2::T) where T<:Number
    for i = 1:half
        @inbounds x[i] *= n1
    end
    for i = half+1:ns
        @inbounds x[i] *= n2
    end
    return x
end
# out of place normalize from x to y
function normalize!(y::AbstractVector{T}, x::AbstractVector{T}, half::Int, ns::Int, n1::T, n2::T) where T<:Number
    for i = 1:half
        @inbounds y[i] = n1*x[i]
    end
    for i = half+1:ns
        @inbounds y[i] = n2*x[i]
    end
    return y
end
function normalize!(y::AbstractArray{T}, iy::Int, incy::Int, x::AbstractVector{T}, half::Int, ns::Int, n1::T, n2::T) where T<:Number
    for i = 1:half
        @inbounds y[iy + (i-1)*incy] = n1*x[i]
    end
    for i = half+1:ns
        @inbounds y[iy + (i-1)*incy] = n2*x[i]
    end
    return y
end
function normalize!(y::AbstractVector{T}, x::AbstractArray{T}, ix::Int, incx::Int, half::Int, ns::Int, n1::T, n2::T) where T<:Number
    for i = 1:half
        @inbounds y[i] = n1*x[ix + (i-1)*incx]
    end
    for i = half+1:ns
        @inbounds y[i] = n2*x[ix + (i-1)*incx]
    end
    return y
end

# predict and update lifting steps inplace on x, forward and backward
# half: half of the length under consideration
# For predict: writes to range 1:half, reads from 1:2*half
# For update : writes to range half+1:2*half, reads from 1:2*half
for step_type in (WT.PredictStep, WT.UpdateStep)
@eval begin
function lift!(x::AbstractVector{T}, half::Int,
                param::WT.LSStepParam{T}, steptype::$step_type) where T<:Number
    lhsr, irange, rhsr, rhsis = getliftranges(half, length(param), param.shift, steptype)
    coefs = param.coef
    # left boundary
    lift_perboundary!(x, half, coefs, lhsr, rhsis, steptype)
    # main loop
    lift_inbounds!(x, coefs, irange, rhsis)
    # right boundary
    lift_perboundary!(x, half, coefs, rhsr, rhsis, steptype)
    return x
end
end # eval begin
end # for

function irlimits(half::Int, nc::Int, shift::Int)
    # conditions for every element i in irange to be in bounds
    # 1 <= i <= half
    # 1 <= i+1-1-shift <= half
    # 1 <= i+nc-1-shift <= half
    return (max(shift+1, 1-nc+shift),
            min(half+1+shift-nc, half+shift))
end

function getliftranges(half::Int, nc::Int, shift::Int, steptype::WT.UpdateStep)
    # define index shift rhsis
    rhsis = -shift-half
    irmin, irmax = irlimits(half, nc, shift)
    if irmin > half || irmax < 1
        irange = 1:0  # empty
    else
        irmin = max(irmin, 1)
        irmax = min(irmax, half)
        irange = irmin + half:irmax + half
    end
    # periodic boundary
    if length(irange)==0
        lhsr = 1 + half:half + half
        rhsr = 1 + half:0 + half
    else
        lhsr = 1 + half:irmin - 1 + half
        rhsr = irmax + 1 + half:half + half
    end
    return (lhsr, irange, rhsr, rhsis)
end

function getliftranges(half::Int, nc::Int, shift::Int, steptype::WT.PredictStep)
    # define index shift rhsis
    rhsis = -shift+half
    irmin, irmax = irlimits(half, nc, shift)
    if irmin > half || irmax < 1
        irange = 1:0  # empty
    else
        irmin = max(irmin, 1)
        irmax = min(irmax, half)
        irange = irmin:irmax
    end
    # periodic boundary
    if length(irange)==0
        lhsr = 1:half
        rhsr = 1:0
    else
        lhsr = 1:irmin - 1
        rhsr = irmax + 1:half
    end
    return (lhsr, irange, rhsr, rhsis)
end

# periodic boundary
for (step_type, puxind) in ((WT.PredictStep, :(mod1(i+k-1+rhsis-half,half)+half)),
                            (WT.UpdateStep,  :(mod1(i+k-1+rhsis,half))) )
@eval begin
function lift_perboundary!(x::AbstractVector{T}, half::Int,
                            c::Vector{T}, irange::AbstractRange, rhsis::Int, ::$step_type) where T<:Number
    nc = length(c)
    for i in irange
        for k in 1:nc
            @inbounds x[i] += c[k]*x[$puxind]
        end
    end
    return x
end
end # eval begin
end # for

# main lift loop
function lift_inbounds!(x::AbstractVector{T}, c::Vector{T}, irange::AbstractRange, rhsis::Int) where T<:Number
    nc = length(c)
    if nc == 1  # hard code the most common cases (1, 2, 3) for speed
        c1 = c[1]
        for i in irange
            @inbounds x[i] += c1*x[i+rhsis]
        end
    elseif nc == 2
        c1,c2 = c[1],c[2]
        rhsisp1 = rhsis+1
        for i in irange
            @inbounds x[i] += c1*x[i+rhsis] + c2*x[i+rhsisp1]
        end
    elseif nc == 3
        c1,c2,c3 = c[1],c[2],c[3]
        rhsisp1 = rhsis+1
        rhsisp2 = rhsis+2
        for i = irange
            @inbounds x[i] += c1*x[i+rhsis] + c2*x[i+rhsisp1] + c3*x[i+rhsisp2]
        end
    else
        for i in irange
            for k in 0:nc-1
                @inbounds x[i] += c[k]*x[i+k+rhsis]
            end
        end
    end
    return x
end


text:
export
    Switch,
    LinearSwitch,
    If,
    choose

"""
    abstract type Switch{N, I, K, O} <: Det{K, O}

`Switch` defines an sfunc that represents choosing between multiple incoming (parent) Sfuncs
based on a test.  A subtype of `Switch` must provide a `choose` function that takes the switch
and an `i` and returns an integer between `1` and `N`.  This is an index into a 'parent array'.

`K` must be a flat tuple type consisting of `I` and `N` occurrences of `O`:  for example,
if I is `Int`

```
    K = extend_tuple_type(Tuple{Int}, NTuple{N, O})
```

If the subtype'd sfunc is not in the Scruff.SFuncs module, the system must
`import Scruff.SFuncs: choose`.

# Additional supported operators
- `support`
- `support_quality`
- `compute_pi`
- `send_lambda`

# Type parameters
- `N`: the count of *sfuncs* from which to choose
- `I`: the type of the second argument of the `choose` function defined for the `Switch`
- `K`: the input type of the `Switch`; see above
- `O`: the output type of the `Switch`

See also:  [`choose`](@ref), [`extend_tuple_type`](@ref)
"""
abstract type Switch{N, I, K, O} <: Det{K, O} end

"""
`choose` interface.  For every subtype of `Switch`, an implementation of this method must
be created, whose first parameter is the subtype, and the second parameter is of type `I`
for the parameter type in `Switch`.

For example, the definition of the `If` *sfunc* is as follows, where `choose` returns
either index `1` or index `2`.

```
struct If{O} <: Switch{2, Bool, Tuple{Bool, O, O}, O} end
choose(::If, b::Bool) = b ? 1 : 2
```
"""
function choose end

# apply(sw::Switch, i, hs...) = hs[choose(sw,i)] # implemented because subtype of `Det`
apply(sw::Switch, i, hs...) = hs[choose(sw,i)] # implemented because subtype of `Det`

struct LinearSwitch{N, K, O} <: Switch{N, Int, K, O}
    n :: Int
    function LinearSwitch(N, O)
        K = extend_tuple_type(Tuple{Int}, NTuple{N, O})
        new{N, K, O}(N)
    end
end

choose(::LinearSwitch, i) = i

struct If{O} <: Switch{2, Bool, Tuple{Bool, O, O}, O} end

choose(::If, b) = b ? 1 : 2

# switch overrides the functions for det that use the Cartesian product of the parents.
# Since only one parent is ever active at a time, we don't need to compute the Cartesian product.
# This leads to exponential time savings in the number of parents.
@impl begin
    struct SwitchSupport end
    function support(sf::Switch{N,I,K,O},
                    parranges::NTuple{M,Vector},
                    size::Integer,
                    curr::Vector{<:O}) where {I,K,O,N,M}
        result = Vector{output_type(sf)}()
        for i in parranges[1]
            h = choose(sf, i)
            append!(result, parranges[h + 1])
        end
        return unique(result)
    end
end

@impl begin
    struct SwitchSupportQuality end

    function support_quality(sf::LinearSwitch, parranges)
        ivals = map(i -> choose(sf, i), parranges[1])
        if all(j -> j in ivals, 1:sf.n)
            return :CompleteSupport
        else
            return :BestEffortSupport
        end
    end

    function support_quality(::Union{<:Support,Nothing}, sf::If, parranges)
        ivals = map(i -> choose(sf, i), parranges[1])
        if all(j -> j in ivals, 1:2)
            return :CompleteSupport
        else
            return :BestEffortSupport
        end
    end
end

# We currently don't overwrite bounded_probs.
# bounded_probs will be replaced with an operation that produces factors.

@impl begin
    struct SwitchComputePi end

    function compute_pi(sf::Switch{N,I,K,O},
                     range::__OptVec{<:O},
                     parranges::NTuple{M,Vector},
                     incoming_pis::Tuple)::Dist{<:O} where {M,N,I,K,O}
        result = zeros(Float64, length(range))
        ipis = incoming_pis[1]
        for (i,ival) in enumerate(parranges[1])
            ipi = cpdf(ipis, (), parranges[1][i])
            h = choose(sf, ival)
            hpis = incoming_pis[h+1]
            hrange = parranges[h+1]
            for (j,jval) in enumerate(range)
                k = indexin([jval], hrange)[1]
                # the range of the parent can be a subset of the range of the switch
                hpi = isnothing(k) ? 0.0 : cpdf(hpis, (), parranges[h+1][k])
                result[j] += ipi * hpi
            end
        end
        return Cat(range, result)
    end
end

@impl begin
    struct SwitchSendLambda end

    function send_lambda(sf::Switch{N,I,K,O},
                       lambda::Score{<:O},
                       range::__OptVec{<:O},
                       parranges::NTuple{M,Vector},
                       incoming_pis::Tuple,
                       parent_ix::Integer)::Score where {M,N,I,K,O}

        # This helper function computes, for a particular choice, the sum of
        # pi times lambda values.
        function compute1(h)
            hrange = parranges[h+1]
            hpis = incoming_pis[h+1]
            tot = 0.0
            for (j,jval) in enumerate(hrange)
                # We need to keep track of indices properly.
                # We cannot assume that the range of the switch and the range passed into
                # send_lambda are the same, even though they are of the same type O.
                # In particular, choices might have different ranges from each other,
                # perhaps restricted or in a different order.
                k = indexin([jval], range)[1]
                tot += cpdf(hpis, (), parranges[h+1][j]) * get_score(lambda, range[k])
            end
            return tot
        end

        # We need to make sure the message is correctly typed to the output type of the appropriate parent
        # Need to make sure the target parent range is a Vector{T} rather than a Vector{Any}
        T = typeof(parranges[parent_ix][1])
        target_parrange :: Vector{T} = parranges[parent_ix]
        if parent_ix == 1
            # We send a lambda message to the selector parent
            # For a given value of the selector, the lambda value is the
            # probability of the child lambda for the corresponding choice.
            # This is equal to the sum, over all values in the range of the
            # choice, of the incoming_pi * lambda for that value.
            return SoftScore(target_parrange,
                map(i -> compute1(choose(sf, i)), target_parrange))

        else
            # The lambda message to a choice consists of two components:
            # A choice-specific component equal to the pi of the I value that leads to that choice, times lambda
            # A constant from all the other choices, equal to the sum, over all choices, of the pi of the I value for the choice
            # times the sum of pi * lambda values for that choice.
            con = 0.0
            spec = nothing
            ipis = incoming_pis[1]
            for (i,ival) in enumerate(parranges[1])
                h = choose(sf, ival)
                if h+1 != parent_ix
                    con += cpdf(ipis, (), ival)  * compute1(h)
                else
                    # lambda is in the order of range, but hrange might be a subset and in a different order
                    hrange = parranges[h+1]
                    qs = []
                    for (j,jval) in enumerate(hrange)
                        # k = indexin([jval], range)[1]
                        # push!(qs, lambda[k])
                        push!(qs, get_score(lambda, jval))
                    end
                    spec = cpdf(ipis, (), ival) .* qs
                end
            end

            return SoftScore(target_parrange, con .+ spec)
        end
    end
end


text:
module YaoTargetQASM

using MLStyle
using OpenQASM
using OpenQASM.Types
using OpenQASM.Tools
using YaoCompiler
using YaoLocations
using Expronicon
using Configurations
using CompilerPluginTools
using RBNF: Token
using YaoCompiler.Intrinsics
using YaoCompiler: YaoCompileTarget
using Core: CodeInfo
using YaoLocations: plain
using OpenQASM.Types: UGate

export OpenQASMTarget, RegInfo, code_qasm, QASM, @qasm_str

@option struct OpenQASMTarget <: YaoCompileTarget
    version::VersionNumber = v"2.0"
    toplevel::Bool = true
    inline_intrinsics::Bool = false
    include_qelib1::Bool = false
end

module QASM

Base.@pure @noinline sin(x) = Base.Math.sin(x)
Base.@pure @noinline cos(x) = Base.Math.cos(x)
Base.@pure @noinline tan(x) = Base.Math.tan(x)
Base.@pure @noinline exp(x) = Base.Math.exp(x)
Base.@pure @noinline log(x) = Base.Math.log(x)
Base.@pure @noinline sqrt(x) = Base.Math.sqrt(x)

end

function code_qasm(op::Type{<:Operation}; kw...)
    target = OpenQASMTarget(;kw...)
    return YaoCompiler.compile(target, Intrinsics.apply, Tuple{AnyReg, op})
end

include("validate.jl")
include("intrinsic.jl")
include("reginfo.jl")
include("codegen.jl")
include("frontend.jl")

end


text:
using JuMP, ForwardDiff, Ipopt
using EAGO
using IntervalArithmetic
using DelimitedFiles
using Statistics
using LatinHypercubeSampling

function Ammonia_oxidation_rate(x)
    c_NH = x[1]
    c_O = x[2]
    r_AOmax = 0.67 # mg N-NH+ /g VSS_AO min
    K_OAO = 0.3 # mg O2/L
    K_SAO = 0.34 # mg N-NH4+/L
    K_IAO = 6200 # mg N-NH4+/L
    r_AO = r_AOmax*(c_NH/(K_SAO + c_NH + c_NH^2/K_IAO))*(c_O/(K_OAO + c_O)); # mg N-NH+ /g VSS_AO min
    return r_AO
end

##
#Data_LHC = randomLHC(100000,2)
#Data_input = scaleLHC(Data_LHC,[(0.0,40.0),(0.0,9.1)])
#Data_output = mapslices(Ammonia_oxidation_rate,Data_input; dims=2)
#out_r_AO = [Data_input Data_output]
#writedlm("TrainingData.csv",out_r_AO,',')

Data_LHC1 = randomLHC(10000,2)
Data_input1 = scaleLHC(Data_LHC1,[(0.0,4.0),(0.0,1.0)])
Data_output1 = mapslices(Ammonia_oxidation_rate, Data_input1; dims=2)
out_r_AO1 = [Data_input1 Data_output1]

Data_LHC2 = randomLHC(10000,2)
Data_input2 = scaleLHC(Data_LHC2,[(0.0,40.0),(0.0,9.1)])
Data_output2 = mapslices(Ammonia_oxidation_rate, Data_input2; dims=2)
out_r_AO2 = [Data_input2 Data_output2]

out_r_AO = [out_r_AO1; out_r_AO2]

writedlm("TrainingData2.csv",out_r_AO,',')

max_stars_repo_name: PSORLab/RobustHybridModels

max_issues_repo_name: PSORLab/RobustHybridModels

max_forks_repo_name: PSORLab/RobustHybridModels


text:
module EGraphs

include("../docstrings.jl")

using RuntimeGeneratedFunctions
RuntimeGeneratedFunctions.init(@__MODULE__)

# import ..Rule
# import ..getrhsfun

using ..TermInterface
using ..Util
using ..Rules
import ..@log

include("enode.jl")
export ENode

include("abstractanalysis.jl")
export AbstractAnalysis

include("eclass.jl")
export EClass
export hasdata
export getdata
export setdata!

include("egg.jl")
export find
export geteclass
export arity
export EGraph
export merge!
export in_same_class
export addexpr!
export addanalysis!
export rebuild!
export prune!

include("analysis.jl")
export analyze!

include("ematch.jl")
include("Schedulers/Schedulers.jl")
export Schedulers
using .Schedulers

include("saturation_reason.jl")
export ReportReasons
include("saturation_report.jl")
include("saturation_params.jl")
export SaturationParams

include("saturation_search.jl")
include("saturation_apply.jl")

include("saturation.jl")
export saturate!
include("equality.jl")
export areequal
export @areequal
export @areequalg

include("extraction.jl")
export extract!
export ExtractionAnalysis
export astsize
export astsize_inv
export @extract

end


text:
function density_gradient!(∇ρ, point, wfn)

    ρ    = 0.0
   ∇ρ   .= 0.0
   ϕ, ∇ϕ = eval_mo1(wfn, point)

   occ = wfn.occ

   # Run again over orbitals
   @inbounds for imo in eachindex(occ)

      @views begin
         n    = occ[imo]
         nϕ   = n   *  ϕ[imo]
          ρ  += nϕ  *  ϕ[imo]
         ∇ρ .+= nϕ .* ∇ϕ[imo,:]
      end
   end
   ∇ρ .*= 2 # + c.c.
   ∇ρ, ρ
end

function eval_mo1(wfn, point) #derivative=1)

   nmo    = wfn.nmo
   natm   = wfn.natm
   coords = wfn.coords
   ngroup = wfn.ngroup
   nuexp  = wfn.nuexp
   rcutte = wfn.rcutte
   oexp   = wfn.oexp
   nzexp  = wfn.nzexp
   ityp   = wfn.ityp
   nlm    = wfn.nlm
   coef   = wfn.coef

    ϕ   = zeros(nmo)
   ∇ϕ   = zeros(nmo,3)
    G1D = @MVector zeros(3) # Gaussian 1D functions x, y, z
   ∇G1D = @MVector zeros(3)

   @views begin
      @inbounds for ic in 1:natm

         # Atomic coordinates of this center
         xcoor = point - coords[ic,:]
         r² = sum( abs2.(xcoor) ) # ≡ dot(xcoor, xcoor)
         # Loop over different shell in this atom
         @inbounds for m in 1:ngroup[ic]

            k = nuexp[ic,m,1]
            # Skip to compute this primitive if distance is too big.
            r² > rcutte[ic,m]^2 && continue

            α = -oexp[k]
            𝟐α = 2α
            # All primitives in a shell share the same exponent.
            aexp = exp(α * r²)
            # Loop over the different primitives in this shell.
            @inbounds for jj in 1:nzexp[ic,m]

               # "i" is the original index of the primitive in the WFN.
               i = nuexp[ic,m,jj]
               itip = ityp[i]
               # Integer coefficients.
               it = nlm[itip,:]
               @inbounds for j in 1:3

                  n = it[j]
                  # r⃗ = ∑ⱼ³ rⱼ e⃗ⱼ # cartesian components
                  rⱼ = xcoor[j]
                  # G1D = r^n exp(α * r²)
                  if n == 0
                     # ∇G1D = 2αr * G1D
                     ∇G1D[j] = 𝟐α * rⱼ
                      G1D[j] = 1.0
                  elseif n == 1
                     ∇G1D[j] = 1 + 𝟐α * rⱼ * rⱼ
                      G1D[j] = rⱼ
                  elseif n == 2
                     rⱼ²   = rⱼ * rⱼ
                     ∇G1D[j] = rⱼ * ( 2.0 + 𝟐α * rⱼ² )
                      G1D[j] = rⱼ²
                  elseif n == 3
                     rⱼ²   = rⱼ  * rⱼ
                     ∇G1D[j] = rⱼ² * ( 3.0 + 𝟐α * rⱼ² )
                      G1D[j] = rⱼ  * rⱼ²
                  elseif n == 4
                     rⱼ²   = rⱼ  * rⱼ
                     ∇G1D[j] = rⱼ² * rⱼ * ( 4.0 + 𝟐α  * rⱼ² )
                      G1D[j] = rⱼ² * rⱼ²
                  elseif n == 5
                     rⱼ²   = rⱼ  * rⱼ
                     ∇G1D[j] = rⱼ² * rⱼ² * ( 5.0 + 𝟐α * rⱼ² )
                      G1D[j] = rⱼ² * rⱼ² * rⱼ
                  end
               end

               # exp(-αr₁²) * exp(-αr₂²) * exp(-αr₃²) = exp(-αr²) ≡ aexp
                 G3D =  prod(G1D) * aexp # Gxyz = Gx * Gy * Gz
               ∇xG3D = ∇G1D[1] *  G1D[2] *  G1D[3] * aexp
               ∇yG3D =  G1D[1] * ∇G1D[2] *  G1D[3] * aexp
               ∇zG3D =  G1D[1] *  G1D[2] * ∇G1D[3] * aexp

               # AO -> MO basis; run over orbitals
               @inbounds for imo in 1:nmo

                  c = coef[imo,i]
                   ϕ[imo]   += c *   G3D
                  ∇ϕ[imo,1] += c * ∇xG3D
                  ∇ϕ[imo,2] += c * ∇yG3D
                  ∇ϕ[imo,3] += c * ∇zG3D
               end
            end
         end
      end
   end
   ϕ, ∇ϕ
end


text:
struct TrainValidTest <: Split
    at::NTuple{2, Float64}

    TrainValidTest(at::NTuple{2, Real} = (0.5, 0.25)) = new(at)
end

function load(
    split::TrainValidTest,
    dataset::Dataset{N, F, T},
) where {N<:Name, F<:Format, T<:Task}

    at = split.at
    train = load(dataset, :train; dopostprocess = false)
    obsdim = samplesdim(F)

    if !hassubset(N, :valid) && !hassubset(N, :test)
        train, valid, test = data_split(dataset, train, obsdim, at)
    elseif !hassubset(N, :test)
        train, test = data_split(dataset, train, obsdim, sum(at))
        valid = load(dataset, :valid; dopostprocess = false)
    elseif !hassubset(N, :valid)
        train, valid = data_split(dataset, train, obsdim, sum(at))
        test = load(dataset, :test; dopostprocess = false)
    else
        valid = load(dataset, :valid; dopostprocess = false)
        test = load(dataset, :test; dopostprocess = false)
    end

    train = postprocess(dataset, train)
    valid = postprocess(dataset, valid)
    test = postprocess(dataset, test)

    return train, valid, test
end

load(dataset::Dataset) = load(TrainValidTest(), dataset)

max_stars_repo_name: VaclavMacha/DatasetProvider

max_issues_repo_name: VaclavMacha/DatasetProvider

max_forks_repo_name: VaclavMacha/DatasetProvider


text:
#############################################################################
#Copyright (c) 2016 Charles Le Losq
#
#The MIT License (MIT)
#
#Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the #Software without restriction, including without limitation the rights to use, copy, #modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, #and to permit persons to whom the Software is furnished to do so, subject to the #following conditions:
#
#The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
#
#THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, #INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
#
# baseline.jl contains several "baseline" functions. It is directly dependent on functions.jl.
#
#############################################################################

using JuMP
using Ipopt

function baseline(x::Array{Float64},y::Array{Float64},roi::Array{Float64},basetype::AbstractString,p::Array{Float64})
    # First we grab the good roi
    interest_index::Array{Int64} = find(roi[1,1] .<= x[:,1] .<= roi[1,2])
    if size(roi)[1] > 1
        for i = 2:size(roi)[1]
            interest_index = vcat(interest_index,  find(roi[i,1] .<= x[:,1] .<= roi[i,2]))
        end
    end
    interest_x = x[interest_index,1]
    interest_y = y[interest_index,1]

    if basetype == "poly"
        # The model for fitting baseline to roi signal
        mod = Model(solver=IpoptSolver(print_level=0))
        n::Int = size(p)[1] # number of coefficients
        m::Int = size(interest_x)[1] # number of data

        @defVar(mod,p_val[i=1:n])
        setValue(p_val[i=1:n], p[i])
        @setNLObjective(mod,Min,sum{( sum{p_val[i]*interest_x[j]^(i-1), i=1:n} - interest_y[j])^2, j=1:m})
        status = solve(mod)
        println("Solver status: ", status)
        best_p::Vector{Float64} = getValue(p_val)
        y_calc::Array{Float64,2} = poly(best_p,x)
        return y[:,1] - y_calc, y_calc# To be continued... Fitting procedure to add there.
    else
        error("Not implemented, choose between poly and [to come]")
    end
end


text:
module TestTypesConversions
using Test
using EcologicalNetworks

N = BipartiteNetwork([true false; true true], [:a, :b], [:A, :B])
@test_throws AssertionError has_interaction(N, :A, :B)
M = convert(UnipartiteNetwork, N)
@test has_interaction(M, :a, :A)
@test !has_interaction(M, :a, :B)
@test has_interaction(M, :b, :A)
@test has_interaction(M, :b, :B)
@test !has_interaction(M, :A, :B)

N = BipartiteQuantitativeNetwork([4 0 ; 2 1], [:a, :b], [:A, :B])
M = convert(UnipartiteQuantitativeNetwork, N)
@test has_interaction(M, :a, :A)
@test !has_interaction(M, :a, :B)
@test has_interaction(M, :b, :A)
@test has_interaction(M, :b, :B)
@test !has_interaction(M, :A, :B)

N = BipartiteProbabilisticNetwork([0.3 0 ; 0.1 0.9], [:a, :b], [:A, :B])
M = convert(UnipartiteProbabilisticNetwork, N)
@test has_interaction(M, :a, :A)
@test !has_interaction(M, :a, :B)
@test has_interaction(M, :b, :A)
@test has_interaction(M, :b, :B)
@test !has_interaction(M, :A, :B)
L = convert(AbstractUnipartiteNetwork, N)
@test has_interaction(L, :a, :A)
@test typeof(L) <: AbstractUnipartiteNetwork
@test typeof(L) <: UnipartiteProbabilisticNetwork

N = BipartiteQuantitativeNetwork([4 0 ; 2 1], [:a, :b], [:A, :B])
M = convert(BipartiteNetwork, N)
@test has_interaction(M, :a, :A)
@test !has_interaction(M, :a, :B)
@test has_interaction(M, :b, :A)
@test has_interaction(M, :b, :B)
L = convert(BinaryNetwork, N)
@test has_interaction(L, :a, :A)
@test typeof(L) <: BipartiteNetwork

N = UnipartiteQuantitativeNetwork([4 0 ; 2 1], [:a, :b])
M = convert(UnipartiteNetwork, N)
@test has_interaction(M, :a, :a)
@test !has_interaction(M, :a, :b)
@test has_interaction(M, :b, :a)
@test has_interaction(M, :b, :b)

N = BipartiteProbabilisticNetwork([1.0 1.0; 1.0 1.0])
M = convert(UnipartiteProbabilisticNetwork, N)
@test typeof(M) <: UnipartiteProbabilisticNetwork
@test M[3,1] == 0.0
@test M[1,4] == 1.0

N = BipartiteNetwork([true true; true true])
M = convert(UnipartiteNetwork, N)
@test typeof(M) <: UnipartiteNetwork

N = BipartiteProbabilisticNetwork([0.2 0.4; 0.4 0.2])
@test typeof(N<0.2) <: BipartiteNetwork
@test links(N<0.2) == 0
@test typeof(N<=0.2) <: BipartiteNetwork
@test links(N<=0.2) == 2
@test typeof(N>0.4) <: BipartiteNetwork
@test links(N>0.4) == 0
@test typeof(N>=0.4) <: BipartiteNetwork
@test links(N>=0.4) == 2

N = BipartiteQuantitativeNetwork([0.2 0.4; 0.4 0.2])
@test typeof(N<0.2) <: BipartiteNetwork
@test links(N<0.2) == 0
@test typeof(N<=0.2) <: BipartiteNetwork
@test links(N<=0.2) == 2
@test typeof(N>0.4) <: BipartiteNetwork
@test links(N>0.4) == 0
@test typeof(N>=0.4) <: BipartiteNetwork
@test links(N>=0.4) == 2

N = UnipartiteProbabilisticNetwork([0.2 0.4; 0.4 0.2])
@test typeof(N<0.2) <: UnipartiteNetwork
@test links(N<0.2) == 0
@test typeof(N<=0.2) <: UnipartiteNetwork
@test links(N<=0.2) == 2
@test typeof(N>0.4) <: UnipartiteNetwork
@test links(N>0.4) == 0
@test typeof(N>=0.4) <: UnipartiteNetwork
@test links(N>=0.4) == 2

N = UnipartiteQuantitativeNetwork([0.2 0.4; 0.4 0.2])
@test typeof(N<0.2) <: UnipartiteNetwork
@test links(N<0.2) == 0
@test typeof(N<=0.2) <: UnipartiteNetwork
@test links(N<=0.2) == 2
@test typeof(N>0.4) <: UnipartiteNetwork
@test links(N>0.4) == 0
@test typeof(N>=0.4) <: UnipartiteNetwork
@test links(N>=0.4) == 2

end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Types and structures for the text backend.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export TextFormat, Highlighter

"""
    TextFormat

# Fields

* `up_right_corner`: Character in the up right corner.
* `up_left_corner`: Character in the up left corner.
* `bottom_left_corner`: Character in the bottom left corner.
* `bottom_right_corner`: Character in the bottom right corner.
* `up_intersection`: Character in the intersection of lines in the up part.
* `left_intersection`: Character in the intersection of lines in the left part.
* `right_intersection`: Character in the intersection of lines in the right
                        part.
* `middle_intersection`: Character in the intersection of lines in the middle of
                         the table.
* `bottom_intersection`: Character in the intersection of the lines in the
                         bottom part.
* `column`: Character in a vertical line inside the table.
* `left_border`: Character used as the left border.
* `right_border`: Character used as the right border.
* `row`: Character in a horizontal line inside the table.
* `hlines`: Horizontal lines that must be drawn by default.
* `vlines`: Vertical lines that must be drawn by default.

# Pre-defined formats

The following pre-defined formats are available: `unicode` (**default**),
`mysql`, `compact`, `markdown`, `simple`, `ascii_rounded`, and `ascii_dots`.

"""
@kwdef struct TextFormat
    up_right_corner::Char                 = '┐'
    up_left_corner::Char                  = '┌'
    bottom_left_corner::Char              = '└'
    bottom_right_corner::Char             = '┘'
    up_intersection::Char                 = '┬'
    left_intersection::Char               = '├'
    right_intersection::Char              = '┤'
    middle_intersection::Char             = '┼'
    bottom_intersection::Char             = '┴'
    column::Char                          = '│'
    row::Char                             = '─'
    hlines::Vector{Symbol}                = [:begin,:header,:end]
    vlines::Union{Symbol,Vector{Symbol}}  = :all
end

"""
    Highlighter

Defines the default highlighter of a table when using the text backend.

# Fields

* `f`: Function with the signature `f(data,i,j)` in which should return `true`
       if the element `(i,j)` in `data` must be highlighter, or `false`
       otherwise.
* `fd`: Function with the signature `f(h,data,i,j)` in which `h` is the
        highlighter. This function must return the `Crayon` to be applied to the
        cell that must be highlighted.
* `crayon`: The `Crayon` to be applied to the highlighted cell if the default
            `fd` is used.

# Remarks

This structure can be constructed using three helpers:

    Highlighter(f::Function; kwargs...)

where it will construct a `Crayon` using the keywords in `kwargs` and apply it
to the highlighted cell,

    Highlighter(f::Function, crayon::Crayon)

where it will apply the `crayon` to the highlighted cell, and

    Highlighter(f::Function, fd::Function)

where it will apply the `Crayon` returned by the function `fd` to the
highlighted cell.

"""
@kwdef struct Highlighter
    f::Function
    fd::Function = (h,data,i,j)->h.crayon

    # Private
    crayon::Crayon = Crayon()
end

Highlighter(f; kwargs...) = Highlighter(f = f, crayon = Crayon(;kwargs...))
Highlighter(f, crayon::Crayon) = Highlighter(f = f, crayon = crayon)
Highlighter(f::Function, fd::Function) = Highlighter(f = f, fd = fd)

"""
    Display

Store the information of the display and the current cursor position. Notice
that this is not the real cursor position with respect to the display, but with
respect to the point in which the table is printed.

# Fields

* `size`: Display size.
* `row`: Current row.
* `col`: Current column.
* `has_color`: Indicates if the display has color support.
* `cont_char`: The character that indicates the line is cropped.
* `cont_space_char`: Space character to be printed before `cont_char`.

"""
@kwdef mutable struct Display
    size::Tuple{Int,Int}  = (-1,-1)
    row::Int              = 1
    col::Int              = 0
    has_color::Bool       = false
    cont_char::Char       = '⋯'
    cont_space_char::Char = ' '

    # Buffer that will store the current line.
    buf_line = IOBuffer()
end

################################################################################
#                                  Constants
################################################################################

# Crayon used to reset all the styling.
const _default_crayon = Crayon()
const _reset_crayon   = Crayon(reset = true)


text:
# Autogenerated wrapper script for Gdbm_jll for armv7l-linux-musleabihf
export gdbm_dump, gdbm_load, gdbmtool, libgdbm

using Libiconv_jll
using Readline_jll
JLLWrappers.@generate_wrapper_header("Gdbm")
JLLWrappers.@declare_executable_product(gdbm_dump)
JLLWrappers.@declare_executable_product(gdbm_load)
JLLWrappers.@declare_executable_product(gdbmtool)
JLLWrappers.@declare_library_product(libgdbm, "libgdbm.so.6")
function __init__()
    JLLWrappers.@generate_init_header(Libiconv_jll, Readline_jll)
    JLLWrappers.@init_executable_product(
        gdbm_dump,
        "bin/gdbm_dump",
    )

    JLLWrappers.@init_executable_product(
        gdbm_load,
        "bin/gdbm_load",
    )

    JLLWrappers.@init_executable_product(
        gdbmtool,
        "bin/gdbmtool",
    )

    JLLWrappers.@init_library_product(
        libgdbm,
        "lib/libgdbm.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module Sabermetrics

export basehit, avg, obp, slg, ops, iso, rc, lazyrc, rc27, babip, woba, wraa, wrc, era, whip, fip, k9, bb9, h9, hr9

# --- consts ---

const DIGITS = 3

# --- 共通function ---

function sabrround(num)
    return round(num, RoundNearestTiesAway, digits=DIGITS)
end

# --- バッティングfunction ---

function basehit(hits, double, triple, hr)
    return (hits - double - triple - hr)
end

function avg(hits, ab)
    return sabrround(hits / ab)
end

function obp(hits, bb, hbp, ab, sf)
    return sabrround((hits + bb + hbp) / (ab + bb + hbp + sf))
end

function slg(hits, double, triple, hr, ab)
    return sabrround((basehit(hits, double, triple, hr) + 2 * double + 3 * triple + 4 * hr) / ab)
end

function ops(hits, double, triple, hr, bb, hbp, ab, sf)
    return sabrround(
        obp(hits, bb, hbp, ab, sf) + slg(hits, double, triple, hr, ab)
    )
end

function iso(hits, double, triple, hr, ab)
    return sabrround(
        slg(hits, double, triple, hr, ab)- avg(hits, ab)
    )
end

function rc(hits, double, triple, hr, bb, ab)
    tb = basehit(hits, double, triple, hr) + 2 * double + 3 * triple + 4 * hr
    return sabrround(((hits + bb) * tb) / (ab + bb))
end

function lazyrc(hits, double, triple, hr, bb, hbp, ab, sf)
    tb = basehit(hits, double, triple, hr) + 2 * double + 3 * triple + 4 * hr
    return sabrround(obp(hits, bb, hbp, ab, sf) * tb)
end

function rc27(hits, double, triple, hr, bb, ab, cs, sh, sf, dp)
    return sabrround(rc(hits, double, triple, hr, bb, ab) * 27 / (ab - hits + cs + sh + sf + dp))
end

function babip(hits, hr, ab, k, sf)
    return sabrround((hits - hr) / (ab - k - hr + sf))
end

function woba(hits, double, triple, hr, ab, bb, hbp, sf)
    return sabrround(
        ((bb + hbp) * 0.7 + 0.9 * basehit(hits, double, triple, hr) + 1.3 * double + 1.6 * triple + 2 * hr) / (ab + bb + hbp + sf)
    )
end

function wraa(hits, double, triple, hr, ab, pa, bb, hbp, sf, lwoba)
    return sabrround(((woba(hits, double, triple, hr, ab, bb, hbp, sf) - lwoba) / 1.2) * pa)
end

function wrc(hits, double, triple, hr, ab, pa, bb, hbp, sf, lwoba, lr, lpa)
    return sabrround(
        (((woba(hits, double, triple, hr, ab, bb, hbp, sf) - lwoba) / 1.2) + (lr / lpa)) * pa
    )
end

# --- ピッチングfunction ---

function era(er, ip)
    return sabrround((er / ip) * 9)
end

function whip(hits, bb, ip)
    return sabrround((hits + bb) / ip)
end

function fip(hr, bb, ibb, hbp, k, ip)
    return sabrround(((13 * hr + 3 * (bb - ibb + hbp) - 2 * k) / ip) + 3)
end

function k9(k, ip)
    return sabrround((k / ip) * 9)
end

function bb9(bb, ip)
    return sabrround((bb / ip) * 9)
end

function h9(h, ip)
    return sabrround((h / ip) * 9)
end

function hr9(hr, ip)
    return sabrround((hr / ip) * 9)
end

end


text:
"""
    COVID_Tracking_Clean

This module is meant to clean the data from the COVID Tracking Project.
"""
module COVID_Tracking_Clean
using Diana: Client, GraphQLClient,
             # HTTP
             HTTP.request
using JSON3: JSON3
using LibPQ: Connection, execute, load!,
             # TimeZones
             ZonedDateTime,
             TimeZones.utc_tz,
             # Dates
             DateTime,
             Dates.format,
             Dates.DateFormat,
             unix2datetime,
             now,
             astimezone,
             Interval
using Parameters: @unpack
using TableOperations: select
using CSV: CSV, File,
    # Tables
    Tables.schema,
    Tables.rowtable,
    # DataFrames
    DataFrame,
    categorical!,
    dropmissing!,
    rename!,
    by
using Base64: base64decode, base64encode
import Base: show, summary
# Constants
const COVID_TRACKING_DT = DateFormat("y/m/d HH:MM Z")

for (root, dirs, files) in walkdir(joinpath(@__DIR__))
    for file in files
        isequal("COVID_Tracking_Clean.jl", file) || include(joinpath(root, file))
    end
end

"""
    main()

Performs the action to refresh the data.
Latest data will be written to data/latest.tsv
"""
function main(deploy::Bool = true)
    ENV["PGHOST"] = get(ENV, "PGHOST", "localhost")
    ENV["PGPORT"] = get(ENV, "PGPORT", "5432")
    ENV["PAT"] = get(ENV, "PAT", "")

    opt = Opt("Nosferican",
              ENV["PAT"],
              host = ENV["PGHOST"],
              port = parse(Int, ENV["PGPORT"]))
    execute(opt.conn,
            """
            CREATE EXTENSION IF NOT EXISTS btree_gist;
            DROP TABLE IF EXISTS daily CASCADE;
            CREATE TABLE IF NOT EXISTS daily (
              state char(2) NOT NULL,
              checkts timestamp NOT NULL,
              positive int,
              negative int,
              pending int,
              hospitalized_currently int,
              hospitalized_cumulative int,
              icu_currently int,
              icu_cumulative int,
              ventilation_currently int,
              ventilation_cumulative int,
              recovered int,
              death int,
              EXCLUDE USING gist (state WITH =, checkts WITH =)
            );
            DROP TABLE IF EXISTS qc CASCADE;
            CREATE TABLE IF NOT EXISTS qc (
              state char(2) NOT NULL,
              during timestamp NOT NULL,
              data_quality_grade varchar(2) NOT NULL
            );
            """
            )
    shas = find_shas(opt.pat)
    qc = vcat((get_tbl(opt.pat, sha) for sha in shas)..., cols = :union)
    qc[!,:checkTimeEt] = astimezone.(qc.checkTimeEt, utc_tz)

    # vscodedisplay(data)

    daily = states_daily()

    execute(opt.conn, "BEGIN;")
    load!(daily, opt.conn, "INSERT INTO daily VALUES($(join(("\$$i" for i in 1:13), ',')));")
    execute(opt.conn, "COMMIT;")

    execute(opt.conn, "BEGIN;")
    load!(qc, opt.conn, "INSERT INTO qc VALUES($(join(("\$$i" for i in 1:3), ',')));")
    execute(opt.conn, "COMMIT;")

    output = execute(
        opt.conn,
        """
        DROP SEQUENCE IF EXISTS serial CASCADE;
        CREATE TEMP SEQUENCE serial START WITH 1;
        SELECT
          nextval('serial');
        CREATE MATERIALIZED VIEW IF NOT EXISTS data_quality_grades AS (
            WITH A AS (
              SELECT
                DISTINCT ON (state, during) state,
                during,
                data_quality_grade
              FROM qc
              ORDER BY
                state ASC,
                during ASC
            ),
            B AS (
              SELECT
                state,
                during,
                data_quality_grade,
                LAG(data_quality_grade) OVER (
                  PARTITION BY state
                  ORDER BY
                    during ASC
                ) AS prev
              FROM A
            ),
            C AS (
              SELECT
                state,
                during,
                data_quality_grade,
                CASE
                  WHEN prev IS null
                  OR prev <> data_quality_grade THEN nextval('serial')
                  ELSE currval('serial')
                END AS grp
              FROM B
            ),
            D AS (
              SELECT
                state,
                MIN(during) AS start_ts,
                MAX(during) AS end_ts,
                data_quality_grade
              FROM C
              GROUP BY
                state,
                data_quality_grade,
                grp
            ),
            E AS (
              SELECT
                *,
                LEAD(start_ts) OVER (
                  PARTITION BY state
                  ORDER BY
                    start_ts
                ) AS until_ts
              FROM D
            )
            SELECT
              state,
              CASE
                WHEN until_ts IS NOT NULL THEN tsrange(start_ts, until_ts)
                WHEN start_ts = end_ts THEN tsrange(start_ts, date_trunc('second', NOW() :: timestamp))
                ELSE tsrange(start_ts, end_ts)
              END AS during,
              data_quality_grade
            FROM E
            ORDER BY
              state ASC,
              during ASC
          );
        SELECT
          A.*,
          B.data_quality_grade
        FROM public.daily A
        LEFT JOIN public.data_quality_grades B ON A.state = B.state
          AND B.during @> A.checkts
        ORDER BY
          state ASC,
          checkts ASC;
        """) |>
        rowtable
    content = base64encode(take!(CSV.write(IOBuffer(), output, delim = '\t')));
    if deploy
        restful(opt.pat, "repos/uva-bi-sdad/COVID_Tracking_Clean/contents/data/daily.tsv",
                method = "PUT",
                params = Dict("message" => "Updating data/states.tsv at $(now(utc_tz))",
                              "content" => content,
                              "sha" => find_shas(opt.pat, "MDEwOlJlcG9zaXRvcnkyNTIwNDUwOTQ=", "data", "daily.tsv", true)))
    end
    true
end

export main
end


text:
p = "non_existing_1234_4321"
@test !ispath(p)

@test_throws(AbstractString, lintfile(p))
@test_throws(AbstractString, lintpkg(p))

# Lint package with full path
msgs = lintpkg(joinpath(Pkg.dir("Lint"), "test", "FakePackage"))
@test isempty(msgs)


text:
#=
fft:
- Julia version:
- Author: ymocquar
- Date: 2019-11-15
=#
# function that reverse the order of the pos lowest bits

function _reverse_num(num, pos)
    result = 0
    pos_m1 = pos - 1
    for i = 0:pos_m1
        if (num & (1 << i)) != 0
            result |= 1 << (pos_m1 - i)
        end
    end
    return result
end

"""
    PrepareFftBig( size_fft::Unsigned, [T=BigFloat])

Immutable structure to operate fft transform,
x is the type of non transformed data also called signal.

# Arguments :
- `size_fft::Integer` : Number of values, must be a power of two
- `[T=BigFloat | x::T ]` : type of the values

# Implementation
- size_fft : size of the signal
- tab_permut : permutation
- root_one : size order roots of one
- root_one_conj : conjugate of root_one

# struct parameters
- T : type of Float
- NUMDIMS : Numbers of dimensions
- DIMS : tuple of dimensions

"""

struct PrepareFftBig{T,NUMDIMS,DIMS}
    size_fft::Any
    tab_permut::Any
    root_one::Any
    root_one_conj::Any
    function PrepareFftBig(
        size_fft::NTuple{N,Integer},
        x::T;
        numdims = 2,
        dims::NTuple{N,Integer} = (1,),
    ) where {T<:AbstractFloat,N}
        tab_permut = Vector{Vector{Int64}}(undef, N)
        root_one = Vector{Vector{Complex{T}}}(undef, N)
        root_one_conj = Vector{Vector{Complex{T}}}(undef, N)
        for (inddim, dim) in enumerate(dims)
            sfft = size_fft[inddim]
            @assert prevpow(2, sfft) == sfft "size_fft[$inddim]=$sfft is not a power of 2"
            power = convert(Int64, log2(sfft))
            prec = precision(T)
            setprecision(prec + 32) do
                tab_permut[inddim] = tperm = Vector{Int64}(undef, sfft)
                root_one[inddim] = rone = Vector{Complex{T}}(undef, sfft)
                root_one_conj[inddim] = rone_conj = Vector{Complex{T}}(undef, sfft)
                for i = 1:sfft
                    tperm[i] = _reverse_num(i - 1, power) + 1
                    rone[i] = round(
                        exp(one(T) * 2big(pi) * im * (i - 1) / sfft),
                        digits = prec + 16,
                        base = 2,
                    )
                    rone_conj[i] = conj(rone[i])
                end
            end
        end
        return new{T,numdims,dims}(size_fft, tab_permut, root_one, root_one_conj)
    end
end
function PrepareFftBig(size_fft, type::DataType; kwargs...)
    return PrepareFftBig(size_fft, one(type); kwargs...)
end
function PrepareFftBig(size_fft; kwargs...)
    return PrepareFftBig(size_fft, one(BigFloat); kwargs...)
end
function PrepareFftBig(s::Integer, x::T; kwargs...) where {T<:AbstractFloat}
    return PrepareFftBig((s,), x; kwargs...)
end

function fftbig!(
    par::PrepareFftBig{T,NUMDIMS,DIMS},
    signal;
    flag_inv = false,
) where {T,NUMDIMS,DIMS}
    for (inddim, dim) in enumerate(DIMS)
        s = size(signal, dim)
        @assert par.size_fft[inddim] == s "size of dim $dim is $s while $(par.sizefft[inddim]) is waiting "
        s_div2 = div(s, 2)
        len = s
        n_len = len >> 1
        nb_r = 1
        rootO = flag_inv ? par.root_one[inddim] : par.root_one_conj[inddim]
        perm = par.tab_permut[inddim]
        #    prec= precision(real(rootO[1]))
        #    setprecision(prec+32) do
        while n_len != 0
            start = 1
            suite = start + n_len
            for i = 1:nb_r
                deb = 1
                for j = deb:nb_r:s_div2
                    if NUMDIMS == 1
                        signal[start], signal[suite] = (signal[start] + signal[suite]),
                        (signal[start] - signal[suite]) * rootO[j]
                    else
                        s_start = selectdim(signal, dim, start)
                        s_suite = selectdim(signal, dim, suite)
                        s_sum = s_start + s_suite
                        s_diff = (s_start - s_suite) * rootO[j]
                        s_start .= s_sum
                        s_suite .= s_diff
                    end
                    # if NDIMS == 1
                    #     signal[start], signal[suite] = (signal[start] + signal[suite]),
                    #     (signal[start] - signal[suite])*rootO[j]
                    # elseif NUMDIM == 1
                    #     signal[start, :], signal[suite, :] = (signal[start, :] + signal[suite, :]),
                    #     (signal[start, :] - signal[suite, :])*rootO[j]
                    # else
                    #     signal[:,start], signal[:,suite] = (signal[:,start] + signal[:,suite]),
                    #     (signal[:,start] - signal[:,suite])*rootO[j]
                    # end
                    start += 1
                    suite += 1
                end
                start = suite
                suite = start + n_len
            end
            len = n_len
            n_len >>= 1
            nb_r <<= 1
        end
        signal .= if NUMDIMS == 1
            signal[perm]
        else
            selectdim(signal, dim, perm)
        end

        if flag_inv
            signal ./= s
        end
    end
    return signal
end
function fftbig(
    par::PrepareFftBig{T,NUMDIMS,DIMS},
    signal;
    flag_inv = false,
) where {T,NUMDIMS,DIMS}
    fl = flag_inv
    return fftbig!(
        par::PrepareFftBig,
        copy(convert(Array{Complex{T}}, signal)),
        flag_inv = fl,
    )
end

# fftgen(_::Any, t::AbstractArray{Complex{Float64}}) = fft(t, (1,))
# fftgen!(_::Any, t::AbstractArray{Complex{Float64}}) = fft!(t, (1,))
# fftgen(_::Any, t::AbstractArray{Float64}) = fft(t, (1,))
# fftgenall(_::Any, t::AbstractArray{Complex{Float64}}) = fft(t, ntuple(x->x,ndims(t)))
# fftgenall!(_::Any, t::AbstractArray{Complex{Float64}}) = fft(t, ntuple(x->x,ndims(t)))
fftgenall(_::Any, t::AbstractArray{Float64}) = fft(t, ntuple(x -> x, ndims(t)))

# function fftgen(
#     _::PrepareFftBig{T, NUMDIMS, DIMS},
#     t::AbstractArray{Complex{Float64}}
# ) where {T, NUMDIMS, DIMS}
#     return fft(t, DIMS)
# end
function fftgen!(
    _::PrepareFftBig{T,NUMDIMS,DIMS},
    t::AbstractArray{Complex{Float64}},
) where {T,NUMDIMS,DIMS}
    return fft!(t, DIMS)
end
function fftgen(
    _::PrepareFftBig{T,NUMDIMS,DIMS},
    t::AbstractArray{Float64},
) where {T,NUMDIMS,DIMS}
    return fft(t, DIMS)
end
fftgenall(p::PrepareFftBig, t::AbstractArray{Float64}) = fftgen(p, t)
# fftgenall(p::PrepareFftBig,t::AbstractArray{Complex{Float64}})=fftgen(p, t)
# fftgenall!(p::PrepareFftBig,t::AbstractArray{Complex{Float64}})=fftgen!(p, t)

fftgen(p::PrepareFftBig, t::AbstractArray{T}) where {T<:AbstractFloat} = fftbig(p, t)
fftgen(p::PrepareFftBig, t::AbstractArray{Complex{T}}) where {T<:AbstractFloat} =
    fftbig(p, t)
fftgen!(p::PrepareFftBig, t::AbstractArray{Complex{T}}) where {T<:AbstractFloat} =
    fftbig!(p, t)
fftgenall(p::PrepareFftBig, t::AbstractArray{T}) where {T} = fftgen(p, t)
fftgenall(p::PrepareFftBig, t::AbstractArray{Complex{T}}) where {T} = fftgen(p, t)
fftgenall!(p::PrepareFftBig, t::AbstractArray{Complex{T}}) where {T} = fftgen!(p, t)

# ifftgen(_::Any, t::AbstractArray{Complex{Float64}}) = ifft(t, (1,))
# ifftgen!(_::Any, t::AbstractArray{Complex{Float64}}) = ifft!(t, (1,))
# ifftgen(_::Any, t::AbstractArray{Float64}) = ifft(t, (1,))
ifftgenall(_::Any, t::AbstractArray{Complex{Float64}}) = ifft(t, ntuple(x -> x, ndims(t)))
# ifftgenall!(_::Any, t::AbstractArray{Complex{Float64}}) = ifft(t, ntuple(x->x,ndims(t)))
# ifftgenall(_::Any, t::AbstractArray{Float64}) = ifft(t, ntuple(x->x,ndims(t)))

function ifftgen(
    _::PrepareFftBig{T,NUMDIMS,DIMS},
    t::AbstractArray{Complex{Float64}},
) where {T,NUMDIMS,DIMS}
    return ifft(t, DIMS)
end
function ifftgen!(
    _::PrepareFftBig{T,NUMDIMS,DIMS},
    t::AbstractArray{Complex{Float64}},
) where {T,NUMDIMS,DIMS}
    return ifft!(t, DIMS)
end
# function ifftgen(
#     _::PrepareFftBig{T, NUMDIMS, DIMS},
#     t::AbstractArray{Float64}
# ) where {T, NUMDIMS, DIMS}
#     return ifft(t, DIMS)
# end
# ifftgenall(p::PrepareFftBig,t::AbstractArray{Float64})=ifftgen(p, t)
ifftgenall(p::PrepareFftBig, t::AbstractArray{Complex{Float64}}) = ifftgen(p, t)
# ifftgenall!(p::PrepareFftBig,t::AbstractArray{Complex{Float64}})=ifftgen!(p, t)

# ifftgen(p::PrepareFftBig, t::AbstractArray{T}) where {T}  = fftbig(p, t, flag_inv = true)
ifftgen(p::PrepareFftBig, t::AbstractArray{Complex{T}}) where {T} =
    fftbig(p, t, flag_inv = true)
ifftgen!(p::PrepareFftBig, t::AbstractArray{Complex{T}}) where {T} =
    fftbig!(p, t, flag_inv = true)
# ifftgenall(p::PrepareFftBig,t::AbstractArray{T}) where {T} =ifftgen(p, t)
ifftgenall(p::PrepareFftBig, t::AbstractArray{Complex{T}}) where {T} = ifftgen(p, t)
ifftgenall!(p::PrepareFftBig, t::AbstractArray{Complex{T}}) where {T} = ifftgen!(p, t)


text:
module LinearAlgebra
  files = ["Commons","HouseHolders"]
  for f in files
    include("$f.jl")
  end
end

max_stars_repo_name: ajermaky/CustomMLLibrary

max_issues_repo_name: ajermaky/CustomMLLibrary

max_forks_repo_name: ajermaky/CustomMLLibrary


text:
using libm
srand(1)

X = exp(randn(10_000_000))

macro testsum(name,fn,X)
    quote
        function test(X)
            s = 0.0
            for x in X
                s += $fn(x)
            end
            s
        end

        test($X)
        test($X)

        println($name)
        @profile test($X)
        Profile.print()
    end
end

# syslog(x::Float64) = Core.Intrinsics.nan_dom_err(ccall((:log,:libm),Float64,(Float64,),x), x)

# @testsum "Openlibm:" log X
# @testsum "System:" syslog X
# @testsum "Julia fdlibm:" libm.log X
@testsum "Julia tang:" libm.log_tang X
# @testsum "fma:" libm.log_fma X
# @testsum "muladd:" libm.log_muladd X


text:
using DeepViz
using Test

@testset "DeepViz.jl" begin
    # Write your own tests here.
end


text:
# the channel example using an event-based approach
#

using DiscreteEvents, Printf, Random

mutable struct Server
  id::Int64
  name::AbstractString
  input::Channel
  output::Channel
  op     # operation to take
  token  # current token

  Server(id, name, input, output, op) = new(id, name, input, output, op, nothing)
end

function take(S::Server)
    if isready(S.input)
        S.token = take!(S.input)
        @printf("%5.2f: %s %d took token %d\n", tau(), S.name, S.id, S.token)
        event!(SF(put, S), after, rand())         # call put after some time
    else
        event!(SF(take, S), SF(isready, S.input)) # call again if input is ready
    end
end

function put(S::Server)
    put!(S.output, S.op(S.id, S.token))
    S.token = nothing
    take(S)
end

reset!(𝐶)
Random.seed!(123)

ch1 = Channel(32)  # create two channels
ch2 = Channel(32)

si = shuffle(1:8)
for i in 1:2:8
    take(Server(si[i], "foo", ch1, ch2, +))
    take(Server(si[i+1], "bar", ch2, ch1, *))
end

put!(ch1, 1) # put first token into channel 1

run!(𝐶, 10)


text:
module RAM

const DIR = @__DIR__
const SAVEDIR = abspath(joinpath(DIR, "..", "checkpoints"))

include("model.jl")
include("train.jl")
include("visualization.jl")
include("exploration.jl")

end # module

max_stars_repo_name: ilkerkesen/RAM

max_issues_repo_name: ilkerkesen/RAM

max_forks_repo_name: ilkerkesen/RAM


text:
@testset "diagnostics" begin
    @testset "bfmi" begin
        idata = load_arviz_data("centered_eight")
        @test bfmi(idata) == ArviZ.arviz.bfmi(idata)
        rng = Random.MersenneTwister(42)
        arr = randn(rng, 4, 100)
        @test bfmi(arr) == ArviZ.arviz.bfmi(arr)
    end

    @testset "ess" begin
        rng = Random.MersenneTwister(42)
        arr = randn(rng, 4, 100)
        @test ess(arr) == ArviZ.arviz.ess(arr)
        @test ess((x = arr,)) isa ArviZ.Dataset
        @test ess((x = arr,)).x.values == ArviZ.arviz.ess(Dict(:x => arr)).x.values
    end

    @testset "mcse" begin
        rng = Random.MersenneTwister(42)
        arr = randn(rng, 4, 100)
        @test mcse(arr) == ArviZ.arviz.mcse(arr)
        @test mcse((x = arr,)) isa ArviZ.Dataset
        @test mcse((x = arr,)).x.values == ArviZ.arviz.mcse(Dict(:x => arr)).x.values
    end

    @testset "rhat" begin
        rng = Random.MersenneTwister(42)
        arr = randn(rng, 4, 100)
        @test rhat(arr) == ArviZ.arviz.rhat(arr)
        @test rhat((x = arr,)) isa ArviZ.Dataset
        @test rhat((x = arr,)).x.values == ArviZ.arviz.rhat(Dict(:x => arr)).x.values
    end

    @testset "geweke" begin
        rng = Random.MersenneTwister(42)
        arr = randn(rng, 1000)
        @test geweke(arr, 0.1, 0.5, 10) == ArviZ.arviz.geweke(arr, 0.1, 0.5, 10)
    end
end


text:
using Random, Plots, Distributions, Statistics

# Example solution for exercise 8.8 using Julia 1.4.1
#
# My solution has similar shape with the book, but different start state value
# under the greedy policy. I am not sure where goes wrong, probably in the
# reward calculation? But my results are similar to all the other people's
# results which I found online (see below reference implementation). So just
# take my solution as one of the references, don't treat it absolutely correct.
#
# By @burmecia at GitHub, 15 May 2020
#
# Other reference implementation:
#
# 1. https://github.com/ShangtongZhang/reinforcement-learning-an-introduction/blob/master/chapter08/trajectory_sampling.py
# 2. https://github.com/JuliaReinforcementLearning/ReinforcementLearningAnIntroduction.jl/blob/b5c718a891a4b3db4fae177b8b33ca506df1ecea/notebooks/Chapter08_Trajectory_Sampling.ipynb
# 3. https://github.com/enakai00/rl_book_solutions/blob/master/Chapter08/Exercise_8_8_Solution.ipynb

stat_len = 10000
T = 200000

ε = 0.1
act_len = 2
term_prob = 0.1
tick_num = 100

# make transition matrix, (state, action) -> (next_state, reward)
# element is b tuples: [(next_state, reward),...]
function make_transition_matrix(b)
    trans = [
        map(s -> (s, randn()), rand(1:stat_len, b))
        for s = 1:stat_len, a = 1:act_len
    ]
    # terminal state, transit to itself with zero reward
    trans[stat_len, 1] = [(stat_len, 0) for _ = 1:b]
    trans[stat_len, 2] = [(stat_len, 0) for _ = 1:b]
    trans
end

# take action, observe the next state and reward
function next_state(s, a, trans)
    if rand() < term_prob
        return stat_len, 0
    end
    sample(trans[s, a])
end

# evaluate start state value under greedy policy
function start_state_value(Q, trans)
    n = 200
    returns = zeros(Float64, n)

    for i = 1:n
        s = 1
        while s != stat_len
            a = argmax(Q[s, :])
            s, r = next_state(s, a, trans)
            returns[i] += r
        end
    end

    mean(returns)
end

function run_uniform(b, trans)
    Q = zeros(Float64, stat_len, act_len)
    interval = T ÷ tick_num
    vπ_s0 = zeros(Float64, tick_num + 1)

    for t = 0:T
        s = t % stat_len + 1
        a = (t % (stat_len * act_len)) ÷ stat_len + 1
        next_states = trans[s, a]

        Q[s, a] =
            (1 - term_prob) *
            mean(map(ns -> (ns[2] + maximum(Q[ns[1], :])), next_states))

        if t % interval == 0
            vπ_s0[t ÷ interval + 1] = start_state_value(Q, trans)
        end
    end

    vπ_s0
end

function run_on_policy(b, trans)
    Q = zeros(Float64, stat_len, act_len)
    interval = T ÷ tick_num
    vπ_s0 = zeros(Float64, tick_num + 1)
    s = 1

    for t = 0:T
        if rand() < ε
            a = rand(1:act_len)
        else
            a = argmax(Q[s, :])
        end
        next_states = trans[s, a]

        Q[s, a] =
            (1 - term_prob) *
            mean(map(ns -> (ns[2] + maximum(Q[ns[1], :])), next_states))

        s, r = next_state(s, a, trans)
        if s == stat_len
            s = 1
        end

        if t % interval == 0
            vπ_s0[t ÷ interval + 1] = start_state_value(Q, trans)
        end
    end

    vπ_s0
end

function experiment(b_list, p, subplot)
    task_cnt = 100
    interval = T ÷ tick_num
    x = 0:interval:T

    for b in b_list
        uniform = zeros(Float64, task_cnt, length(x))
        on_policy = zeros(Float64, task_cnt, length(x))

        println("start subplot $(subplot), b=$(b)")
        for i = 1:task_cnt
            if i % 10 == 0
                println("  task $(i)")
            end
            trans = make_transition_matrix(b)
            uniform[i, :] = run_uniform(b, trans)
            on_policy[i, :] = run_on_policy(b, trans)
        end

        uniform = mean(uniform, dims = 1)
        on_policy = mean(on_policy, dims = 1)

        plot!(
            p,
            x,
            uniform[1, :],
            subplot = subplot,
            title = "$(stat_len) states",
            label = "b=$(b), uniform",
        )
        plot!(
            p,
            x,
            on_policy[1, :],
            subplot = subplot,
            label = "b=$(b), on-policy",
        )
    end
end

function main()
    global stat_len, T

    p = plot(
        xlabel = "Computation time, in expected updates",
        ylabel = "Start state value",
        legend = :best,
        size = (600, 800),
        layout = (2, 1),
    )

    stat_len, T = 1000, 20000
    experiment([1, 3, 10], p, 1)
    stat_len, T = 10000, 200000
    experiment([1, 3], p, 2)

    display(p)
end

# uncomment the main function below to start learning
#main()

max_stars_repo_path: cheatsheet/Solutions/Chapter 8/ex8_8.jl

max_issues_repo_path: cheatsheet/Solutions/Chapter 8/ex8_8.jl

max_forks_repo_path: cheatsheet/Solutions/Chapter 8/ex8_8.jl


text:
"""
    dominates(p1::Array, p2::Array, fns)

    Return true if each element in p1 is not less than the corresponding element in p2 and at least one element in p1 is bigger than the corresponding element in p2.

# Arguments
- `p1::Array`: Numeric array of n elements.
- `p2::Array`: Numeric array of n elements.
- `fns::Array{Function, 1}`: Vector of functions with elements of maximum or minimum.

# Examples
```julia-repl
julia> dominates([1,2,3], [1,2,1], makeminmax([maximum, maximum, maximum]))
true
julia> dominates([0,0,0,0], [1,0,0,0], makeminmax([maximum, maximum, maximum]))
false
```

# References
Deb, Kalyanmoy, et al. "A fast elitist non-dominated sorting genetic algorithm for multi-objective optimization: NSGA-II."
International conference on parallel problem solving from nature. Springer, Berlin, Heidelberg, 2000.
"""
function dominates(p1::Array, p2::Array, fns::Array{Function, 1})::Bool
    n = length(p1)
    notworse = count(i -> if fns[i] == maximum p1[i] < p2[i] else p1[i] > p2[i] end, 1:n)
    better   = count(i -> if fns[i] == maximum p1[i] > p2[i] else p1[i] < p2[i] end , 1:n)
    return (notworse == 0) && (better > 0)
end

"""
    ndsranks(data, fns)

    Sort multidimensional data usin non-dominated sorting algorithm.

# Arguments
- `data::DataFrame`: DataFrame of variables.
- `fns::Array{Function, 1}`: Vector of functions with elements of maximum or minimum.

# References
Deb, Kalyanmoy, et al. "A fast elitist non-dominated sorting genetic algorithm for multi-objective optimization: NSGA-II."
International conference on parallel problem solving from nature. Springer, Berlin, Heidelberg, 2000.
"""
function ndsranks(data::DataFrame, fns::Array{Function, 1})::Array{Int}

    #mat = convert(Matrix, data)
    mat = Matrix(data)
    return ndsranks(mat, fns)

end

"""
    ndsranks(data, fns)

    Sort multidimensional data using non-dominated sorting algorithm.

# Arguments
- `data::Matrix`: n x k matrix of observations where n is number of observations and k is number of variables.
- `fns::Array{Function, 1}`: Vector of functions with elements of maximum or minimum.

# Examples
```julia-repl

```

# References
Deb, Kalyanmoy, et al. "A fast elitist non-dominated sorting genetic algorithm for multi-objective optimization: NSGA-II."
International conference on parallel problem solving from nature. Springer, Berlin, Heidelberg, 2000.
"""
function ndsranks(data::Matrix, fns::Array{Function, 1})::Array{Int64}

    n, _ = size(data)

    ranks = zeros(Int64, n)

    mat = convert(Matrix, data)

    @inbounds for i in 1:n
        @inbounds for j in 1:n
            if i != j
                if dominates(mat[i,:], mat[j,:], fns)
                    ranks[i] += 1
                end
            end
        end
    end

    return ranks
end

"""
    nds(data, fns)

    Sort multidimensional data using non-dominated sorting algorithm.

# Arguments
- `data::DataFrame`: n x k desicion matrix with n cases and k criteria.
- `fns::Array{Function, 1}`: Vector of functions with elements of maximum or minimum.

# Output
- `::NDSResult`: NDSResult object that holds multiple outputs including ranks and best index.

# Examples
```julia-repl
 julia> cases = [
               1.0 2.0 3.0;
               2.0 1.0 3.0;
               1.0 3.0 2.0;
               4.0 5.0 6.0
           ];

julia> nd = makeDecisionMatrix(cases)
4×3 DataFrame
 Row │ Crt1     Crt2     Crt3
     │ Float64  Float64  Float64
─────┼───────────────────────────
   1 │     1.0      2.0      3.0
   2 │     2.0      1.0      3.0
   3 │     1.0      3.0      2.0
   4 │     4.0      5.0      6.0

julia> result = nds(nd, makeminmax([maximum, maximum, maximum]));

julia> result.ranks
4-element Array{Int64,1}:
 0
 0
 0
 3

julia> result.bestIndex
4
```

# References
Deb, Kalyanmoy, et al. "A fast elitist non-dominated sorting genetic algorithm for multi-objective optimization: NSGA-II."
International conference on parallel problem solving from nature. Springer, Berlin, Heidelberg, 2000.
"""
function nds(data::DataFrame, fns::Array{Function, 1})::NDSResult

    ranks = ndsranks(data, fns)

    bestIndex = sortperm(ranks) |> last

    result = NDSResult(
        ranks,
        bestIndex
    )

    return result
end

max_stars_repo_name: jbytecode/jMCDM

max_issues_repo_name: jbytecode/jMCDM

max_forks_repo_name: jbytecode/jMCDM