text:
using LinearAlgebra: Adjoint, Transpose
using Base: StridedVector, StridedMatrix, StridedVecOrMat
AdjOrTransSparseMatrix1DVBC{W, Tv, Ti} = Union{SparseMatrix1DVBC{W, Tv, Ti}, Adjoint{<:Any,<:SparseMatrix1DVBC{W, Tv, Ti}}, Transpose{<:Any, <:SparseMatrix1DVBC{W, Tv, Ti}}}
const DenseMatrixUnion = Union{StridedMatrix, LowerTriangular, UnitLowerTriangular, UpperTriangular, UnitUpperTriangular, BitMatrix}
const AdjOrTransDenseMatrix = Union{DenseMatrixUnion,Adjoint{<:Any,<:DenseMatrixUnion},Transpose{<:Any,<:DenseMatrixUnion}}
const DenseInputVector = Union{StridedVector, BitVector}
const DenseInputVecOrMat = Union{AdjOrTransDenseMatrix, DenseInputVector}

function LinearAlgebra.mul!(y::StridedVector, A::SparseMatrix1DVBC{W, Tv, Ti}, x::StridedVector, α::Number, β::Number) where {W, Tv<:SIMD.VecTypes, Ti}
    Δw = fld(DEFAULT_SIMD_SIZE, sizeof(eltype(y)))
    return _mul!(y, A, x, α, β, Val(Δw))
end
@generated function _mul!(y::StridedVector, A::SparseMatrix1DVBC{W, Tv, Ti}, x::StridedVector, α::Number, β::Number, ::Val{Δw}) where {W, Tv<:SIMD.VecTypes, Ti, Δw}
    if Δw == 1
        ws = (1:W...,)
    else
        ws = (1, Δw : Δw : (W + Δw - 1)...,)
    end

    function stripe_body(safe, w)
        if false && safe #TODO how do I zero out the other entries efficiently?
            thk = quote
                tmp = convert(Vec{$w, eltype(y)}, vload(Vec{$w, eltype(x)}, x, j))
            end
        else
            thk = quote
                tmp = convert(Vec{$w, eltype(y)}, Vec{$w, eltype(x)}(($(map(Δj -> :($Δj < w ? x[j + $Δj] : zero(eltype(x))), 0:w-1)...),)))
            end
        end
        thk = quote
            $thk
            q = A_ofs[l]
            for Q = A_pos[l]:(A_pos[l + 1] - 1)
                y[A_idx[Q]] += sum(convert(Vec{$w, eltype(y)}, vload(Vec{$w, Tv}, A_val, q)) * tmp)
                q += w
            end
        end

        return thk
    end

    thunk = quote
        @inbounds begin
            size(A, 1) == size(y, 1) || throw(DimensionMismatch())
            size(A, 2) == size(x, 1) || throw(DimensionMismatch())
            (m, n) = size(A)

            yα = convert(eltype(y), α)

            if β != 1
                β != 0 ? rmul!(y, β) : fill!(y, zero(eltype(y)))
            end

            Φ_spl = A.Φ.spl
            A_pos = A.pos
            A_idx = A.idx
            A_ofs = A.ofs
            A_val = A.val
            L = length(A.Φ)
            L_safe = L
            while L_safe >= 1 && n + 1 - Φ_spl[L_safe] < $(Δw * cld(W, Δw)) L_safe -= 1 end
            for l = 1:L_safe
                j = Φ_spl[l]
                w = Φ_spl[l + 1] - j
                $(le_nest(w->stripe_body(true, w), :w, ws))
            end
            for l = L_safe + 1:L
                j = Φ_spl[l]
                w = Φ_spl[l + 1] - j
                $(le_nest(w->stripe_body(false, w), :w, ws))
            end
            return y
        end
    end

    if eltype(y) <: SIMD.FloatingTypes
        thunk = quote
            @fastmath $thunk
        end
    end

    return thunk
end

function LinearAlgebra.mul!(y::StridedVector, adjA::Union{Adjoint{<:Any,<:SparseMatrix1DVBC{W, Tv, Ti}}, Transpose{<:Any, <:SparseMatrix1DVBC{W, Tv, Ti}}}, x::StridedVector, α::Number, β::Number) where {W, Tv<:SIMD.VecTypes, Ti}
    Δw = fld(DEFAULT_SIMD_SIZE, sizeof(eltype(y)))
    _mul!(y, adjA, x, α, β, Val(Δw))
end

@generated function _1DVBR_mul!(y::StridedVector, x::StridedVector, α::Number, β::Number, l, j, w, A_ofs, A_pos, A_idx, A_val, ::Val{Δw}, ::Val{W}, ::Val{Tv}, ::Val{Ti}) where {Δw, W, Tv, Ti}
    if Δw == 1
        ws = (1:W...,)
    else
        ws = (1, Δw : Δw : (W + Δw - 1)...,)
    end

    function stripe_body(safe, w)
        thk = quote
            tmp = Vec{$w, eltype(y)}(zero(eltype(y)))
            q = A_ofs[l]
            for Q = A_pos[l]:(A_pos[l + 1] - 1)
                tmp += convert(Vec{$w, eltype(y)}, vload(Vec{$w, Tv}, A_val, q)) * convert(eltype(y), x[A_idx[Q]])
                q += w
            end
        end
        if safe
            thk = quote
                $thk
                vstore(tmp, y, j)
            end
        else
            thk = quote
                $thk
                for Δj = 1:w
                    y[j + Δj - 1] = tmp[Δj]
                end
            end
        end

        return thk
    end

    thunk = quote
        @inbounds begin
            $(le_nest(w->stripe_body(false, w), :w, ws))
        end
    end

    if eltype(y) <: SIMD.FloatingTypes
        thunk = :(@fastmath $thunk)
    end

    return :(Base.@_inline_meta; $thunk)
end

@inline function _mul!(y::StridedVector, adjA::Union{Adjoint{<:Any,<:SparseMatrix1DVBC{W, Tv, Ti}}, Transpose{<:Any, <:SparseMatrix1DVBC{W, Tv, Ti}}}, x::StridedVector, α::Number, β::Number, Δw) where {W, Tv<:SIMD.VecTypes, Ti}
    @inbounds begin
        A = adjA.parent
        size(A, 2) == size(y, 1) || throw(DimensionMismatch())
        size(A, 1) == size(x, 1) || throw(DimensionMismatch())
        (m, n) = size(A)

        yα = convert(eltype(y), α)

        if β != 1
            β != 0 ? rmul!(y, β) : fill!(y, zero(eltype(y)))
        end

        Φ_spl = A.Φ.spl
        A_pos = A.pos
        A_idx = A.idx
        A_ofs = A.ofs
        A_val = A.val
        L = length(A.Φ)
        #=
        L_safe = L
        while L_safe >= 1 && n + 1 - Φ_spl[L_safe] < $(Δw * cld(W, Δw)) L_safe -= 1 end
        for l = 1:L_safe
            j = Φ_spl[l]
            w = Φ_spl[l + 1] - j
            $(le_nest(w->stripe_body(true, w), :w, ws))
        end
        for l = L_safe + 1:L
            j = Φ_spl[l]
            w = Φ_spl[l + 1] - j
            $(le_nest(w->stripe_body(false, w), :w, ws))
        end
        =#
        l′ = Atomic{Int}(1)
        @threads for t = 1:nthreads()
            while (l = atomic_add!(l′, 1)) <= L
                j = Φ_spl[l]
                w = Φ_spl[l + 1] - j

                _1DVBR_mul!(y, x, α, β, l, j, w, A_ofs, A_pos, A_idx, A_val, Δw, Val(W), Val(Tv), Val(Ti))
            end
        end
        return y
    end
end

Base.:*(adjA::AdjOrTransSparseMatrix1DVBC, x::StridedVector{Tx}) where {Tx} =
    (T = Base.promote_op(LinearAlgebra.matprod, eltype(adjA), Tx); mul!(similar(x, T, size(adjA, 1)), adjA, x, true, false))
Base.:*(adjA::AdjOrTransSparseMatrix1DVBC, B::AdjOrTransDenseMatrix) =
    (T = Base.promote_op(LinearAlgebra.matprod, eltype(adjA), eltype(B)); mul!(similar(B, T, (size(adjA, 1), size(B, 2))), adjA, B, true, false))


text:
using NLsolve

# set moving cost and time horizon vectors (first 12 elements are for type 1; last 12 elements are for type 2)
MC = -[-2.654795205, -3.144665335, -6.574142317, -7.064012447, -8.476504970, -8.966375099, -7.628074092, -8.117944222, -8.527064851, -9.016934981, -9.141404761, -9.631274891, -2.910429570, -3.400299700, -6.829776683, -7.319646812, -8.732139335, -9.222009465, -7.883708458, -8.373578588, -8.782699217, -9.272569347, -9.397039127, -9.886909257]
horzns = 65 .- [18, 18, 39, 39, 39, 39, 25, 25, 40, 40, 55, 55, 18, 18, 39, 39, 39, 39, 25, 25, 40, 40, 55, 55]
incs = [34530.18, 34530.18, 41364.76, 41364.76, 41364.76, 41364.76, 27142.27, 27142.27, 44071.81, 44071.81, 42667.49, 42667.49, 34530.18, 34530.18, 41364.76, 41364.76, 41364.76, 41364.76, 27142.27, 27142.27, 44071.81, 44071.81, 42667.49, 42667.49]

function executor(MCvec,horzns)
    # Present value formula
    function f!(F,x,a)
        # a = [beta, tau, T, MC]
        F[1] = sum( ( a[1]^(t-a[2]) )*log(x[1]) for t in a[2]:a[3] )-a[4]
    end

    beta=0.9
    tau=1
    for i=1:length(MCvec)
        g!(F,x) = f!(F,x,[beta, tau, horzns[i], MCvec[i]]) # create closure
        res = nlsolve(g!, [1.5])
        output = round.(res.zero[1]; digits=4)
        println("MC = ",MC[i],"; T = ",horzns[i],"; tau = ",tau,";wage ratio = ",output,";wage pct = ",round(100*(output-1);digits=1))
    end
end

executor(MC,horzns)

function executorPV(MCvec,horzns,incs)
    # Present value formula
    function f!(F,x,a)
        # a = [beta, tau, T, MC]
        F[1] = sum( ( a[1]^(t-a[2]) )*log(x[1]) for t in a[2]:a[3] )-a[4]
    end

    beta=0.9
    tau=1
    for i=1:length(MCvec)
        g!(F,x) = f!(F,x,[beta, tau, horzns[i], MCvec[i]]) # create closure
        res = nlsolve(g!, [1.5])
        output = round.(res.zero[1]; digits=4)
        println("MC = ",MC[i],"; T = ",horzns[i],"; tau = ",tau,";wage ratio = ",output,";PV = ",round( sum((beta^(t-1)) * ((res.zero[1]-1)*incs[i]) for t=1:horzns[i]); digits=0) )
    end
end

executorPV(MC,horzns,incs)

println("std dev of local amenities: ",sum((0.9^(t-1)) * (.077*incs[1]) for t=1:20) )
println("range of local amenities: ",sum((0.9^(t-1)) * (.302*incs[1]) for t=1:20) )
println("birth state bonus: ",sum((0.9^(t-1)) * (.189*incs[1]) for t=1:20) )

ext: jl

max_stars_repo_name: tyleransom/frictions-JHR

max_issues_repo_name: tyleransom/frictions-JHR

max_forks_repo_name: tyleransom/frictions-JHR


text:
using Documenter
using Markdown

# Types and functions to generate a Markdown table with links to package badges etc.
struct PackageDefinition
    name :: String
    url :: String
    docs :: Vector{Pair{String, String}} # type => URL
    buildbadges :: Vector{Pair{String, String}} # badge => URL
end

function markdown(p::PackageDefinition)
    row = Any[]
    push!(row, Markdown.Link(p.name, p.url))
    push!(row, [
        Markdown.Link(
            [Markdown.Image("https://img.shields.io/badge/docs-$(ver)-blue.svg", "$(ver)")],
            url
        )
        for (ver, url) in p.docs
    ])
    push!(row, [
        Markdown.Link(
            [Markdown.Image(image, "")],
            url
        )
        for (image, url) in p.buildbadges
    ])
end

function package_table_markdown(packages)
    titles = map(["Package", "Documentation", "Build Status"]) do s
        Markdown.Bold(s)
    end
    table = Markdown.Table([titles], [:l, :c, :c])
    for p in packages
        push!(table.rows, markdown(p))
    end
    Markdown.MD(table)
end

# Build the docs
makedocs(
    sitename = "JuliaDocs",
    pages = [
        "Home" => "index.md",
    ],
    format = Documenter.HTML(
        edit_link = "source",
        assets = ["assets/favicon.ico"],
        analytics = "UA-136089579-1",
    ),
)

deploydocs(
    repo = "github.com/JuliaDocs/juliadocs.github.io.git",
    branch = "master",
    devbranch = "source",
    versions = nothing,
    push_preview = true,
)


text:
import Base: *, sort
using Random, Statistics, DataFrames, CSV
using BioSequences, Combinatorics

# Overload the concatenation operator to combine a sequence and
# a single base, i.e. dna"AGCGTGC" * DNA_T
function *(seq::LongDNASeq, base::DNA)
    seq * LongDNASeq([base])
end

# Overload sort to order DNA codes by degeneracy:
#   A,G,C,T,M,R,W,S,Y,K,V,H,D,B,N
function sort(seq::LongSequence{DNAAlphabet{4}}; rev::Bool=false)
    sort(convert(Vector{DNA}, seq), by=degeneracy, rev=rev)
end

