text:
module t

using BioSequences

function composition(x::BioSequences.EachKmerIterator{T,S}) where {T,S}
    if x.step == 1
        return composition_nostep(T, x.seq)
    else
        return composition_step(x)
    end
end

############## Everything below here is internals ##########################

# Array lookups are faster to determine the two-bit value of a nucleotide
# than using trailing_zeros - this shaves off about 15% time.
const kmerof = [0x04, 0x00, 0x01, 0x04,
                0x02, 0x04, 0x04, 0x04,
                0x03, 0x04, 0x04, 0x04,
                0x04, 0x04, 0x04, 0x04]

# We define this function by itself, so that `composition` type-specializes for
# K > 8 for Dict and Array, respectively.
increment!(collection::Array, T, index) = @inbounds collection[index] += 1

function increment!(collection::Dict, T, index)
    key = reinterpret(T, index-1)
    # Defining "value" and using that in the setindex! instead of using += is faster
    get!(collection, key, 0)
    collection[key] = value + 1
    return nothing
end

function convert_to_composition(counts, T)
    if counts isa Dict
        return Composition{T}(counts)
    else
        dictionary = Dict{T,Int}()
        for i in eachindex(counts)
            @inbounds c = counts[i]
            if c > 0
                dictionary[reinterpret(T, i-1)] = c
            end
        end
        return Composition{T}(dictionary)
    end
end

function composition_nostep(::Type{T}, x::BioSequence{A}) where {T<:Kmer, A<:NucAlphs}
    BioSequences.checkkmer(T)
    kmer = typemin(UInt64)

    # When the unfilled bases reaches 0, we have observed K unambiguous nucleotides.
    # and we can increment the Kmercounter. Reset unfilled if we see an
    # ambiguous nucleotide.
    unfilled = BioSequences.kmersize(T) - 1
    mask = UInt64(1) << (2 * BioSequences.kmersize(T)) - UInt64(1)

    if BioSequences.kmersize(T) ≤ 8
        counts = zeros(Int, 4^BioSequences.kmersize(T))
    else
        counts = Dict{T, Int}()
    end

    for i in 1:length(x)
        nuc = BioSequences.inbounds_getindex(x, i)
        @inbounds value = kmerof[reinterpret(UInt8, nuc) + 1]
        if value == 0x04 # ambiguous nucleotide
            unfilled = BioSequences.kmersize(T) - 1
        elseif unfilled == 0
            kmer |= value
            increment!(counts, T, kmer & mask + 1)
            kmer <<= 2
        else
            kmer |= value
            unfilled -= 1
            kmer <<= 2
        end
    end
    return convert_to_composition(counts, T)
end

function comp_kmers(x::BioSequences.OneStepKmerIterator{T}) where {T}
    counts = zeros(Int, 4^BioSequences.kmersize(T))
    for (pos, kmer) in x
        @inbounds counts[reinterpret(Int, kmer) + 1] += 1
    end
    convert_to_composition(counts, T)
end

function composition_nostep2(::Type{T}, x::BioSequence{A}) where {T<:Kmer, A<:NucAlphs}
    BioSequences.checkkmer(T)
    kmer = typemin(UInt64)

    # When the unfilled bases reaches 0, we have observed K unambiguous nucleotides.
    # and we can increment the Kmercounter. Reset unfilled if we see an
    # ambiguous nucleotide.
    unfilled = BioSequences.kmersize(T) - 1
    mask = UInt64(1) << (2 * BioSequences.kmersize(T)) - UInt64(1)

    if BioSequences.kmersize(T) ≤ 8
        counts = zeros(Int, 4^BioSequences.kmersize(T))
    else
        counts = Dict{T, Int}()
    end

    for i in 1:length(x)
        nuc = BioSequences.inbounds_getindex(x, i)
        if isambiguous(nuc)
            unfilled = BioSequences.kmersize(T) - 1
        elseif unfilled == 0
            kmer |= trailing_zeros(nuc)
            increment!(counts, T, kmer & mask + 1)
            kmer <<= 2
        else
            kmer |= trailing_zeros(nuc)
            unfilled -= 1
            kmer <<= 2
        end
    end
    return convert_to_composition(counts, T)
end

function composition_step(x::BioSequences.EachKmerIterator{T}) where {T}
    BioSequences.checkkmer(T)
    kmer = typemin(UInt64)

    # When the unfilled bases reaches 0, we have observed K unambiguous nucleotides.
    # and we can increment the Kmercounter. Reset unfilled if we see an
    # ambiguous nucleotide.
    unfilled = BioSequences.kmersize(T) - 1
    mask = UInt64(1) << (2 * BioSequences.kmersize(T)) - UInt64(1)
    jumpsize = max(1, x.step - BioSequences.kmersize(T) + 1)

    if BioSequences.kmersize(T) ≤ 8
        counts = zeros(Int, 4^BioSequences.kmersize(T))
    else
        counts = Dict{T, Int}()
    end

    frame = 0
    i = x.start
    while i <= length(x.seq)
        nuc = BioSequences.inbounds_getindex(x.seq, i)
        @inbounds value = kmerof[reinterpret(UInt8, nuc) + 1]
        if value == 0x04 # ambiguous nucleotide
            unfilled = BioSequences.kmersize(T) - 1
            continue
        elseif unfilled == 0
            kmer |= value
            if frame == 0
                increment!(counts, T, kmer & mask + 1)
            end
            kmer <<= 2
        else
            unfilled -= jumpsize
            kmer |= value
            kmer <<= 2
        end
        i += jumpsize
        unfilled = jumpsize - 1
        frame += jumpsize
        if frame == x.step
            frame = 0
        end
    end
    return convert_to_composition(counts, T)
end

function extract_kmer_impl(seq, from, k)
    kmer::UInt64 = 0
    isok = true
    for i in 1:k
        nt = inbounds_getindex(seq, from+i-1)
        kmer = kmer << 2 | trailing_zeros(nt)
        isok &= iscertain(nt)
    end
    return kmer, isok
end

####################### Benchmarking ########################
# The new implementation is slower

# The new implementation is faster in all cases, but the difference
# is only appreciable for long sequences, especially at low K and higher stepsizes
# In the best cases, the new implementation is ~6x faster on my laptop.
# In more typical cases, it's 3-4 times faster.

using BenchmarkTools

function benchmark()
    rand_100_4bit = randdnaseq(100)
    rand_100_2bit = BioSequence{DNAAlphabet{2}}(rand_100_4bit)
    rand_10k_4bit = randdnaseq(10000)
    rand_10k_2bit = BioSequence{DNAAlphabet{2}}(rand_10k_4bit)
    rand_1m_4bit = randdnaseq(1000000)
    rand_1m_2bit = BioSequence{DNAAlphabet{2}}(rand_1m_4bit)

    # rand_100_4bit, rand_100_2bit,
    for sequence in (rand_10k_4bit, rand_10k_2bit,
                     rand_1m_4bit, rand_1m_2bit)
        for k in [4, 20] # One which fits in array, one which doesn't
            for step in [1,] # k is worst case, highvalue (e.g. 25) is best case, 1 is usual case
                kmerit = each(Kmer{DNA, k}, sequence, step)
                bits = typeof(sequence).parameters[1].parameters[1]
                println(length(sequence), " nt, ", bits, " bit, k = ", k, ", step = ", step)

                print("Old: ")
                x = BioSequences.composition(kmerit)
                @btime BioSequences.composition($kmerit);

                print("New: ")
                x = composition(kmerit)
                @btime composition($kmerit);
                println("")
            end
        end
    end
end

end # module

max_stars_repo_name: jakobnissen/biobag

max_issues_repo_name: jakobnissen/biobag

max_forks_repo_name: jakobnissen/biobag


text:
abstract type Controller end

mutable struct PID{T} <: Controller
    P::T
    I::T
    D::T

    eqcid::Int64
    goal::T
    integratedError::T

    function PID(eqcid::Int64, goal::T; P::T = zero(T), I::T = zero(T), D::T = zero(T)) where T
        @assert !(P == 0 && I == 0 && D == 0)

        new{T}(P, I, D, eqcid, goal, 0)
    end
end

@inline stateError(mechanism, eqc, goal, N) = goal - minimalCoordinates(eqc, mechanism, K = N)[1]

function control!(mechanism, pid::PID)
    e0 = stateError(mechanism, geteqconstraint(mechanism, pid.eqcid), pid.goal, 1)
    e1 = stateError(mechanism, geteqconstraint(mechanism, pid.eqcid), pid.goal, 2)

    pid.integratedError += e1 * mechanism.Δt
    differentialError = (e1 - e0) / mechanism.Δt
    perror = e1

    # return pid.P*perror + pid.I*pid.integratedError + pid.D*differentialError
    τ = SVector{3,Float64}(1, 0, 0) * (pid.P * perror + pid.I * pid.integratedError + pid.D * differentialError)
    setForce!(mechanism, mechanism.bodies[1], τ = τ)
end

# abstract type Controller end

# mutable struct PID{T,N} <: Controller
#     P::T
#     I::T
#     D::T

#     eqcids::SVector{N,Int64}
#     goal::SVector{N,T}
#     integratedError::SVector{N,T}

#     function PID(eqcids::AbstractVector{Int64},goal::AbstractVector{T}; P::T=zero(T), I::T=zero(T), D::T=zero(T)) where T
#         @assert !(P==0 && I==0 && D==0)

#         new{T}(P,I,D,eqcids,goal,zeros(T,length(goal)))
#     end

#     function PID(eqcids::Int64,goal::T; P::T=zero(T), I::T=zero(T), D::T=zero(T)) where T
#         @assert !(P==0 && I==0 && D==0)

#         new{T}(P,I,D,[eqcids],[goal],[0])
#     end
# end

# @inline stateError(mechanism,eqc,goal,N) = goal - minimalCoordinates(eqc,mechanism,K=N)

# function control!(mechanism,pid::PID)
#     e0 = stateError(mechanism,geteqconstraint(mechanism,pid.eqcid),pid.goal,1)
#     e1 = stateError(mechanism,geteqconstraint(mechanism,pid.eqcid),pid.goal,2)

#     pid.integratedError += e1*mechanism.Δt
#     differentialError = (e1-e0)/mechanism.Δt
#     perror = e1

#     # F, τ =
#     # return pid.P*perror + pid.I*pid.integratedError + pid.D*differentialError
#     τ = SVector{3,Float64}(1,0,0)*(pid.P*perror + pid.I*pid.integratedError + pid.D*differentialError)
#     setForce!(mechanism,mechanism.bodies[1],τ=τ)
# end

ext: jl


text:
struct vdWParam <: EoSParam
    Tc::SingleParam{Float64}
    pc::SingleParam{Float64}
    Mw::SingleParam{Float64}
    a::PairParam{Float64}
    b::PairParam{Float64}

end

abstract type vdWModel <: ABCubicModel end
@newmodel vdW vdWModel vdWParam

export vdW
function vdW(components::Array{String,1}; userlocations::Array{String,1}=String[], verbose=false,idealmodel=BasicIdeal)
    params = getparams(components, ["properties/critical.csv", "properties/molarmass.csv","SAFT/PCSAFT/PCSAFT_unlike.csv"]; userlocations=userlocations, verbose=verbose)
    Mw = params["Mw"]
    k = params["k"]
    _pc = params["pc"]
    _Tc = params["Tc"]
    pc = _pc.values
    Tc = _Tc.values

    a = epsilon_LorentzBerthelot(SingleParam(params["pc"], @. 27/64*R̄^2*Tc^2/pc), k)
    b = sigma_LorentzBerthelot(SingleParam(params["pc"], @. 1/8*R̄*Tc/pc))

    packagedparams = vdWParam(_Tc,_pc,Mw,a,b)
    model = vdW(packagedparams,idealmodel)
    return model
end

function ab_consts(::Type{<:vdWModel})
    Ωa =  27/64
    Ωb =  1/8
    return Ωa,Ωb
end

function cubic_ab(model::vdWModel,T,z=SA[1.0],Σn=sum(z))
    inv2Σn = (one(T)/Σn)^2
    a = model.params.a.values
    b = model.params.b.values
    ā = dot(z,Symmetric(a),z)*inv2Σn
    b̄ = dot(z,Symmetric(b),z)*inv2Σn
    return ā,b̄
end

function cubic_abp(model::vdWModel, V, T, z=@SVector [1.0])
    n = ∑(z)
    v = V/n
    a,b = cubic_ab(model,T,z,n)
    p = R̄*T/(v-b) - a/(v*v)
    return a,b,p
end

function cubic_poly(model::vdWModel,p,T,z)
    x = z/sum(z)
    a,b = cubic_ab(model,T,x)
    _1 = one(a+b)
    RT⁻¹ = 1/(R̄*T)
    A = a*p*RT⁻¹*RT⁻¹
    B = b*p*RT⁻¹
    return [-A*B, A, -B-_1, _1]
end

function a_res(model::vdWModel, V, T, z)
    n = sum(z)
    ā,b̄ = cubic_ab(model,T,z,n)
    RT⁻¹ = 1/(R̄*T)
    ρ = n/V
    return -log(1-b̄*ρ) - ā*ρ*RT⁻¹
    # -log(1-n*b/B)-a*n/(R̄*T*B)
    return -log(1-n*b̄/V) - ā*n/(R̄*T*V)
    #return -log(V-n*b̄) - ā*n/(R̄*T*V) + log(V)