function is_palindrome(randomer)
    return occursin(randomer,reverse_complement(randomer))
end

"""
    isvalid(seq, sites)

Check if any of the sites appear in the sequence.
"""
function isvalid(seq, sites)
    for site in sites
        if occursin(site, seq)
            return false
        end
    end
    return true
end

"""
    find_valid_randomer(prefix, bases, sites)

Randomly select bases to add the prefix, ensuring no site
appears in the randomer.
"""
function find_valid_randomer(prefix, bases, sites)
    n = length(bases)
    pre_n = length(prefix)
    randomer = prefix * (dna"-" ^ n)
    for i = 1:n
        candidates = shuffle(bases[i])
        for cand in candidates
            randomer[pre_n+i] = cand
            if isvalid(randomer[1:pre_n+i], sites)
                break
            else
                randomer[pre_n+i] = DNA_Gap
            end
        end
        if randomer[pre_n+i] == DNA_Gap
            break
        end
    end
    return randomer
end

"""
    degeneracy(base::DNA[, uselog2=true])

Calculate the degeneracy of a single base. The degeneracy
is the number of bases in the code, so
    degeneracy(DNA_N) == 4
and
    degeneracy(DNA_A) == 1

If `uselog2`, the degeneracy is log2 transformed. By convention,
the degneracy of `DNA_Gap` is zero in both the log2-transformed
and the untransformed case.
"""
function degeneracy(base::DNA; uselog2=true)
    deg = 0.0
    if base == DNA_N
        deg = 4.0
    elseif base == DNA_A || base == DNA_C || base == DNA_G || base == DNA_T
        deg = 1.0
    elseif base == DNA_B || base == DNA_D || base == DNA_H || base == DNA_V
        deg = 3.0
    elseif (base == DNA_K || base == DNA_Y || base == DNA_S || base == DNA_W
            || base == DNA_R || base == DNA_M)
        deg = 2.0
    else # DNA_Gap
        deg = 0.0
    end

    if uselog2
        if base == DNA_Gap
            deg = 2^-100
        end
        return log2(deg)
    else
        return deg
    end
end

"""
    degeneracy(seq::LongSequence{DNAAlphabet{4}}; uselog2=true)

Calculate the degeneracy of a DNA sequence.
"""
function degeneracy(seq::LongSequence{DNAAlphabet{4}}; uselog2=true)
    if uselog2
        return sum(degeneracy.(seq, uselog2=true))
    else
        return prod(degeneracy.(seq ;uselog2=false))
    end
end

## Read in the list of sequences you want to exclude in the designed primers.
# Function finds every unique hexamer in the list of sequences and stores it as an array
#. For example, use list of rRNA and tRNA sequences. This funciton takes care of orientation. We end up using the reverse complement since we are finding primers that won't hybridize
## Example
#blocking_sites=read_blocking_sites(CSV.read("./SMU_UA159_rRNA_tRNA.csv"))
function read_invalid_sites(data; randomerlen=6)
    blocking_sites=[];
    for i =1:nrow(data)
        seq=reverse_complement(LongDNASeq(string(data[i,:Sequence])));
        for j=1:length(seq)-(randomerlen-1)
            if isvalid(seq[j:j+(randomerlen-1)],blocking_sites)
                push!(blocking_sites,seq[j:j+(randomerlen-1)])
            end
        end
    end
    return blocking_sites
end

function edit_distance(seq1::LongSequence{DNAAlphabet{4}},seq2::LongSequence{DNAAlphabet{4}})
    m=length(seq1);
    n=length(seq2);
    d=zeros(m+1,n+1)
    for i = 1:m
        d[i+1,1]=i
    end
    for j = 1:n
        d[1,j+1]=j
    end
    for j = 1:n
        for i=1:m
            substitution_cost=1
            if iscompatible(seq1[i],seq2[j])
                substitution_cost=0;
            end
            d[i+1,j+1]= min(d[i,j+1]+1,d[i+1,j]+1,d[i,j]+substitution_cost)
        end
    end
    return d[m+1,n+1]
end

function hamming_distance(seq1::LongSequence{DNAAlphabet{4}},seq2::LongSequence{DNAAlphabet{4}})
    dist_counter=0
    if  length(seq1) != length(seq2)
        error("Sequences must be the same length.")
    end
    n=length(seq1)
    for i=1:n
        if !iscompatible(seq1[i],seq2[i])
            dist_counter+=1
        end
    end
    return dist_counter
end
function reward_function_factory(barcode_pool)
    function reward_function(barcode;distance=hamming_distance)
        sum_distances=0
        min_distance=length(barcode)
        for code in barcode_pool
            dist=distance(barcode,code)
            sum_distances+=dist
            min_distance=min(min_distance,dist)

        end
        if is_palindrome(barcode)
            sum_distances=0;
            min_distance=0;
        end
        return mean(sum_distances), min_distance
    end
    return reward_function
end

function generate_available_ns(base_randomer::LongSequence{DNAAlphabet{4}})
    randomerlen=length(base_randomer)
    basedict=Dict(DNA_A => dna"A", DNA_T => dna"T", DNA_C => dna"C", DNA_G => dna"G",DNA_R =>dna"AG",DNA_Y=>dna"CT",DNA_S=>dna"GC",
    DNA_W=>dna"AT",DNA_K=>dna"GT",DNA_M=>dna"AC",DNA_B=>dna"CGT",DNA_D=>dna"AGT",DNA_H=>dna"ACT",DNA_V=>dna"ACG",DNA_N=>dna"ACGT",DNA_Gap=>dna"-")

    available_ns=[dna"-" for _ in 1:randomerlen]
    for i=1:randomerlen
        available_ns[i]=basedict[base_randomer[i]]
    end
    return available_ns
end

function generate_randomer_sample(available_ns)
    randomerlen=length(available_ns);
    randomer=dna"-"^randomerlen
    for i =1:randomerlen
        randomer[i]=shuffle(available_ns[i])[1]
    end
    return randomer
end

function simulate_random(prefix, bases, sites; nsims=1000, horizon=length(bases),reward=reward_function,kwargs...)
    score = zeros(nsims)
    min_distance=zeros(nsims)
    true_horizon = min(horizon, length(bases)) # don't go beyond the randomer length
    for i = 1:nsims
        score[i],min_distance[i] = reward(find_valid_randomer(prefix, bases[1:true_horizon], sites))
    end
    return mean(score),mean(min_distance)
end

function oligo_rollout(bases, sites; simulate=simulate_random, kwargs...)
    n = length(bases)
    randomer = dna"-" ^ n

    # Find the longest restriction site and use this as the horizon.
    max_len = map(length, sites) |> maximum
    horizon = max_len

    for i = 1:n
        best_base = DNA_Gap
        best_min_distance = 0
        best_score=0
        start = max(i - max_len, 1) # only need to look back max_len when checking the oligos
        for base in bases[i]
            if ~isvalid(randomer[start:i-1] * base, sites)
                # adding this bases creates a restriction site; skip it
                continue
            end
            score,min_distance = simulate(randomer[start:i-1] * base, bases[i+1:end], sites; horizon=horizon, kwargs...)

            # Update use this base if:
            #   1. the mean randomer degeneracy is higher than the previous best, OR
            #   2. the randomer degeneracy is tied with the previous best, but the new base has
            #      higher individual degeneracy
            if min_distance > best_min_distance || (min_distance== best_min_distance && score > best_score)
                best_min_distance=min_distance
                best_score= score
                best_base = base
            end
        end

        if best_base == DNA_Gap
            # The sequence has terminated; there are no bases that can be added without making
            # a valid barcode
            break
        else
            randomer[i] = best_base
        end
    end

    return randomer
end

function BarcodeRL(sites;nbarcodes=16,barcodelen=8,kwargs...)
    barcode_pool=[];
    base_barcode_ns=dna"N"^barcodelen
    available_ns=generate_available_ns(base_barcode_ns)
    for i= 1:nbarcodes
        reward_function=reward_function_factory(barcode_pool)
        barcode=oligo_rollout(available_ns,sites;reward=reward_function, kwargs...)
        push!(barcode_pool,barcode)
        push!(sites,barcode)
        push!(sites,reverse_complement(barcode))
    end
    reward_function=reward_function_factory(barcode_pool[1:end-1])
    score,min_distance=reward_function(barcode_pool[end];distance=hamming_distance)
    return barcode_pool,score,min_distance
end

function generate_initial_population(nbarcodes,barcode_form,population_size)
    barcodelen=length(barcode_form)
    initial_population=[[dna"-"^barcodelen for i =1:nbarcodes] for j=1:population_size]
    available_ns=generate_available_ns(barcode_form)
    for j=1:population_size
        for k=1:nbarcodes
            for i=1:barcodelen
                initial_population[j][k][i]=shuffle(available_ns[i])[1]
            end
        end
    end
    return initial_population

end

function fitness(barcodes,sites,x_combs)
    n=length(barcodes)
    fits=zeros(length(x_combs))
    for i=1:length(x_combs)
        fits[i]=hamming_distance(barcodes[x_combs[i][1]],barcodes[x_combs[i][2]])
    end
    fitness=minimum(fits)
    for i=1:length(barcodes)
        if ~isvalid(barcodes[i],vcat(sites,barcodes[union(1:i-1,i+1:n)]))
            fitness=0
            break
        end
    end
    return fitness
end

function mutation(barcodes,sites,barcode_form;mutation_prob=0.25,kwargs...)
    new_barcodes=deepcopy(barcodes)
    n=length(new_barcodes)
    mutation_sites=findall(y-> y==1,rand(n).<mutation_prob)
    available_ns=generate_available_ns(barcode_form)
    order=shuffle(mutation_sites)
    for i in order
        barcode_pool=new_barcodes[union(1:i-1,i+1:n)]
        reward_function=reward_function_factory(barcode_pool)
        new_barcodes[i]=oligo_rollout(available_ns,sites;reward=reward_function,kwargs...)
    end
    return new_barcodes
end
#barcodes=[dna"ATGCGT",dna"ATGGGG",dna"GCACGT",dna"ATTGAT"]
#sites=read_invalid_sites(CSV.read("./Invalid_Sequences_Test.csv",DataFrame);randomerlen=6)
#test=mutation(barcodes,sites)

function pair_off(candidates)
    npairs=Int(floor(length(candidates)/2))
    pairs=[zeros(Int64,2) for i =1:npairs]
    cands=shuffle(candidates)
    for i =1:npairs
        pairs[i][1]=cands[2*i-1]
        pairs[i][2]=cands[2*i]
    end
    return pairs

end
function recombination(population;recombination_prob=0.25,kwargs...)
    new_population=deepcopy(population)
    n_barcodes=length(population[1])
    n_chromosomes=length(population);
    recombination_sites=findall(y-> y==1, rand(n_chromosomes).< recombination_prob)
    r_pairs=pair_off(recombination_sites)
    for i=1:length(r_pairs)
        r_index=Int(ceil(rand()*n_barcodes))
        new_population[r_pairs[i][1]][r_index:end],new_population[r_pairs[i][2]][r_index:end]=new_population[r_pairs[i][2]][r_index:end],new_population[r_pairs[i][1]][r_index:end];
    end

    return new_population[recombination_sites]
end
#Test
#population=[[dna"AAAA",dna"AAAA",dna"AAAA",dna"AAAA"],[dna"GGGG",dna"GGGG",dna"GGGG",dna"GGGG"]]
#recombination(population;recombination_prob=1)

function rank_probability(n_chromosomes;Pc=0.25,kwargs...)
    probs=zeros(n_chromosomes)
    cumul_prob=zeros(n_chromosomes)
    for i=1:n_chromosomes-1
        probs[i]=Pc*(1-Pc)^(i-1)
        cumul_prob[i]=sum(probs[1:i])
    end
    probs[end]=(1-Pc)^(n_chromosomes-1)
    cumul_prob[end]=sum(probs)
    return cumul_prob
end

function selection(population,fitnesses;kwargs...)
    n_chromosomes=length(population)
    ranks=sortperm(fitnesses,rev=true)
    cumul_prob=rank_probability(n_chromosomes;kwargs...)
    sample=rand()
    cutoff=minimum(findall(y-> y>=sample,cumul_prob));
    idx=ranks[cutoff]
    return population[idx]
end
#Test
#population=[[dna"AAAA",dna"AAAA",dna"AAAA",dna"AAAA"],[dna"GGGG",dna"GGGG",dna"GGGG",dna"GGGG"]]
#fitnesses=[20, 10]
#selection(population,fitnesses)

#results=zeros(Int64,100000)
#for i=1:100000
#    results[i]=selection(population,fitnesses)
#end
#x=length(findall(y->y==1,results))
#t=length(findall(y->y==2,results))
#thy=length(findall(y->y==30,results))
#thy/100000

function GA_BarcodeRL(sites;nbarcodes=16,barcode_form=dna"NNNNNNNN",generations=50,population_size=30,kwargs...)
    x_comb=collect(combinations(1:nbarcodes,2));
    population=generate_initial_population(nbarcodes,barcode_form,population_size)
    for i=1:generations
        for j=1:population_size
            population=push!(population,mutation(population[j],sites,barcode_form;kwargs...))
        end
        recomb_pop=recombination(population;kwargs...)
        for j =1:population_size-1
            recomb_pop=vcat(recomb_pop,recombination(population;kwargs...))
        end
        population=vcat(population,recomb_pop)
        fitnesses=fitness.(population,(sites,),(x_comb,))
        new_pop_indicies=zeros(Int64,population_size)
        v_indicies=collect(1:length(population))
        for i=1:population_size
            new_pop_cand=selection(population[v_indicies],fitnesses[v_indicies];kwargs...)
            new_pop_indicies[i]=findall(y-> y == new_pop_cand,population)[1]
            v_indicies=filter(y->y != new_pop_indicies[i],v_indicies)
        end
        population=population[new_pop_indicies]
        println("Generation: $i")
    end

    final_fitness=fitness.(population,(sites,),(x_comb,))
    max_fitness,max_pop_idx=findmax(final_fitness)
    return population[max_pop_idx], max_fitness
end

sites=read_invalid_sites(CSV.read("./Invalid_Sequences_Test.csv",DataFrame);randomerlen=8)

barcodes,best_fitness=GA_BarcodeRL(sites;nbarcodes=16,barcode_form=dna"NNNNNNNN",generations=10,population_size=30,nsims=100,Pc=0.25,recombination_prob=0.25,mutation_prob=0.25)

function benchmark_GA_barcodeRL()
    num_barcodes=repeat([4 ,8 ,12, 16, 20, 24, 28, 32, 36, 40, 44, 48],3)
    nexps=length(num_barcodes)
    time=zeros(nexps)
    min_edit_distance=zeros(nexps)
    sites=read_invalid_sites(CSV.read("./Invalid_Sequences_Test.csv",DataFrame);randomerlen=8)
    barcodes=["-" for i=1:nexps]
    for i =1:nexps
        data= @timed GA_BarcodeRL(sites;nbarcodes=num_barcodes[i],barcode_form=dna"NNNNNNNN",generations=10,population_size=50,nsims=100,Pc=0.25,recombination_prob=0.25,mutation_prob=0.25)
        time[i]=data.time
        min_edit_distance[i]=data.value[2]
        barcodes[i]=join(convert.(String,data.value[1]),",")
    end
    DF=DataFrame(num_barcodes=num_barcodes,barcodes=barcodes,min_edit_distance=min_edit_distance,time=time)
    return DF
end

data=benchmark_GA_barcodeRL()
CSV.write("GA_BarcodeRL_benchmark_replicates.csv",data)


text:
list    = include(Pkg.dir("Colors/src/names_data.jl"))

palette = ASCIIString[]

for (k,v) in list
    push!(palette, k)
end

const default_colors = ["dodgerblue", "red1", "goldenrod", "wheat3", "olive", "burlywood"]


text:
import Base.rand

export BallInf

"""
    BallInf{N<:Real} <: AbstractHyperrectangle{N}

Type that represents a ball in the infinity norm.

### Fields

- `center` -- center of the ball as a real vector
- `radius` -- radius of the ball as a real scalar (``≥ 0``)

### Notes

Mathematically, a ball in the infinity norm is defined as the set

```math
\\mathcal{B}_∞^n(c, r) = \\{ x ∈ \\mathbb{R}^n : ‖ x - c ‖_∞ ≤ r \\},
```
where ``c ∈ \\mathbb{R}^n`` is its center and ``r ∈ \\mathbb{R}_+`` its radius.
Here ``‖ ⋅ ‖_∞`` denotes the infinity norm, defined as
``‖ x ‖_∞ = \\max\\limits_{i=1,…,n} \\vert x_i \\vert`` for any
``x ∈ \\mathbb{R}^n``.

### Examples

Create the two-dimensional unit ball and compute its support function along the
positive ``x=y`` direction:

```jldoctest
julia> B = BallInf(zeros(2), 1.0)
BallInf{Float64}([0.0, 0.0], 1.0)
julia> dim(B)
2
julia> ρ([1., 1.], B)
2.0
```
"""
struct BallInf{N<:Real} <: AbstractHyperrectangle{N}
    center::Vector{N}
    radius::N

    # default constructor with domain constraint for radius
    function BallInf{N}(center::Vector{N}, radius::N) where {N<:Real}
        @assert radius >= zero(N) "radius must not be negative"
        return new{N}(center, radius)
    end
end

# convenience constructor without type parameter
BallInf(center::Vector{N}, radius::N) where {N<:Real} =
    BallInf{N}(center, radius)

# --- AbstractHyperrectangle interface functions ---

"""
    radius_hyperrectangle(B::BallInf{N}, i::Int)::N where {N<:Real}

Return the box radius of a infinity norm ball in a given dimension.

### Input

- `B` -- infinity norm ball
- `i` -- dimension of interest

### Output

The box radius of the ball in the infinity norm in the given dimension.
"""
function radius_hyperrectangle(B::BallInf{N}, i::Int)::N where {N<:Real}
    return B.radius
end

"""
    radius_hyperrectangle(B::BallInf{N})::Vector{N} where {N<:Real}

Return the box radius of a infinity norm ball, which is the same in every
dimension.

### Input

- `B` -- infinity norm ball

### Output

The box radius of the ball in the infinity norm.
"""
function radius_hyperrectangle(B::BallInf{N})::Vector{N} where {N<:Real}
    return fill(B.radius, dim(B))
end

# --- AbstractCentrallySymmetric interface functions ---

"""
    center(B::BallInf{N})::Vector{N} where {N<:Real}

Return the center of a ball in the infinity norm.

### Input

- `B` -- ball in the infinity norm

### Output

The center of the ball in the infinity norm.
"""
function center(B::BallInf{N})::Vector{N} where {N<:Real}
    return B.center
end

# --- LazySet interface functions ---

"""
    radius(B::BallInf, [p]::Real=Inf)::Real

Return the radius of a ball in the infinity norm.

### Input

- `B` -- ball in the infinity norm
- `p` -- (optional, default: `Inf`) norm

### Output

A real number representing the radius.

### Notes

The radius is defined as the radius of the enclosing ball of the given
``p``-norm of minimal volume with the same center.
"""
function radius(B::BallInf, p::Real=Inf)::Real
    return (p == Inf) ? B.radius : norm(fill(B.radius, dim(B)), p)
end

"""
    rand(::Type{BallInf}; [N]::Type{<:Real}=Float64, [dim]::Int=2,
         [rng]::AbstractRNG=GLOBAL_RNG, [seed]::Union{Int, Nothing}=nothing
        )::BallInf{N}

Create a random ball in the infinity norm.

### Input

- `BallInf` -- type for dispatch
- `N`       -- (optional, default: `Float64`) numeric type
- `dim`     -- (optional, default: 2) dimension
- `rng`     -- (optional, default: `GLOBAL_RNG`) random number generator
- `seed`    -- (optional, default: `nothing`) seed for reseeding

### Output

A random ball in the infinity norm.

### Algorithm

All numbers are normally distributed with mean 0 and standard deviation 1.
Additionally, the radius is nonnegative.
"""
function rand(::Type{BallInf};
              N::Type{<:Real}=Float64,
              dim::Int=2,
              rng::AbstractRNG=GLOBAL_RNG,
              seed::Union{Int, Nothing}=nothing
             )::BallInf{N}
    rng = reseed(rng, seed)
    center = randn(rng, N, dim)
    radius = abs(randn(rng, N))
    return BallInf(center, radius)
end

"""
    translate(B::BallInf{N}, v::AbstractVector{N}) where {N<:Real}

Translate (i.e., shift) a ball in the infinity norm by a given vector.

### Input

- `B` -- ball in the infinity norm
- `v` -- translation vector

### Output

A translated ball in the infinity norm.

### Algorithm

We add the vector to the center of the ball.
"""
function translate(B::BallInf{N}, v::AbstractVector{N}) where {N<:Real}
    @assert length(v) == dim(B) "cannot translate a $(dim(B))-dimensional " *
                                "set by a $(length(v))-dimensional vector"
    return BallInf(center(B) + v, B.radius)
end


text:
export SparseRadialMapComponent, component, construct, clearcoeff!,
       extract_offcoeff, extract_diagcoeff, evaluate, negative_likelihood, off_diagonal,
       set_id

#### Structure for the Nx-th component SparseRadialMapComponent of the lower triangular map U

# For now RadialMapComponent = Σi=1,k ui(k)

# The last array for ξ, σ and a will be of size p+2, p+2, p+3 versus p, p, p+1
# for the components i=1,...,Nx-1

struct SparseRadialMapComponent
        Nx::Int64
        p::Array{Int64,1}
        activedim::Array{Int64, 1}
        ξ::Array{Array{Float64,1}, 1}
        σ::Array{Array{Float64,1}, 1}
        coeff::Array{Array{Float64,1}, 1}

        # Cache for evaluation of function, gradient and hessian on basis

        # Inner constructor
        function SparseRadialMapComponent(Nx::Int64, p::Array{Int64,1}, activedim::Array{Int64,1}, ξ::Array{Array{Float64,1},1}, σ::Array{Array{Float64,1},1}, coeff::Array{Array{Float64,1},1})

        @assert size(activedim,1) <= Nx "Size of activedim should be smaller than the order of SparseRadialMapComponent"
        @assert size(ξ,1)==Nx "Size of ξ doesn't match the order of SparseRadialMapComponent"
        @assert size(σ,1)==Nx "Size of σ doesn't match the order of SparseRadialMapComponent"
        @assert size(coeff,1)==Nx "Size of coeff doesn't match the order of SparseRadialMapComponent"

        # activedim = filter(x-> x != -1, p)
                return new(Nx, p, activedim, ξ, σ, coeff)
        end
end

function SparseRadialMapComponent(Nx::Int64, p::Array{Int64,1})
        @assert Nx>=1 "Nx must be >0"
        activedim = Int64[]
        @inbounds for i=1:Nx
                if p[i] != -1
                        push!(activedim, copy(i))
                end
        end

        if Nx==1
                if p[1]==-1
                return SparseRadialMapComponent(Nx, p, activedim, [Float64[]], [Float64[]],  [Float64[]])
                end
                if p[1]==0
                return SparseRadialMapComponent(Nx, p, activedim, [Float64[]], [Float64[]],  [zeros(2)])
                end
                if p[1]>0
                return SparseRadialMapComponent(Nx, p, activedim, [zeros(p[1]+2)], [ones(p[1]+2)], [zeros(p[1]+3)])
                end

        else # Nx>1
                ξ = Array{Float64,1}[]
                σ = Array{Float64,1}[]
                coeff =  Array{Float64,1}[]
                # Fill Off Diagonal components
                for pj in p[1:end-1]
                if pj==-1
                        # Null component
                        push!(ξ, Float64[])
                        push!(σ, Float64[])
                        push!(coeff,Float64[])
                elseif pj==0
                        # Linear function
                        push!(ξ, Float64[])
                        push!(σ, Float64[])
                        push!(coeff,zeros(1))
                else    # Linear function + p rbf functions
                        push!(ξ, zeros(pj))
                        push!(σ, ones(pj))
                        push!(coeff,zeros(pj+1))
                end
                end

                # Fill Diagonal Component
                if p[end]==-1
                        # Null component
                        push!(ξ, Float64[])
                        push!(σ, Float64[])
                        push!(coeff,Float64[])
                elseif p[end]==0
                        # Affine function
                        push!(ξ, Float64[])
                        push!(σ, Float64[])
                        push!(coeff,zeros(2))
                else
                        #Constant + p+2 ψ functions
                        push!(ξ, zeros(p[end]+2))
                        push!(σ, ones(p[end]+2))
                        push!(coeff,zeros(p[end]+3))
                end
                return SparseRadialMapComponent(Nx, p, activedim, ξ, σ, coeff)
        end
end

SparseRadialMapComponent(Nx::Int64, p::Int64) = SparseRadialMapComponent(Nx, fill(p, Nx))

function show(io::IO, C::SparseRadialMapComponent)
    println(io,"Sparse radial map component of dimension "*string(C.Nx)*" and order p = $(C.p)")
end

function set_id(C::SparseRadialMapComponent)
        Nx = C.Nx
if Nx==1
        C.p[1] = -1
        C.ξ[1] = Float64[]
        C.σ[1] = Float64[]
        C.coeff[1] = Float64[]
else
        C.p .= rmul!(ones(Nx), -1.0)

        @inbounds for i=1:Nx
        C.ξ[i] = Float64[]
        C.σ[i] = Float64[]
        C.coeff[i] = Float64[]
        end
end
end

# Transform vector of coefficients to a form SparseRadialMapComponent.coeff
function modifycoeff!(C::SparseRadialMapComponent, A::Array{Float64,1})
        @assert size(A,1)==size(extractcoeff(C),1) "A doesn't match the size of C"
        @get C (Nx, p)
                count = 0
        if Nx==1
                C.coeff[1] .= A
                count += size(A,1)
                # C.coeff[1] .= A
        else
                # count will store the last position of each assignment

                # Off-diagonal components
                @inbounds for idx=1:Nx-1
                pidx = p[idx]
                if pidx==-1
                        #No assignment
                elseif pidx==0
                        C.coeff[idx] .= [A[count+1]]
                        count += 1
                else
                        tmp_coeffi  = view(C.coeff,idx)[1]
                        tmp_coeffi .= view(A, count+1:count+(pidx+1))
                        count += deepcopy(pidx)+1
                end
                end

                # Diagonal components
                pidx = p[Nx]
                if pidx==-1
                        #No assignment, set to the identity without coefficient
                elseif pidx==0
                        # Affine function with two coefficients
                        tmp_coeffk  = view(C.coeff,Nx)[1]
                        tmp_coeffk .= view(A, count+1:count+2)

                        count +=2
                else
                        # Constant + p+2 ψ functions
                        tmp_coeffk  = view(C.coeff,Nx)[1]
                        tmp_coeffk .= view(A, count+1:count+pidx+3)
                        count += deepcopy(pidx)+3
                        # C.coeff[Nx] .= view(A, (Nx-1)*(p+1)+1:Nx*(p+1)+2)
                        # C.coeff[Nx] .= A[(Nx-1)*(p+1)+1:Nx*(p+1)+2]
                end
        end
        @assert count == size(A,1) "Error in the assignment of the coefficients"