end

cubic_zc(::vdWModel) = 3/8


text:
include("agent/DQN.jl")
# include("agent/Online.jl")

# Settings
include("agent/settings/DQN.jl")


text:
__precompile__()

module DPMMSubClustersStreaming

using Distributed
using DistributedArrays
using StatsBase
using Distributions
using SpecialFunctions
using LinearAlgebra
using Clustering
using KernelFunctions

import Random:seed!

include("ds.jl")

#Priors:
include("priors/niw.jl")
include("priors/multinomial_prior.jl")
include("priors/compact_mnm_prior.jl")

#Distributions:
include("distributions/mv_gaussian.jl")
include("distributions/multinomial_dist.jl")
include("distributions/compact_mnm_dist.jl")

include("utils.jl")
include("shared_actions.jl")
include("local_clusters_actions.jl")
include("global_params.jl")
include("dp-parallel-sampling.jl")
include("data_generators.jl")

export generate_gaussian_data, generate_mnmm_data, dp_parallel_sampling, dp_parallel, run_model_from_checkpoint, save_model, calculate_posterior, fit, get_labels_histogram,run_model_streaming,dp_parallel_streaming,get_labels,predict

end # module


text:
##################################
# Operador de cruce  por media ponderada
##############################
using Random

"""
    MediaPonderada((C1, C2)
Operador de cruce por ponderación
Genera dos hijos
"""
function MediaPonderada(C1, C2)
    a1 = rand() # valor entre cero y 1
    a2 = rand() # valor entre cero y 1
    return [a1*C1 + (1-a1)*C2, a2*C1 + (1-a2)*C2]
end


text:
# X and Z components in fT
# conducting 100 Ωm
## model
zfixed   = [-1e6,   0, 40, 60]
rho      = [1e12,   20, 10, 1000]

ross_primary = [25.7847     -18.7313]
ross_secondary = [6.7652     -8.0762
    4.7404     -6.6033
    3.7871     -5.7704
    2.9434     -4.9302
    2.0732      -3.932
    1.2823     -2.8452
   0.66508     -1.7978
   0.29382    -0.99918
   0.11343    -0.49729
  0.038092    -0.22124
  0.011105   -0.088032
0.0028947   -0.032018
0.00070238   -0.010946
0.00015894  -0.0036284
1.0855e-05  -0.0010876]

max_stars_repo_name: GeoscienceAustralia/HiQGA

max_issues_repo_name: GeoscienceAustralia/HiQGA


text:
#! /usr/bin/julia

# Rosetta Code, Vampire number

function divisors{T<:Integer}(n::T)
    !isprime(n) || return [one(T), n]
    d = [one(T)]
    for (k, v) in factor(n)
        e = T[k^i for i in 1:v]
        append!(d, vec([i*j for i in d, j in e]))
    end
    sort(d)
end

function vampirefangs{T<:Integer}(n::T)
    fangs = T[]
    isvampire = false
    vdcnt = ndigits(n)
    fdcnt = vdcnt>>1
    iseven(vdcnt) || return (isvampire, fangs)
    !isprime(n) || return (isvampire, fangs)
    vdigs = sort(digits(n))
    d = divisors(n)
    len = length(d)
    len = iseven(len) ? len>>1 : len>>1 + 1
    for f in d[1:len]
        ndigits(f) == fdcnt || continue
        g = div(n, f)
        f%10!=0 || g%10!=0 || continue
        sort([digits(f), digits(g)]) == vdigs || continue
        isvampire = true
        append!(fangs, [f, g])
    end
    if isvampire
        fangs = reshape(fangs, (2,length(fangs)>>1))'
    end
    return (isvampire, fangs)
end

function showvampire{T<:Integer}(i::T, n::T, fangs::Array{T,2})
    s = @sprintf "%6d  %14d %s\n" i n join(fangs[1,:], "\u00d7")
    for i in 2:size(fangs)[1]
        s *= " "^23*join(fangs[i,:], "\u00d7")*"\n"
    end
    return s
end

vgoal = 25
vcnt = 0
dcnt = 0
println("Finding the first ", vgoal, " vampire numbers.")
println("     N         Vampire Fangs")
while vcnt < vgoal
    dcnt += 2
    for i in (10^(dcnt-1)):(10^dcnt-1)
        (isvampire, fangs) = vampirefangs(i)
        isvampire || continue
        vcnt += 1
        print(showvampire(vcnt, i, fangs))
        vcnt < vgoal || break
    end
end

test = [16758243290880, 24959017348650, 14593825548650]
println()
println("Checking a few numbers.")
println("     N         Vampire Fangs")
for (i, v) in enumerate(test)
    (isvampire, fangs) = vampirefangs(v)
    if isvampire
        print(showvampire(i, v, fangs))
    else
        println(@sprintf "%6d  %14d is not a vampire" i v)
    end
end

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
"""
spy(x::Range, y::Range, z::AbstractSparseArray)

Visualizes big sparse matrices.
Usage:
```julia
N = 200_000
x = sprand(Float64, N, N, (3(10^6)) / (N*N));
spy(x)
# or if you want to specify the range of x and y:
spy(0..1, 0..1, x)
```
## Attributes
$(ATTRIBUTES)
"""
@recipe(Spy, x, y, z) do scene
    Attributes(
        marker = automatic,
        markersize = automatic,
        colormap = theme(scene, :colormap),
        colorrange = automatic,
        framecolor = :black,
        framesize = 1,
        inspectable = theme(scene, :inspectable)
    )
end

function convert_arguments(::Type{<: Spy}, x::SparseArrays.AbstractSparseArray)
    (0..size(x, 1), 0..size(x, 2), x)
end
function convert_arguments(::Type{<: Spy}, x, y, z::SparseArrays.AbstractSparseArray)
    (x, y, z)
end

function calculated_attributes!(::Type{<: Spy}, plot)
end

function plot!(p::Spy)
    rect = lift(p.x, p.y) do x, y
        xe = extrema(x)
        ye = extrema(y)
        FRect2D((xe[1], ye[1]), (xe[2] - xe[1], ye[2] - ye[1]))
    end
    # TODO FastPixel isn't accepting marker size in data coordinates
    # but instead in pixel - so we need to fix that in GLMakie for consistency
    # and make this nicer when redoing unit support
    markersize = lift(p.markersize, rect, p.z) do msize, rect, z
        if msize === automatic
            widths(rect) ./ Vec2f0(size(z))
        else
            msize
        end
    end
    # TODO correctly align marker
    xycol = lift(rect, p.z, markersize) do rect, z, markersize
        x, y, color = SparseArrays.findnz(z)
        points = map(x, y) do x, y
            (((Point2f0(x, y) .- 1) ./ Point2f0(size(z) .- 1)) .*
            widths(rect) .+ minimum(rect))
        end
        points, convert(Vector{Float32}, color)
    end

    replace_automatic!(p, :colorrange) do
        lift(xycol) do (xy, col)
            extrema_nan(col)
        end
    end

    marker = lift(p.marker) do x
        return x === automatic ? FastPixel() : x
    end

    scatter!(
        p,
        lift(first, xycol), color = lift(last, xycol),
        marker = marker, markersize = markersize, colorrange = p.colorrange,
        inspectable = p.inspectable
    )

    lines!(p, rect, color = p.framecolor, linewidth = p.framesize, inspectable = p.inspectable)
end


text:
using OpenSMILES
using LightGraphs
using SimpleWeightedGraphs
using Test

@testset "Read Next Element" begin
    @test OpenSMILES.tryparseelement( "ScC", 1, OpenSMILES.bracket ) == ("Sc", 3)
    @test OpenSMILES.tryparseelement( "ScAs", 1, OpenSMILES.bracket ) == ("Sc", 3)
    @test OpenSMILES.tryparseelement( "ScAs", 3, OpenSMILES.bracket ) == ("As", 5)
    @test OpenSMILES.tryparseelement( "Sc", 1, OpenSMILES.bracket ) == ("Sc", 3)
    @test OpenSMILES.tryparseelement( "CO", 1, OpenSMILES.bracket ) == ("C", 2)
    @test OpenSMILES.tryparseelement( "C", 1, OpenSMILES.bracket ) == ("C", 2)
    @test OpenSMILES.tryparseelement( "s", 1, OpenSMILES.aromatics ) == ("s", 2)
    @test OpenSMILES.tryparseelement( "Oc", 1, [OpenSMILES.aliphatics; OpenSMILES.aromatics] ) == ("O", 2)
end

@testset "Read Next Numeric" begin
    @test OpenSMILES.tryparseint16("123UrC", 1) == (123, 4)
    @test OpenSMILES.tryparseint16("UrC", 1) == (nothing, 1)
    @test OpenSMILES.tryparseint16("12CH4", 1) == (12, 3)
    @test OpenSMILES.tryparseint16("12CH4", 5) == (4, 6)
end

@testset "Read Next Charge" begin
    @test OpenSMILES.tryparsecharge("+", 1) == (1, 2)
    @test OpenSMILES.tryparsecharge("-", 1) == (-1, 2)
    @test OpenSMILES.tryparsecharge("+++", 1) == (3, 4)
    @test OpenSMILES.tryparsecharge("---", 1) == (-3, 4)
    @test OpenSMILES.tryparsecharge("+2", 1) == (2, 3)
    @test OpenSMILES.tryparsecharge("-2", 1) == (-2, 3)
    @test OpenSMILES.tryparsecharge("Na+Cl-", 3) == ( 1, 4)
    @test OpenSMILES.tryparsecharge("Na+Cl-", 6) == (-1, 7)
end

@testset "Brackets Known" begin
    @test OpenSMILES.ParseBracket("22NaH") == OpenSMILES.Element("Na", 22, false, Int16[], 1, 0, 0)
    @test OpenSMILES.ParseBracket("O-") == OpenSMILES.Element("O", nothing, false, Int16[], 0, 0, -1)
    @test OpenSMILES.ParseBracket("CH4") == OpenSMILES.Element("C", nothing, false, Int16[], 4, 0, 0)
    @test OpenSMILES.ParseBracket("CH3-") == OpenSMILES.Element("C", nothing, false, Int16[], 3, 0, -1)
end

@testset "Brackets Equivalences" begin
    @test OpenSMILES.ParseBracket("CH2--") == OpenSMILES.ParseBracket("CH2-2")
    @test OpenSMILES.ParseBracket("CH1---") == OpenSMILES.ParseBracket("CH1-3")
    @test OpenSMILES.ParseBracket("Fe++") == OpenSMILES.ParseBracket("Fe+2")
    @test OpenSMILES.ParseBracket("Fe+++") == OpenSMILES.ParseBracket("Fe+3")
end

@testset "ParseOpenSMILES Checking Implicit Hydrogens on Simple Molecules" begin
    _, Data = OpenSMILES.ParseSMILES("CCCC(C)C")
    checkhydrogens = Dict( OpenSMILES.countitems( OpenSMILES.H.( Data ) ) )
    @test checkhydrogens[1] == 1
    @test checkhydrogens[2] == 2
    @test checkhydrogens[3] == 3

    #Check a Ring
    _, Data = OpenSMILES.ParseSMILES("C1CCC(C)C1")
    checkhydrogens = Dict( OpenSMILES.countitems( OpenSMILES.H.( Data ) ) )
    @test checkhydrogens[1] == 1
    @test checkhydrogens[2] == 4
    @test checkhydrogens[3] == 1

    # Atoms with multiple valence states (here, S)
    _, Data = OpenSMILES.ParseSMILES("OS(=O)(=O)O")
    @test OpenSMILES.EmpiricalFormula( Data ) == "H2O4S"
end

@testset "Chirality" begin
    # Chirality is not fully supported, but we shouldn't error
    # Tetrahedral
    for smiles in ("N[C@](Br)(O)C", "N[C@TH1](Br)(O)C", "N[C@@](Br)(O)C", "N[C@TH2](Br)(O)C")
        _, Data = OpenSMILES.ParseSMILES(smiles)
        @test OpenSMILES.EmpiricalFormula( Data ) == "BrC2H6NO"
    end
    # Allenal
    for smiles in ("NC(Br)=[C@]=C(O)C", "NC(Br)=[C@AL1]=C(O)C", "NC(Br)=[C@@]=C(O)C", "NC(Br)=[C@AL2]=C(O)C")
        _, Data = OpenSMILES.ParseSMILES(smiles)
        @test OpenSMILES.EmpiricalFormula( Data ) == "BrC4H6NO"
    end
    # Square-planar
    for smiles in ("F[Po@SP1](Cl)(Br)I", "F[Po@SP2](Cl)(Br)I", "F[Po@SP3](Cl)(Br)I")
        _, Data = OpenSMILES.ParseSMILES(smiles)
        @test OpenSMILES.EmpiricalFormula( Data ) == "BrClFIPo"
    end
    # Trigonal-bipyramidal
    for smiles in ("S[As@TB1](F)(Cl)(Br)N", "S[As@TB2](Br)(Cl)(F)N", "S[As@TB5](F)(N)(Cl)Br", "F[As@TB10](S)(Cl)(N)Br",
                   "F[As@TB15](Cl)(S)(Br)N", "Br[As@TB20](Cl)(S)(F)N")
        _, Data = OpenSMILES.ParseSMILES(smiles)
        @test OpenSMILES.EmpiricalFormula( Data ) == "AsBrClFH3NS"
    end
    # Octahedral
    for smiles in ("C[Co@](F)(Cl)(Br)(I)S", "F[Co@@](S)(I)(C)(Cl)Br", "S[Co@OH5](F)(I)(Cl)(C)Br",
                   "Br[Co@OH9](C)(S)(Cl)(F)I", "Br[Co@OH12](Cl)(I)(F)(S)C", "Cl[Co@OH15](C)(Br)(F)(I)S",
                   "Cl[Co@OH19](C)(I)(F)(S)Br", "I[Co@OH27](Cl)(Br)(F)(S)C")
        _, Data = OpenSMILES.ParseSMILES(smiles)
        @test OpenSMILES.EmpiricalFormula( Data ) == "BrCClCoFH4IS"
    end

    # Hydrogens as bracket atoms are equivalent to Hcount hydrogens
    strh = "N[C@H](O)C"                 # hcount
    strb = "N[C@]([H])(O)C"             # bracket
    _, Data = OpenSMILES.ParseSMILES(strh)
    @test OpenSMILES.EmpiricalFormula( Data ) == "C2H7NO"
    _, Data = OpenSMILES.ParseSMILES(strb)
    @test OpenSMILES.EmpiricalFormula( Data ) == "C2H7NO"