end

function clearcoeff!(C::SparseRadialMapComponent)
        @inbounds for i in C.activedim
                fill!(C.coeff[i], 0.0)
        end
end

extractcoeff(C::SparseRadialMapComponent) = vcat(C.coeff...)

function extract_offcoeff(C::SparseRadialMapComponent)
        if C.Nx == 1
                return Float64[]
        else
                return  vcat(C.coeff[1:C.Nx-1]...)
        end
end

extract_diagcoeff(C::SparseRadialMapComponent) = C.coeff[end]

size(C::SparseRadialMapComponent) = (C.Nx, C.p)

# Evaluate the map SparseRadialMapComponent at z = (x,...,x)
function (C::SparseRadialMapComponent)(z::T) where {T<:Real}
        @get C (Nx, p)

        out = uk(p[Nx], C.ξ[Nx], C.σ[Nx], C.coeff[Nx])(z)
        if Nx>1
                @inbounds for i in intersect(1:Nx-1, C.activedim)
                        out += ui(p[i], C.ξ[i], C.σ[i], C.coeff[i])(z)
                end
        end
        return out
end

# Evaluate the map SparseRadialMapComponent at z = (z1,...,zNx)
# Optimized with views
function (C::SparseRadialMapComponent)(z)
        @get C (Nx, p)
        @assert size(z,1)<=Nx  "The vector z has more components than RadialMapComponent"
        if p[Nx]<1
                out = uk(p[Nx], C.ξ[Nx], C.σ[Nx], C.coeff[Nx])(z[Nx])
        else
                out = uk(p[Nx], view(C.ξ,Nx)[1], view(C.σ,Nx)[1], view(C.coeff,Nx)[1])(z[Nx])
        end

        if Nx>1
                @inbounds for idx in intersect(1:Nx-1, C.activedim)
                        pidx = p[idx]

                        if pidx==-1
                        # ui(z) = 0.0 can skip this
                        elseif pidx==0
                        out += ui(pidx, C.ξ[idx], C.σ[idx], C.coeff[idx])(z[idx])
                        else
                        out += ui(pidx, view(C.ξ,idx)[1], view(C.σ,idx)[1], view(C.coeff,idx)[1])(z[idx])
                        end
                end
        end
        return out
end

# Evaluate the map C at z = (z1,...,zk)
# Optimized with views
function off_diagonal(C::SparseRadialMapComponent, z)
        @get C (Nx, p)
        @assert size(z,1)<=Nx  "The vector z has more components than RadialMapComponent"
        out = 0.0
        if Nx>1
        @inbounds for idx in intersect(1:Nx-1, C.activedim)
                        pidx = p[idx]
                        if pidx==-1
                                # No computation
                        elseif pidx==0
                        out += ui(pidx, C.ξ[idx], C.σ[idx], C.coeff[idx])(z[idx])
                        else
                        out += ui(pidx, view(C.ξ,idx)[1], view(C.σ,idx)[1], view(C.coeff,idx)[1])(z[idx])
                        end
                end
        end
        return out
end

function D!(C::SparseRadialMapComponent, z::AbstractVector{Float64})
        @get C (Nx, p)
        @assert size(z,1) == Nx  "The vector z has more components than the SparseRadialMapComponent"
        return D!(uk(p[Nx], C.ξ[Nx], C.σ[Nx], C.coeff[Nx]), z[Nx])
end

D(C::SparseRadialMapComponent) =  z-> D!(C, z)

# function negative_likelihood(C::RadialMapComponent, X::AbstractMatrix{Float64}, λ::Float64, δ::Float64)
#         @get C (Nx, p)
#         J = 0.0
#         NxX, Ne = size(X)
#
#         @assert NxX == Nx "Wrong dimension of the ensemble matrix `X`"
#         ∂kC = D(C)
#         @inbounds for i=1:Ne
#                 col = view(X,:,i)
#                 # Quadratic term
#                 J += 0.5*C(col)^2
#                 # Log barrier term
#                 J -= log(∂kC(col))
#         end
#         J *= 1/Ne
#         return J
# end

function negative_likelihood(C::SparseRadialMapComponent, X::AbstractMatrix{Float64}, λ::Float64, δ::Float64)
        @get C (Nx, p)
        J = 0.0
        NxX, Ne = size(X)
        @assert NxX == Nx "Wrong dimension of the ensemble matrix `X`"

        ψ_off, ψ_diag, dψ_diag = compute_weights(C, X)
        n_off = size(ψ_off,2)
        # ∂kC = D(C)
        x_diag = extract_diagcoeff(C)

        ## Quadratic term
        if n_off>0
                x_off = extract_offcoeff(C)

                J = 0.5*norm(ψ_off*x_off + ψ_diag*x_diag[2:end] .+ x_diag[1])^2
        else
                J = 0.5*norm(ψ_diag*x_diag[2:end] .+ x_diag[1])^2
        end

        ##  Log barrier term

        # Penalty on the log barrier temr
        δlog = reshape(δ*sum(dψ_diag, dims = 2), Ne)
        for i=1:Ne
                J -= log(dot(view(dψ_diag,i,:), x_diag[2:end]) + δlog[i])
        end
        J *= 1/Ne

        # Penalty term
        # if n_off>0
        #         J += δ*dot(x_off, sum(ψ_off'*ψ_off, dims = 1))
        # end
        #
        # J += δ*dot(x_diag[2:end], sum(ψ_diag'*ψ_diag, dims = 1))
        # @inbounds for i=1:Ne
        #         col = view(X,:,i)
        #         # Quadratic term
        #         J += 0.5*(C(col))^2
        #         # Log barrier term
        #         J -= log(∂kC(col))
        # end
        return J
end

function component(C::SparseRadialMapComponent, idx::Int64)
        @get C (Nx, p)
        @assert Nx>= idx "Can't access this component idx>Nx"
        if p[idx]<1
                if idx==Nx
                        return uk(p[idx], C.ξ[idx], C.σ[idx], C.coeff[idx])
                else
                        return ui(p[idx], C.ξ[idx], C.σ[idx], C.coeff[idx])
                end
        else
                if idx==Nx
                        return uk(p[idx], C.ξ[idx], C.σ[idx], C.coeff[idx])
                else
                        return ui(p[idx], C.ξ[idx], C.σ[idx], C.coeff[idx])
                end
        end
end


text:
include("../01_perturbation_model/perturbation_model.jl")

function sqrt_epsilon_fun(ω, z, t, z1, z2, zmax, dN, z0, d_z, t_r, t_d, tau_r, tau_d)
    N_model = N_model_fun(z, t, z1, z2, zmax, dN, z0, d_z, t_r, t_d, tau_r, tau_d)
    ωₚ = 2pi*Ne_to_fpe(N_model)
    if ωₚ > ω
        return NaN
    end
    return sqrt(1 - ωₚ^2/ω^2)
end

function sqrt_eps_integral(z_step, ω, t, z1, z2, zmax, dN, z0, d_z, t_r, t_d, tau_r, tau_d)
    z = z1
    S = sqrt_epsilon_fun(ω, z, t, z1, z2, zmax, dN, z0, d_z, t_r, t_d, tau_r, tau_d)/2
    while true
        z = z + z_step
        s = sqrt_epsilon_fun(ω, z, t, z1, z2, zmax, dN, z0, d_z, t_r, t_d, tau_r, tau_d)
        if isnan(s)
            S = S - sqrt_epsilon_fun(ω, z-z_step, t, z1, z2, zmax, dN, z0, d_z, t_r, t_d, tau_r, tau_d)/2
            break
        end
        S = S + s
    end
    return S*z_step
end

# sqrt_eps_integral(1000., 45e5, 5., z1, z2, zmax, dN, z0, d_z, t_r, t_d, tau_r, tau_d)

function phase_variations(ω, t1, t2, z_step, z1, z2, zmax, dN, z0, d_z, t_r, t_d, tau_r, tau_d)
    c = 299792458.
    phase2 = 2ω/c * sqrt_eps_integral(z_step, ω, t2, z1, z2, zmax, dN, z0, d_z, t_r, t_d, tau_r, tau_d)
    phase1 = 2ω/c * sqrt_eps_integral(z_step, ω, t1, z1, z2, zmax, dN, z0, d_z, t_r, t_d, tau_r, tau_d)
    return phase2 - phase1
end

# f_axe = 4e6:1000.:5e6
# ω_axe = 2pi*f_axe
# t_axe = 0.0:0.1:5.
# z_step = 10.
# DPhi = zeros(length(t_axe),length(f_axe))

# for i=2:length(t_axe)
#     println(i)
#     for j=1:length(f_axe)
#         # println(i, " ",j)
#         DPhi[i,j] = phase_variations(ω_axe[j], t_axe[i-1], t_axe[i], z_step, z1, z2, zmax, dN, z0, d_z, t_r, t_d, tau_r, tau_d)
#     end
# end

# using Plots

# i=12; plot(f_axe/1e6, DPhi[i,:])
# title(string("t=",t_axe[i]))
# ylims((0,15))


text:
using Documenter
using cNODE

makedocs(
    sitename = "cNODE.jl",
    modules = [cNODE],
    pages = [
            "Home" => "index.md",
            "Generate Data" => "generator.md",
            "Load Data" => "loader.md",
            "Use cNODE" => "trainer.md"],
    format = Documenter.HTML()
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    repo = "github.com/michel-mata/cNODE.jl.git",
    devbranch = "main",
    target = "build",
    versions = ["stable" => "v^"]
)


text:
const ID = Union{Integer, UUID}
const NameOrId = Union{AStr, ID}

"""
    get_user(::Forge[, name_or_id::Union{$AStr, Integer, UUID}])

Get the currently authenticated user, or a user by name or ID.
"""
@endpoint get_user()
@endpoint get_user(name_or_id::NameOrId)

"""
    get_users(::Forge)

Get all users.
"""
@endpoint get_users()

"""
    update_user(::Forge[, id::Union{Integer, UUID}]; kwargs...)

Update the currently authenticated user, or a user by ID.
"""
@endpoint update_user()
@endpoint update_user(id::ID)

"""
    create_user(::Forge; kwargs...)

Create a new user.
"""
@endpoint create_user()

"""
    delete_user(::Forge, id::Union{Integer, UUID})

Delete a user by ID.
"""
@endpoint delete_user(id::ID)

"""
    get_user_repos(::Forge[, name_or_id::Union{$AStr, Integer, UUID}])

Get the currently authenticated user's repositories, or those of a user by name or ID.
"""
@endpoint get_user_repos()
@endpoint get_user_repos(name_or_id::NameOrId)

"""
    get_repo(::Forge, owner_repo::$AStr)
    get_repo(::Forge, owner::$AStr, repo::$AStr)
    get_repo(::Forge, id::Integer)
    get_repo(::Forge, id::UUID)
    get_repo(::Forge, owner::$AStr, subgroup::$AStr, repo::$AStr)

Get a repository by owner and name or ID.
"""
@endpoint get_repo(owner_repo::AStr)
@endpoint get_repo(owner::AStr, repo::AStr)
@endpoint get_repo(id::Integer)
@endpoint get_repo(id::UUID)
@endpoint get_repo(owner::AStr, subgroup::AStr, repo::AStr)

"""
    create_repo(::Forge, owner::$AStr)
    create_repo(::Forge)
Create a repository.
If using GitHub and you want to create a repository in an organization, pass the organization name as argument.
"""
@endpoint create_repo(owner::AStr, repo::AStr)
@endpoint create_repo(owner::NameOrId)
@endpoint create_repo()

"""
    get_branch(::Forge, owner::AbstractString, repo::AbstractString, branch::AbstractString)

Get a branch from a repository.
"""
@endpoint get_branch(owner::AStr, repo::AStr, branch::AStr)

"""
    get_branches(::Forge, owner::AbstractString, repo::AbstractString)

Get all branches for a repository.
"""
@endpoint get_branches(owner::AStr, repo::AStr)

"""
    delete_branch(::Forge, owner::AbstractString, repo::AbstractString, branch::AbstractString)

Delete a branch from a repository.
"""
@endpoint delete_branch(owner::AStr, repo::AStr, branch::AStr)

"""
    get_file_contents(
        ::Forge,
        owner::$AStr,
        repo::$AStr,
        path::$AStr,
    )
    get_file_contents(f::Forge, id::Integer, path::$AStr)

Get a file from a repository.
"""
@endpoint get_file_contents(owner::AStr, repo::AStr, path::AStr)
@endpoint get_file_contents(id::Integer, path::AStr)

"""
    get_pull_request(::Forge, owner::$AStr, repo::$AStr, number::Integer)
    get_pull_request(::Forge, project::Integer, number::Integer)

Get a specific pull request.
"""
@endpoint get_pull_request(owner::AStr, repo::AStr, number::Integer)
@endpoint get_pull_request(project::Integer, number::Integer)

"""
    get_pull_requests(::Forge, owner::$AStr, repo::$AStr)
    get_pull_requests(::Forge, project::Integer)

List a repository's pull requests.
"""
@endpoint get_pull_requests(owner::AStr, repo::AStr)
@endpoint get_pull_requests(repo::Integer)

"""
    create_pull_request(::Forge, owner::$AStr, repo::$AStr; kwargs...)
    create_pull_request(::Forge, project::Integer; kwargs...)

Create a pull request.
"""
@endpoint create_pull_request(owner::AStr, repo::AStr)
@endpoint create_pull_request(project::Integer)

"""
    update_pull_request(::Forge, owner::$AStr, repo::$AStr, number::Integer; kwargs...)
    update_pull_request(::Forge, project::Integer, number::Integer; kwargs...)

Update a pull request.
"""
@endpoint update_pull_request(owner::AStr, repo::AStr, number::Integer)
@endpoint update_pull_request(project::Integer, number::Integer)

"""
    subscribe_to_pull_request(::Forge, project::Integer, pull_request_id::Integer)

Subscribe to notifications for the supplied project and pull request ID.
"""
@endpoint subscribe_to_pull_request(project::Integer, pull_request_id::Integer)

"""
    unsubscribe_from_pull_request(::Forge, project::Integer, pull_request_id::Integer)

Unsubscribe from notifications for the supplied project and pull request ID.
"""
@endpoint unsubscribe_from_pull_request(project::Integer, pull_request_id::Integer)

"""
    get_commit(::Forge, owner::$AStr, repo::$AStr, ref::$AStr)
    get_commit(::Forge, project::Integer, ref::$AStr)

Get a commit from a repository.
"""
@endpoint get_commit(owner::AStr, repo::AStr, ref::AStr)
@endpoint get_commit(project::Integer, ref::AStr)

"""
    is_collaborator(
        ::Forge,
        owner::$AStr,
        repo::$AStr,
        name_or_id::Union{$AStr, Integer, UUID},
    )

Check whether or not a user is a collaborator on a repository.
"""
@endpoint is_collaborator(owner::AStr, repo::AStr)
@endpoint is_collaborator(owner::AStr, repo::AStr, name_or_id::NameOrId)

"""
    is_member(::Forge, org::$AStr, name_or_id::Union{$AStr, Integer, UUID})

Check whether or not a user is a member of an organization.
"""
@endpoint is_member(org::AStr, name_or_id::NameOrId)

"""
    groups(::Forge)
Search for groups in GitLab.
"""
@endpoint groups()

"""
    get_tags(::Forge, owner::$AStr, repo::$AStr)
    get_tags(::Forge, project::Integer)

Get a list of tags from a repository.
"""
@endpoint get_tags(owner::AStr, repo::AStr)
@endpoint get_tags(project::Integer)

"""
    list_pull_request_comments(
        ::Forge, owner::$AStr, repo::$AStr, pull_request_id::Integer;
        kwargs...
    )
    list_pull_request_comments(
        ::Forge, project::Integer, pull_request_id::Integer;
        kwargs...
    )

Get a list of comments for a given pull request.
"""
@endpoint list_pull_request_comments(owner::AStr, repo::AStr, pull_request_id::Integer)
@endpoint list_pull_request_comments(project::Integer, pull_request_id::Integer)

"""
    get_pull_request_comment(
        ::Forge, owner::$AStr, repo::$AStr, comment_id::Integer
    )
    get_pull_request_comment(
        ::Forge, project::Integer, pull_request_id::Integer, comment_id::Integer
    )

Get a single comment for a given pull request and comment ID.
"""
@endpoint get_pull_request_comment(owner::AStr, repo::AStr, comment_id::Integer)
@endpoint get_pull_request_comment(
    project::Integer, pull_request_id::Integer, comment_id::Integer
)

"""
    create_pull_request_comment(
        ::Forge, owner::$AStr, repo::$AStr, pull_request_id::Integer;
        kwargs...
    )
    create_pull_request_comment(
        ::Forge, project::Integer, pull_request_id::Integer;
        kwargs...
    )

Create a new comment for a given pull request.
"""
@endpoint create_pull_request_comment(owner::AStr, repo::AStr, pull_request_id::Integer)
@endpoint create_pull_request_comment(project::Integer, pull_request_id::Integer)

"""
    update_pull_request_comment(
        ::Forge, owner::$AStr, repo::$AStr, comment_id::Integer;
        kwargs...
    )
    update_pull_request_comment(
        ::Forge, project::Integer, pull_request_id::Integer, comment_id::Integer;
        kwargs...
    )

Update the body for a given pull request comment.
"""
@endpoint update_pull_request_comment(owner::AStr, repo::AStr, comment_id::Integer)
@endpoint update_pull_request_comment(
    project::Integer, pull_request_id::Integer, comment_id::Integer
)

"""
    delete_pull_request_comment(
        ::Forge, owner::$AStr, repo::$AStr, comment_id::Integer
    )
    delete_pull_request_comment(
        ::Forge, project::Integer, pull_request_id::Integer, comment_id::Integer
    )

Delete a given pull request comment.
"""
@endpoint delete_pull_request_comment(owner::AStr, repo::AStr, comment_id::Integer)
@endpoint delete_pull_request_comment(
    project::Integer, pull_request_id::Integer, comment_id::Integer
)

"""
    list_pipeline_schedules(::Forge, project::Integer; kwargs...)
    list_pipeline_schedules(::Forge, owner::$AStr, repo::$AStr; kwargs...)

List all existing pipeline schedules for a project.
"""
@endpoint list_pipeline_schedules(project::Integer)
@endpoint list_pipeline_schedules(owner::AStr, repo::AStr)


text:
using StatsBase, DataFrames, Flux
include("2048_fn.jl")

using Flux
model = Flux.Chain(
  Dense(14*16+4, 128, relu),
  Dense(128, 64, relu),
  Dense(64, 32, relu),
  Dense(32, 1, relu));

loss(x,y) = sum((model(x) .- y).^2)
opt = ADAM(params(model))

ϵ = 1
@time for i=1:1000
  global ϵ
  ϵ -= 0.0001
  simulate_one_game_and_train!(model, opt, ϵ)
  if i % 100 == 0
    println("epsilon: $ϵ")
    init_grid = init_game()
    res = simulate_game!(init_grid, model, 0)
    display(res[2])
    assess_model(model)
  end
end

using JLD2
@save "no_neg_model.jld2" model

# play one game to see how it goes
function assess_model(model)
  init_grid = init_game()
  @time (init_grid, fnl, seq, cart, tf, reward) = simulate_game!(init_grid, model, 0.5)

  # obtain the discounted rewards
  y = discounted_rewards(reward, 0.99)[1:end-1];

  # score every move using  model
  updated_model = [Tracker.data(model(grid_move_to_feature(sim_seq(init_grid, seq[1:i+1], cart[1:i+1], tf[1:i+1]), seq[i])))[1] for i=1:length(seq)-1]
  #plot(cumsum(reward))
  plot(updated_model, label="model")
  plot!(y, label="discounted rewards")
end

assess_model(model)


text:
include("Constraint.jl")
include("Variables.jl")
include("Solution.jl")

include("Change.jl")

include("Model.jl")

include("Preprocessing.jl")

mutable struct Compteur
	ind::Int
end

function Arc_Consistency(model::Model, setOfConstraint::Union{Set{Constraint}, Nothing} = nothing; verbose::Bool = false)

	stop = false

	if setOfConstraint == nothing
		setOfConstraint = copy(model.constraints)
	end
	dictChanges = Dict{Variable, Change}(var => Change(var) for var in model.varsInter)

	for g in 1:model.g
		for w in 1:model.w
			dictChanges[model.X[g, w]] = Change(model.X[g, w])
		end
	end

	while !isempty(setOfConstraint) && !stop

		verbose && println("N° de contraintes : $(length(setOfConstraint))")

		constr = pop!(setOfConstraint)

		modifiedVariables, changes, stop = filtrage!(constr)

		for change in changes
			union!(dictChanges[change.var].added, change.added)
			union!(dictChanges[change.var].removed, change.removed)

			dictChanges[change.var].cardAdded += change.cardAdded
			dictChanges[change.var].cardRemoved += change.cardRemoved

			dictChanges[change.var].fixed |= change.fixed

			stop |= change.cardAdded < 0
			stop |= change.cardRemoved > 0
		end

		for var in modifiedVariables
			union!(setOfConstraint, var.linkedConstraint)
		end

	end
	#=
	if modelNotFixed(model)
		println("WIP : Branch-And-Bound")
	elseif isEmpty(model)
		return solution
	else
		solution = createSolution(model)
	end

	return solution
	=#

	vars = append!([model.X[g, w] for g in 1:model.g for w in 1:model.w], model.varsInter)

	for var in vars
		if isFixed(var) && !var.isFixed
			fix!(var, dictChanges[var])
		end
	end

	filter!(var-> !(var.isFixed), model.varsNotFixed)

	# !stop && println("Yo : ", map(x->x.cardAdded, values(dictChanges)))

	return dictChanges, stop
end

#TODO
#
# Je doids gérer le cas d'erreur -> SI on tombre sur une solution non faisable ca fait quoi
# Faire les fonctions unforced, returnParent, Solution(model)

function branch(model = ModelTest();setOfConstraint::Union{Nothing, Set{Constraint}} = nothing, nbAppel::Int = 1, verbose::Bool = false, debug::Real = Inf, compteur::Compteur = Compteur(0))

	# println("N° : ", compteur.ind)
	filter!(var-> !(var.isFixed), model.varsNotFixed)

	changes = Dict{Variable, Change}()
	changes, stop = Arc_Consistency(model)

	if stop
		returnParent!(changes, model)
		return nothing, false
	else
	@assert compteur.ind <= debug "Mwahahaahahahahhahahahaha"

		# println("Profondeur : $nbAppel")

		mini = Inf
		varToTest = nothing
		indVar = 1

		nbVarToTest = length(model.varsNotFixed)

		if nbVarToTest > 0

			while indVar <= nbVarToTest && mini > 1
				var = model.varsNotFixed[indVar]
				if length(var.upperBound) - length(var.lowerBound) < mini
					mini = length(var.upperBound) - length(var.lowerBound)
					varToTest = var
				end
				indVar += 1
			end

			valueToTest = union(setdiff(varToTest.upperBound, varToTest.lowerBound), [nothing])
			nbValueToTest = length(valueToTest)
			sol = nothing
			stop = false
			indValue = 1

			# println("_______________")
			# println("Variable : ", varToTest)
			# println("Valeurs : ", valueToTest)
			# println("_______________")

			while !stop && indValue <= nbValueToTest
				value = valueToTest[indValue]

				if value != nothing || varToTest.cardinalInf == length(varToTest.lowerBound)
					if value == nothing || value in varToTest.upperBound
						#println()
						#println("On force $varToTest à $value")

						changeForce = forced!(varToTest, value, model)
						# println("$varToTest : $(varToTest.lowerBound)")

						# println("	"^(nbAppel-1), "lowerBound : ", varToTest.lowerBound)
						# println("	"^(nbAppel-1), "upperBound : ", varToTest.upperBound)

						compteur.ind += 1
						sol, stop = branch(model, setOfConstraint = copy(varToTest.linkedConstraint), nbAppel = nbAppel + 1, compteur = compteur, debug = debug)
						#println()
						#println("On déforce $varToTest pour $value")

						# println("	"^(nbAppel-1), "Après")
						# println("	"^(nbAppel-1), "lowerBound : ", varToTest.lowerBound)
						# println("	"^(nbAppel-1), "upperBound : ", varToTest.upperBound)
						# println()
						# println("	"^(nbAppel-1), sol, " - ", stop)

						unforced!(varToTest, changeForce, model)

					end
				end
				indValue += 1
			end
			returnParent!(changes, model)
			return sol, stop
		else
			return Solution(model), true
		end
	end

end

function modelFixed()
	model = main()
	model.X[1, 1].lowerBound = [1, 2, 3]
	model.X[2, 1].lowerBound = [4, 5, 6]
	model.X[3, 1].lowerBound = [7, 8, 9]
	model.X[4, 1].lowerBound = [10, 11, 12]

	model.X[1, 1].upperBound = [1, 2, 3]
	model.X[2, 1].upperBound = [4, 5, 6]
	model.X[3, 1].upperBound = [7, 8, 9]
	model.X[4, 1].upperBound = [10, 11, 12]

	for x in model.X[:, 1]
		x.isFixed = true
	end
	model.X[4, 1].isFixed = false

	model.varsInter[1].lowerBound = [1, 2, 3, 4, 5, 6]
	model.varsInter[2].lowerBound = [1, 2, 3, 4, 5, 6, 7, 8, 9]
	model.varsInter[3].lowerBound = [1, 2, 3, 4, 5, 6 ,7, 8, 9, 10, 11 ,12]

	model.varsInter[1].upperBound = [1, 2, 3, 4, 5, 6]
	model.varsInter[2].upperBound = [1, 2, 3, 4, 5, 6, 7, 8, 9]
	model.varsInter[3].upperBound = [1, 2, 3, 4, 5, 6 ,7, 8, 9, 10, 11 ,12]

	for f in model.varsInter
		f.isFixed = true
	end

	return model
end

function debug(ind::Int)
	model = ModelTest(6, 3, 2)

	try
		test = branch(model, debug = ind)
	catch y
		println(y)
	end
	return model
end

function main(p = 4, g = 2, w = 2)
	model = ModelTest(p, g, w)

	preprocessing(model)
	sol, stop = branch(model)

	return sol, stop
end


text:
## LICENSE
##   Copyright (c) 2018 GEOPHYSTECH LLC
##
##   Licensed under the Apache License, Version 2.0 (the "License");
##   you may not use this file except in compliance with the License.
##   You may obtain a copy of the License at
##
##       http://www.apache.org/licenses/LICENSE-2.0
##
##   Unless required by applicable law or agreed to in writing, software
##   distributed under the License is distributed on an "AS IS" BASIS,
##   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
##   See the License for the specific language governing permissions and
##   limitations under the License.

## initial release by Andrey Stepnov, email: a.stepnov@geophsytech.ru

## VS30 data point
"""
Mutable type for fill vs30 grid (Array).

  `Point_vs30(lon,lat,vs30)`

Latitude and longitude assumes degrees for WGS84 ellipsoid. `vs30` in meters per second. All fields are ::Float64.
"""
mutable struct Point_vs30
  lon::Float64
  lat::Float64
  vs30::Float64