end

@testset "Ring identifiers" begin
    # spiro[5.5]undecane, https://pubchem.ncbi.nlm.nih.gov/#query=spiro%5B5.5%5Dundecane
    g1, elms1 = ParseSMILES("C12(CCCCC1)CCCCC2")     # single-digit ring identifiers
    g2, elms2 = ParseSMILES("C%012(CCCCC%01)CCCCC2") # two-digit ring identifiers
    @test g1 == g2 && elms1 == elms2
end

@testset "ParseOpenSMILES Empirical Formulas of Complicated Molecules" begin
    #Phenol
    g, Data = OpenSMILES.ParseSMILES("Oc1ccccc1")
    @test OpenSMILES.EmpiricalFormula( Data ) == "C6H6O"
    for i = 1:6
        @test SimpleWeightedGraphs.has_edge(g, i+1, mod1(i-1, 6)+1)   # +1 is because O is the first atom
        @test SimpleWeightedGraphs.has_edge(g, i+1, mod1(i+1, 6)+1)   #              "
    end
    #Anthracene
    _, Data = OpenSMILES.ParseSMILES("C1=CC=C2C=C3C=CC=CC3=CC2=C1")
    @test OpenSMILES.EmpiricalFormula( Data ) == "C14H10"
    #Tryptophan
    _, Data = OpenSMILES.ParseSMILES("C1=CC=C2C(=C1)C(=CN2)CC(C(=O)O)N")
    @test OpenSMILES.EmpiricalFormula( Data ) == "C11H12N2O2"
    #Lysergic Acid Diethylamide
    _, Data = OpenSMILES.ParseSMILES("CCN(CC)C(=O)C1CN(C2CC3=CNC4=CC=CC(=C34)C2=C1)C")
    @test OpenSMILES.EmpiricalFormula( Data ) == "C20H25N3O"
    #Firefly Luciferin
    _, Data = OpenSMILES.ParseSMILES("O=C(O)[CH]1N=C(SC1)c2sc3cc(O)ccc3n2")
    @test OpenSMILES.EmpiricalFormula( Data ) == "C11H8N2O3S2"
end

@testset "instantiate_hydrogens" begin
    g, atomnodes = OpenSMILES.ParseSMILES("O")
    gH, atomnodesH = instantiate_hydrogens(g, atomnodes)
    # Ensure the input is not modified
    @test SimpleWeightedGraphs.nv(g) == 1
    @test length(atomnodes) == 1 && atomnodes[1].symbol == "O"
    @test OpenSMILES.EmpiricalFormula(atomnodes) == "H2O"
    # Check the instantiated version
    @test SimpleWeightedGraphs.nv(gH) == 3
    @test length(atomnodesH) == 3
    @test atomnodesH[1].symbol == "O"
    @test atomnodesH[2].symbol == "H"
    @test atomnodesH[3].symbol == "H"
    @test SimpleWeightedGraphs.has_edge(gH, 1, 2) && SimpleWeightedGraphs.has_edge(gH, 1, 3) && !SimpleWeightedGraphs.has_edge(gH, 2, 3)
    @test OpenSMILES.EmpiricalFormula(atomnodesH) == "H2O"

    g, atomnodes = OpenSMILES.ParseSMILES("[H]O[H]")
    gH, atomnodesH = instantiate_hydrogens(g, atomnodes)
    @test g == gH && atomnodes == atomnodesH
    @test g !== gH && atomnodes !== atomnodesH    # it's a copy, not the same object
end


text:
const SVF = MOI.SingleVariable
const SAF = MOI.ScalarAffineFunction{Float64}

# indices
const VI = MOI.VariableIndex
const CI = MOI.ConstraintIndex

# sets
const BOUNDS = Union{
    MOI.EqualTo{Float64},
    MOI.GreaterThan{Float64},
    MOI.LessThan{Float64},
    MOI.Interval{Float64}
}

const VAR_TYPES = Union{
    MOI.ZeroOne,
    MOI.Integer
}

"""
Optimizer struct
"""
mutable struct Optimizer <: MOI.AbstractOptimizer
    inner::Union{CoM, Nothing}
    variable_info::Vector{Variable}
    # which variable index, (:leq,:geq,:eq,:Int,:Bin), and lower and upper bound
    var_constraints::Vector{Tuple{Int64,Symbol,Int64,Int64}}
    status::MOI.TerminationStatusCode
    options::SolverOptions
end

include("variables.jl")
include("constraints.jl")
include("objective.jl")
include("results.jl")

MOI.get(::Optimizer, ::MOI.SolverName) = "ConstraintSolver"

"""
Optimizer struct constructor
"""
function Optimizer(;options...)
    com = CS.init()
    options = combine_options(options)
    return Optimizer(
        com,
        [],
        [],
        MOI.OPTIMIZE_NOT_CALLED,
        options
    )
end

"""
    MOI.is_empty(model::Optimizer)
"""
function MOI.is_empty(model::Optimizer)
    return isempty(model.variable_info) &&
           isempty(model.var_constraints)
end

"""
    MOI.empty!(model::Optimizer)
"""
function MOI.empty!(model::Optimizer)
    model.inner = CS.init()
    empty!(model.variable_info)
    empty!(model.var_constraints)
    model.status = MOI.OPTIMIZE_NOT_CALLED
    # !important => don't remove the options
end

"""
Copy constructor for the optimizer
"""
MOIU.supports_default_copy_to(model::Optimizer, copy_names::Bool) = !copy_names
function MOI.copy_to(model::Optimizer, src::MOI.ModelLike; kws...)
    return MOI.Utilities.automatic_copy_to(model, src; kws...)
end

"""
    MOI.optimize!(model::Optimizer)
"""
function MOI.optimize!(model::Optimizer)
    # check if every variable has bounds and is an Integer
    check_var_bounds(model)

    set_pvals!(model)
    set_constraint_hashes!(model)

    status = solve!(model)
    set_status!(model, status)
end


text:
module ControlSystemIdentification

using DSP, LinearAlgebra, Statistics, Random, Optim, ControlSystems, FillArrays, Parameters, TotalLeastSquares, RecipesBase, FFTW, Roots, MonteCarloMeasurements
import Optim: minimizer, Options

export StateSpaceNoise, pem, simulation_errors, prediction_errors, predict, simulate, noise_model
export getARXregressor, getARregressor, find_na, arx, ar, bodeconfidence, tls, wtls_estimator, plr
export FRD, tfest, coherence, coherenceplot, simplot, simplot!, predplot, predplot!

include("utils.jl")
include("types.jl")
include("pem.jl")
include("arx.jl")
include("frd.jl")

function predict(sys, y, u, x0=zeros(sys.nx))
	model = SysFilter(sys, copy(x0))
	yh = [model(yt,ut) for (yt,ut) in observations(y,u)]
	oftype(y,yh)
end
predict(sys::ControlSystems.TransferFunction, args...) = predict(ss(sys), args...)

"""
	yh = predict(ar::TransferFunction, y)

Predict AR model
"""
function predict(G::ControlSystems.TransferFunction, y)
	_,a,_ = params(G)
	yr,A = getARregressor(y,length(a))
	yh = A*a
	oftype(y,yh)
end

function simulate(sys, u, x0=zeros(sys.nx))
	model = SysFilter(sys, copy(x0))
	yh = map(observations(u,u)) do (ut,_)
		model(ut)
	end
	oftype(u,yh)
end
simulate(sys::ControlSystems.TransferFunction, args...) = simulate(ss(sys), args...)

@userplot Simplot
"""
simplot(sys, y, u, x0=zeros(sys.nx); ploty=true)
Plot system simulation and measured output to compare them.
`ploty` determines whether or not to plot the measured signal
"""
simplot
@recipe function simplot(p::Simplot; ploty=true)
	sys,y,u = p.args[1:3]
	y = oftype(randn(2,2), y)
	u = oftype(randn(2,2), u)
	x0 = length(p.args) > 3 ? p.args[4] : zeros(sys.nx)
	yh = simulate(sys,u, x0)
	xguide --> "Time [s]"
	yguide --> "Output"
	t = range(0, step=sys.Ts, length=length(y))
	err = nrmse(y,yh)
	ploty && @series begin
		label --> "y"
		t,y'
	end
	@series begin
		label --> "sim fit :$(round.(err, digits=2))%"
		t,yh'
	end
	nothing
end

@userplot Predplot
"""
predplot(sys, y, u, x0=zeros(sys.nx); ploty=true)
Plot system simulation and measured output to compare them.
`ploty` determines whether or not to plot the measured signal
"""
predplot
@recipe function predplot(p::Predplot; ploty=true)
	sys,y,u = p.args[1:3]
	y = oftype(randn(2,2), y)
	u = oftype(randn(2,2), u)
	x0 = length(p.args) > 3 ? p.args[4] : zeros(sys.nx)
	yh = predict(sys,y,u, x0)
	xguide --> "Time [s]"
	yguide --> "Output"
	t = range(0, step=sys.Ts, length=length(y))
	err = nrmse(y,yh)
	ploty && @series begin
		label --> "y"
		t,y'
	end
	@series begin
		label --> "pred fit :$(round.(err, digits=2))%"
		t,yh'
	end
	nothing
end

function ControlSystems.lsim(sys::StateSpaceNoise, u; x0=zeros(sys.nx))
	simulate(sys, u, x0)
end

end # module


text:
#CLIENT

#HTTP.request sends a HTTP Request Message and returns a Response Message.

r = HTTP.request("GET", "http://httpbin.org/ip"; verbose=3)
println(r.status)
println(String(r.body))

#HTTP.open sends a HTTP Request Message and opens an IO stream from which the Response can be read.
HTTP.open(:GET, "https://tinyurl.com/bach-cello-suite-1-ogg") do http
    open(`vlc -q --play-and-exit --intf dummy -`, "w") do vlc
        write(vlc, http)
    end
end

#SERVERS

#Using HTTP.Servers.listen:
#The server will start listening on 127.0.0.1:8081 by default.

using HTTP

HTTP.listen() do http::HTTP.Stream
    @show http.message
    @show HTTP.header(http, "Content-Type")
    while !eof(http)
        println("body data: ", String(readavailable(http)))
    end
    HTTP.setstatus(http, 404)
    HTTP.setheader(http, "Foo-Header" => "bar")
    HTTP.startwrite(http)
    write(http, "response body")
    write(http, "more response body")
end

#Using HTTP.Handlers.serve:

using HTTP

HTTP.serve() do request::HTTP.Request
   @show request
   @show request.method
   @show HTTP.header(request, "Content-Type")
   @show HTTP.payload(request)
   try
       return HTTP.Response("Hello")
   catch e
       return HTTP.Response(404, "Error: $e")
   end
end

#WebSocket Examples
@async HTTP.WebSockets.listen("127.0.0.1", UInt16(8081)) do ws
    while !eof(ws)
        data = readavailable(ws)
        write(ws, data)
    end
end

HTTP.WebSockets.open("ws://127.0.0.1:8081") do ws
    write(ws, "Hello")
    x = readavailable(ws)
    @show x
    println(String(x))
end;
x = UInt8[0x48, 0x65, 0x6c, 0x6c, 0x6f]
#Output: Hello

#=Custom HTTP Layer Examples
Notes:
There is no enforcement of a "well-defined" stack, you can insert a layer anywhere in the stack even if it logically does not make sense
When creating a custom layer, you need to create a request(), see below for an example
Custom layers is only implemented with the "low-level" request() calls, not the "convenience" functions such as HTTP.get(), HTTP.put(), etc.
module TestRequest=#
        import HTTP: Layer, request, Response

        abstract type TestLayer{Next <: Layer} <: Layer{Next} end
        export TestLayer, request

        function request(::Type{TestLayer{Next}}, io::IO, req, body; kw...)::Response where Next
                println("Insert your custom layer logic here!")
                return request(Next, io, req, body; kw...)
        end
end

using HTTP
using ..TestRequest

custom_stack = insert(stack(), StreamLayer, TestLayer)

result = request(custom_stack, "GET", "https://httpbin.org/ip")

# Insert your custom layer logic here!