end
"""
Mutable type for fill vs30 grid (Array) with `Dl` data. `Dl` is the top depth to the layer whose S-wave velocity is `l` (in `[m/s]`) at the site. See Morikawa Fujiwara 2013 for further reading.

  `Point_vs30(lon,lat,vs30,dl)`

Latitude and longitude assumes degrees for WGS84 ellipsoid. `vs30` in meters per second. `dl` in kilometers. All fields are ::Float64.
"""
mutable struct Point_vs30_dl
  lon::Float64
  lat::Float64
  vs30::Float64
  dl::Float64
end
## Output PGA data point
"""
Mutable type for output PGA data from GMPE modeling funtions

  Fields:
```
  lon   :: Float64
  lat   :: Float64
  pga   :: Float64
```
Latitude and longitude assumes degrees for WGS84 ellipsoid. `pga` is Acceleration of gravity in percent (%g) rounded to ggg.gg.
"""
mutable struct Point_pga_out
  lon::Float64
  lat::Float64
  pga::Float64 #Acceleration of gravity  in percent (%g) rounded to ggg.gg
end
## Output PGV data point
"""
Mutable type for output PGV data from GMPE modeling funtions

  Fields:
```
  lon   :: Float64
  lat   :: Float64
  pgv   :: Float64
```
Latitude and longitude assumes degrees for WGS84 ellipsoid. `pgv` is [cm/s].
"""
mutable struct Point_pgv_out
  lon::Float64
  lat::Float64
  pgv::Float64 #Acceleration of gravity in percent rounded to ggg.gg
end
## Output PSA data point
"""
Mutable type for output PSA data

  Fields:
```
  lon   :: Float64
  lat   :: Float64
  psa   :: Float64
```
Latitude and longitude assumes degrees for WGS84 ellipsoid. `psa` is damped pseudo-spectral acceleration (%g) in percent rounded to ggg.gg.
"""
mutable struct Point_psa_out
  lon::Float64
  lat::Float64
  psa::Float64 # damped pseudo-spectral acceleration of gravity in percent rounded to ggg.gg (%g)
end
## earthquake location data
"""
Mutable type for earthquake location data.

  Earthquake(lat,lon,depth,local_mag,moment_mag)

Latitude and longitude assumes degrees for WGS84 ellipsoid. Depth in km.
Mw=0 in case of moment magnitude is not specified.
All fields are ::Float64.
"""
mutable struct Earthquake
  lon::Float64
  lat::Float64
  depth::Float64
  local_mag::Float64
  moment_mag::Float64
end
Earthquake(x,y,z,k) = Earthquake(x,y,z,k,0) # Mw usually not ready right after earthquake
## AS2008 GMPE parameters
mutable struct Params_as2008
  a1::Float64
  a2::Float64
  a3::Float64
  a4::Float64
  a5::Float64
  a8::Float64
  a10::Float64
  a18::Float64
  b::Float64
  c::Float64
  c1::Float64
  c4::Float64
  n::Float64
  vlin::Float64
  v1::Float64
  vs30_1100::Float64
  ground_motion_type::String
end
## Si-Midorikawa 1999 GMPE parameters
mutable struct Params_simidorikawa1999
  a::Float64
  h::Float64
  d1::Float64
  d2::Float64
  d3::Float64
  e_::Float64
  k::Float64
  S1::Float64
  S2::Float64
  S3::Float64
  ground_motion_type::String
end
## Morikawa-Fujiwara 2013 GMPE parameters
mutable struct Params_mf2013
  Mw0::Float64
  a::Float64
  b::Float64
  c::Float64
  d::Float64
  e::Float64
  sigma::Float64
  pd::Float64
  Dlmin::Float64
  D0::Float64
  ps::Float64
  Vsmax::Float64
  V0::Float64
  gamma::Float64
  ASID::Bool # Anomalous Seismic Intensity Distribution
  ground_motion_type::String
end


text:
using QuantumOpticsBase
using QuantumOpticsBase: check_samebases, check_multiplicable

import OrdinaryDiffEq, DiffEqCallbacks

function recast! end

"""
    integrate(tspan, df::Function, x0::Vector{ComplexF64},
            state::T, dstate::T, fout::Function; kwargs...)

Integrate using OrdinaryDiffEq
"""
function integrate(tspan, df, x0,
            state, dstate, fout;
            alg = OrdinaryDiffEq.DP5(),
            steady_state = false, tol = 1e-3, save_everystep = false, saveat=tspan,
            callback = nothing, kwargs...) where {T,X}

    function df_(dx, x, p, t)
        recast!(state,x)
        recast!(dstate,dx)
        df(t, state, dstate)
        recast!(dx,dstate)
    end
    function fout_(x, t, integrator)
        recast!(state,x)
        fout(t, state)
    end

    tType = float(eltype(tspan))
    out_type = pure_inference(fout, Tuple{tType,typeof(state)})

    out = DiffEqCallbacks.SavedValues(tType,out_type)

    scb = DiffEqCallbacks.SavingCallback(fout_,out,saveat=saveat,
                                         save_everystep=save_everystep,
                                         save_start = false)

    prob = OrdinaryDiffEq.ODEProblem{true}(df_, x0,(convert(tType, tspan[1]),convert(tType, tspan[end])))

    if steady_state
        affect! = function (integrator)
            !save_everystep && scb.affect!(integrator,true)
            OrdinaryDiffEq.terminate!(integrator)
        end
        _cb = OrdinaryDiffEq.DiscreteCallback(
                                SteadyStateCondtion(copy(state),tol,state),
                                affect!;
                                save_positions = (false,false))
        cb = OrdinaryDiffEq.CallbackSet(_cb,scb)
    else
        cb = scb
    end

    full_cb = OrdinaryDiffEq.CallbackSet(callback,cb)

    sol = OrdinaryDiffEq.solve(
                prob,
                alg;
                reltol = 1.0e-6,
                abstol = 1.0e-8,
                save_everystep = false, save_start = false,
                save_end = false,
                callback=full_cb, kwargs...)
    out.t,out.saveval
end

function integrate(tspan, df, x0,
            state, dstate, ::Nothing; kwargs...)
    function fout(t, state)
        copy(state)
    end
    integrate(tspan, df, x0, state, dstate, fout; kwargs...)
end

function check_psi0(psi0::Ket{B, T} where B<:Basis where T<:Array)
    if !isa(psi0, Ket{B, Array{T, 1}} where B<:Basis where T<:Complex)
        error("Method only works for complex-valued initual wave functions.")
    end
end

function check_psi0(psi0)
end

struct SteadyStateCondtion{T,T2,T3}
    rho0::T
    tol::T2
    state::T3
end
function (c::SteadyStateCondtion)(rho,t,integrator)
    timeevolution.recast!(rho,c.state)
    dt = integrator.dt
    drho = tracedistance(c.rho0, c.state)
    c.rho0.data[:] = c.state.data
    drho/dt < c.tol
end

const QO_CHECKS = Ref(true)
"""
    @skiptimechecks

Macro to skip checks during time-dependent problems.
Useful for `timeevolution.master_dynamic` and similar functions.
"""
macro skiptimechecks(ex)
    return quote
        QO_CHECKS.x = false
        local val = $(esc(ex))
        QO_CHECKS.x = true
        val
    end
end

Base.@pure pure_inference(fout,T) = Core.Compiler.return_type(fout, T)


text:
using DiffEqFlux, DiffEqSensitivity, Flux, OrdinaryDiffEq, Zygote, Optim, NLopt, BlackBoxOptim, Test #using Plots
using MultistartOptimization, NLopt

function lotka_volterra(du,u,p,t)
  x, y = u
  α, β, δ, γ = p
  du[1] = dx = (α - β*y)x
  du[2] = dy = (δ*x - γ)y
end
p = [2.2, 1.0, 2.0, 0.4]
u0 = [1.0,1.0]
prob = ODEProblem(lotka_volterra,u0,(0.0,10.0),p)

# Reverse-mode

function predict_rd(p)
  Array(solve(prob,Tsit5(),p=p,saveat=0.1,reltol=1e-4))
end
function loss_rd(p)
  sum(abs2,x-1 for x in predict_rd(p))
end

loss_rd(p)

grads = Zygote.gradient(loss_rd, p)
@test !iszero(grads[1])

cb = function (p,l)
  display(l)
  false
end

# Display the ODE with the current parameter values.
loss1 = loss_rd(p)
pmin = DiffEqFlux.sciml_train(loss_rd, p, ADAM(0.1), cb = cb, maxiters = 100)
loss2 = loss_rd(pmin.minimizer)
@test 10loss2 < loss1

pmin = DiffEqFlux.sciml_train(loss_rd, p, BFGS(initial_stepnorm = 0.01), cb = cb)
loss2 = loss_rd(pmin.minimizer)
@test 10loss2 < loss1
pmin = DiffEqFlux.sciml_train(loss_rd, p, Fminbox(BFGS(initial_stepnorm = 0.01)), lower_bounds = [0.0 for i in 1:4], upper_bounds = [5.0 for i in 1:4], cb = cb)
loss2 = loss_rd(pmin.minimizer)
@test 10loss2 < loss1

pmin = DiffEqFlux.sciml_train(loss_rd, p, maxiters=100, lower_bounds = [0.0 for i in 1:4], upper_bounds = [5.0 for i in 1:4], cb = cb)
loss2 = loss_rd(pmin.minimizer)
@test 10loss2 < loss1

pmin = DiffEqFlux.sciml_train(loss_rd, p, TikTak(100),
                              local_method = NLopt.LN_BOBYQA,
                              maxiters=100,
                              lower_bounds = [0.0 for i in 1:4], upper_bounds = [5.0 for i in 1:4], cb = cb)
loss2 = loss_rd(pmin.minimizer)
@test 10loss2 < loss1

#=
using QuadDIRECT
pmin = DiffEqFlux.sciml_train(loss_rd, p, QuadDIRECT(),
                              splits = (),
                              maxiters=100,
                              lower_bounds = [0.0 for i in 1:4], upper_bounds = [5.0 for i in 1:4], cb = cb)
loss2 = loss_rd(pmin.minimizer)
@test 10loss2 < loss1
=#

# Forward-mode, R^n -> R^m layer

p = [2.2, 1.0, 2.0, 0.4]
function predict_fd(p)
  vec(Array(solve(prob,Tsit5(),p=p,saveat=0.0:0.1:1.0,reltol=1e-4,sensealg=ForwardDiffSensitivity())))
end
loss_fd(p) = sum(abs2,x-1 for x in predict_fd(p))
loss_fd(p)

grads = Zygote.gradient(loss_fd, p)
@test !iszero(grads[1])

opt = ADAM(0.1)
cb = function (p,l)
  display(l)
  #display(plot(solve(remake(prob,p=p),Tsit5(),saveat=0.1),ylim=(0,6)))
  false
end

# Display the ODE with the current parameter values.
loss1 = loss_fd(p)
pmin = DiffEqFlux.sciml_train(loss_fd, p, opt, cb = cb, maxiters = 100)
loss2 = loss_fd(pmin.minimizer)
@test 10loss2 < loss1

pmin = DiffEqFlux.sciml_train(loss_fd, p, BFGS(initial_stepnorm = 0.01), cb = cb)
loss2 = loss_fd(pmin.minimizer)
@test 10loss2 < loss1

pmin = DiffEqFlux.sciml_train(loss_fd, p, TikTak(100),
                              local_method = NLopt.LN_BOBYQA,
                              maxiters=100,
                              lower_bounds = [0.0 for i in 1:4], upper_bounds = [5.0 for i in 1:4], cb = cb)
loss2 = loss_fd(pmin.minimizer)
@test 10loss2 < loss1

#=
pmin = DiffEqFlux.sciml_train(loss_fd, p, QuadDIRECT(),
                              splits = (),
                              maxiters=100,
                              lower_bounds = [0.0 for i in 1:4], upper_bounds = [5.0 for i in 1:4], cb = cb)
loss2 = loss_fd(pmin.minimizer)
@test 10loss2 < loss1
=#

# Adjoint sensitivity
p = [2.2, 1.0, 2.0, 0.4]
ps = Flux.params(p)
function predict_adjoint(p)
    vec(Array(solve(prob,Tsit5(),u0=eltype(p).(prob.u0),p=p,saveat=0.1,reltol=1e-4)))
end
loss_reduction(sol) = sum(abs2,x-1 for x in vec(sol))
loss_adjoint(p) = loss_reduction(predict_adjoint(p))
loss_adjoint(p)

grads = Zygote.gradient(loss_adjoint, p)
@test !iszero(grads[1])

opt = ADAM(0.1)
cb = function (p,l)
  display(l)
  #display(plot(solve(remake(prob,p=p),Tsit5(),saveat=0.1),ylim=(0,6)))
  false
end

# Display the ODE with the current parameter values.
loss1 = loss_adjoint(p)
pmin = DiffEqFlux.sciml_train(loss_adjoint, p, opt, cb = cb, maxiters = 100)
loss2 = loss_adjoint(pmin.minimizer)
@test 10loss2 < loss1

pmin = DiffEqFlux.sciml_train(loss_adjoint, p, BFGS(initial_stepnorm = 0.01), cb = cb)
loss2 = loss_adjoint(pmin.minimizer)
@test 10loss2 < loss1

pmin = DiffEqFlux.sciml_train(loss_adjoint, p, Fminbox(BFGS(initial_stepnorm = 0.01)), lower_bounds = [0.0 for i in 1:4], upper_bounds = [5.0 for i in 1:4], cb = cb)
loss2 = loss_adjoint(pmin.minimizer)
@test 10loss2 < loss1