# HTTP.Messages.Response:
# """
# HTTP/1.1 200 OK
# Access-Control-Allow-Credentials: true
# Access-Control-Allow-Origin: *
# Content-Type: application/json
# Date: Fri, 30 Aug 2019 14:13:17 GMT
# Referrer-Policy: no-referrer-when-downgrade
# Server: nginx
# X-Content-Type-Options: nosniff
# X-Frame-Options: DENY
# X-XSS-Protection: 1; mode=block
# Content-Length: 45
# Connection: keep-alive

# {
#   "origin": "--Redacted--"
# }
# """


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

using Markdown
using InteractiveUtils

# ╔═╡ c436330a-fced-11ea-132d-85f5197a22b8
using Pkg, DrWatson

# ╔═╡ c43683be-fced-11ea-0575-f78e4551118a
begin
	#@quickactivate "StatisticalRethinkingStan"
	using StanSample
	using StatisticalRethinking
end

# ╔═╡ 1e3eaf90-fced-11ea-12c7-9f994f24a65d
md"## Clip-05-13s.jl"

# ╔═╡ 4bd1668a-81df-11eb-2096-ebf30ab5cd8a
begin
	df = CSV.read(sr_datadir("WaffleDivorce.csv"), DataFrame);
	scale!(df, [:Marriage, :MedianAgeMarriage, :Divorce])
end;

# ╔═╡ 97b5bff4-81df-11eb-35f8-47a6d481ed6e
stan5_3 = "
data {
  int N;
  vector[N] divorce_s;
  vector[N] marriage_s;
  vector[N] medianagemarriage_s;
}
parameters {
  real a;
  real bA;
  real bM;
  real<lower=0> sigma;
}
model {
  vector[N] mu = a + + bA * medianagemarriage_s + bM * marriage_s;
  a ~ normal( 0 , 0.2 );
  bA ~ normal( 0 , 0.5 );
  bM ~ normal( 0 , 0.5 );
  sigma ~ exponential( 1 );
  divorce_s ~ normal( mu , sigma );
}
";

# ╔═╡ 74ff3772-81e1-11eb-00f2-11ca8eb479ea
# Rethinking results
rethinking_results = "
	   mean   sd  5.5% 94.5%
a      0.00 0.10 -0.16  0.16
bM    -0.07 0.15 -0.31  0.18
bA    -0.61 0.15 -0.85 -0.37
sigma  0.79 0.08  0.66  0.91
";

# ╔═╡ 97b5fcee-81df-11eb-153a-1d3c25d4ae12
begin
	m5_3s = SampleModel("m5.3", stan5_3);
	m5_3_data = Dict(
	  "N" => size(df, 1),
	  "divorce_s" => df[:, :Divorce_s],
	  "marriage_s" => df[:, :Marriage_s],
	  "medianagemarriage_s" => df[:, :MedianAgeMarriage_s]
	)
	rc5_3s = stan_sample(m5_3s, data=m5_3_data);
	if success(rc5_3s)

		post5_3s_df = read_samples(m5_3s, :dataframe)
		PRECIS(post5_3s_df)
	end
end

# ╔═╡ c4470928-fced-11ea-14d1-df38dd138e54
if success(rc5_3s)
	begin
		title = "Divorce rate vs. Marriage rate" *
			"\nshowing predicted and hpd range"
		plotbounds(
			df, :Marriage, :Divorce,
			post5_3s_df, [:a, :bM, :sigma];
			title=title,
			colors=[:lightgrey, :darkgrey],
			bounds=[:predicted, :hpdi]
		)
	end
end

# ╔═╡ c447c624-fced-11ea-23e3-297d5e57104e
md"## End of clip-05-13s.jl"

# ╔═╡ Cell order:
# ╟─1e3eaf90-fced-11ea-12c7-9f994f24a65d
# ╠═c436330a-fced-11ea-132d-85f5197a22b8
# ╠═c43683be-fced-11ea-0575-f78e4551118a
# ╠═4bd1668a-81df-11eb-2096-ebf30ab5cd8a
# ╠═97b5bff4-81df-11eb-35f8-47a6d481ed6e
# ╠═74ff3772-81e1-11eb-00f2-11ca8eb479ea
# ╠═97b5fcee-81df-11eb-153a-1d3c25d4ae12
# ╠═c4470928-fced-11ea-14d1-df38dd138e54
# ╟─c447c624-fced-11ea-23e3-297d5e57104e


text:
using LightGraphs, MetaGraphs
using Combinatorics: combinations, powerset

macro arrow_str(str)
    symbols = ['*', 'o', '>', '<', '-']
    symDict = Dict('o' => :circle,
                   '>' => :arrow,
                   '<' => :arrow,
                   '-' => :tail,
                   '*' => :star)

    @assert length(str)==3
    @assert str[1] ∈ symbols
    @assert str[2] == '-'
    @assert str[3] ∈ symbols

    return (symDict[str[1]], symDict[str[3]])
end

"""
    has_marks(dg, v1, v2, t::Tuple{Symbol, Symbol}

test if the edge between node v1 and v2 has the edge markers given by the tuple t (use
the arrow macro to simplify use)

Example:
has_marks(dg, 1, 2, arrow"o->")
"""
function has_marks(dg, v1, v2, t::Tuple{Symbol, Symbol})
    if t[2]!=:star
        if t[1]!=:star
            result = get_prop(dg, v2, v1, :mark)==t[1] && get_prop(dg, v1, v2, :mark)==t[2]
        else
            result = get_prop(dg, v1, v2, :mark)==t[2]
        end
    else
        result = get_prop(dg, v2, v1, :mark)==t[1]
    end

    return result
end

"""
    set_marks!(dg, v1, v2, t::Tuple{Symbol, Symbol})

set edge marks between node v1 and v2.

Example:
set_marks!(dg, 1, 2, arrow"*->")
"""
function set_marks!(dg, v1, v2, t::Tuple{Symbol, Symbol})
    if t[1]!=:star
        set_prop!(dg, v2, v1, :mark, t[1])
    end

    if t[2]!=:star
        set_prop!(dg, v1, v2, :mark, t[2])
    end
end

"""
    is_collider(dg, v1, v2, v3)

check if egde v1, v2 and v3 form a collider
"""
function is_collider(dg, v1, v2, v3)
    return has_marks(dg, v1, v2, arrow"*->") && has_marks(dg, v2, v3, arrow"<-*")
end

"""
    is_parent(dg, v1, v2)

check if v1 is a parent of v2
"""
function is_parent(dg, v1, v2)
    return has_edge(dg, v1, v2) && has_marks(dg, v1, v2, arrow"-->")
end

"""
    is_triangle(dg, v1, v2, v3)

check if v1, v2 and v3 form a triangle
"""
function is_triangle(dg, v1, v2, v3)
    return isadjacent(dg, v1, v2) && isadjacent(dg, v2, v3) && isadjacent(dg, v3, v1)
end

"""
    is_discriminating_path(dg, path)

check if `path` is a discriminating path
"""
function is_discriminating_path(dg, path)
    # a discriminating path consists of at least four edges
    if length(path)<4
        return false
    end

    triples = collect(zip(path[1:end-2], path[2:end-1], path[3:end]))[1:end-1]
    colliders = map(t->is_collider(dg, t...), triples)
    parents = map(t->is_parent(dg, t[2], path[end]), triples)
    if all(colliders) && all(parents)
        return true
    else
        return false
    end
end

"""
    is_uncovered_circle_path(dg, path)

check if `path` is an uncovered circle path
"""
function is_uncovered_circle_path(dg, path)
    if length(path)<4
        return false
    end

    edges = collect(zip(path[1:end-1], path[2:end]))
    triples = collect(zip(path[1:end-2], path[3:end]))

    unshielded = map(t->!isadjacent(dg, t[1], t[2]), triples)
    circles = map(e->has_marks(dg, e[1], e[2], arrow"o-o"), edges)

    return all(unshielded) && all(circles)
end

"""
    is_uncovered_PD_path(dg, path)

check if `path` is an uncovered potentially directed path
"""
function is_uncovered_PD_path(dg, path)
    if length(path)<4
        return false
    end

    edges = collect(zip(path[1:end-1], path[2:end]))
    triples = collect(zip(path[1:end-2], path[3:end]))
    unshielded = map(t->!isadjacent(dg, t[1], t[2]), triples)
    directions = map(e->(!has_marks(dg, e[1], e[2], arrow"<-*") &&
                         !has_marks(dg, e[1], e[2], arrow"*--")), edges)

    return all(unshielded) && all(directions)
end

"""
    fcialg(n::V, I, par...; augmented=true, verbose=false, kwargs...)

Perform the FCI algorithm for a set of `n` variables using the test

    I(u, v, [s1, ..., sn], par...)

Returns the PAG as a MetaDiGraph
"""
function fcialg(n::V, I, par...; augmented=true, verbose=false, kwargs...) where {V<:Integer}

    # Step F1 and F2
    g, S = skeleton(n, I, par...; kwargs...)

    # Apply R0 once
    Z = unshielded(g, S)
    dg = MetaDiGraph(g) # use g to keep track of unoriented edges

    # construct initial PAG
    for e in edges(dg)
        set_prop!(dg, e, :mark, :circle)
    end

    for (u, v, w) in Z
        # is this check actually needed?
        if has_edge(dg, (u, v))
            set_marks!(dg, u, v, arrow"*->")
        end
        if has_edge(dg, (v, w))
            set_marks!(dg, v, w, arrow"<-*")
        end
    end

    # find possible d-separation sets (MAGs are trickier than DAGs...)
    pdsep = Dict()

    for v in vertices(g)
        pdsep[v] = Set{Int64}()
        for w in vertices(g)
            if w == v
                continue
            end
            paths = yen_k_shortest_paths(g, v, w, LightGraphs.weights(g), 100).paths
            for path in paths
                if length(path)==2
                    push!(pdsep[v], w)
                else
                    triples = zip(path[1:end-2], path[2:end-1], path[3:end])
                    tconds = map(t->is_triangle(dg, t...) || is_collider(dg, t...), triples)
                    if all(tconds)
                        push!(pdsep[v], w)
                    end
                end
            end
        end
    end

    # need to collect edges here since graph could change while looping
    for e in collect(edges(g))
        v = src(e)
        w = dst(e)
        if !(has_edge(g, v, w) || has_edge(g, v, w))
            # edge has already been removed
            continue
        end
        seps1 = Set(powerset([d for d in pdsep[v] if d!=w]))
        seps2 = Set(powerset([d for d in pdsep[w] if d!=v]))
        seps = union(seps1, seps2)

        for s in seps
            if I(src(e), dst(e), s, par...; kwargs...)
                @debug "Found hidden dsep set: $(src(e)) - $(dst(e)) given $(s)"
                rem_edge!(g, e)
                if !(e in keys(S))
                    S[e] = s
                end
                break
            end
        end
    end

    #  step F3
    Z = unshielded(g, S)
    dg = MetaDiGraph(g)

    for e in edges(dg)
        set_marks!(dg, src(e), dst(e), arrow"o-o")
    end

    for (u, v, w) in Z
        verbose && println("R0 with ($(u), $(v), $(w))")
        if has_edge(dg, (u, v))
            set_marks!(dg, u, v, arrow"*->")
        end
        if has_edge(dg, (v, w))
            set_marks!(dg, v, w, arrow"<-*")
        end
    end

    # main loop for rules R1 to R4
    # loop is repeated until none of them apply
    loop = true
    while loop
        loop = false
        for e in edges(dg)
            (α, β) = Tuple(e)

            for γ in inneighbors(dg, β)
                if γ==α
                    continue
                end

                # R1
                if (!has_edge(dg, α, γ) &&
                    has_marks(dg, α, β, arrow"*->") &&
                    has_marks(dg, β, γ, arrow"o-*"))

                    set_marks!(dg, β, γ, arrow"-->")
                    loop = true
                    verbose && println("R1 with $(α)-$(β)-$(γ)")
                end

                # R2
                if (has_edge(dg, α, γ) &&
                    has_marks(dg, α, γ, arrow"*-o") &&
                    ((has_marks(dg, α, β, arrow"-->") && has_marks(dg, β, γ, arrow"o-*")) ||
                     (has_marks(dg, α, β, arrow"o->") && has_marks(dg, β, γ, arrow"-->"))))

                    set_marks!(dg, α, γ, arrow"*->")
                    loop = true
                    verbose && println("R1 with $(α)-$(β)-$(γ)")
                end

                #R3
                if !isadjacent(dg, α, γ)
                    for θ in inneighbors(dg, γ)
                        if (θ ∈ inneighbors(dg, α) &&
                            θ ∈ inneighbors(dg, β) &&
                            has_marks(dg, α, β, arrow"*->") &&
                            has_marks(dg, β, γ, arrow"<-*") &&
                            has_marks(dg, α, θ, arrow"*-o") &&
                            has_marks(dg, γ, θ, arrow"o-*") &&
                            has_marks(dg, θ, β, arrow"*-o"))
                            set_marks!(dg, θ, β, arrow"*->")
                            loop = true
                            verbose && println("R3 with $(α)-$(β)-$(γ)")
                        end
                    end
                end

            end

            # R4
            for x in vertices(dg)
                paths = yen_k_shortest_paths(g, x, α, LightGraphs.weights(g), 100).paths
                for path in paths
                    if (is_discriminating_path(dg, path) &&
                        has_marks(dg, path[end-1], path[end], arrow"o-*"))
                        if (haskey(S, Edge(path[1], path[end])) &&
                            path[end-1] ∈ S[Edge(path[1], path[end])])
                            set_marks!(dg, path[end-1], path[end], arrow"-->")
                        else
                            set_marks!(dg, path[end-1], path[end], arrow"<->")
                            set_marks!(dg, path[end-2], path[end-1], arrow"<->")
                        end
                        loop=true
                        verbose && println("R4 with $(α) and $(path)")
                    end
                end
            end
        end

    end

    if !augmented
        return dg
    end

    # rules R5 to R10

    loop = true
    while loop
        loop = false
        for e in edges(dg)
            (α, β) = Tuple(e)

            # R5
            if has_marks(dg, α, β, arrow"o-o")
                paths = yen_k_shortest_paths(g, α, β, LightGraphs.weights(g), 100).paths

                for path in paths
                    if (is_uncovered_circle_path(dg, path) &&
                        !isadjacent(dg, path[1], path[end-1]) &&
                        !isadjacent(dg, path[2], path[end]))
                        verbose && println("R5: $(α)-$(β) with $(path)")
                        set_marks!(dg, α, β, arrow"---")
                        for (e1, e2) in zip(path[1:end-1], path[2:end])
                            set_marks!(dg, e1, e2, arrow"---")
                        end
                        loop = true
                        break
                    end
                end
            end

            for γ in inneighbors(dg, β)

                if γ == α
                    continue
                end

                #R6
                if has_marks(dg, α, β, arrow"---") && has_marks(dg, β, γ, arrow"o-*")
                    set_marks!(dg, β, γ, arrow"--*")
                    loop = true
                    verbose && println("R6: $(α)-$(β)-$(γ)")
                end

                # R7
                if(!isadjacent(dg, α, γ) &&
                   has_marks(dg, α, β, arrow"--o") &&
                   has_marks(dg, β, γ, arrow"o-*"))
                    set_marks!(dg, β, γ, arrow"--*")
                    loop = true
                    verbose && println("R7: $(α)-$(β)-$(γ)")
                end

                # R8
                if(has_edge(dg, α, γ) &&
                   has_marks(dg, α, γ, arrow"o->") &&
                   (has_marks(dg, α, β, arrow"-->") || has_marks(dg, α, β, arrow"--o")) &&
                   has_marks(dg, β, γ, arrow"-->"))
                    set_marks!(dg, α, γ, arrow"-->")
                    loop = true
                    verbose && println("R8: $(α)-$(β)-$(γ)")
                end
            end

            # R9
            if has_marks(dg, α, β, arrow"o->")
                paths = yen_k_shortest_paths(g, α, β, LightGraphs.weights(g), 100).paths
                for path in paths
                    if (length(path)>3 &&
                        is_uncovered_PD_path(dg, path) &&
                        !isadjacent(dg, path[2], path[end]))
                        set_marks!(dg, α, β, arrow"-->")
                        loop = true
                        verbose && println("R9: $(α)-$(β) with $(path)")
                        break
                    end
                end
            end

            #R10
            if has_marks(dg, α, β, arrow"o->")
                for (γ,θ) in combinations(inneighbors(dg, β), 2)
                    if(θ==α || γ==α)
                        continue
                    end
                    if(has_marks(dg, γ, β, arrow"-->") &&
                       has_marks(dg, β, θ, arrow"<--"))

                        p1 = yen_k_shortest_paths(g, α, γ, LightGraphs.weights(g), 100).paths
                        p2 = yen_k_shortest_paths(g, α, θ, LightGraphs.weights(g), 100).paths

                        for path1 in p1
                            if is_uncovered_PD_path(dg, path1)
                                for path2 in p2
                                    if is_uncovered_PD_path(dg, path2)
                                        μ = path1[2]
                                        ω = path2[2]
                                        if(μ != ω && !isadjacent(dg, μ, ω))
                                            set_marks!(dg, α, β, arrow"-->")
                                            loop=true
                                            verbose && println("R10: $(α)-$(β) with $(p1) and $(p2)")
                                            break
                                        end
                                    end
                                end
                            end
                        end
                    end
                end
            end
        end

    end
    dg