opt = Opt(:LD_MMA, 4)
pmin = DiffEqFlux.sciml_train(loss_adjoint, p, opt)
loss2 = loss_adjoint(pmin.minimizer)
@test 10loss2 < loss1

pmin = DiffEqFlux.sciml_train(loss_adjoint, p, TikTak(100),
                              local_method = NLopt.LN_BOBYQA,
                              maxiters=100,
                              lower_bounds = [0.0 for i in 1:4], upper_bounds = [5.0 for i in 1:4], cb = cb)
loss2 = loss_adjoint(pmin.minimizer)
@test 10loss2 < loss1

#=
pmin = DiffEqFlux.sciml_train(loss_adjoint, p, QuadDIRECT(),
                              splits = (),
                              maxiters=100,
                              lower_bounds = [0.0 for i in 1:4], upper_bounds = [5.0 for i in 1:4], cb = cb)
loss2 = loss_adjoint(pmin.minimizer)
@test 10loss2 < loss1
=#

function lotka_volterra2(u,p,t)
  x, y = u
  α, β, δ, γ = p
  dx = (α - β*y)x
  dy = (δ*x - γ)y
  [dx,dy]
end
p = [2.2, 1.0, 2.0, 0.4]
u0 = [1.0,1.0]
prob = ODEProblem{false}(lotka_volterra2,u0,(0.0,10.0),p)
function predict_adjoint(p)
    vec(Array(solve(prob,Tsit5(),p=p,saveat=0.1,reltol=1e-4)))
end
loss_reduction(sol) = sum(abs2,x-1 for x in vec(sol))
loss_adjoint(p) = loss_reduction(predict_adjoint(p))

pmin = DiffEqFlux.sciml_train(loss_adjoint, p, Newton())
loss2 = loss_adjoint(pmin.minimizer)
@test 10loss2 < loss1

pmin = DiffEqFlux.sciml_train(loss_adjoint, p, NewtonTrustRegion())
loss2 = loss_adjoint(pmin.minimizer)
@test 10loss2 < loss1

pmin = DiffEqFlux.sciml_train(loss_adjoint, p, Optim.KrylovTrustRegion(), maxiters = 100)
loss2 = loss_adjoint(pmin.minimizer)
@test 10loss2 < loss1


text:
module TestPromise

using ConcurrentUtils
using Random: randperm
using Test

function test_serial()
    p = Promise()
    @test try_race_fetch(p) == Err(NotSetError())
    put!(p, 1)
    @test fetch(p) == 1
    @test wait(p) === nothing
    @test try_race_fetch(p) == Ok{Any}(1)
    @test try_race_put_with!(error, p) == Err{Any}(1)
    @test race_put_with!(error, p) == 1
    @test try_race_put!(p, 2) == Err{Any}(1)
    @test_throws OccupiedError{Any}(1) put!(p, 2)
end

function test_serial_try_race_fetch_or()
    p = Promise()
    @test try_race_put_with!(() -> 1, p) == Ok{Any}(1)
    @test try_race_put_with!(error, p) == Err{Any}(1)
    @test fetch(p) == 1
end

function test_serial_race_fetch_or()
    p = Promise()
    @test race_put_with!(() -> 1, p) == 1
    @test race_put_with!(error, p) == 1
    @test fetch(p) == 1
end

function test_serial_try_race_put()
    p = Promise()
    @test try_race_put!(p, 1) == Ok{Any}(1)
    @test try_race_put!(p, 2) == Err{Any}(1)
end

function check_concurrent_put_fetch(ntasks, ntries)
    promises = [[Promise{Int}() for _ in 1:ntasks] for _ in 1:ntries]
    shufflers = [randperm(ntasks) for _ in 1:ntries]
    outputs = [[Threads.Atomic{Int}(0) for _ in 1:ntasks] for _ in 1:ntries]
    @sync begin
        for itask in 1:ntasks
            Threads.@spawn for itry in 1:ntries
                put!(promises[itry][itask], itask)
                jtask = shufflers[itry][itask]
                x = fetch(promises[itry][jtask])
                Threads.atomic_add!(outputs[itry][jtask], x)
            end
        end
    end
    return [outputs[itry][itask][] for itask in 1:ntasks, itry in 1:ntries]
end

function test_concurrent_put_fetch(; ntries = 128)
    @testset for ntasks in [Threads.nthreads(), 64 * Threads.nthreads()]
        actual = check_concurrent_put_fetch(ntasks, ntries)
        desired = [itask for itask in 1:ntasks, _ in 1:ntries]
        @test actual == desired
    end
end

end  # module


text:
# Tree data type used for a routing tree.
type Tree
    value::Any
    children::Array{Tree,1}
end
Tree(value::Any) = Tree(value, Tree[])

function ischild(element, tree::Tree)
    for child in tree.children
        if isequal(child.value, element)
            return true
        end
    end
    return false
end

function getchild(tree::Tree, element)
    for child in tree.children
        if isequal(child.value, element)
            return child
        end
    end
end

function insert!(tree::Tree, values::Array)
    if isempty(values)
        tree
    elseif ischild(values[1], tree)
        t = getchild(tree, values[1])
        insert!(t, values[2:end])
    else
        t = Tree(values[1])
        insert!(t, values[2:end])
        push!(tree.children, t)
    end
end

abstract SearchStrategy

type DFS <: SearchStrategy end

type PRUNE end

function search(tree::Tree, pred::Function)
    action = pred(tree.value)
    if action == true
        tree.value
    elseif action == PRUNE
        nothing
    else
        for child in tree.children
            val = search(child, pred)
            if val != nothing
                return val
            end
        end
        nothing
    end
end


text:
include("testfunctions.jl")
using DiferentialEvolution
using PlotlyJS

#Plots
len=500
x = LinRange(-600,600,len)
y = LinRange(-600,600,len)
X = ones(len)*x'
Y = y*ones(len)'
Z = Eggholder.(X,Y)
Eggholder(-2.5153506020027042e64, 1.79281509087981e64)
@time critic1 = DE_JADE(Eggholder, 500, 25, [400,400], [600,600], 0.2, 1/20)
@time critic2 = DE_dither(Eggholder, 1000, 500, [-512,-512], [512,512],0.9,0.5)
p_init = contour(x =x, y= y,z = Z)
crit = scatter(x=[critic[2][1]],y=[critic[2][2]], mode="markers",color=2)
plot([p_init,crit])


text:
import Base.transpose

export translate, transpose, randomnotes, subdivision
export velocities, positions, pitches, durations

velocities(notes::Notes) = [Int(x.velocity) for x in notes]
positions(notes::Notes) = [Int(x.position) for x in notes]
pitches(notes::Notes) = [Int(x.pitch) for x in notes]
durations(notes::Notes) = [Int(x.duration) for x in notes]

"""
    randomnotes(n::Int, tpq = 960)
Generate some random notes that start sequentially.
"""
function randomnotes(n::Int, tpq = 960, randchannel = false)
    notes = Note[]
    prevpos = 0
    durran = 1:tpq
    posran = 0:4*tpq
    for i in 1:n
        if randchannel
            note = Note(rand(UInt8), rand(0:0x7F), prevpos + rand(posran), rand(durran), rand(0:127))
        else
            note = Note(rand(UInt8), rand(0:127), prevpos + rand(posran), rand(durran))
        end
        push!(notes, note)
        prevpos = note.position
    end
    return Notes(notes, tpq)
end

"""
    translate(notes, ticks)
Translate the `notes` for the given amount of `ticks`.
"""
translate(notes::Notes, ticks) = Notes(translate(notes.notes, ticks), notes.tpq)
translate(notes::Vector{N}, ticks) where {N<:AbstractNote} =
[Note(n.pitch, n.velocity, n.position + ticks, n.duration, n.channel) for n in notes]

"""
    transpose(notes, semitones)
Transpose the `notes` for the given amount of `semitones`.
"""
transpose(notes::Notes, semitones) = Notes(transpose(notes.notes, semitones), notes.tpq)
transpose(notes::Vector{N}, semitones) where {N<:AbstractNote} =
[Note(n.pitch + semitones, n.velocity, n.position, n.duration, n.channel) for n in notes]

"""
    subdivision(n::Int, tpq)
Return how many ticks is the duration of
the subdivision of a 4/4-bar into `n` equal parts, assuming the ticks per quarter
note are `tpq`.

For example, for sixteenth notes you would do `subdivision(16, tpq)`, for
eigth-note triplets `subdivision(12, tpq)` and so on.
"""
subdivision(n::Int, tpq)::Int = (4*tpq)/n

function timesort(notes::Notes)
    issorted(notes, by = x -> x.position) && return notes
    return Notes(sort(notes.notes, by = x -> x.position), notes.tpq)
end


text:
using Random, Test
using ModelStatisticsLH

## --------------   Setup

function split_sname(statName :: Symbol)
    return split(string(statName), "_");
end

function ModelStatisticsLH.groups_from_name(m :: StatsCollection{:g},
    statName :: Symbol)

    v = split_sname(statName);
    return Symbol(v[2])
end

function ModelStatisticsLH.statname_from_name(m :: StatsCollection{:g},
    statName :: Symbol)

    v = split_sname(statName);
    return Symbol(v[1])
end

function scollection_test()
    rng = MersenneTwister(434);

    @testset "StatsCollection" begin
        T = :g;
        sc = StatsCollection{T}();
        @test groups_from_name(sc, :workTime_xy) == :xy
        @test statname_from_name(sc, :workTime_g) == :workTime

        TMS = :tms;
        ms = make_test_model_stats(TMS, 5, rng);
        add_mstats!(sc, ms);
        @test has_mstats(sc, TMS)
        ms2 = get_mstats(sc, TMS);
        @test sort(collect(var_names(ms))) == sort(collect(var_names(ms2)))

        TMS3 = (:tms,:three);
        ms3 = make_test_model_stats(TMS3, 4, rng);
        add_mstats!(sc, ms3);
        @test has_mstats(sc, TMS3);
        ms4 = make_test_model_stats(TMS3, 3, rng);
        replace_mstats!(sc, ms4);
        ms4a = get_mstats(sc, TMS3);
        @test sort(collect(var_names(ms4))) == sort(collect(var_names(ms4a)))

        grpV = get_groups(sc);
        @test length(grpV) == 2

        delete_mstats!(sc, TMS);
        @test !has_mstats(sc, TMS)
	end
end

function get_stats_test()
    rng = MersenneTwister(43);
    @testset "Access stats directly" begin
        T = :g;
        sc = make_test_scollection(T, rng);
        grpV = collect(get_groups(sc));
        grp = grpV[2];

        @test get_stats(sc, grp, :notThere; defaultValue = 1) == 1;
        ms = get_mstats(sc, grp);
        for varName in var_names(ms)
            @test has_variable(sc, grp, varName);
            x = get_stats(sc, grp, varName);
            @test isequal(x, getproperty(ms, varName))

            # Test syntactic sugar
            statName = Symbol("$(varName)_$(grp)");
            @test isequal(x, get_stats(sc, statName))

            if isa(x, AbstractVector)
                idx = length(x);
            elseif isa(x, AbstractMatrix)
                idx = [1, size(x, 2)];
            else
                idx = nothing;
            end
            if !isnothing(idx)
                x2 = getindex(sc, statName, idx...);
                @test isequal(x[idx...], x2)
            end
        end
    end
end

function dot_notation_test()
    rng = MersenneTwister(434);
    @testset "Dot notation" begin
        varName = :x1;
        varMeta = VarInfo(varName, "x1", Matrix{Float64}; size = (4,3));
        xValues = rand(rng, size(varMeta)...);
        @test size(xValues) == (4,3)
        mMeta = Dict{Symbol, Any}([:test => true]);
        grp = :xy;
        statName = :x1_xy;
        ms = ModelStats{grp}(mMeta);
        add_variable!(ms, varMeta, xValues);

        T = :g;
        sc = StatsCollection{T}();
        add_mstats!(sc, ms);

        x = get_stats(sc, statName);
        @test x == xValues
        x2 = sc.x1_xy;
        @test x == x2

        # setstats!
        set_stats!(sc, :x1_xy, x .+ 1);
        @test all(sc.x1_xy .== x .+ 1);
        sc.x1_xy = xValues;
        @test all(sc.x1_xy .== xValues);

        # getindex
        sc.x1_xy .= xValues;
        @test sc.x1_xy[2,3] == xValues[2,3]
        @test sc.x1_xy[3] == xValues[3]
        @test sc.x1_xy[:,2] == xValues[:,2]

        # setstats!
        newValues = xValues .+ 0.5;
        set_stats!(sc, grp, varName, newValues);
        x3 = get_stats(sc, statName);
        @test isapprox(newValues, x3)

        # setindex!
        sc.x1_xy .= xValues;
        sc.x1_xy[3,2] = 17.0;
        @test sc.x1_xy[3,2] == 17.0
        @test sc.x1_xy[2,2] == xValues[2,2];
        sc.x1_xy[:,2] .= 17.0;
        @test all(sc.x1_xy[:,2] .== 17.0)
        @test sc.x1_xy[:,1] == xValues[:,1];
    end
end

@testset "StatsCollection" begin
    scollection_test();
    get_stats_test();
    dot_notation_test();
end

# ------------


text:
#===============================================================================

    SFCModelDef.jl

    struct containing SFC model definition

    Author: Tom Price
    Date:   June 2019

===============================================================================#

import DataStructures.OrderedDict