end

function fcialg(t, p::Float64, test::typeof(gausscitest); kwargs...)
    @assert Tables.istable(t)

    c = Tables.columns(t)
    sch = Tables.schema(t)
    n = length(sch.names)

    X = reduce(hcat, map(c->t[c], 1:n))
    N = size(X,1)
    C = Statistics.cor(X)
    return fcialg(n, gausscitest, (C,N), quantile(Normal(), 1-p/2); kwargs...)
end

function fcialg(t, p::Float64, test::typeof(cmitest); kwargs...)
    @assert Tables.istable(t)
    @assert all(t->t==Float64, Tables.schema(t).types)

    c = Tables.columns(t)
    sch = Tables.schema(t)
    n = length(sch.names)

    return fcialg(n, cmitest, c, p; kwargs...)
end


text:
# ---
# title: 189. Rotate Array
# id: problem189
# author: Indigo
# date: 2021-01-17
# difficulty: Medium
# categories: Array
# link: <https://leetcode.com/problems/rotate-array/description/>
# hidden: true
# ---
#
# Given an array, rotate the array to the right by _k_ steps, where  _k_  is
# non-negative.
#
# **Follow up:**
#
#   * Try to come up as many solutions as you can, there are at least 3 different ways to solve this problem.
#   * Could you do it in-place with O(1) extra space?
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [1,2,3,4,5,6,7], k = 3
#     Output: [5,6,7,1,2,3,4]
#     Explanation:
#     rotate 1 steps to the right: [7,1,2,3,4,5,6]
#     rotate 2 steps to the right: [6,7,1,2,3,4,5]
#     rotate 3 steps to the right: [5,6,7,1,2,3,4]
#
#
# **Example 2:**
#
#
#
#     Input: nums = [-1,-100,3,99], k = 2
#     Output: [3,99,-1,-100]
#     Explanation:
#     rotate 1 steps to the right: [99,-1,-100,3]
#     rotate 2 steps to the right: [3,99,-1,-100]
#
#
#
#
# **Constraints:**
#
#   * `1 <= nums.length <= 2 * 104`
#   * `-231 <= nums[i] <= 231 - 1`
#   * `0 <= k <= 105`
#
#
## @lc code=start
using LeetCode

function rotate_arr(nums::Vector{Int}, k::Int)
    k %= length(nums)
    reverse!(nums)
    reverse!(@view nums[1:k])
    reverse!(@view nums[k+1:end])
    nums
end
## @lc code=end


text:
function nlp_mi_005_010(
    optimizer,
    objective_tol,
    primal_tol,
    dual_tol,
    termination_target = TERMINATION_TARGET_LOCAL,
    primal_target = PRIMAL_TARGET_LOCAL,
)
    # Test Goals:
    # - function /

    model = Model(optimizer)

    @variable(model, x >= 0, Int)
    @variable(model, y >= 0)

    @objective(model, Min, x + y)
    @NLconstraint(model, y >= 1 / (x + 0.1) - 0.5)
    @NLconstraint(model, x >= y^(-2) - 0.5)
    @NLconstraint(model, 4 / (x + y + 0.1) >= 1)

    optimize!(model)

    check_status(model, FEASIBLE_PROBLEM, termination_target, primal_target)
    check_objective(model, 1.8164965727459055, tol = objective_tol)
    return check_solution([x, y], [1, 0.816496581496872], tol = primal_tol)
end


text:
using Colors, FixedPointNumbers, JLD2
using Test
using ColorTypes: eltype_default