mutable struct SFCModelDef ## non-linear model, one feature (`x`)
    ## included in SFC_MODEL_BASE
    name            ::SFCModelName
    linear          ::Bool
    _x_strs         ::Vector{String}
    X_strs          ::Vector{String}
    coef_strs       ::Vector{String}
    coef_cnstrnts   ::AbstractVector ## assume all linear
    func_init_coefs ::Function
    pred_strs       ::OrderedDict{Symbol,String}
    ## added by `add*!`` functions
    func_pred_strs  ::OrderedDict{Symbol,String}
    funcs_pred      ::OrderedDict{Symbol,Function}
    func_fit_str    ::String
    func_fit        ::Function
end

max_stars_repo_name: MakerButt/chaipcr

max_issues_repo_name: MakerButt/chaipcr

max_forks_repo_name: MakerButt/chaipcr


text:
#==========================================================================================#

# TYPE

mutable struct FrölichMelly <: TwoStageModel

    first_stage::OneStageModel
    second_stage::OLS
    pscore::Vector{Float64}
    eweights::ProbabilityWeights{Float64, Float64, Vector{Float64}}

    FrölichMelly() = new()
end

#==========================================================================================#

# FIRST STAGE

function first_stage(::Type{FrölichMelly}, MM::Type{<:OneStageModel}, MD::Microdata; kwargs...)

    FSD                    = Microdata(MD)
    FSD.mapping[:response] = MD.mapping[:instrument]

    pop!(FSD.mapping, :treatment)
    pop!(FSD.mapping, :instrument)

    fit(MM, FSD; kwargs...)
end

#==========================================================================================#

# ESTIMATION

function fit(
        ::Type{FrölichMelly},
        MM::Type{<:OneStageModel},
        MD::Microdata;
        novar::Bool = false,
        kwargs...
    )

    m = first_stage(FrölichMelly, MM, MD; novar = novar)
    fit(FrölichMelly, m, MD; novar = novar, kwargs...)
end

function fit(
        ::Type{FrölichMelly},
        MM::OneStageModel,
        MD::Microdata;
        novar::Bool = false,
        trim::AbstractFloat = 0.0,
    )

    w = getweights(MD)
    d = getvector(MD, :treatment)
    z = getvector(MD, :instrument)
    π = predict(MM)
    v = [(2.0 * di - 1.0) * (zi - πi) / (πi * (1.0 - πi)) for (di, zi, πi) in zip(d, z, π)]

    v[((trim .> π) .| (1.0 - trim .< π))] .= 0.0

    SSD                   = Microdata(MD)
    SSD.mapping[:control] = asgn(MD.model, (MD.mapping[:treatment], InterceptTerm{true}()))

    obj              = FrölichMelly()
    obj.first_stage  = MM
    obj.second_stage = OLS(SSD)
    obj.pscore       = π
    obj.eweights     = ProbabilityWeights(v)

    _fit!(second_stage(obj), reweight(w, obj.eweights))
    novar || _vcov!(obj, getcorr(obj), w)

    return obj
end

#==========================================================================================#

# SCORE (MOMENT CONDITIONS)

function score(obj::FrölichMelly)
    return lmul!(Diagonal(obj.eweights), score(second_stage(obj)))
end

# EXPECTED JACOBIAN OF SCORE × NUMBER OF OBSERVATIONS

jacobian(obj::FrölichMelly, ::UnitWeights) = jacobian(second_stage(obj), obj.eweights)

function jacobian(obj::FrölichMelly, w::AbstractWeights)
    jacobian(second_stage(obj), reweight(w, obj.eweights))
end

# EXPECTED JACOBIAN OF SCORE W.R.T. FIRST-STAGE PARAMETERS × NUMBER OF OBSERVATIONS

function crossjacobian(obj::FrölichMelly, ::UnitWeights)

    d = getvector(obj, :treatment)
    z = getvector(obj, :instrument)
    π = obj.pscore
    D = [- (2.0 * di - 1.0) * (zi / abs2(πi) + (1.0 - zi) / abs2(1.0 - πi))
         for (di, zi, πi) in zip(d, z, π)]

    D[iszero.(obj.eweights)] .= 0.0

    g₁ = jacobexp(obj.first_stage)
    g₂ = score(second_stage(obj))

    return g₂' * lmul!(Diagonal(D), g₁)
end

function crossjacobian(obj::FrölichMelly, w::AbstractWeights)

    d = getvector(obj, :treatment)
    z = getvector(obj, :instrument)
    π = obj.pscore
    D = [- wi * (2.0 * di - 1.0) * (zi / abs2(πi) + (1.0 - zi) / abs2(1.0 - πi))
         for (di, zi, πi, wi) in zip(d, z, π, w)]

    D[iszero.(obj.eweights)] .= 0.0

    g₁ = jacobexp(obj.first_stage)
    g₂ = score(second_stage(obj))

    return g₂' * lmul!(Diagonal(D), g₁)
end

#==========================================================================================#

# FITTED VALUES

predict(obj::FrölichMelly, MD::Microdata) = predict(second_stage(obj), MD)

# DERIVATIVE OF FITTED VALUES

jacobexp(obj::FrölichMelly) = jacobexp(second_stage(obj))

#==========================================================================================#

# UTILITIES

coefnames(obj::FrölichMelly) = coefnames(second_stage(obj))

max_stars_repo_name: lbittarello/Microeconometrics

max_issues_repo_name: lbittarello/Microeconometrics

max_forks_repo_name: lbittarello/Microeconometrics


text:
VERSION >= v"0.4.0-dev+6521" && __precompile__(true)

using Base.Dates, TimeSeries, DataFrames

module Quandl

using Base.Dates, TimeSeries, DataFrames, Requests, CSV

export quandlget,
       quandl,
       quandlsearch,
       set_auth_token,
       interactivequandl

include("api.jl")
include("timearray.jl")
include("dataframe.jl")

# Create empty auth token if none exists
if !isfile(joinpath(dirname(@__FILE__),"../token/auth_token"))
    set_auth_token("")
end

end


text:
using Base.Test, MathProgComplex, OPFInstances

## Beware, rank relaxation value is for problem without objective constant term

@testset "ACOPFs - real formulation - order 1" begin

    testsCases = OrderedSet([("WB2", 1),
                            ("LMBM3", 1),
                            ("WB5", 1),
                            ("case6ww", 1),
                            ("case9mod", 1),
                            ("case9", 1),
                            ("case14", 1)])
                            # ("case24_ieee_rts", 1),
                            # ("case_ieee30", 1),
                            # ("case30pwl", 1),
                            # ("case30", 1)])

    testfolder = joinpath(Pkg.dir("MathProgComplex"), "Mosek_runs", "tests", "sos_example6")

    @testset "instance $instance, (rmeqs=$rmeqs)" for (instance, order) in testsCases, rmeqs in Set([false])

        problem_c, point = import_from_dat(getinstancepath("Matpower", "QCQP", instance))
        problem = pb_cplx2real(problem_c)

        relax_ctx = set_relaxation(problem; hierarchykind=:Real,
                                            d = order,
                                            params = Dict(:opt_outlev=>0,
                                                          :opt_solver=>testsolver))

        primobj, dualobj = run_hierarchy(problem, relax_ctx);

        ε_rel = OPFsols[(instance, 1)].rel_opt_gap
        ε_abs = max(primobj, dualobj) * ε_rel

        @show primobj, dualobj,  OPFsols[(instance, 1)].primal_solvalue, ε_abs, instance

        @test primobj ≈ OPFsols[(instance, 1)].primal_solvalue atol=ε_abs
        @test dualobj ≈ primobj atol=ε_abs
    end
end


text:
import FEM.extload
import FEM.assembleK!
import FEM.intf
import FEM.assemble_intf
import FEM.updatedofs!

facts("FEM.FEProblem") do

# Setup a FEProblem
mesh = GeoMesh()

nodes = [GeoNode2(1, [0, 0]), GeoNode2(2, [1, 1]), GeoNode2(3, [1, 2]), GeoNode2(4, [0, 1])]

push!(mesh, nodes)

push!(mesh, GeoTrig(1, [1, 2, 3]))
push!(mesh, GeoTrig(2, [1, 2, 4]))

bottom_set = NodeSet("y0", [1])
top_set = NodeSet("x0", [2, 3])
push!(mesh, bottom_set)
push!(mesh, top_set)

element_set = ElementSet("all", [1, 2])
push!(mesh, element_set)

bcs =  Any[DirichletBC("0.1", [FEM.Du, FEM.Dv], mesh.node_sets["x0"])]
loads =  Any[NodeLoad("10e5", [FEM.Dv], mesh.node_sets["y0"])]

# Element set
mat_section = MaterialSection(LinearIsotropic(250e9, 0.3))
push!(mat_section, mesh.element_sets["all"])

ele_section = ElementSection(LinTrig)
push!(ele_section, mesh.element_sets["all"])

fp = create_feproblem("cook_example_quad", mesh, [ele_section], [mat_section], bcs, loads)

context("FEM.FEProblem") do
    K = FEM.create_sparse_structure(fp::FEProblem)
    colptrs = FEM.get_colptrs(K, fp::FEProblem)

    load = extload(fp)

    int_f = assemble_intf(fp)
    assembleK!(K, fp, colptrs)
    du = K \ (load - int_f)

    updatedofs!(fp, du)

    # Verify force balance
    load = extload(fp)
    int_f = assemble_intf(fp)
    @fact norm(load - int_f) / norm(load) => roughly(0.0, 10.0^(-10))

end # context

end # facts


text:
# AUTO GENERATED FILE - DO NOT EDIT

export tabs

"""
    tabs(;kwargs...)
    tabs(children::Any;kwargs...)
    tabs(children_maker::Function;kwargs...)

A Tabs component.
Switch between different views. For more information, see: https://mantine.dev/core/tabs/
Keyword arguments:
- `children` (a list of or a singular dash component, string or number; optional): <Tab /> components only
- `id` (String; optional): The ID of this component, used to identify dash components in callbacks
- `active` (Real; optional): Index of active tab, overrides internal state
- `className` (String; optional): Often used with CSS to style elements with common properties
- `color` (a value equal to: "dark", "gray", "red", "pink", "grape", "violet", "indigo", "blue", "cyan", "teal", "green", "lime", "yellow", "orange"; optional): Active tab color from theme.colors
- `grow` (Bool; optional): True if tabs should take all available space
- `orientation` (a value equal to: "horizontal", "vertical"; optional): Controls tab orientation
- `position` (a value equal to: "right", "center", "left", "apart"; optional): Tab controls position
- `style` (Dict; optional): Inline style override
- `tabPadding` (a value equal to: "xs", "sm", "md", "lg", "xl"; optional): Controls tab content padding-top
- `variant` (a value equal to: "default", "outline", "pills"; optional): Controls appearance
"""
function tabs(; kwargs...)
        available_props = Symbol[:children, :id, :active, :className, :color, :grow, :orientation, :position, :style, :tabPadding, :variant]
        wild_props = Symbol[]
        return Component("tabs", "Tabs", "dash_mantine_components", available_props, wild_props; kwargs...)
end

tabs(children::Any; kwargs...) = tabs(;kwargs..., children = children)
tabs(children_maker::Function; kwargs...) = tabs(children_maker(); kwargs...)


text:
@testset "LineString" begin
    ps = [SVector(0,0), SVector(0,1), SVector(1,1), SVector(1,0)]
    ls = LineString(ps)

    @test length(ls) == 3
    @test ls[2] == Line(SVector(0,1), SVector(1,1))

    @test boundingbox(ls) == BoundingBox(0,0,typemin(Int),1,1,typemax(Int))
    @test isclosed(ls) == false
    @test issimple(ls) == false

    @test convert2d(ls) == ls
    @test convert3d(ls) == LineString([SVector(0,0,0), SVector(0,1,0), SVector(1,1,0), SVector(1,0,0)])
    @test convert3d(ls, 1) == LineString([SVector(0,0,1), SVector(0,1,1), SVector(1,1,1), SVector(1,0,1)])
    @test convert3d(LineString([SVector(0,0,1), SVector(0,1,1), SVector(1,1,1), SVector(1,0,1)]), 2) == LineString([SVector(0,0,2), SVector(0,1,2), SVector(1,1,2), SVector(1,0,2)])

    push!(ps, SVector(0,0))
    @test length(ls) == 4
    @test isclosed(ls) == true
    @test issimple(ls) == true

    @test issimple(LineString([SVector(0,0), SVector(0,1), SVector(1,0), SVector(1,1), SVector(0,0)])) == false

    @testset "Winding number" begin
        ls_clockwise = LineString([SVector(0.0,0.0), SVector(0.0,1.0), SVector(1.0,1.0), SVector(1.0,0.0), SVector(0.0,0.0)])
        ls_anticlockwise = LineString([SVector(0.0,0.0), SVector(1.0,0.0), SVector(1.0,1.0), SVector(0.0,1.0), SVector(0.0,0.0)])

        p_centre = SVector(0.5, 0.5)
        ps_outer_ring = SVector{2,Float64}[[-1.0,-1.0],[0.0,-1.0],[0.5,-1.0],[1.0,-1.0],
                                           [2.0,-1.0],[2.0,0.0],[2.0,0.5],[2.0,1.0],
                                           [2.0,2.0],[1.0,2.0],[0.5,2.0],[0.0,2.0],
                                           [-1.0,2.0],[-1.0,1.0],[-1.0,0.5],[-1.0,0.0]]

        @test winding_number(p_centre, ls_clockwise) == 1
        @test winding_number(p_centre, ls_anticlockwise) == -1

        for p in ps_outer_ring
            @test winding_number(p, ls_clockwise) == 0
            @test winding_number(p, ls_anticlockwise) == 0
        end
    end
end