@testset "Conversion" begin
    r8(x) = reinterpret(N0f8, x)

    # Promotions
    a, b = promote(RGB(1,0,0), Gray(0.8))
    @test isa(a, RGB{Float64}) && isa(b, RGB{Float64})
    a, b = promote(RGBA(1,0,0), Gray(0.8))
    @test isa(a, RGBA{Float64}) && isa(b, RGBA{Float64})
    a, b = promote(RGBA(1,0,0), GrayA(0.8))
    @test isa(a, RGBA{Float64}) && isa(b, RGBA{Float64})
    a, b = promote(RGB(1,0,0), GrayA(0.8))
    @test isa(a, RGBA{Float64}) && isa(b, RGBA{Float64})
    a, b = promote(RGB(1,0,0), AGray(0.8))
    @test isa(a, ARGB{Float64}) && isa(b, ARGB{Float64})

    # Color parsing
    redN0f8 = parse(Colorant, "red")
    @test colorant"red" == redN0f8
    @test isa(redN0f8, RGB{N0f8})
    @test redN0f8 == RGB(1,0,0)
    redF64 = convert(RGB{Float64}, redN0f8)
    @test parse(RGB{Float64}, "red") === RGB{Float64}(1,0,0)
    @test isa(parse(HSV, "blue"), HSV)
    @test parse(Colorant, "rgb(55,217,127)") === RGB{N0f8}(r8(0x37),r8(0xd9),r8(0x7f))
    @test colorant"rgb(55,217,127)" === RGB{N0f8}(r8(0x37),r8(0xd9),r8(0x7f))
    @test parse(Colorant, "rgba(55,217,127,0.5)") === RGBA{N0f8}(r8(0x37),r8(0xd9),r8(0x7f),0.5)
    @test parse(Colorant, "rgb(55,217,127)") === RGB{N0f8}(r8(0x37),r8(0xd9),r8(0x7f))
    @test parse(Colorant, "rgba(55,217,127,0.5)") === RGBA{N0f8}(r8(0x37),r8(0xd9),r8(0x7f),0.5)
    @test parse(Colorant, "hsl(120, 100%, 50%)") === HSL{Float32}(120,1.0,.5)
    @test colorant"hsl(120, 100%, 50%)" === HSL{Float32}(120,1.0,.5)
    @test parse(RGB{N0f8}, "hsl(120, 100%, 50%)") === convert(RGB{N0f8}, HSL{Float32}(120,1.0,.5))
    @test_throws ErrorException  parse(Colorant, "hsl(120, 100, 50)")
    @test parse(Colorant, "#D0FF58") === RGB(r8(0xD0),r8(0xFF),r8(0x58))
    @test parse(Colorant, "#FB0") === RGB(r8(0xFF),r8(0xBB),r8(0x00))

    @test parse(Colorant, :red) === colorant"red"
    @test parse(Colorant, colorant"red") === colorant"red"

    @test hex(RGB(1,0.5,0)) == "FF8000"
    @test hex(RGBA(1,0.5,0,0.25)) == "40FF8000"

    # srgb_compand / invert_srgb_compand
    @test Colors.srgb_compand(0.5) ≈ 0.7353569830524494 atol=eps()
    @test Colors.invert_srgb_compand(0.7353569830524494) ≈ 0.5 atol=eps()
    # issue #351
    l_pow_x_y() = for i=1:1000; (i/1000)^(1/2.4) end
    l_exp_y_log_x() = for i=1:1000; exp(1/2.4*log(i/1000)) end
    l_pow_x_y(); t_pow_x_y = @elapsed l_pow_x_y()
    l_exp_y_log_x(); t_exp_y_log_x = @elapsed l_exp_y_log_x()
    if t_exp_y_log_x > t_pow_x_y
        @warn "Optimization in `[invert_]srgb_compand()` may have the opposite effect."
    end

    fractional_types = (RGB, BGR, RGB1, RGB4)  # types that support Fractional

    red24 = reinterpret(RGB24, 0x00ff0000)
    red32 = reinterpret(ARGB32, 0xffff0000)
    for T in (Float64, Float32, N0f8)
        c = RGB(one(T), zero(T), zero(T))
        @test eltype(c) == T
        c64 = convert(RGB{Float64}, c)
        @test typeof(c64) == RGB{Float64}
        @test c64 == redF64
        cr = convert(RGB{T}, redF64)
        @test cr == c
    end
    @test RGB(1,0,0) == redF64
    @test RGB(convert(UInt8, 1),0,0) == redF64
    @test RGB(convert(UInt8, 1),convert(UInt8, 0),convert(UInt8, 0)) == redF64
    @test convert(RGB, red24) == redF64

    @test convert(Gray{N0f8}, Gray{N0f8}(0.1)) == Gray{N0f8}(0.1)
    @test convert(Gray{N0f8}, Gray(0.1))     == Gray{N0f8}(0.1)
    @test convert(Gray{N0f8}, Gray24(0.1))   == Gray{N0f8}(0.1)
    @test convert(Gray24, Gray{N0f8}(0.1))   == Gray24(0.1)

    @test convert(RGB{N0f8}, Gray{N0f8}(0.1)) == RGB{N0f8}(0.1,0.1,0.1)
    @test convert(RGB{N0f8}, Gray24(0.1))   == RGB{N0f8}(0.1,0.1,0.1)

    for Cto in ColorTypes.parametric3
        for Cfrom in ColorTypes.parametric3
            for Tto in (Float32, Float64)
                for Tfrom in (Float32, Float64)
                    c = convert(Cfrom{Tfrom}, redF64)
                    @test typeof(c) == Cfrom{Tfrom}
                    c1 = convert(Cto, c)
                    @test eltype(c1) == Tfrom
                    c2 = convert(Cto{Tto}, c)
                    @test typeof(c2) == Cto{Tto}
                end
            end
        end
    end
    for Cto in ColorTypes.parametric3
        @test typeof(convert(Cto, red24)) == Cto{eltype_default(Cto)}
        @test typeof(convert(Cto{Float64}, red24)) == Cto{Float64}
    end

    # Test conversion from Normed types
    for Cto in ColorTypes.parametric3
        for Cfrom in fractional_types
            for Tto in (Float32, Float64)
                for Tfrom in (N0f8, N6f10, N4f12, N2f14, N0f16)
                    c = convert(Cfrom{Tfrom}, redF64)
                    @test typeof(c) == Cfrom{Tfrom}
                    if !(eltype_default(Cto) <: FixedPoint)
                        c1 = convert(Cto, c)
                        @test eltype(c1) == eltype_default(Cto)
                    end
                    c2 = convert(Cto{Tto}, c)
                    @test typeof(c2) == Cto{Tto}
                end
            end
        end
    end

    # Test conversion to Normed types
    for Cto in fractional_types
        for Cfrom in ColorTypes.parametric3
            for Tto in (N0f8, N6f10, N4f12, N2f14, N0f16)
                for Tfrom in (Float32, Float64)
                    c = convert(Cfrom{Tfrom}, redF64)
                    @test typeof(c) == Cfrom{Tfrom}
                    c2 = convert(Cto{Tto}, c)
                    @test typeof(c2) == Cto{Tto}
                end
            end
        end
    end

    ac = RGBA(redF64)

    @test convert(RGB, ac) == RGB(1,0,0)
    @test convert(RGB{N0f8}, ac) == RGB{N0f8}(1,0,0)
    @test convert(RGBA{N0f8}, ac) == RGBA{N0f8}(1,0,0,1)
    @test convert(HSVA, ac) == HSVA{Float64}(convert(HSV, redF64), 1.0)
    @test convert(HSVA{Float32}, ac) == HSVA{Float32}(convert(HSV{Float32}, redF64), 1.0f0)
    @test convert(RGBA, redF64) == ac

    @test convert(ARGB32, ac) == reinterpret(ARGB32, 0xffff0000)
    @test convert(ARGB32, RGB(1,0,0)) == reinterpret(ARGB32, 0xffff0000)
    ac32 = convert(ARGB32, RGB(1,0,0), 0.5)
    @test ac32 == reinterpret(ARGB32, 0x80ff0000)
    @test convert(ARGB32, ac32) === ac32
    @test convert(ARGB32, BGRA(1,0,0,0.5)) == reinterpret(ARGB32, 0x80ff0000)
    @test reinterpret(UInt32, convert(ARGB32, ac)) == 0xffff0000
    @test convert(RGB24, RGB(r8(0xff),r8(0x00),r8(0x00))) == reinterpret(RGB24, 0x00ff0000)
    @test reinterpret(UInt32, convert(RGB24, RGB(r8(0xff),r8(0x00),r8(0x00)))) == 0x00ff0000
    redhsv = convert(HSV, redF64)
    @test convert(RGB24, red24) === red24
    @test convert(RGB24, redhsv) == reinterpret(RGB24, 0x00ff0000)
    @test_throws ArgumentError convert(RGB24,  RGB(0, 1.1, 0))
    @test_throws ArgumentError convert(ARGB32, RGBA(0, 1.1, 0, 0.8))
    @test_throws ArgumentError convert(ARGB32, RGBA(0, 0.8, 0, 1.1))

    @test convert(RGB{N0f8}, red24) == RGB{N0f8}(1,0,0)
    @test convert(RGBA{N0f8}, red32) == RGBA{N0f8}(1,0,0,1)
    @test convert(HSVA{Float64}, red32) == HSVA{Float64}(360, 1, 1, 1)

    @test_throws MethodError AlphaColor(RGB(1,0,0), r8(0xff))

    # whitepoint conversions
    @test isa(convert(XYZ, convert(Lab, redF64), Colors.WP_DEFAULT), XYZ{Float64})
    @test isa(convert(XYZ{Float32}, convert(Lab, redF64), Colors.WP_DEFAULT), XYZ{Float32})
    @test isa(convert(XYZ, convert(Luv, redF64), Colors.WP_DEFAULT), XYZ{Float64})
    @test isa(convert(XYZ{Float32}, convert(Luv, redF64), Colors.WP_DEFAULT), XYZ{Float32})
    @test isa(convert(Lab, convert(XYZ, redF64), Colors.WP_DEFAULT), Lab{Float64})
    @test isa(convert(Lab{Float32}, convert(XYZ, redF64), Colors.WP_DEFAULT), Lab{Float32})
    @test isa(convert(Luv, convert(XYZ, redF64), Colors.WP_DEFAULT), Luv{Float64})
    @test isa(convert(Luv{Float32}, convert(XYZ, redF64), Colors.WP_DEFAULT), Luv{Float32})

    # Test vector space operations
    @test LMS{Float64}(0.125,0.5,0.0)+LMS{Float64}(0.2,0.7,0.4) ≈ LMS{Float64}(0.325,1.2,0.4) atol=91eps()
    @test 3LMS{Float64}(0.125,0.5,0.03) ≈ LMS{Float64}(0.375,1.5,0.09) atol=91eps()

    @test XYZ{Float64}(0.125,0.5,0.0)+XYZ{Float64}(0.2,0.7,0.4) ≈ XYZ{Float64}(0.325,1.2,0.4) atol=91eps()
    @test 3XYZ{Float64}(0.125,0.5,0.03) ≈ XYZ{Float64}(0.375,1.5,0.09) atol=91eps()

    #59
    @test Colors.xyz_to_uv(XYZ{Float64}(0.0, 0.0, 0.0)) === (0.0, 0.0)
    @test Colors.xyz_to_uv(XYZ{Float64}(0.0, 1.0, 0.0)) === (0.0, 0.6)
    @test Colors.xyz_to_uv(XYZ{Float64}(0.0, 1.0, 1.0)) === (0.0, 0.5)
    @test Colors.xyz_to_uv(XYZ{Float64}(1.0, 0.0, 1.0)) === (1.0, 0.0)
    @test Colors.xyz_to_uv(XYZ{Float64}(1.0, 0.0, 0.0)) === (4.0, 0.0)

    # ColorTypes.jl issue #40
    @test convert(HSL, RGB{N0f8}(0.678, 0.847, 0.902)) ≈ HSL{Float32}(194.73685f0,0.5327105f0,0.7901961f0) atol=100eps(Float32)

    # YIQ
    @test convert(YIQ, RGB(1,0,0)) == YIQ{Float32}(0.299, 0.595716, 0.211456)
    @test convert(YIQ, RGB(0,1,0)) == YIQ{Float32}(0.587, -0.274453, -0.522591)
    @test convert(YIQ, RGB(0,0,1)) == YIQ{Float32}(0.114, -0.321263, 0.311135)
    @test convert(RGB, YIQ(1.0,0.0,0.0)) == RGB(1,1,1)
    v = 0.5957
    @test convert(RGB, YIQ(0.0,1.0,0.0)) == RGB(0.9563*v,0,0)
    v = -0.5226
    @test convert(RGB, YIQ(0.0,0.0,-1.0)) == RGB(0,-0.6474*v,0)

    # Gray
    c = Gray{N0f16}(0.8)
    @test convert(RGB, c) == RGB{N0f16}(0.8,0.8,0.8)
    @test convert(RGB{Float32}, c) == RGB{Float32}(0.8,0.8,0.8)

    for C in ColorTypes.parametric3
        c = convert(C, RGB(1,1,1))
        @test gray(convert(Gray, c)) ≈ 1 atol=0.01
        @test gray(convert(Gray{Float64}, c)) ≈ 1 atol=0.01
    end

    # More AbstractRGB
    r4 = RGB4(1,0,0)
    @test convert(RGB, r4) == RGB(1,0,0)
    @test convert(RGB{N0f8}, r4) == RGB{N0f8}(1,0,0)
    @test convert(RGB4{N0f8}, r4) == RGB4{N0f8}(1,0,0)
    @test convert(RGB4{Float32}, r4) == RGB4{Float32}(1,0,0)
    @test convert(BGR{Float32}, r4) == BGR{Float32}(1,0,0)

    # Test accuracy of conversion
    csconv = jldopen(joinpath(dirname(@__FILE__), "test_conversions.jld2")) do file
        read(file, "csconv")
    end

    function convcompare(from, to, eps; showfailure::Bool=false)
        errmax = 0.0
        for i = 1:length(from)
            t = to[i]
            f = convert(typeof(t), from[i])
            diff = abs(comp1(t)-comp1(f)) + abs(comp2(t)-comp2(f)) + abs(comp3(t)-comp3(f))
            mag = abs(comp1(t)+comp1(f)) + abs(comp2(t)+comp2(f)) + abs(comp3(t)+comp3(f))
            if showfailure && diff>eps*mag
                original = from[i]
                @show original f t
            end
            errmax = max(errmax, diff/mag)
        end
        errmax > eps && @warn("Error on conversions from $(eltype(from)) to $(eltype(to)), relative error = $errmax")
        errmax <= eps
    end

    for i = 1:length(csconv)
        f, t = csconv[i]
        if !convcompare(f, t, 1e-3)
            println("  index $i")
        end
    end

    # Images issue #382
    @test convert(Gray, RGBA(1,1,1,1)) == Gray(N0f8(1))

    # https://github.com/timholy/Images.jl/pull/445#issuecomment-189866806
    @test convert(Gray, RGB{N0f8}(0.145,0.145,0.145)) == Gray{N0f8}(0.145)

    # Issue #257
    c = RGB{Float16}(0.9473,0.962,0.9766)
    hsi = convert(HSI, c)
    @test hsi.i > 0.96 && hsi.h ≈ 210
end


text:
export numstring

"""
`numstring(x, n=2)`

Left pads a number (`x`) with zeros to a total of `n` charaters.

## Example
```julia-repl
julia> numstring(11, 2)
"11"

julia> numstring(11, 3)
"011"

julia> numstring(11, 4)
"0011"
```
"""
numstring(x, n=2) = string(10^n + x)[2:end]


text:
abstract type AbstractApproximator end

function construct_layer_array(node_array, act::AbstractVector)
    layer_array = []
    for i in 2:length(node_array)
        node_prev = node_array[i-1]
        node = node_array[i]
        _act = i == length(node_array) ? Flux.identity : act[i-1]
        push!(layer_array, Dense(node_prev, node, _act))
    end
    model = Chain(layer_array...)
    return model
end

"""
    construct_layer_array(node_array, act; act_terminal=Flux.identity)

A convenient way to generate a layer array.

# Example
node_array = [1, 2, 3]
act = Flux.relu
layer_array = PCApprox.construct_layer_array(node_array, act)
model = Chain(layer_array...)
"""
function construct_layer_array(node_array, act; act_terminal=Flux.identity)
    l = length(node_array)
    act_vec = vcat(repeat([act], l-2), [act_terminal])
    construct_layer_array(node_array, act_vec)
end

function number_of_parameters(approximator::AbstractApproximator)
    Flux.params(approximator) |> Map(length) |> sum
end

include("FNN.jl")
include("parametrised_convex_approximators/parametrised_convex_approximators.jl")


text:
using Base.Sys: cpu_info
using Base.Threads
using ArgParse
using Dates
using JLD
using Random
using CellListsBenchmarks

# --- Arguments ---
s = ArgParseSettings()
@add_arg_table s begin
    "benchmark"
        help = "Choose one of : $(keys(benchmark_functions))"
        arg_type = String
        required = true
    "--ns"
        help = "String that is evaluates to `Vector{Int}` type. For example: [1, 10, 100:50:200...]"
        arg_type = String
        required = true
    "-d"
        help = "Dimension. Positive integer."
        arg_type = Int
        required = true
    "-r"
        help = "Radius. Positive float."
        arg_type = Float64
        required = true
    "--seed"
        default = 1
        arg_type = Int
    "--iterations", "-i"
        default = 1
        arg_type = Int
    "--seconds", "-s"
        default = 1.0
        arg_type = Float64
    "--dir"
        default = "output"
        arg_type = AbstractString
    "--float32"
        help = "Use Float32 insteads of Float64."
        action = :store_true
end
args = parse_args(s)

# --- Parameters ---
d = args["d"]
r = args["r"]
ns = Vector{Int}(eval(Meta.parse(args["ns"])))
seed = args["seed"]
iterations = args["iterations"]
seconds = args["seconds"]
directory = args["dir"]
benchmark = args["benchmark"]
float32 = args["float32"]

# --- Validate parameters ---
length(ns) < 1 && throw(DomainError(""))
d < 1 && throw(DomainError(""))
r ≤ 0 && throw(DomainError(""))
if float32
    r = Float32(r)
end

# --- Trials ---
@info versioninfo()
benchmarks = [Benchmark(benchmark, n, d, r, seed, iterations, seconds) for n in ns]

# --- IO ---
if !ispath(directory)
    mkpath(directory)
end

filepath = joinpath(directory, "n$(ns[end])-d$(d)-r$(r).jld")
@info "Saving results to $(filepath)"
JLD.save(filepath, "benchmarks", benchmarks)

max_stars_repo_name: jaantollander/cell-lists-benchmarks

max_issues_repo_name: jaantollander/cell-lists-benchmarks

max_forks_repo_name: jaantollander/cell-lists-benchmarks


text:
"""
    solve_orbital_equation!(Pj, okw, method, tol)

Solve the orbital equation for the `j`th orbital, by solving the
eigenproblem of the [`OrthogonalKrylovWrapper`](@ref) operator `okw`,
and store the result in the radial orbital vector of expansion
coefficients `Pj`. The solution is computed using `method`; valid
choices are

1. `:arnoldi`, which uses ArnoldiMethod.jl. Really small grid
   spacings/large grid extents can be problematic with this method,
   due to high condition numbers.

2. `:lobpcg`, which uses IterativeSolvers.jl. Too coarse grids can
   cause never-ending loops or errors from the Cholesky decomposition,
   complaining about non-positive definite subproblems.

Both methods are controlled by the stopping tolerance `tol`.
"""
function solve_orbital_equation!(Pj, okw, method, tol;
                                 facttol=√(eps(real(eltype(Pj)))),
                                 io=stdout, verbosity=0,
                                 kwargs...)
    verbosity > 0 && println(io, "Improving orbital using $method")

    update!(okw)

    ϕ = if method == :arnoldi
        # It would be preferrable if the Arnoldi state could
        # be preserved between iterations, pending
        # https://github.com/haampie/ArnoldiMethod.jl/issues/91
        schur,history = partialschur(okw, nev=1, tol=tol, which=SR())
        verbosity > 2 && println(io,"Orbital improvement: ", history)

        view(schur.Q, :, 1)
    elseif method==:lobpcg
        res = lobpcg(okw, okw.S, false, reshape(Pj, :, 1), 1,
                     tol=tol, C=okw.C)
        verbosity > 2 && show(io, res)

        view(res.X, :, 1)
    else
        throw(ArgumentError("Unknown diagonalization method $(method)"))
    end

    copyto!(Pj, ϕ)
end


text:
using Distributed

addprocs(2)

@everywhere using MonteCarlo, Distributed

@everywhere function simulate(mu)
    # Connect other workers
    pool = filter(x -> x != myid(), workers())
    println(myid(), " => ", pool)
    connect(pool, wait = true)

    println("Creating model")
    model = HubbardModelAttractive(4, 2, U = 1.0, mu = mu)

    println("Creating scheduler")
    # scheduler = SimpleScheduler(
    #     LocalSweep(3), GlobalFlip(),
    #     LocalSweep(3), GlobalShuffle(),
    #     LocalSweep(3), ReplicaExchange(pool[1]),
    #     LocalSweep(3), GlobalFlip(),
    #     LocalSweep(3), GlobalShuffle(),
    #     LocalSweep(3), ReplicaPull()
    # )
    scheduler = AdaptiveScheduler(
        (
            LocalSweep(3), Adaptive(),
            LocalSweep(3), ReplicaPull(),
            LocalSweep(3), Adaptive(),
            LocalSweep(3), ReplicaExchange(pool[1])
        ),
        (GlobalFlip(), GlobalShuffle())
    )

    println("Creating dqmc")
    dqmc = DQMC(
        model,
        beta=5.0, recorder=Discarder, scheduler = scheduler,
        thermalization=10000, sweeps=10000
    )

    dqmc[:G]    = greens_measurement(dqmc, model)
    dqmc[:E]    = total_energy(dqmc, model)
    dqmc[:Occs] = occupation(dqmc, model)
    dqmc[:CDC]  = charge_density_correlation(dqmc, model)
    dqmc[:Mx]   = magnetization(dqmc, model, :x)
    dqmc[:My]   = magnetization(dqmc, model, :y)
    dqmc[:Mz]   = magnetization(dqmc, model, :z)
    dqmc[:SDCx] = spin_density_correlation(dqmc, model, :x)
    dqmc[:SDCy] = spin_density_correlation(dqmc, model, :y)
    dqmc[:SDCz] = spin_density_correlation(dqmc, model, :z)
    dqmc[:PC]   = pairing_correlation(dqmc, model, K=4)

    println("Running Simulation")
    run!(dqmc, verbose=!true)

    dqmc
end

begin
    @everywhere MonteCarlo.enable_benchmarks()
    @time dqmcs = pmap(simulate, [1.0, 0.0])
    @everywhere MonteCarlo.reset_timer!()
    @time dqmcs = pmap(simulate, [1.0, 0.0])
    MonteCarlo.show_statistics(stdout, dqmcs[1].scheduler)
    MonteCarlo.show_statistics(stdout, dqmcs[2].scheduler)
end

TOs = pmap(_ -> MonteCarlo.TimerOutputs.DEFAULT_TIMER, workers());


text:
# ---
# title: 458. Poor Pigs
# id: problem458
# author: Indigo
# date: 2021-02-02
# difficulty: Hard
# categories: Math
# link: <https://leetcode.com/problems/poor-pigs/description/>
# hidden: true
# ---
#
# There are `buckets` buckets of liquid, where **exactly one** of the buckets is
# poisonous. To figure out which one is poisonous, you feed some number of
# (poor) pigs the liquid to see whether they will die or not. Unfortunately, you
# only have `minutesToTest` minutes to determine which bucket is poisonous.
#
# You can feed the pigs according to these steps:
#
#   1. Choose some live pigs to feed.
#   2. For each pig, choose which buckets to feed it. The pig will consume all the chosen buckets simultaneously and will take no time.
#   3. Wait for `minutesToDie` minutes. You may **not** feed any other pigs during this time.
#   4. After `minutesToDie` minutes have passed, any pigs that have been fed the poisonous bucket will die, and all others will survive.
#   5. Repeat this process until you run out of time.
#
# Given `buckets`, `minutesToDie`, and `minutesToTest`, return _the **minimum**
# number of pigs needed to figure out which bucket is poisonous within the
# allotted time_.
#
#
#
# **Example 1:**
#
#
#
#     Input: buckets = 1000, minutesToDie = 15, minutesToTest = 60
#     Output: 5
#
#
# **Example 2:**
#
#
#
#     Input: buckets = 4, minutesToDie = 15, minutesToTest = 15
#     Output: 2
#
#
# **Example 3:**
#
#
#
#     Input: buckets = 4, minutesToDie = 15, minutesToTest = 30
#     Output: 2
#
#
#
#
# **Constraints:**
#
#   * `1 <= buckets <= 1000`
#   * `1 <= minutesToDie <= minutesToTest <= 100`
#
#
## @lc code=start
using LeetCode

function poor_pigs(buckets::Int, minutesToDie::Int, minutesToTest::Int)
    state = minutesToTest ÷ minutesToDie + 1
    return ceil(Int, log(buckets) / log(state))
end
## @lc code=end


text:
using SmartShockFinder

datatyp = "NNRH"
n_dims = 1
η = 0.002                # learning rate
β = 0.009                # regularization paramater
number_epochs = 100
Sb = 500                 # batch size
L = 10                   # early stopping parameter
hidden_units_1 = 10
hidden_units_2 = 10
hidden_units_3 = 10
hidden_units_4 = 10
hidden_units_5 = 10

traindata = joinpath(@__DIR__, "traindata1dNNRH.h5")
download("https://gist.github.com/JuliaOd/f4d10153e539c026daf7076f087dc937/raw/traindata1dNNRH.h5", traindata)
validdata = joinpath(@__DIR__, "validdata1dNNRH.h5")
download("https://gist.github.com/JuliaOd/f4d10153e539c026daf7076f087dc937/raw/validdata1dNNRH.h5", validdata)

train_network1d(datatyp, n_dims, η, β, number_epochs, Sb, L,
                hidden_units_1, hidden_units_2, hidden_units_3,
                hidden_units_4, hidden_units_5)


text:
# get hold of all the predefined math constants in Julia
import Base.MathConstants: π, pi, ℯ, e, γ, eulergamma, catalan, φ, golden
Base.@irrational twoπ  6.28318530717958647692528676655900576839433879  big(2.0)*pi


text:
# -------------------------------------------------------------------------------
# Progress bar and State
# -------------------------------------------------------------------------------
mutable struct ProgressBar{P<:ProgressMeter.Progress, T<:Real}
    bar::P
    L0::T
    L::T
    L_primal::T
end

struct State{S, D<:Dict, T, A}
    seed::S
    dict::D
    time_init::T
    coordinates::A
end

function ProgStateInit(solver::S,
                       model::M,
                       data::D,
                       scores::AbstractVector;
                       kwargs...) where {S<:AbstractSolver, M<:AbstractModel, D<:AbstractData}

    msg  = "$(nameof(M)) $(nameof(D)) loss - $(nameof(S)) solver: "
    bar  = ProgressMeter.Progress(solver.maxiter, 0.1, msg)
    L    = objective(model, data, values(kwargs)..., scores)
    vals = (values(kwargs)..., time = 0, L = L)
    L_primal = Inf
    if isa(data, Dual)
        L_primal = primal_objective(model, data, values(kwargs)..., scores)
        vals = (; vals..., L_primal = L_primal)
    end
    dict = Dict{Union{Symbol, Int64}, Any}(:initial => vals)

    if S <: Coordinate
        coordinates = Array{Int64}(undef, solver.maxiter, 2)
    else
        coordinates = Int64[]
    end
    return State(solver.seed, dict, time(), coordinates), ProgressBar(bar, L, L, L_primal)
end

ProgStateInit(solver::General, model::AbstractModel, tm::Real; kwargs...) =
    State(solver.seed, Dict{Union{Symbol, Int64}, Any}(:optimal => (values(kwargs)..., time = tm)), time(), Int64[])

function update!(state::State,
                 progress::ProgressBar,
                 solver::S,
                 model::AbstractModel,
                 data::AbstractData,
                 iter::Integer,
                 scores::AbstractVector;
                 k::Int = 0,
                 l::Int = 0,
                 gap::Real = Inf,
                 kwargs...) where {S<:AbstractSolver}

    vars = values(kwargs)
    condition_1 = iter == solver.maxiter
    condition_2 = iter in solver.iters
    condition_3 = mod(iter, ceil(Int, solver.maxiter/10)) == 0 && solver.verbose

    if S <: Coordinate
        state.coordinates[iter, 1] = k
        state.coordinates[iter, 2] = l
    end

    if condition_1 || condition_2 || condition_3
        L = objective(model, data, vars..., scores)
        progress.L = L
        if isa(data, Dual)
            L_primal = primal_objective(model, data, vars..., scores)
            progress.L_primal = L_primal
        end
    end

    if solver.verbose
        vals = [(:L0, progress.L0), (:L, progress.L)]
        if isa(data, Dual)
            push!(vals, (:L_primal, progress.L_primal))
            push!(vals, (:gap, dualitygap(progress)))
        end
        ProgressMeter.next!(progress.bar; showvalues = vals)
    end

    if condition_2
        key = iter
        if isa(data, Dual)
            state.dict[key] = (vars..., time = time() - state.time_init, L = L, L_primal = L_primal)
        else
            state.dict[key] = (vars..., time = time() - state.time_init, L = L)
        end
    end
    if condition_1
        key = :optimal
        if isa(data, Dual)
            state.dict[key] = (vars..., time = time() - state.time_init, L = L, L_primal = L_primal)
        else
            state.dict[key] = (vars..., time = time() - state.time_init, L = L)
        end
    end
end

dualitygap(progress::ProgressBar) = progress.L_primal - progress.L

# -------------------------------------------------------------------------------
# Gradient descent utilities
# -------------------------------------------------------------------------------
function minimize!(solver::AbstractSolver, x, Δ)
    Optimise.apply!(solver.optimizer, x, Δ)
    x .-= Δ
end

function maximize!(solver::AbstractSolver, x, Δ)
    Optimise.apply!(solver.optimizer, x, Δ)
    x .+= Δ
end

# -------------------------------------------------------------------------------
# Coordinate descent utilities
# -------------------------------------------------------------------------------
solution(a::Real, b::Real, Δlb::Real, Δub::Real) = min(max(Δlb, - b/a), Δub)

function update!(best::BestUpdate, k, l, Δ, L, vars)
    if L >= best.L
        best.k = k
        best.l = l
        best.Δ = Δ
        best.L = L
        best.vars = vars
    end
end

function select_rule(model::AbstractModel, data::Dual{<:DTrain}, k, args...)
    best = BestUpdate(1, 2, 0.0, -Inf, (αk = 0.0, αl = 0.0))

    for l in 1:data.n
        l == k && continue

        if k <= data.nα && l <= data.nα
            rule_αα!(model, data, best, k, l, args...)
        elseif k <= data.nα && l > data.nα
            rule_αβ!(model, data, best, k, l, args...)
        elseif k > data.nα && l <= data.nα
            rule_αβ!(model, data, best, l, k, args...)
        else
            rule_ββ!(model, data, best, k, l, args...)
        end
    end
    return best
end

function scores!(data::Dual{<:DTrain}, best::BestUpdate, s)
    if best.k <= data.nα && best.l > data.nα
        s .+= best.Δ*(data.K[:, best.k] + data.K[:, best.l])
    else
        s .+= best.Δ*(data.K[:, best.k] - data.K[:, best.l])
    end
end

find_βmax(βsort, β, k) = βsort[1] != β[k] ? βsort[1] : βsort[2]

function find_βmax(βsort, β, k, l)
    if βsort[1] ∉ [β[k], β[l]]
        return βsort[1]
    elseif βsort[2] ∉ [β[k], β[l]]
        return βsort[2]
    else
        return βsort[3]
    end
end

function βsorted!(data::Dual{<:DTrain}, best::BestUpdate, β, βsort)
    if haskey(best.vars, :βk)
        deleteat!(βsort, searchsortedfirst(βsort, β[best.k - data.nα]; rev = true))
        insert!(βsort, searchsortedfirst(βsort, best.vars.βk; rev = true), best.vars.βk)
    end
    if haskey(best.vars, :βl)
        deleteat!(βsort, searchsortedfirst(βsort, β[best.l - data.nα]; rev = true))
        insert!(βsort, searchsortedfirst(βsort, best.vars.βl; rev = true), best.vars.βl)
    end
end

# -------------------------------------------------------------------------------
# Objective from named tuples
# -------------------------------------------------------------------------------
objective(model::AbstractModel, data::Primal, solution::NamedTuple) =
    objective(model, data, solution.w, solution.t)

objective(model::AbstractPatMat, data::Dual{<:DTrain}, solution::NamedTuple) =
    objective(model, data, solution.α, solution.β, solution.δ)

objective(model::AbstractTopPushK, data::Dual{<:DTrain}, solution::NamedTuple) =
    objective(model, data, solution.α, solution.β)

# -------------------------------------------------------------------------------
# Exact thresholds
# -------------------------------------------------------------------------------
exact_threshold(model::PatMat, data::Primal, s) =
    any(isnan.(s)) ? NaN : quantile(s, 1 - model.τ)

exact_threshold(model::PatMatNP, data::Primal, s) =
    any(isnan.(s)) ? NaN : quantile(s[data.ind_neg], 1 - model.τ)

function exact_threshold(model::AbstractTopPushK, data::Primal, s)
    K = getK(model, data)
    if K >= length(data.type.ind_neg)
        return mean(s[data.ind_neg])
    else
        return mean(partialsort(s[data.ind_neg], 1:K, rev = true))
    end
end

exact_threshold(model::TopPush, data::Primal, s) =
    maximum(s[data.ind_neg])

exact_threshold(model::AbstractModel, data::Primal, w, T) =
    exact_threshold(model, data, scores(model, data, w))

exact_threshold(model::PatMat, data::Dual{<:Union{DTrain, DValidation}}, s) =
    any(isnan.(s)) ? NaN : quantile(s, 1 - model.τ)

exact_threshold(model::PatMatNP, data::Dual{<:Union{DTrain, DValidation}}, s) =
    any(isnan.(s)) ? NaN : quantile(s[data.type.ind_neg], 1 - model.τ)

function exact_threshold(model::AbstractTopPushK, data::Dual{<:Union{DTrain, DValidation}}, s)
    K = getK(model, data)
    if K >= length(data.type.ind_neg)
        return mean(s[data.type.ind_neg])
    else
        return mean(partialsort(s[data.type.ind_neg], 1:K, rev = true))
    end
end

exact_threshold(model::TopPush, data::Dual{<:Union{DTrain, DValidation}}, s) =
    maximum(s[data.type.ind_neg])

exact_threshold(model::AbstractModel, data::Dual{<:Union{DTrain, DValidation}}, α, β) =
    exact_threshold(model, data, scores(model, data, α, β))


text:
@testset "Testing the BBox structure" begin
    X = [
        0.512363  0.523803  0.610511;
        0.207648  0.54746   0.21826;
        0.472262  0.919383  0.313127;
        0.812478  0.223824  0.593296;
        0.109389  0.569246  0.625016
    ]
    Y = [
        0.551976   0.105557   0.733297
        0.0419843  0.0205705  0.545029
        0.237319   0.844252   0.217371
        0.759811   0.638246   0.404687
    ]

    bx = BBox(X)
    by = BBox(Y)

    # constructor
    @test_throws DomainError BBox([1 2 3 4])

    # minimum
    @test get_min(bx) == [0.109389, 0.223824, 0.21826]
    @test get_min(by) == [0.0419843, 0.0205705, 0.217371]

    # maximum
    @test get_max(bx) == [0.812478, 0.919383, 0.625016]
    @test get_max(by) == [0.759811, 0.844252, 0.733297]

    # center
    @test get_ctr(bx) ≈ [0.4609335, 0.5716035, 0.421638]
    @test get_ctr(by) ≈ [0.40089765, 0.43241125, 0.475334]

    # edge
    @test get_edg(bx) ≈ [0.703089, 0.695559, 0.406756]
    @test get_edg(by) ≈ [0.7178267, 0.8236815, 0.515926]
end

max_stars_repo_name: marcbakry/HierarchicalMatrices

max_issues_repo_name: marcbakry/HierarchicalMatrices

max_forks_repo_name: marcbakry/HierarchicalMatrices


text:
#=
Copyright (c) 2015, Intel Corporation
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
- Redistributions of source code must retain the above copyright notice,
  this list of conditions and the following disclaimer.
- Redistributions in binary form must reproduce the above copyright notice,
  this list of conditions and the following disclaimer in the documentation
  and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
THE POSSIBILITY OF SUCH DAMAGE.
=#

module BroadcastTest
using ParallelAccelerator
using CompilerTools

#ParallelAccelerator.DomainIR.set_debug_level(3)
#ParallelAccelerator.ParallelIR.set_debug_level(4)
#ParallelAccelerator.CGen.set_debug_level(3)
#ParallelAccelerator.set_debug_level(3)
#CompilerTools.LambdaHandling.set_debug_level(3)

@acc function twoargs(A, B)
    broadcast(*, A, B)
end

@acc function threeargs(A, B, C)
    broadcast(*, A, B, C)
end

function test1()
    A = [1, 2, 3]
    B = [4, 5, 6]
    twoargs(A, B) == @noacc twoargs(A, B)
end

function test2()
    A = [1, 2, 3]
    B = 3
    twoargs(A, B) == @noacc twoargs(A, B)
end

function test3()
    A = [1, 2]
    B = [1 2 3; 4 5 6]
    C = twoargs(A, B)
    println("test3 C = ", C)
    C == @noacc twoargs(A, B)
end

function test4()
    A = [1 2 3]
    B = [1 2 3; 4 5 6]
    twoargs(A, B) == @noacc twoargs(A, B)
end

function test5()
    A = [1 2 3]
    B = [1 2 3; 4 5 6]
    C = 1
    threeargs(A, B, C) == @noacc threeargs(A, B, C)
end

function test6()
    A = [1 2 3]
    B = [1 2 3; 4 5 6]
    C = [1, 2]
    threeargs(A, B, C) == @noacc threeargs(A, B, C)
end

end

using Base.Test

println("Testing broadcast...")
@test BroadcastTest.test1()
@test BroadcastTest.test2()
@test BroadcastTest.test3()
@test BroadcastTest.test4()
#@test BroadcastTest.test5()
#@test BroadcastTest.test6()
println("Done testing broadcast.")


text:
using Coconet
using Flux
using BSON
using Plots
using MusicTransformer
using Dates

function train_model(;
    model = nothing,
    opt = nothing,
    losses = nothing,

    head_count = 8,
    latent_size = 32,
    word_size = 3,
    N = 6,

    sequence_length = 32,
    batch_size = 5,
    epochs = 1,
    set_length = nothing,
    max_regions = nothing,

    scramble_rate = 0.1,
    scramble_increase = 0.1,
    scramble_max = 0.9,
    graduation_rate = 0.5,

    dataset_name = "Tonnetz_Bach",
)

    difficulty_level = 1
    difficulties = [scramble_rate:scramble_increase:scramble_max...]

    pos_data = nothing

    ## Get model
    epoch_start = 1

    if model == nothing
        @info "Building model..."
        model = Transformer(N, head_count, word_size, latent_size)
        opt = Flux.ADAM()
    end

    model = model |> gpu
    ps = params(model)

    ## Data
    @info "Wrangling data"
    path = "datasets"
    dataset = BSON.load(joinpath(path, string(dataset_name, ".bson")))

    if set_length != nothing
        dataset[:train] = dataset[:train][1:set_length]
        dataset[:test] = dataset[:test][1:set_length]
    end

    ## Loss
    function loss(query_data, mask, y)
        l = y .- model(query_data, pos_data, mask=mask)
        l = l .* l
        sum(l)
    end

    function test_loss(set::Symbol)
        L_test = 0 # Total loss
        i = 0 # counter for the progress bar

        test = Coconet.get_batch(
            dataset,
            sequence_length,
            batch_size,
            set = set,
            max_regions = max_regions,
        )

        for batch in test
            L_batch = 0 # Loss per batch
            for song in batch
                L_song = 0 # loss per song
                song = gpu(song)
                for region in song
                    # data
                    mask = gpu(MusicTransformer.get_mask(region, scramble_rate))

                    query_data = region ./ (mask[1,1,:,:,:] .+ 1)
                    key_data = query_data
                    y = region
                    pos_data = gpu(get_positional_encoding(region))

                    # Loss per region
                    L_song += loss(query_data, mask, y)
                end

                # Average of song losses across batch
                L_batch += L_song

                # Progress bar
                if i % 4 == 0
                    if i % 16 == 0
                        print(' ')
                    end
                    print('/')
                end
                i += 1
            end

            # Losses across test set
            L_test += L_batch
        end

        print(string("> ", L_test, "\n"))
        return L_test
    end

    function ratio()
        @info string("Current Ratio")
        print(string("Train: ", losses[length(losses)][1]/losses[1][1], "\n"))
        print(string("Test: ", losses[length(losses)][2]/losses[1][2], "\n"))
    end

    ## Plot
    function save_to_plot(epoch)
        train = [loss[1] for loss in losses] ./ losses[1][1]
        test = [loss[2] for loss in losses] ./ losses[1][2]
        x = 0:1:(length(losses)-1)

        p = plot(x, train, label="Train")
        plot!(x, test, label="Test")

        ylims!(0, max(train..., test..., 1))
        xlims!(0, length(losses) - 1)

        xlabel!("Epochs")
        p = ylabel!("percent of loss_0")

        p = plot(
            p,
            size = (1280,720),
            xticks = x,
            yticks = 0:0.1:1,
        )

        savefig(p, string("checkpoint/loss_curve"))
        return p
    end

    # Get initial loss
    @info string("Initial Loss:\n")
    if losses == nothing
        print("Train: ")
        tr_l = test_loss(:train)

        print("Test: ")
        te_l = test_loss(:test)

        losses = [(tr_l, te_l)]
    else
        print(string("Train: ", losses[1][1],"\n"))
        print(string("Test: ", losses[1][2],"\n"))
    end
    ratio()
    best_loss = losses[length(losses)][2]
    last_improvement = best_loss / losses[1][2]

    # Train!
    for epoch = epoch_start:(epoch_start+epochs - 1)
        i = 0 # counter for the progress bar
        @info string("Epoch ", epoch)
        train = Coconet.get_batch(
            dataset,
            sequence_length,
            batch_size,
            set = :train,
            max_regions = max_regions,
        )

        L_train = 0
        batch_num = 1

        # Epoch
        for batch in train
            print(string("Batch ", batch_num," of ",length(train)))

            ∇_batch = [] # Gradient per batch
            L_batch = 0 # Loss per batch
            for song in batch
                ∇_song = [] # Gradient per song
                L_song = 0 # Loss per song
                song = gpu(song)
                for region in song
                    # data
                    mask = gpu(MusicTransformer.get_mask(region, scramble_rate))

                    query_data = region ./ (mask[1,1,:,:,:] .+ 1)
                    key_data = query_data
                    y = region
                    pos_data = gpu(get_positional_encoding(region))

                    L_region = 0

                    #time
                    t = now()

                    # Calculate gradient over region
                    ∇_region = gradient(ps) do
                        L_region = loss(query_data, mask, y)
                    end

                    #time
                    print(string("\n", Dates.canonicalize(now() - t), "\n"))

                    push!(∇_song, cpu(∇_region))
                    L_song += L_region
                end

                # Store song gradients
                push!(∇_batch, ∇_song)
                L_batch += L_song

                # Progress bar
                print('.')
            end
            print("\n")

            ∇_batch = gpu(∇_batch)
            # Combine gradients
            ∇_batch = sum([sum(∇_song) for ∇_song in ∇_batch])
            L_train += L_batch

            # Apply gradients (a.k.a. Train!)
            Flux.update!(opt, ps, ∇_batch)
            batch_num += 1
        end

        # Get loss
        L_test = test_loss(:test)
        push!(losses, (L_train, L_test))

        # Print current state
        ratio()

        #FUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUU
        if isnan(L_test)
            print("FUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUU")
            break
        end

        # Make it harder if it's doing too well
        train_ratio = losses[length(losses)][2] / losses[1][2]
        improvement = 1 - (train_ratio / last_improvement)

        if epoch > 0 && epoch % 2 == 0
            if difficulty_level < length(difficulties)
                difficulty_level += 1
                last_improvement = train_ratio
                @info string("Increasing scramble rate to ", difficulties[difficulty_level])
            end
        elseif improvement < 0
            last_improvment = train_ratio
            @warn string("Looks like they're struggling...")
        end

        # Save checkpoint
        if losses[length(losses)][2] < best_loss
            @info "Saving model"
            best_loss = losses[length(losses)][2]
            m = model |> cpu
            o = opt |> cpu
            BSON.@save string("checkpoint/my_model.bson") model=m opt=o epoch_start=epoch+1 losses=losses scramble_rate
        end

        # Plot losses
        save_to_plot(epoch)
    end
end

max_stars_repo_name: MacKenzieHnC/Coconet

max_issues_repo_name: MacKenzieHnC/Coconet

max_forks_repo_name: MacKenzieHnC/Coconet