text:
# From https://discourse.julialang.org/t/aliases-for-union-t-nothing-and-union-t-missing/15402/4
const Maybe{T} = Union{T,Nothing}

# The following default values are picked from `<QE source>/Modules/read_namelists.f90`
"""
    ControlNamelist <: Namelist

Represent the `CONTROL` namelist of `cp.x`.
"""
@with_kw struct ControlNamelist <: Namelist
    calculation::String = "cp"
    title::String = "MD Simulation"
    verbosity::String = "low"
    isave::Int = 100
    restart_mode::String = "restart"
    nstep::Int = 50
    iprint::Int = 10
    tstress::Bool = false
    tprnfor::Bool = false
    dt::Float64 = 1.0
    outdir::String = "./"
    saverho::Bool = true
    prefix::String = "cp"
    ndr::Int = 50
    ndw::Int = 50
    tabps::Bool = false
    max_seconds::Float64 = 1e7
    etot_conv_thr::Float64 = 1e-4
    forc_conv_thr::Float64 = 1e-3
    ekin_conv_thr::Float64 = 1e-6
    disk_io::String = "default"
    memory::String = "default"
    pseudo_dir::String = raw"$HOME/espresso/pseudo/"
    tefield::Bool = false
    # These checks are from https://github.com/QEF/q-e/blob/4132a64/Modules/read_namelists.f90#L1282-L1369.
    @assert(
        calculation in
        ("cp", "scf", "nscf", "relax", "vc-relax", "vc-cp", "cp-wf", "vc-cp-wf")
    )
    @assert verbosity in ("high", "low", "debug", "medium", "default", "minimal")
    @assert isave >= 1
    @assert nstep >= 0
    @assert iprint >= 1
    @assert dt >= 0
    @assert ndr >= 50
    @assert ndw <= 0 || ndw >= 50
    @assert max_seconds >= 0
    @assert etot_conv_thr >= 0
    @assert forc_conv_thr >= 0
    @assert ekin_conv_thr >= 0
    @assert memory in ("small", "default", "large")
end # struct ControlNamelist

"""
    SystemNamelist <: Namelist

Represent the `SYSTEM` namelist of `cp.x`.
"""
@with_kw struct SystemNamelist <: Namelist
    ibrav::Int = -1
    celldm::Vector{Maybe{Float64}} = zeros(6)  # Must specify
    A::Float64 = 0.0
    B::Float64 = 0.0
    C::Float64 = 0.0
    cosAB::Float64 = 0.0
    cosAC::Float64 = 0.0
    cosBC::Float64 = 0.0
    nat::Int = 0
    ntyp::Int = 0
    nbnd::Int = 0
    tot_charge::Float64 = 0.0
    tot_magnetization::Float64 = -1.0
    ecutwfc::Float64 = 0.0
    ecutrho::Float64 = 0.0
    nr1::Int = 0
    nr2::Int = 0
    nr3::Int = 0
    nr1s::Int = 0
    nr2s::Int = 0
    nr3s::Int = 0
    nr1b::Int = 0
    nr2b::Int = 0
    nr3b::Int = 0
    occupations::String = "fixed"
    degauss::Float64 = 0.0
    smearing::String = "gaussian"
    nspin::Int = 1
    ecfixed::Float64 = 0.0
    qcutz::Float64 = 0.0
    q2sigma::Float64 = 0.1  # The default value in QE's source code is 0.01
    input_dft::String = "none"
    exx_fraction::Float64 = 0.25
    lda_plus_u::Bool = false
    Hubbard_U::Vector{Union{Nothing,Float64}} = []
    vdw_corr::String = "none"
    london_s6::Float64 = 0.75
    london_rcut::Float64 = 200.0
    ts_vdw::Bool = false
    ts_vdw_econv_thr::Float64 = 1e-6
    ts_vdw_isolated::Bool = false
    assume_isolated::String = "none"
    # These checks are from https://github.com/QEF/q-e/blob/4132a64/Modules/read_namelists.f90#L1378-L1499.
    @assert ibrav in union(-1:1:14, (-3, -5, -9, -12, -13))
    @assert(
        if ibrav in -1:0
            true  # Skip the check, cannot use `nothing`
        else
            celldm[1] != 0 || A != 0  # Cannot use `iszero` to compare now!
        end,
        "invalid lattice parameters (`celldm` $celldm or `A` $A)!"
    )
    @assert(
        if ibrav == 14
            length(celldm) == 6
        elseif ibrav in (5, -5, 12, 13)
            4 <= length(celldm) <= 6
        elseif ibrav in (4, 6, 7, 8, 9, -9, 10, 11)
            3 <= length(celldm) <= 6
        elseif ibrav == -13  # `-13` is new from QE 6.4
            5 <= length(celldm) <= 6
        else
            1 <= length(celldm) <= 6
        end,
        "`celldm` must have length between 1 to 6! See `ibrav`'s doc!"
    )
    @assert nat >= 0
    @assert(0 <= ntyp <= 10, "`ntyp` $ntyp is either less than zero or too large!")
    @assert nspin in (1, 2, 4)
    @assert ecutwfc >= 0
    @assert ecutrho >= 0
    @assert(iszero(degauss), "`degauss` is not used in CP!")
    @assert ecfixed >= 0
    @assert qcutz >= 0
    @assert q2sigma >= 0
end # struct SystemNamelist

"""
    ElectronsNamelist <: Namelist

Represent the `ELECTRONS` namelist of `cp.x`.
"""
@with_kw struct ElectronsNamelist <: Namelist
    electron_maxstep::Int = 100
    electron_dynamics::String = "none"
    conv_thr::Float64 = 1e-6
    niter_cg_restart::Int = 20
    efield::Float64 = 0.0
    epol::Int = 3
    emass::Float64 = 400.0
    emass_cutoff::Float64 = 2.5
    orthogonalization::String = "ortho"
    ortho_eps::Float64 = 1e-8
    ortho_max::Int = 20
    ortho_para::Int = 0
    electron_damping::Float64 = 0.1
    electron_velocities::String = "default"
    electron_temperature::String = "not_controlled"
    ekincw::Float64 = 0.001
    fnosee::Float64 = 1.0
    startingwfc::String = "random"
    tcg::Bool = false
    maxiter::Int = 100
    passop::Float64 = 0.3
    n_inner::Int = 2
    ninter_cold_restart::Int = 1
    lambda_cold::Float64 = 0.03
    grease::Float64 = 1.0
    ampre::Float64 = 0.0
    # These checks are from https://github.com/QEF/q-e/blob/4132a64/Modules/read_namelists.f90#L1508-L1543.
    @assert electron_dynamics in ("none", "sd", "damp", "verlet", "cg")  # Different from code
    @assert emass > 0
    @assert emass_cutoff > 0
    @assert orthogonalization in ("ortho", "Gram-Schmidt")
    @assert ortho_eps > 0
    @assert ortho_max >= 1
    @assert electron_velocities in ("zero", "default", "change_step")  # New in 6.4
    @assert electron_temperature in ("nose", "rescaling", "not_controlled")
    @assert ekincw > 0
    @assert fnosee > 0
    @assert startingwfc in ("atomic", "random")
end # struct ElectronsNamelist

"""
    IonsNamelist <: Namelist

Represent the `IONS` namelist of `cp.x`.

Input this namelist only if `calculation` is `"cp"`, `"relax"`, `"vc-relax"`, `"vc-cp"`, `"cp-wf"`, or `"vc-cp-wf"`.
"""
@with_kw struct IonsNamelist <: Namelist
    ion_dynamics::String = "none"
    ion_positions::String = "default"
    ion_velocities::String = "default"
    ion_damping::Float64 = 0.2
    ion_radius::Vector{Union{Nothing,Float64}} = [0.5]  # The default value in QE's source code is just one 0.5
    iesr::Int = 1
    ion_nstepe::Int = 1
    remove_rigid_rot::Bool = false
    ion_temperature::String = "not_controlled"
    tempw::Float64 = 300.0
    fnosep::Float64 = 1.0
    tolp::Float64 = 100.0
    nhpcl::Int = 1
    nhptyp::Int = 0
    nhgrp::Vector{Union{Nothing,Int}} = [0]
    fnhscl::Vector{Union{Nothing,Float64}} = zeros(1)
    ndega::Int = 0
    tranp::Vector{Union{Nothing,Bool}} = falses(1)
    amprp::Vector{Union{Nothing,Float64}} = zeros(1)
    greasp::Float64 = 1.0
    # These checks are from https://github.com/QEF/q-e/blob/4132a64/Modules/read_namelists.f90#L1552-L1585.
    @assert ion_dynamics in ("none", "sd", "cg", "damp", "verlet")
    @assert ion_positions in ("default", "from_input")
    @assert ion_velocities in ("default", "change_step", "random", "from_input", "zero")
    @assert ion_nstepe > 0
    @assert ion_temperature in ("nose", "rescaling", "not_controlled")
    @assert tempw > 0
    @assert fnosep > 0
    @assert 0 <= nhpcl <= 4
end # struct IonsNamelist

"""
    CellNamelist <: Namelist

Represent the `CELL` namelist of `cp.x`.

Input this namelist only if `calculation` is `"vc-relax"`, `"vc-cp"`, or `"vc-cp-wf"`.
"""
@with_kw struct CellNamelist <: Namelist
    cell_parameters::String = "default"
    cell_dynamics::String = "none"
    cell_velocities::String = "default"
    cell_damping::Float64 = 0.1
    press::Float64 = 0.0
    wmass::Float64 = 0.0
    cell_factor::Float64 = 1.2
    cell_temperature::String = "not_controlled"
    temph::Float64 = 0.0
    fnoseh::Float64 = 1.0
    greash::Float64 = 1.0
    cell_dofree::String = "all"
    # These checks are from https://github.com/QEF/q-e/blob/4132a64/Modules/read_namelists.f90#L1596-L1625.
    @assert cell_parameters in ("default", "from_input")
    @assert cell_dynamics in ("none", "sd", "damp-pr", "pr")
    @assert cell_velocities in ("zero", "default")
    @assert wmass >= 0
    @assert cell_temperature in ("nose", "rescaling", "not_controlled")
    @assert(
        cell_dofree in (
            "all",
            "x",
            "y",
            "z",
            "xy",
            "xz",
            "yz",
            "xyz",
            "shape",
            "2Dxy",
            "2Dshape",
            "volume",
        )
    )
end # struct CellNamelist

"""
    PressAiNamelist <: Namelist

Represent the `PRESS_AI` namelist of `cp.x`.

Input this namelist only when `tabps` is `true`.
"""
@with_kw struct PressAiNamelist <: Namelist
    abivol::Bool = false
    abisur::Bool = false
    P_ext::Float64 = 0.0
    pvar::Bool = false
    P_in::Float64 = 0.0
    P_fin::Float64 = 0.0
    Surf_t::Float64 = 0.0
    rho_thr::Float64 = 0.0
    dthr::Float64 = 0.0
end # struct PressAiNamelist

"""
    WannierNamelist <: Namelist

Represent the `WANNIER` namelist of `cp.x`.

Input this namelist only if `calculation` is `"cp-wf"` or `"vc-cp-wf"`.
"""
@with_kw struct WannierNamelist <: Namelist
    wf_efield::Bool = false
    wf_switch::Bool = false
    sw_len::Int = 1
    efx0::Float64 = 0.0
    efy0::Float64 = 0.0
    efz0::Float64 = 0.0
    efx1::Float64 = 0.0
    efy1::Float64 = 0.0
    efz1::Float64 = 0.0
    wfsd::Int = 1
    wfdt::Float64 = 5.0
    maxwfdt::Float64 = 0.3
    nit::Int = 10
    nsd::Int = 10
    wf_q::Float64 = 1500.0
    wf_friction::Float64 = 0.3
    nsteps::Int = 20
    tolw::Float64 = 1e-8
    adapt::Bool = true
    calwf::Int = 3
    nwf::Int = 0
    wffort::Int = 40
    writev::Bool = false
    exx_neigh::Int = 60
    exx_dis_cutoff::Float64 = 8.0
    exx_poisson_eps::Float64 = 1e-6
    exx_ps_rcut_self::Float64 = 6.0
    exx_ps_rcut_pair::Float64 = 5.0
    exx_me_rcut_self::Float64 = 10.0
    exx_me_rcut_pair::Float64 = 7.0
    # These checks are from https://github.com/QEF/q-e/blob/4132a64/Modules/read_namelists.f90#L1634-L1650.
    @assert 1 <= calwf <= 5
    @assert 1 <= wfsd <= 3
end # struct WannierNamelist


text:
using Dates
using Printf

t1 = now()

sleep(3)

t2 = now()

etime = t2 - t1
outtime = ["temp"]
outtime[1] = string(etime)

et = replace(outtime[1], " milliseconds"=>"")
et = parse(Float64, et)

# et_s = string( Int(floor(et / 10^3)))
# et_m = string( Int(floor(et / 10^3 / 60)))
# et_h = string( Int(floor(et / 10^3 / 60 /60)))

et_h = floor(et / 10^3 / 60 /60)
et = et - et_h * (10^3 * 60 * 60)
et_m = floor(et / 10^3 / 60 )
et = et - et_m * (10^3 * 60)
et_s = floor(et / 10^3)
et = et - et_s * (10^3)

et_h = @sprintf( "%02d", et_h)
et_m = @sprintf( "%02d", et_m)
et_s = @sprintf( "%02d", et_s)
et = et_h*":"*et_m*":"*et_s

println("\n")
println("\n")
println("---------------------------- \n")
println("  fin program  ")
println("\n")
println("  elapse time (hh:mm:ss) = " * et)
println("\n")
println("---------------------------- \n")

ext: jl

max_stars_repo_name: hide-dog/test

max_issues_repo_name: hide-dog/test

max_forks_repo_name: hide-dog/test


text:
@testset "plotting_range" begin
    @testset "types" begin
        @test UnicodePlots.plotting_range(0, 1) === (0.0, 1.0)
        @test UnicodePlots.plotting_range(0.0, 1) === (0.0, 1.0)
        @test UnicodePlots.plotting_range(0, 1f0) === (0.0, 1.0)
        @test UnicodePlots.plotting_range(0x0, 0x1) === (0.0, 1.0)
    end

    @test UnicodePlots.plotting_range(0.0001, 0.002) === (0.0, 0.01)
    @test UnicodePlots.plotting_range(0.001,  0.02)  === (0.0, 0.1)
    @test UnicodePlots.plotting_range(0.01,   0.2)   === (0.0, 1.0)
    @test UnicodePlots.plotting_range(0.1,    2.)    === (0.0, 10.0)
    @test UnicodePlots.plotting_range(0, 2) === (0.0, 10.0)
    @test UnicodePlots.plotting_range(0, 5) === (0.0, 10.0)
end

@testset "plotting_range_narrow" begin
    @testset "types" begin
        @test UnicodePlots.plotting_range_narrow(0, 1) === (0.0, 1.0)
        @test UnicodePlots.plotting_range_narrow(0.0, 1) === (0.0, 1.0)
        @test UnicodePlots.plotting_range_narrow(0, 1f0) === (0.0, 1.0)
        @test UnicodePlots.plotting_range_narrow(0x0, 0x1) === (0.0, 1.0)
    end

    @test UnicodePlots.plotting_range_narrow(0.0001, 0.002) === (0.0, 0.002)
    @test UnicodePlots.plotting_range_narrow(0.001,  0.02)  === (0.0, 0.02)
    @test UnicodePlots.plotting_range_narrow(0.01,   0.2)   === (0.0, 0.2)
    @test UnicodePlots.plotting_range_narrow(0.1,    2.)    === (0.0, 2.0)
    @test UnicodePlots.plotting_range_narrow(0, 2) === (0.0, 2.0)
    @test UnicodePlots.plotting_range_narrow(0, 5) === (0.0, 5.0)
end

@testset "extend_limits" begin
    @test UnicodePlots.extend_limits([1,2,3,4], [0.1,2]) === (0.1, 2.)
    @test UnicodePlots.extend_limits([1,2,3,4], [0,1.1]) === (0., 1.1)
    @test UnicodePlots.extend_limits([1,2,3,4], [2,3]) === (2., 3.)
    @test UnicodePlots.extend_limits([1,2,3,4], [0,0]) === (1., 4.)
    @test UnicodePlots.extend_limits([1,2,3,4], [1,1]) === (0., 2.)
end

@testset "bordermap" begin
    @test length(keys(UnicodePlots.bordermap)) == 8
    @test haskey(UnicodePlots.bordermap, :none)
    @test haskey(UnicodePlots.bordermap, :solid)
    @test haskey(UnicodePlots.bordermap, :corners)
    @test haskey(UnicodePlots.bordermap, :barplot)
    @test haskey(UnicodePlots.bordermap, :bold)
    @test haskey(UnicodePlots.bordermap, :dotted)
    @test haskey(UnicodePlots.bordermap, :dashed)
    @test haskey(UnicodePlots.bordermap, :ascii)
    for (k,v) in UnicodePlots.bordermap
        @test length(keys(v)) == 8
        @test haskey(v, :tl)
        @test haskey(v, :tr)
        @test haskey(v, :bl)
        @test haskey(v, :br)
        @test haskey(v, :r)
        @test haskey(v, :l)
        @test haskey(v, :t)
        @test haskey(v, :b)
    end
end


text:
eltype_or(::Type{G}, or) where G <: (AbstractGeometry{N, T} where N) where T = T
eltype_or(::Type{G}, or) where G <: (AbstractGeometry{N, T} where {N, T}) = or

ndims_or(::Type{G}, or) where G <: (AbstractGeometry{N, T} where T) where N = N
ndims_or(::Type{G}, or) where G <: (AbstractGeometry{N, T} where {N, T}) = or

Base.eltype(T::Type{<:AbstractGeometry}) = eltype_or(T, Any)

function Base.ndims(T::Type{<:AbstractGeometry})
    ndims_or(T, Any)
end

Base.eltype(x::AbstractGeometry{N, T}) where {N, T} = T
Base.ndims(x::AbstractGeometry{N, T}) where {N, T} = N

Base.eltype(::Type{AFG{T}}) where {T} = T
Base.eltype(::Type{FG}) where {FG <: AFG} = eltype(supertype(FG))

using Base.Cartesian

"""
deleteat{N,T,i}(x::NTuple{N,T}, ::Type{Val{i}})

Return copy of x with ith entry ommited.
"""
@generated function deleteat(x::NTuple{N,T}, ::Type{Val{i}}) where {N,T,i}
    (1 <= i <= N) || throw(MethodError(drop_index, (x,Val{i})))
    args = [:(x[$j]) for j in deleteat!([1:N...], i)]
    Expr(:tuple, args...)
end

"""
deleteat{N,T}(x::NTuple{N,T}, i::Int)
"""
@generated deleteat(x::NTuple{N,T}, i::Int) where {N,T} = quote
    # would be nice to eliminate boundscheck
    (1 <= i <= N) || throw(BoundsError(x,i))
    @nif $(N) d->(i==d) d-> deleteat(x, Val{d})
end

"""
best_arg, best_val = argmax(f, iter)

Computes the first element of iterator, on which f takes its maximum.
"""
function argmax(f, iter)
    iter_result = iterate(iter)
    if iter_result === nothing
        throw(ArgumentError("Cannot compute the argmax of an empty collection"))
    end
    best_arg, state = iter_result
    best_val = f(best_arg)
    while true
        iter_result = iterate(iter, state)
        iter_result === nothing && break
        arg, state = iter_result
        val = f(arg)
        if val > best_val
            best_val = val
            best_arg = arg
        end
    end
    best_arg, best_val
end

w_component(::Type{Point{N, T}}) where {N, T} = T(1)
w_component(::Type{Vec{N, T}}) where {N, T} = T(0)

@generated function transform_convert(::Type{T1}, x::T2) where {T1 <: StaticVector, T2 <: StaticVector}
    w = w_component(T1)
    n1 = length(T1)
    n2 = length(T2)
    n1 <= n2 && return :(T1(x))
    tupl = Expr(:tuple)
    ET = eltype(T1)
    for i = 1:n2
        push!(tupl.args, :($ET(x[$i])))
    end
    for i = 1:(n1 - n2 - 1)
        push!(tupl.args, :($ET(0)))
    end
    push!(tupl.args, :($w))
    :(T1($tupl))
end

@inline sqnorm(x) = sum(x.^2)

"""
pinvli{n,m,T}(a::Mat{n,m,T})

Compute pseudo inverse of matrix with linear independent columns.
"""
function pinvli(a::Mat{n,m,T}) where {n,m,T}
    @assert n >= m
    return inv(a'*a)*a'
end


text:
@testset "Position weight matrix" begin
    @testset "PFM" begin
        m = [
            1   2  3
            4   5  6
            7   8  9
            10 11 12
        ]
        pfm = PFM{DNA}(m)
        @test pfm isa PFM{DNA,Int}
        @test convert(Matrix, pfm) == m
        @test all(pfm[i] == m[i] for i in eachindex(pfm))
        @test all(pfm[i,j] == m[i,j] for i in 1:4, j in 1:3)
        @test all(pfm[ACGT[i],j] == m[i,j] for i in 1:4, j in 1:3)
        @test startswith(sprint(show, pfm), string(summary(pfm), ":\n"))
        @test_throws ArgumentError PFM{DNA}(vcat(m, [13 14 15]))

        # indexing
        pfm[5] = 100
        @test pfm[5] == pfm[DNA_A,2] == 100
        pfm[2,3] = 200
        @test pfm[2,3] == pfm[DNA_C,3] == 200
        pfm[DNA_G,1] = 300
        @test pfm[3,1] == pfm[DNA_G,1] == 300
        @test_throws ArgumentError pfm[DNA_N,1]

        # broadcast
        @test log.(pfm) isa PFM{DNA,Float64}
        @test log.(pfm) == PFM{DNA}(log.(m))
        @test pfm .+ 1 isa PFM{DNA,Int}
        @test pfm .+ 1 == PFM{DNA}(m .+ 1)
        @test pfm .+ 1.1 isa PFM{DNA,Float64}
        @test pfm .+ 1.1 == PFM{DNA}(m .+ 1.1)
        @test pfm .+ [0,1,2,3] isa PFM{DNA,Int}
        @test pfm .+ [0,1,2,3] == PFM{DNA}(m .+ [0,1,2,3])

        set = DNAKmer.(split(
        """
        ACG
        ATG
        AGG
        ACC
        """))
        pfm = PFM(set)
        @test pfm == [
            4 0 0  # A
            0 2 1  # C
            0 1 3  # G
            0 1 0  # T
        ]
        @test pfm == PFM(Set(set))
        @test_throws ArgumentError PFM(DNAKmer[])
        @test_throws ArgumentError PFM(["foo"])
        @test_throws ArgumentError PFM([DNAKmer("AA"), RNAKmer("AA")])
        @test_throws ArgumentError PFM([DNAKmer("AA"), DNAKmer("AAA")])
    end

    @testset "PWM" begin
        m = [
             1  2  3
             4  5  6
             7  8  9
            10 11 12
        ]
        pfm = PFM{DNA}(m)
        pwm = PWM(pfm)
        raw = log2.((pfm ./ sum(pfm, dims=1)) ./ fill(1/4, 4))
        @test pwm isa PWM{DNA,Float64}
        @test pwm == raw
        @test PWM(pfm, prior=[0.1, 0.4, 0.4, 0.1]) == log2.((pfm ./ sum(pfm, dims=1)) ./ [0.1, 0.4, 0.4, 0.1])
        @test PWM(pfm) == PWM{DNA}(raw)
        @test maxscore(pwm) ≈ sum(maximum(raw, dims=1))
        @test maxscore(PWM{DNA}(zeros(4, 0))) === 0.0
        @test PWM(pfm .+ 0.1) isa PWM{DNA,Float64}  # pseudo count
        @test_throws ArgumentError PWM{DNA}(hcat(m, [13, 14, 15]))
        @test_throws ArgumentError PWM(pfm, prior=normalize([0,1,2,3], 1))
        @test_throws ArgumentError PWM(pfm, prior=normalize([0,1,2,3], 1).+1e-3)
        @test all(pwm[i] === pwm[i] for i in eachindex(pwm))
        @test all(pwm[i,j] === raw[i,j] for i in 1:4, j in 1:3)
        @test all(pwm[ACGT[i],j] === pwm[i,j] for i in 1:4, j in 1:3)
        @test startswith(sprint(show, pwm), string(summary(pwm), ":\n"))
    end

    @testset "findfirst" begin
        seq = dna"ACGATNATCGCGTANTG"
        data = [
            1.0 0.1 0.2
            0.0 0.2 0.3
            0.1 0.2 0.0
            0.9 0.5 0.2
        ]
        pwm = PWM{DNA}(data)
        @test maxscore(pwm) == 1.8
        @test scoreat(seq, pwm, 1) === 1.2
        @test findfirst(pwm, seq, 1.0) === 1
        @test findfirst(pwm, seq, 1.4) === 4
        @test findfirst(pwm, seq, 1.8) === 7
        @test findfirst(pwm, seq, 2.0) === nothing
        @test_throws ArgumentError findfirst(pwm, RNASequence(seq), 1.0)
    end
end


text:
function genlib(cpucompiler::String, gpucompiler::String)
# Examples: genlib("g++","");
#           genlib("g++","nvcc");

if length(cpucompiler)==0
    error("cpucompiler is empty");
end

str = `$cpucompiler -fPIC -O3 -c commonCore.cpp`;
run(str);
str = `ar rvs commonCore.a commonCore.o`;
run(str);

str = `$cpucompiler -fPIC -O3 -c opuCore.cpp`;
run(str);
str = `ar rvs opuCore.a opuCore.o`;
run(str);

# str = `$cpucompiler -fPIC -O3 -c cpuCore.cpp -fopenmp`;
# run(str);
# str = `ar rvs cpuCore.a cpuCore.o`;
# run(str);

if length(gpucompiler)>0
    str = `$gpucompiler -D_FORCE_INLINES -O3 -c --compiler-options '-fPIC' gpuCore.cu`;
    run(str);
    str = `ar rvs gpuCore.a gpuCore.o `;
    run(str);
end

if Sys.isapple()
#    mv("*.a",pwd() * "/Mac");
    if !isdir("Mac")
        mkdir("Mac")
    end
    if isfile("commonCore.a")
        run(`mv commonCore.a Mac`);
    end
    if isfile("opuCore.a")
        run(`mv opuCore.a Mac`);
    end
    if isfile("cpuCore.a")
        run(`mv cpuCore.a Mac`);
    end
    if isfile("gpuCore.a")
        run(`mv gpuCore.a Mac`);
    end
elseif Sys.isunix()
    if !isdir("Linux")
        mkdir("Linux")
    end
    if isfile("commonCore.a")
        run(`mv commonCore.a Linux`);
    end
    if isfile("opuCore.a")
        run(`mv opuCore.a Linux`);
    end
    if isfile("cpuCore.a")
        run(`mv cpuCore.a Linux`);
    end
    if isfile("gpuCore.a")
        run(`mv gpuCore.a Linux`);
    end
elseif Sys.iswindows()
    if !isdir("Windows")
        mkdir("Windows")
    end
    if isfile("commonCore.a")
        run(`mv commonCore.a Windows`);
    end
    if isfile("opuCore.a")
        run(`mv opuCore.a Windows`);
    end
    if isfile("cpuCore.a")
        run(`mv cpuCore.a Windows`);
    end
    if isfile("gpuCore.a")
        run(`mv gpuCore.a Windows`);
    end
end

end

# delete *.o;
# delete *.a;


text:
module PersonalPlayground

greet() = print("Hello World! Thank you for using this useless package!")

end # module


text:
"""
HDF5.Filters

This module contains the interface for using filters in HDF5.jl.

# Example Usage
```julia
using HDF5
using HDF5.Filters

# Create a new file
fn = tempname()

# Create test data
data = rand(1000, 1000)

# Open temp file for writing
f = h5open(fn, "w")

# Create datasets
dsdeflate = create_dataset(f, "deflate", datatype(data), dataspace(data),
                           chunk=(100, 100), deflate=3)

dsshufdef = create_dataset(f, "shufdef", datatype(data), dataspace(data),
                           chunk=(100, 100), shuffle=true, deflate=3)

dsfiltdef = create_dataset(f, "filtdef", datatype(data), dataspace(data),
                           chunk=(100, 100), filters=Filters.Deflate(3))

dsfiltshufdef = create_dataset(f, "filtshufdef", datatype(data), dataspace(data),
                               chunk=(100, 100), filters=[Filters.Shuffle(), Filters.Deflate(3)])

# Write data
write(dsdeflate, data)
write(dsshufdef, data)
write(dsfiltdef, data)
write(dsfiltshufdef, data)

close(f)
```

## Additonal Examples

See [test/filter.jl](https://github.com/JuliaIO/HDF5.jl/blob/master/test/filter.jl) for further examples.
"""
module Filters

# builtin filters
export Deflate, Shuffle, Fletcher32, Szip, NBit, ScaleOffset

import ..HDF5: Properties, h5doc, API

"""
    Filter

Abstract type to describe HDF5 Filters.
See the Extended Help for information on implementing a new filter.

# Extended Help

## Filter interface

The Filter interface is implemented upon the Filter subtype.

See [`API.h5z_register`](@ref) for details.

### Required Methods to Implement
* `filterid` - registered filter ID
* `filter_func` - implement the actual filter

### Optional Methods to Implement
* `filtername` - defaults to "Unnamed Filter"
* `encoder_present` - defaults to true
* `decoder_present` - defaults to true
* `can_apply_func` - defaults to nothing
* `set_local_func` - defaults to nothing

### Advanced Methods to Implement
* `can_apply_cfunc` - Defaults to wrapping @cfunction around the result of `can_apply_func`
* `set_local_cfunc` - Defaults to wrapping @cfunction around the result of `set_local_func`
* `filter_cfunc` - Defaults to wrapping @cfunction around the result of `filter_func`
* `register_filter` - Defaults to using the above functions to register the filter

Implement the Advanced Methods to avoid @cfunction from generating a runtime closure which may not work on all systems.
"""
abstract type Filter end

"""
    FILTERS

Maps filter id to filter type.
"""
const FILTERS = Dict{API.H5Z_filter_t, Type{<: Filter}}()

"""
    filterid(F) where {F <: Filter}

The internal filter id of a filter of type `F`.
"""
filterid

"""
    encoder_present(::Type{F}) where {F<:Filter}

Can the filter have an encode or compress the data?
Defaults to true.
Returns a Bool. See [`API.h5z_register`](@ref).
"""
encoder_present(::Type{F}) where {F<:Filter} = true

"""
    decoder_present(::Type{F}) where {F<:Filter}

Can the filter decode or decompress the data?
Defaults to true.
Returns a Bool.
See [`API.h5z_register`](@ref)
"""
decoder_present(::Type{F}) where {F<:Filter} = true

"""
    filtername(::Type{F}) where {F<:Filter}

What is the name of a filter?
Defaults to "Unnamed Filter"
Returns a String describing the filter. See [`API.h5z_register`](@ref)
"""
filtername(::Type{F}) where {F<:Filter} = "Unnamed Filter"

"""
    can_apply_func(::Type{F}) where {F<:Filter}

Return a function indicating whether the filter can be applied or `nothing` if no function exists.
The function signature is `func(dcpl_id::API.hid_t, type_id::API.hid_t, space_id::API.hid_t)`.
See [`API.h5z_register`](@ref)
"""
can_apply_func(::Type{F}) where {F<:Filter} = nothing

"""
    can_apply_cfunc(::Type{F}) where {F<:Filter}

Return a C function pointer for the can apply function.
By default, this will return the result of using `@cfunction` on the function
specified by `can_apply_func(F)` or `C_NULL` if `nothing`.

Overriding this will allow `@cfunction` to return a `Ptr{Nothing}` rather
than a `CFunction`` closure which may not work on all systems.
"""
function can_apply_cfunc(::Type{F}) where {F<:Filter}
    func = can_apply_func(F)
    if func === nothing
        return C_NULL
    else
        return @cfunction($func, API.herr_t, (API.hid_t,API.hid_t,API.hid_t))
    end
end

"""
    set_local_func(::Type{F}) where {F<:Filter}

Return a function that sets dataset specific parameters or `nothing` if no function exists.
The function signature is `func(dcpl_id::API.hid_t, type_id::API.hid_t, space_id::API.hid_t)`.
See [`API.h5z_register`](@ref).
"""
set_local_func(::Type{F}) where {F<:Filter} = nothing

"""
    set_local_cfunc(::Type{F}) where {F<:Filter}

Return a C function pointer for the set local function.
By default, this will return the result of using `@cfunction` on the function
specified by `set_local_func(F)` or `C_NULL` if `nothing`.

Overriding this will allow `@cfunction` to return a `Ptr{Nothing}` rather
than a `CFunction`` closure which may not work on all systems.
"""
function set_local_cfunc(::Type{F}) where {F<:Filter}
    func = set_local_func(F)
    if func === nothing
        return C_NULL
    else
        return @cfunction($func, API.herr_t, (API.hid_t,API.hid_t,API.hid_t))
    end
end

"""
    filter_func(::Type{F}) where {F<:Filter}

Returns a function that performs the actual filtering.

See [`API.h5z_register`](@ref)
"""
filter_func(::Type{F}) where {F<:Filter} = nothing

"""
    filter_cfunc(::Type{F}) where {F<:Filter}

Return a C function pointer for the filter function.
By default, this will return the result of using `@cfunction` on the function
specified by `filter_func(F)` or will throw an error if `nothing`.

Overriding this will allow `@cfunction` to return a `Ptr{Nothing}` rather
than a `CFunction`` closure which may not work on all systems.
"""
function filter_cfunc(::Type{F}) where {F<:Filter}
    func = filter_func(F)
    if func === nothing
        error("Filter function for $f must be defined via `filter_func`.")
    end
    c_filter_func = @cfunction($func, Csize_t,
                               (Cuint, Csize_t, Ptr{Cuint}, Csize_t,
                               Ptr{Csize_t}, Ptr{Ptr{Cvoid}}))
    return c_filter_func
end

# Generic implementation of register_filter
"""
    register_filter(::Type{F}) where F <: Filter

Register the filter with the HDF5 library via [`API.h5z_register`](@ref).
Also add F to the FILTERS dictionary.
"""
function register_filter(::Type{F}) where F <: Filter
    id = filterid(F)
    encoder = encoder_present(F)
    decoder = decoder_present(F)
    name = filtername(F)
    can_apply = can_apply_cfunc(F)
    set_local = set_local_cfunc(F)
    func = filter_cfunc(F)
    GC.@preserve name begin
        API.h5z_register(API.H5Z_class_t(
            API.H5Z_CLASS_T_VERS,
            id,
            encoder,
            decoder,
            pointer(name),
            can_apply,
            set_local,
            func
        ))
    end
    FILTERS[id] = F
    return nothing
end

"""
    UnknownFilter(filter_id::API.H5Z_filter_t, flags::Cuint, data::Vector{Cuint}, name::String, config::Cuint)

An unknown filter.
"""
struct UnknownFilter <: Filter
    filter_id::API.H5Z_filter_t
    flags::Cuint
    data::Vector{Cuint}
    name::String
    config::Cuint
end
function UnknownFilter(filter_id, flags, data::Integer...)
    UnknownFilter(filter_id, flags, Cuint[data...], "Unknown Filter with id $filter_id", 0)
end
filterid(filter::UnknownFilter) = filter.filter_id
filtername(filter::UnknownFilter) = filter.name
filtername(::Type{UnknownFilter}) = "Unknown Filter"
encoder_present(::Type{UnknownFilter}) = false
decoder_present(::Type{UnknownFilter}) = false

"""
    ExternalFilter

External filter type. Alias for `UnkownFilter` (see related documentation).
Intended to support arbitrary, unregistered, external filters. Allows the
quick creation of filters using internal/proprietary filters without subtyping
`HDF5.Filters.Filter`.
Users are instead encouraged to define subtypes on `HDF5.Filters.Filter`.
"""
const ExternalFilter = UnknownFilter

"""
    FilterPipeline(plist::DatasetCreateProperties)

The filter pipeline associated with `plist`. Acts like a `AbstractVector{Filter}`,
supporting the following operations:

- `length(pipeline)`: the number of filters.
- `pipeline[i]` to return the `i`th filter.
- `pipeline[FilterType]` to return a filter of type `FilterType`
- `push!(pipline, filter)` to add an extra filter to the pipeline.
- `append!(pipeline, filters)` to add multiple filters to the pipeline.
- `delete!(pipeline, FilterType)` to remove a filter of type `FilterType` from the pipeline.
- `empty!(pipeline)` to remove all filters from the pipeline.
"""
struct FilterPipeline{P<:Properties} <: AbstractVector{Filter}
    plist::P
end

function Base.length(f::FilterPipeline)
    API.h5p_get_nfilters(f.plist)
end
Base.size(f::FilterPipeline) = (length(f),)

function Base.getindex(f::FilterPipeline, i::Integer)
    id = API.h5p_get_filter(f.plist, i-1, C_NULL, C_NULL, C_NULL, 0, C_NULL, C_NULL)
    F = get(FILTERS, id, UnknownFilter)
    return getindex(f, F, i)
end

function Base.getindex(f::FilterPipeline, ::Type{UnknownFilter}, i::Integer, cd_values::Vector{Cuint} = Cuint[])
    flags = Ref{Cuint}()
    cd_nelmts = Ref{Csize_t}(length(cd_values))
    namebuf = Array{UInt8}(undef, 256)
    config = Ref{Cuint}()
    id = API.h5p_get_filter(f.plist, i-1, flags, cd_nelmts, cd_values, length(namebuf), namebuf, config)
    if cd_nelmts[] > length(cd_values)
        resize!(cd_values, cd_nelmts[])
        return getindex(f, UnknownFilter, i, cd_values)
    end
    resize!(namebuf, findfirst(isequal(0), namebuf)-1)
    resize!(cd_values, cd_nelmts[])
    return UnknownFilter(id, flags[], cd_values, String(namebuf), config[])
end

function Base.getindex(f::FilterPipeline, ::Type{F}, i::Integer) where {F<:Filter}
    @assert isbitstype(F)
    ref = Ref{F}()
    GC.@preserve ref begin
        id = API.h5p_get_filter(f.plist, i-1, C_NULL, div(sizeof(F), sizeof(Cuint)), pointer_from_objref(ref), 0, C_NULL, C_NULL)
    end
    @assert id == filterid(F)
    return ref[]
end
function Base.getindex(f::FilterPipeline, ::Type{F}) where {F<:Filter}
    @assert isbitstype(F)
    ref = Ref{F}()
    GC.@preserve ref begin
        API.h5p_get_filter_by_id(f.plist, filterid(F), C_NULL, div(sizeof(F), sizeof(Cuint)), pointer_from_objref(ref), 0, C_NULL, C_NULL)
    end
    return ref[]
end

function Base.empty!(filters::FilterPipeline)
    API.h5p_remove_filter(filters.plist, API.H5Z_FILTER_ALL)
    return filters
end
function Base.delete!(filters::FilterPipeline, ::Type{F}) where {F<:Filter}
    API.h5p_remove_filter(filters.plist, filterid(F))
    return filters
end
function Base.append!(filters::FilterPipeline, extra::Union{AbstractVector{<:Filter}, NTuple{N, Filter} where N})
    for filter in extra
        push!(filters, filter)
    end
    return filters
end
function Base.push!(p::FilterPipeline, f::F) where F <: Filter
    ref = Ref(f)
    GC.@preserve ref begin
        API.h5p_set_filter(p.plist, filterid(F), API.H5Z_FLAG_OPTIONAL, div(sizeof(F), sizeof(Cuint)), pointer_from_objref(ref))
    end
    return p
end
function Base.push!(p::FilterPipeline, f::UnknownFilter)
    GC.@preserve f begin
        API.h5p_set_filter(p.plist, f.filter_id, f.flags, length(f.data), pointer(f.data))
    end
end

# Convert a Filter to an Integer subtype using filterid
function Base.convert(::Type{I}, ::Type{F}) where {I <: Integer, F <: Filter}
    Base.convert(I, filterid(F))
end

include("builtin.jl")
include("filters_midlevel.jl")

end # module


text:
using MadMax
using Test

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

X = [ 2 -1  0 -7
     -1  2 -1 -1
      0 -1  2 -1]

x, a = GaussSeidelEval(X, [0;0;0], N=5)
x, a = GaussJacobiEval(X, [0;0;0], N=5)
x, a = SOREval(X, [5.5;4.5;3.5], 1.1716, N=2)
println(x)


text:
# JulTran.jl

struct NONALLOCATO
    DT::DataType
end

macro declare(A, T, NX)
    return :( $(esc(A)) = NONALLOCATO(Array{$(esc(T)),$(esc(NX))}) )
end

macro Ndeclare(A, T)
    return :( $(esc(A)) = NONALLOCATO($(esc(T))) )
end

#https://stackoverflow.com/questions/39931112/quote-unquote-idiom-in-julia-concatenating-expr-objects
import Base.zero
zero(::Type{T}) where{T<:AbstractString} = string("")

macro declareNtimes(vars...)
    ex = Expr(:block)
    t = vars[1]
    ex.args = [:($(esc(n)) = Base.zero($(esc(t)))) for n in vars[2:end]]
    ex
end

#@declareNtimes Bool a b c

macro Ninit(A, T, V)
    return :( $(esc(A)) = Base.convert($(esc(T)),$(esc(V))) )
end

macro NinitNtimes(vars...)
    ex = Expr(:block)
    t = vars[1]
    ex.args = [:( $(esc(c.args[1])) = Base.convert($(esc(t)), $(esc(c.args[2]))) ) for c in vars[2:end]]
    ex
end

#@NinitNtimes Bool (a,false) (b,true)
#@assert !a && b

macro allocate(A, T, S)
    return :( $(esc(A)) = Base.fill(Base.zero($(esc(T))),$(esc(S))) )
end

macro deallocate(a)
    return :( $(esc(a)) = NONALLOCATO(Base.typeof($(esc(a)))) )
end

macro isallocated(a)
    return :( ! Base.isa($(esc(a)), NONALLOCATO) )
end


text:
using Images

function fitness(img, orig_img)
    ssdn(img, orig_img)
end

max_stars_repo_name: llrt/evolisa-jl

max_issues_repo_name: llrt/evolisa-jl

max_forks_repo_name: llrt/evolisa-jl


text:
module GPFlux

using LinearAlgebra
# using BackwardsLinalg
using Flux
using Flux: @functor, params
import Base: reset
import Distributions: MvNormal

export GaussProcess, MvNormal, negloglik, predict, params, dispatch!, flatten_params

export ConstantMean, SimpleNeuralNetworkMean

export ArdGaussKernel, IsoGaussKernel,
			 ArdRQKernel, IsoRQKernel,
			 ArdLinearKernel, IsoLinearKernel,
			 IsoPeriodKernel

export ProductCompositeKernel, AddCompositeKernel,
 			 SE_mul_PeriodKernel, Lin_mul_LinKernel, Period_mul_LinKernel, SE_mul_LinKernel,
 			 SE_add_PeriodKernel

export NeuralKernelNetwork, Primitive, Linear, Product, allProduct, allSum, f64

export norm2_metric, square_metric, inner_prod_metric

export gradient_check, model_gradient_check, reset

include("gp.jl")
include("mean.jl")
include("kernels/kernels.jl")
include("metrics/metrics.jl")
include("zygote.jl")

end # module

max_issues_repo_name: HamletWantToCode/GPFlux

max_forks_repo_name: HamletWantToCode/GPFlux


text:
using WebIO
using Pkg
Pkg.build("WebIO")
WebIO.install_jupyter_labextension([ARGS[1]])
WebIO.install_jupyter_nbextension([ARGS[1]])

max_stars_repo_name: benbenz/jupyterlab-econ-desktop

max_issues_repo_name: benbenz/jupyterlab-econ-desktop

max_forks_repo_name: benbenz/jupyterlab-econ-desktop


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

using Markdown
using InteractiveUtils

# ╔═╡ 88071c54-619d-11eb-3007-7db2c4a86215
begin
	using MLJ
	using LeastSquaresSVM
	using MLJModels
	using DataFrames
	using StatsPlots
	using BenchmarkTools
end

# ╔═╡ ae6803ac-619f-11eb-0d67-cb1be3400e9c
@load SVC pkg = LIBSVM;

# ╔═╡ b1598d44-619d-11eb-1efa-eb785775d940
gr();

# ╔═╡ cb3ab102-619d-11eb-1339-81cf03974fee
X, y = MLJ.make_moons(250; noise=0.3);

# ╔═╡ d7e2967c-619d-11eb-0d7f-9f0e79b48902
dfCircles = DataFrame(X);

# ╔═╡ e37d9ae0-619d-11eb-2ff0-3394d6447c9f
dfCircles.y = y;

# ╔═╡ 15c33f64-619e-11eb-0745-8fe5b6cf13a9
first(dfCircles, 3)

# ╔═╡ d0318284-619e-11eb-3fa1-59fb348e5d06
d = MLJ.int(y; type=Int);

# ╔═╡ 6db3e5d8-619f-11eb-1512-4d23ffd7adb3
dfCircles.c = d;

# ╔═╡ e7202e7e-619d-11eb-1b4a-e77c27e2d200
@df dfCircles scatter(:x1, :x2, colour=:c)

# ╔═╡ 3ca54cf6-619e-11eb-1a89-d7e83e9de721
function train_lssvm_hp(X, y, train)
	model = LSSVClassifier()

	r1 = range(model, :σ, lower=1e-2, upper=10.0)
	r2 = range(model, :γ, lower=1, upper=200.0)

	self_tuning_model = TunedModel(
		model=model,
		tuning=RandomSearch(),
		resampling=CV(nfolds=5),
		range=[r1, r2],
        measure=accuracy,
        n=1000
	)

	mach = machine(self_tuning_model, X, y)
	fit!(mach, rows=train)

	return mach
end;

# ╔═╡ a17664ba-619f-11eb-317c-79e5a717faca
function train_svm_hp(X, y, train)
	model = SVC()

	r1 = range(model, :cost, lower=1, upper=1000.0)
	r2 = range(model, :gamma, lower=1e-2, upper=10.0)

	self_tuning_model = TunedModel(
		model=model,
		tuning=RandomSearch(),
		resampling=CV(nfolds=5),
		range=[r1, r2],
        measure=accuracy,
        n=1000
	)

	mach = machine(self_tuning_model, X, y)
	fit!(mach, rows=train)

	return mach
end;

# ╔═╡ a6aa28f2-619f-11eb-0f6a-9167af8e877a
Xstd = MLJ.transform(MLJ.fit!(MLJ.machine(Standardizer(), X)), X);

# ╔═╡ e03600a0-619f-11eb-15a2-3f4677c2f723
train, test = MLJ.partition(eachindex(y), 0.8, shuffle=true);

# ╔═╡ 67d60f1a-61a4-11eb-3d38-35c2ccbd2744
@benchmark train_lssvm_hp($Xstd, $y, $train)

# ╔═╡ 70b7286c-61a4-11eb-28d0-4fbb79a7f373
@benchmark train_svm_hp($Xstd, $y, $train)

# ╔═╡ c1e62e68-619f-11eb-32b0-a744377d1b3c
mach1 = train_lssvm_hp(Xstd, y, train);

# ╔═╡ 2792dda8-619e-11eb-0da5-9712e7cb14de
mach2 = train_svm_hp(Xstd, y, train);

# ╔═╡ c4f54008-61a4-11eb-3b9e-994353e3ff5b
function predict_model(y, test, mach)
	results = MLJ.predict(mach, rows=test)
	acc = MLJ.accuracy(results, y[test])

	return acc
end;

# ╔═╡ 3689b474-61a5-11eb-0c60-a5c7c2a8e51e
acc1 = predict_model(y, test, mach1)

# ╔═╡ 4d1edf3c-61a5-11eb-30a9-3dc0a13fcba0
acc2 = predict_model(y, test, mach2)

# ╔═╡ 53c88e0c-61a5-11eb-3581-852491d59ffb

# ╔═╡ f9e794a6-61a3-11eb-316c-21c05947feaa

# ╔═╡ f6c26bfc-61a3-11eb-240c-97f3e8e48d63

# ╔═╡ e1dfbd84-61a3-11eb-2f44-f7146ce3045e

# ╔═╡ dbaacbf2-61a3-11eb-25b6-d33f4d04a74d

# ╔═╡ d449e050-61a3-11eb-1087-29a8e7b2b7dc

# ╔═╡ 03119e10-619e-11eb-28c6-ef029e684719

# ╔═╡ fcdf5636-619d-11eb-3019-953155bd46d3

# ╔═╡ b55b6cc8-619d-11eb-2829-7dbc5071a4e4

# ╔═╡ b544cca2-619d-11eb-30a5-c95dcb2294ea

# ╔═╡ b52c910a-619d-11eb-28af-b52571b41a82

# ╔═╡ b5135dd4-619d-11eb-226b-ad3792f806d6

# ╔═╡ b4fbe730-619d-11eb-294b-d100ca550814

# ╔═╡ b4e0b232-619d-11eb-3807-898664a71425

# ╔═╡ b4c5fc4c-619d-11eb-1567-650d59a0de08

# ╔═╡ b4abcbec-619d-11eb-3a17-b9064e2eedbc

# ╔═╡ b491fa8c-619d-11eb-3863-3ddd7fc0aa54

# ╔═╡ b4768d88-619d-11eb-0500-d3d24f8dbc5e

# ╔═╡ b45e3d8c-619d-11eb-081c-b9b953986714

# ╔═╡ b442d9b4-619d-11eb-36b5-81efc98e22cd

# ╔═╡ b42a464e-619d-11eb-11cb-870135a4f5a9

# ╔═╡ b40e3044-619d-11eb-0f41-8b657055d962

# ╔═╡ b3f4735c-619d-11eb-1610-057c920d8c69

# ╔═╡ b3a6646e-619d-11eb-1f5e-17020f493ab5

# ╔═╡ Cell order:
# ╠═88071c54-619d-11eb-3007-7db2c4a86215
# ╠═ae6803ac-619f-11eb-0d67-cb1be3400e9c
# ╠═b1598d44-619d-11eb-1efa-eb785775d940
# ╠═cb3ab102-619d-11eb-1339-81cf03974fee
# ╠═d7e2967c-619d-11eb-0d7f-9f0e79b48902
# ╠═e37d9ae0-619d-11eb-2ff0-3394d6447c9f
# ╠═15c33f64-619e-11eb-0745-8fe5b6cf13a9
# ╠═d0318284-619e-11eb-3fa1-59fb348e5d06
# ╠═6db3e5d8-619f-11eb-1512-4d23ffd7adb3
# ╠═e7202e7e-619d-11eb-1b4a-e77c27e2d200
# ╠═3ca54cf6-619e-11eb-1a89-d7e83e9de721
# ╠═a17664ba-619f-11eb-317c-79e5a717faca
# ╠═a6aa28f2-619f-11eb-0f6a-9167af8e877a
# ╠═e03600a0-619f-11eb-15a2-3f4677c2f723
# ╠═67d60f1a-61a4-11eb-3d38-35c2ccbd2744
# ╠═70b7286c-61a4-11eb-28d0-4fbb79a7f373
# ╠═c1e62e68-619f-11eb-32b0-a744377d1b3c
# ╠═2792dda8-619e-11eb-0da5-9712e7cb14de
# ╠═c4f54008-61a4-11eb-3b9e-994353e3ff5b
# ╠═3689b474-61a5-11eb-0c60-a5c7c2a8e51e
# ╠═4d1edf3c-61a5-11eb-30a9-3dc0a13fcba0
# ╠═53c88e0c-61a5-11eb-3581-852491d59ffb
# ╠═f9e794a6-61a3-11eb-316c-21c05947feaa
# ╠═f6c26bfc-61a3-11eb-240c-97f3e8e48d63
# ╠═e1dfbd84-61a3-11eb-2f44-f7146ce3045e
# ╠═dbaacbf2-61a3-11eb-25b6-d33f4d04a74d
# ╠═d449e050-61a3-11eb-1087-29a8e7b2b7dc
# ╠═03119e10-619e-11eb-28c6-ef029e684719
# ╠═fcdf5636-619d-11eb-3019-953155bd46d3
# ╠═b55b6cc8-619d-11eb-2829-7dbc5071a4e4
# ╠═b544cca2-619d-11eb-30a5-c95dcb2294ea
# ╠═b52c910a-619d-11eb-28af-b52571b41a82
# ╠═b5135dd4-619d-11eb-226b-ad3792f806d6
# ╠═b4fbe730-619d-11eb-294b-d100ca550814
# ╠═b4e0b232-619d-11eb-3807-898664a71425
# ╠═b4c5fc4c-619d-11eb-1567-650d59a0de08
# ╠═b4abcbec-619d-11eb-3a17-b9064e2eedbc
# ╠═b491fa8c-619d-11eb-3863-3ddd7fc0aa54
# ╠═b4768d88-619d-11eb-0500-d3d24f8dbc5e
# ╠═b45e3d8c-619d-11eb-081c-b9b953986714
# ╠═b442d9b4-619d-11eb-36b5-81efc98e22cd
# ╠═b42a464e-619d-11eb-11cb-870135a4f5a9
# ╠═b40e3044-619d-11eb-0f41-8b657055d962
# ╠═b3f4735c-619d-11eb-1610-057c920d8c69
# ╠═b3a6646e-619d-11eb-1f5e-17020f493ab5

max_stars_repo_name: edwinb-ai/LSSVM-benchmarks

max_issues_repo_name: edwinb-ai/LSSVM-benchmarks

max_forks_repo_name: edwinb-ai/LSSVM-benchmarks


text:
module MPQ

import Base.GMP.MPZ: mpz_t
import ..BigRationals: BigRational

const mpq_t = Ref{BigRational}

gmpq(op::Symbol) = (Symbol(:__gmpq_, op), :libgmp)

for (op) in (:canonicalize, :init, :clear)
    op! = Symbol(op, :!)
    @eval $op!(x::BigRational) = (ccall($(gmpq(op)), Cvoid, (mpq_t,), x); x)
end

for (op, T) in (:set_ui => Culong, :set_si => Clong)
    op! = Symbol(op, :!)
    @eval begin
        $op!(x::BigRational, a, b) = (ccall($(gmpq(op)), Cvoid, (mpq_t, $T, Culong), x, a, b); x)
        $op(a, b) = $op!(BigRational(), a, b)
    end
end

for (op, T) in (:set => mpq_t, :set_z => mpz_t, :set_d => Cdouble, :set_num => mpz_t, :set_den => mpz_t)
    op! = Symbol(op, :!)
    @eval $op!(x::BigRational, a) = (ccall($(gmpq(op)), Cvoid, (mpq_t, $T), x, a); x)
end

for op in (:get_num, :get_den)
    op! = Symbol(op, :!)
    @eval begin
        $op!(x::BigInt, a::BigRational) = (ccall($(gmpq(op)), Cvoid, (mpz_t, mpq_t), x, a); x)::BigInt
        $op(a::BigRational) = $op!(BigInt(), a)
    end
end

for (op) in (:inv, :neg, :abs)
    op! = Symbol(op, :!)
    @eval begin
        $op!(x::BigRational, a::BigRational) = (ccall($(gmpq(op)), Cvoid, (mpq_t, mpq_t), x, a); x)
        $op!(x::BigRational) = $op!(x, x)
    end
end

for op in (:set, :set_z, :set_d, :set_num, :inv, :neg, :abs)
    op! = Symbol(op, :!)
    @eval $op(a) = $op!(BigRational(), a)
end

for (op) in (:add, :sub, :mul, :div)
    op! = Symbol(op, :!)
    @eval begin
        $op!(x::BigRational, a::BigRational, b::BigRational) = (ccall($(gmpq(op)), Cvoid, (mpq_t, mpq_t, mpq_t), x, a, b); x)
        $op(a::BigRational, b::BigRational) = $op!(BigRational(), a, b)
        $op!(x::BigRational, b::BigRational) = $op!(x, x, b)
    end
end

for op in (:mul_2exp, :div_2exp)
    op! = Symbol(op, :!)
    @eval begin
        $op!(x::BigRational, a::BigRational, b) = (ccall($(gmpq(op)), Cvoid, (mpq_t, mpq_t, Culong), x, a, b); x)
        $op(a::BigRational, b) = $op!(BigRational(), a, b)
        $op!(x::BigRational, b) = $op!(x, x, b)
    end
end

cmp(a::BigRational, b::BigRational) = Int(ccall((:__gmpq_cmp, :libgmp), Cint, (mpq_t, mpq_t), a, b))
equal(a::BigRational, b::BigRational) = Bool(ccall((:__gmpq_equal, :libgmp), Cint, (mpq_t, mpq_t), a, b))

cmp_z(a::BigRational, b::BigInt) = Int(ccall((:__gmpq_cmp_z, :libgmp), Cint, (mpq_t, mpz_t), a, b))
cmp_ui(a::BigRational, b, c) = Int(ccall((:__gmpq_cmp_ui, :libgmp), Cint, (mpq_t, Culong, Culong), a, b, c))
cmp_si(a::BigRational, b, c) = Int(ccall((:__gmpq_cmp_si, :libgmp), Cint, (mpq_t, Clong, Culong), a, b, c))

get_d(a::BigRational) = ccall((:__gmpq_get_d, :libgmp), Cdouble, (mpq_t,), a)

get_str!(x, a, b::BigRational) = (ccall((:__gmpq_get_str, :libgmp), Ptr{Cchar}, (Ptr{Cchar}, Cint, mpq_t), x, a, b); x)
set_str!(x::BigRational, a, b) = Int(ccall((:__gmpq_set_str, :libgmp), Cint, (mpq_t, Ptr{UInt8}, Cint), x, a, b))

swap!(x::BigRational, y::BigRational) = ccall((:__gmpq_swap, :libgmp), Cvoid, (mpq_t, mpq_t), x, y)

end # module MPQ


text:
using Plots
x = Array(0:1:50)
y = 88 .+ (22-88)*exp.(-x/14)
scatter(x,y,title = "Theoretical Thermocouple Response", minorticks ="true")
xlabel!("Time (ms)")
ylabel!("Temperature (C)")


text:
using ArgCheck
using LazyArtifacts
using DataStructures: OrderedDict
using DocStringExtensions
################################################################################
##### testdatapath
################################################################################
function testdatapath(args...)
    joinpath(@__DIR__, "..", "test", "data", args...)
end

using .CAPI
using .CAPI: juliatype, EXECUTION_PROVIDERS
export InferenceSession, load_inference

"""
    $TYPEDEF

Represents an infernence session. Should only be created by calling [`load_inference`](@ref).
"""
struct InferenceSession
    api::OrtApi
    execution_provider::Symbol
    #env::OrtEnv
    session::OrtSession
    meminfo::OrtMemoryInfo
    allocator::OrtAllocator
    _input_names::Vector{String}
    _output_names::Vector{String}
end
function Base.show(io::IO, o::InferenceSession)
    print(io,
"""InferenceSession
    input_names:        $(input_names(o))
    output_names:       $(output_names(o))
    execution_provider: $(repr(o.execution_provider))
""")
end

input_names(o::InferenceSession) = o._input_names
output_names(o::InferenceSession) = o._output_names

function input_names(api::OrtApi, session::OrtSession, allocator::OrtAllocator)::Vector{String}
    n = SessionGetInputCount(api, session)
    map(0:n-1) do i
        SessionGetInputName(api, session, Csize_t(i), allocator)
    end
end
function output_names(api::OrtApi, session::OrtSession, allocator::OrtAllocator)::Vector{String}
    n = SessionGetOutputCount(api, session)
    map(0:n-1) do i
        SessionGetOutputName(api, session, Csize_t(i), allocator)
    end
end

"""
    $TYPEDSIGNATURES
"""
function load_inference(path::AbstractString; execution_provider::Symbol=:cpu,
        envname::AbstractString="defaultenv",
    )::InferenceSession
    api = GetApi(;execution_provider)
    env = CreateEnv(api, name=envname)
    if execution_provider === :cpu
        session_options = CreateSessionOptions(api)
    elseif execution_provider === :cuda
        if !(isdefined(@__MODULE__, :CUDA))
            @warn """
            The $(repr(execution_provider)) requires the CUDA.jl package to be available. Try adding `import CUDA` to your code.
            """
        end
        session_options = CreateSessionOptions(api)
        cuda_options = OrtCUDAProviderOptions()
        SessionOptionsAppendExecutionProvider_CUDA(api, session_options, cuda_options)
    else
        error("Unsupported execution_provider $execution_provider")
    end
    session = CreateSession(api, env, path, session_options)
    meminfo = CreateCpuMemoryInfo(api)
    allocator = CreateAllocator(api, session, meminfo)
    _input_names =input_names(api, session, allocator)
    _output_names=output_names(api, session, allocator)
    # TODO Is aliasing supported by ONNX? It will cause bugs, so lets forbid it.
    @check allunique(_input_names)
    @check allunique(_output_names)
    return InferenceSession(api, execution_provider, session, meminfo, allocator,
                            _input_names,
                            _output_names,
                           )
end

function make_input_tensor(o::InferenceSession, inputs, key)
    arr = inputs[keytype(inputs)(key)]
    cstorage = vec(reversedims(arr)::Array)
    CreateTensorWithDataAsOrtValue(o.api, o.meminfo, cstorage, size(arr))
end

function prepare_inputs(o::InferenceSession, inputs)
    names = input_names(o)
    tens = OrtValue[make_input_tensor(o, inputs, key) for key in names]
    names, tens
end

keytype(o) = Base.keytype(o)
keytype(o::NamedTuple) = Symbol

function make_output(o, inputs::NamedTuple, output_names, output_tensors)
    @argcheck length(output_names) == length(output_tensors)
    # TODO
    # we can probably optimize this
    # use unsafe_GetTensorMutableData and afterwards set `val.isalive = false` and
    # `empty!(gchandles)`
    pairs = (Symbol(key) => GetTensorMutableData(o.api, val) for (key, val) in zip(output_names, output_tensors))
    (;pairs...)
end
function make_output(o, inputs::AbstractDict, output_names, output_tensors)
    @argcheck length(output_names) == length(output_tensors)
    ret = OrderedDict{keytype(inputs), AbstractArray}()
    for (key, val) in zip(output_names, output_tensors)
        ret[key] = GetTensorMutableData(o.api, val)
    end
    ret
end

"""
    (o::InferenceSession)(inputs [,output_names])

Run an [`InferenceSession`](@ref) on a collection of inputs. Here `inputs` can either
be a `NamedTuple` or an `AbstractDict`. Optionally `output_names` can be passed.
In this case only the outputs whose name is contained in `output_names` are computed.
"""
function (o::InferenceSession)(
                               inputs,
                               output_names=nothing
                              )
    if output_names === nothing
        output_names = @__MODULE__().output_names(o)
    end
    @argcheck o.execution_provider in EXECUTION_PROVIDERS
    @argcheck eltype(output_names) <: Union{AbstractString, Symbol}
    @argcheck keytype(inputs) <: Union{AbstractString, Symbol}
    expected_input_names = ONNXRunTime.input_names(o)
    for key in keys(inputs)
        if !(String(key) in expected_input_names)
            msg = """
            Invalid input name.
            Expected: $(expected_input_names)
            Got: $(key)
            """
            throw(ArgumentError(msg))
        end
    end
    expected_output_names = ONNXRunTime.output_names(o)
    for name in output_names
        if !(String(name) in expected_output_names)
            msg = """
            Invalid output name.
            Expected: $(expected_output_names)
            Got: $(name)
            """
            throw(ArgumentError(msg))
        end
    end
    inp_names, input_tensors = prepare_inputs(o, inputs)
    run_options    = nothing
    output_tensors = Run(o.api, o.session, run_options, inp_names, input_tensors, output_names)
    make_output(o, inputs, output_names, output_tensors)
end


text:
#=
A lot of visualization forms in GLVisualize are realised in the form of instanced
particles. This is because they can be handled very efficiently by OpenGL.
There are quite a few different ways to feed instances with different attributes.
The main constructor for particles is a tuple of (Primitive, Position), whereas
position can come in all forms and shapes. You can leave away the primitive.
In that case, GLVisualize will fill in some default that is anticipated to make
the most sense for the datatype.
=#

#3D primitives
typealias Primitives3D  Union{AbstractGeometry{3}, AbstractMesh}
#2D primitives AKA sprites, since they are shapes mapped onto a 2D rectangle
typealias Sprites Union{AbstractGeometry{2}, Shape, Char, Type}
typealias AllPrimitives Union{AbstractGeometry, Shape, Char}

"""
We plot simple Geometric primitives as particles with length one.
At some point, this should all be appended to the same particle system to increase
performance.
"""
function _default{G<:GeometryPrimitive{2}}(
        geometry::TOrSignal{G}, s::Style, data::Dict
    )
    data[:offset] = Vec2f0(0)
    _default((geometry, const_lift(x->Point2f0[minimum(x)], geometry)), s, data)
end

"""
Vectors of floats are treated as barplots, so they get a HyperRectangle as
default primitive.
"""
function _default{T<:AbstractFloat}(main::VecTypes{T}, s::Style, data::Dict)
    _default((centered(HyperRectangle{2, Float32}), main), s, data)
end
"""
Matrices of floats are represented as 3D barplots with cubes as primitive
"""
function _default{T<:AbstractFloat}(main::MatTypes{T}, s::Style, data::Dict)
    _default((AABB(Vec3f0(-0.5,-0.5,0), Vec3f0(1.0)), main), s, data)
end
"""
Vectors of n-dimensional points get ndimensional rectangles as default
primitives. (Particles)
"""
function _default{N, T}(main::VecTypes{Point{N, T}}, s::Style, data::Dict)
    @gen_defaults! data begin
        scale = N == 2 ? Vec2f0(30) : Vec3f0(0.03) # for 2D points we assume they're in pixels
    end
    _default((centered(HyperRectangle{N, Float32}), main), s, data)
end

"""
3D matrices of vectors are 3D vector field with a pyramid (arrow) like default
primitive.
"""
function _default{T<:Vec}(main::ArrayTypes{T, 3}, s::Style, data::Dict)
    _default((Pyramid(Point3f0(0,0,-0.5), 1f0, 0.2f0), main), s, data)
end
"""
2D matrices of vectors are 2D vector field with a an unicode arrow as the default
primitive.
"""
function _default{T<:Vec}(main::ArrayTypes{T, 2}, s::Style, data::Dict)
    _default(('⬆', main), s, data)
end

"""
Vectors with `Vec` as element type are treated as vectors of rotations.
The position is assumed to be implicitely on the grid the vector defines (1D,2D,3D grid)
"""
function _default{P<:AllPrimitives, T<:Vec, N}(
        main::Tuple{P, ArrayTypes{T, N}}, s::Style, data::Dict
    )
    primitive, rotation_s = main
    rotation_v = value(rotation_s)
    @gen_defaults! data begin
        color_norm = const_lift(extrema2f0, rotation_s)
        ranges = ntuple(i->linspace(0f0, 1f0, size(rotation_v, i)), N)
    end
    grid = Grid(rotation_v, ranges)

    if N == 1
        scalevec = Vec2f0(step(grid.dims[1]), 1)
    elseif N == 2
        scalevec = Vec2f0(step(grid.dims[1]), step(grid.dims[2]))
    else
        scalevec = Vec3f0(ntuple(i->step(grid.dims[i]), 3)).*Vec3f0(0.4,0.4, 1/value(color_norm)[2]*4)
    end
    if P <: Char # we need to preserve proportion of the glyph
        glyphscale = primitive_scale(primitive)
        glyphscale /= max(glyphscale...)
        scalevec = Vec2f0(scalevec).*glyphscale
        @gen_defaults! data begin # for chars we need to make sure they're centered
            offset = -scalevec/2f0
        end
    end
    @gen_defaults! data begin
        rotation   = const_lift(vec, rotation_s)
        color_map  = default(Vector{RGBA})
        scale      = scalevec
        color      = nothing
    end
    _default((primitive, grid), s, data)
end

"""
arrays of floats with any geometry primitive, will be spaced out on a grid defined
by `ranges` and will use the floating points as the
height for the primitives (`scale_z`)
"""
function _default{P<:AbstractGeometry, T<:AbstractFloat, N}(
        main::Tuple{P, ArrayTypes{T,N}}, s::Style, data::Dict
    )
    primitive, heightfield_s = main
    heightfield = value(heightfield_s)
    @gen_defaults! data begin
        ranges = ntuple(i->linspace(0f0, 1f0, size(heightfield, i)), N)
    end
    grid = Grid(heightfield, ranges)
    @gen_defaults! data begin
        scale = nothing
        scale_x::Float32 = step(grid.dims[1])
        scale_y::Float32 = N==1 ? 1f0 : step(grid.dims[2])
        scale_z = const_lift(vec, heightfield_s)
        color = nothing
        color_map  = color == nothing ? default(Vector{RGBA}) : nothing
        color_norm = color == nothing ? const_lift(extrema2f0, heightfield_s) : nothing
    end
    _default((primitive, grid), s, data)
end

"""
arrays of floats with the sprite primitive type (2D geometry + picture like),
will be spaced out on a grid defined
by `ranges` and will use the floating points as the
z position for the primitives.
"""
function _default{P<:Sprites, T<:AbstractFloat, N}(
        main::Tuple{P, ArrayTypes{T,N}}, s::Style, data::Dict
    )
    primitive, heightfield_s = main
    heightfield = value(heightfield_s)
    @gen_defaults! data begin
        ranges = ntuple(i->linspace(0f0, 1f0, size(heightfield, i)), N)
    end
    grid = Grid(heightfield, ranges)
    @gen_defaults! data begin
        position_z = const_lift(vec, heightfield_s)
        scale      = Vec2f0(step(grid.dims[1]), N>=2 ? step(grid.dims[2]) : 1f0)
        color_map  = default(Vector{RGBA})
        color      = nothing
        color_norm = const_lift(extrema2f0, heightfield_s)
    end
    _default((primitive, grid), s, data)
end
"""
Sprites primitives with a vector of floats are treated as something barplot like
"""
function _default{P<:Sprites, T<:AbstractFloat}(
        main::Tuple{P, VecTypes{T}}, s::Style, data::Dict
    )
    primitive, heightfield_s = main
    heightfield = value(heightfield_s)
    @gen_defaults! data begin
        ranges = linspace(0f0, 1f0, length(heightfield))
    end
    grid = Grid(heightfield, ranges)
    @gen_defaults! data begin
        scale            = nothing
        scale_x::Float32 = step(grid.dims[1])
        scale_y          = heightfield_s
        scale_z::Float32 = 1f0
    end
    _default((primitive, grid), s, data)
end

# There is currently no way to get the two following two signatures
# under one function, which is why we delegate to meshparticle
function _default{Pr <: Primitives3D, P <: Point}(
        p::Tuple{TOrSignal{Pr}, VecTypes{P}}, s::Style, data::Dict
    )
    meshparticle(p, s, data)
end

function _default{Pr <: Primitives3D, G <: Tuple}(
        p::Tuple{TOrSignal{Pr}, G}, s::Style, data::Dict
    )
    @gen_defaults! data begin
        primitive = p[1]
        position         = nothing => TextureBuffer
        position_x       = p[2][1] => TextureBuffer
        position_y       = p[2][2] => TextureBuffer
        position_z       = length(p[2]) > 2 ? p[2][3] : 0f0 => TextureBuffer
        instances        = const_lift(length, position_x)
    end
    meshparticle(p, s, data)
end

function _default{Pr <: Primitives3D, G <: Grid}(
        p::Tuple{TOrSignal{Pr}, G}, s::Style, data::Dict
    )
    meshparticle(p, s, data)
end

# make conversion of mesh signals work. TODO move to GeometryTypes?
function Base.convert{T<:GeometryTypes.HomogenousMesh}(::Type{T}, mesh::Signal)
    map(T, mesh)
end
"""
This is the main function to assemble particles with a GLNormalMesh as a primitive
"""
function meshparticle(p, s, data)
    @gen_defaults! data begin
        primitive::GLNormalMesh = p[1]
        position         = p[2] => TextureBuffer
        position_x       = nothing => TextureBuffer
        position_y       = nothing => TextureBuffer
        position_z       = nothing => TextureBuffer

        scale            = Vec3f0(1) => TextureBuffer
        scale_x          = nothing => TextureBuffer
        scale_y          = nothing => TextureBuffer
        scale_z          = nothing => TextureBuffer

        rotation         = Vec3f0(0,0,1) => TextureBuffer
    end
    inst = _Instances(
        position, position_x, position_y, position_z,
        scale, scale_x, scale_y, scale_z,
        rotation, value(primitive)
    )
    @gen_defaults! data begin
        color_map  = nothing => Texture
        color_norm = nothing
        intensity = nothing => TextureBuffer
        color      = if color_map == nothing
            default(RGBA{Float32}, s)
        else
             nothing
        end => TextureBuffer

        instances   = const_lift(length, position)
        boundingbox = const_lift(GLBoundingBox, inst)
        shader      = GLVisualizeShader(
            "util.vert", "particles.vert", "fragment_output.frag", "standard.frag",
            view=Dict("position_calc"=>position_calc(position, position_x, position_y, position_z, TextureBuffer))
        )
    end
end

"""
This is the most primitive particle system, which uses simple points as primitives.
This is supposed to be the fastest way of displaying particles!
"""
_default{T <: Point}(position::VecTypes{T}, s::style"speed", data::Dict) = @gen_defaults! data begin
    vertex       = position => GLBuffer
    color_map    = nothing  => Vec2f0
    color        = (color_map == nothing ? default(RGBA{Float32}, s) : nothing) => GLBuffer

    color_norm   = nothing  => Vec2f0
    intensity = nothing  => GLBuffer
    point_size   = 2f0
    #boundingbox  = ParticleBoundingBox(position, Vec3f0(1), SimpleRectangle(-point_size/2,-point_size/2, point_size, point_size))
    prerender    = ()->glPointSize(point_size)
    shader       = GLVisualizeShader("fragment_output.frag", "dots.vert", "dots.frag")
    gl_primitive = GL_POINTS
end

"""
returns the Shape for the distancefield algorithm
"""
primitive_shape(::Char) = DISTANCEFIELD
primitive_shape{X}(x::X) = primitive_shape(X)
primitive_shape{T<:Circle}(::Type{T}) = CIRCLE
primitive_shape{T<:SimpleRectangle}(::Type{T}) = RECTANGLE
primitive_shape{T<:HyperRectangle{2}}(::Type{T}) = RECTANGLE
primitive_shape(x::Shape) = x

"""
Extracts the scale from a primitive.
"""
primitive_scale(prim::GeometryPrimitive) = Vec2f0(widths(prim))
primitive_scale(::Shape) = Vec2f0(40)
primitive_scale(c::Char) = Vec2f0(glyph_scale!(c))
primitive_scale(c) = Vec2f0(0.1)

"""
Extracts the offset from a primitive.
"""
#primitive_offset(prim::GeometryPrimitive) = Vec2f0(minimum(prim))

primitive_offset(x, scale::Void) = Vec2f0(0) # default offset
primitive_offset(x, scale) = const_lift(./, scale, -2f0)  # default offset

"""
Extracts the uv offset and width from a primitive.
"""
primitive_uv_offset_width(c::Char) = glyph_uv_width!(c)
primitive_uv_offset_width(x) = Vec4f0(0,0,1,1)

"""
Gets the texture atlas if primitive is a char.
"""
primitive_distancefield(x) = nothing
primitive_distancefield(::Char) = get_texture_atlas().images

"""
Particles with an image as primitive
"""
function _default{Pr <: Images.Image, P <: Point}(
        p::Tuple{TOrSignal{Pr}, VecTypes{P}}, s::Style, data::Dict
    )
    _default((const_lift(img -> img.data, p[1]), p[2]), s, data)
end
function _default{C <: Colorant, P <: Point}(
        p::Tuple{TOrSignal{Matrix{C}}, VecTypes{P}}, s::Style, data::Dict
    )
    data[:image] = p[1] # we don't want this to be overwritten by user
    @gen_defaults! data begin
        scale = map(Vec2f0, const_lift(size, p[1]))
        shape = RECTANGLE
        offset = Vec2f0(0)
    end
    sprites(p, s, data)
end
function _default{C <: AbstractFloat, P <: Point}(
        p::Tuple{TOrSignal{Matrix{C}}, VecTypes{P}}, s::Style, data::Dict
    )
    data[:distancefield] = p[1] # we don't want this to be overwritten by user
    @gen_defaults! data begin
        scale = map(Vec2f0, const_lift(size, p[1]))
        shape = RECTANGLE
        offset = Vec2f0(0)
    end
    sprites(p, s, data)
end

function _default{C <: Colorant, P <: Point}(
        p::Tuple{Vector{Matrix{C}}, VecTypes{P}}, s::Style, data::Dict
    )
    images = p[1]
    isempty(images) && error("Can not display empty vector of images as primitive")
    images = sort(images, by=size)
    sizes = map(size, images)
    scale = Vec2f0(first(sizes))
    if !all(x-> x == sizes[1], sizes) # if differently sized
        # create texture atlas
        maxdims = sum(map(Vec{2, Int}, sizes))
        rectangles = map(x->SimpleRectangle(0,0,x...), sizes)
        rpack = RectanglePacker(SimpleRectangle(0, 0, maxdims...))
        uv_coordinates = [push!(rpack, rect).area for rect in rectangles]
        max_xy = mapreduce(maximum, max, uv_coordinates)
        texture_atlas = Texture(C, tuple(max_xy...))
        for (area, img) in zip(uv_coordinates, images)
            texture_atlas[area] = img #transfer to texture atlas
        end
        scale = map(Vec2f0, map(widths, uv_coordinates))
        data[:uv_offset_width] = map(uv_coordinates) do uv
            mini = minimum(uv) ./ max_xy
            maxi = maximum(uv) ./ max_xy
            Vec4f0(mini..., maxi...)
        end
        images = texture_atlas
    end
    data[:image] = images # we don't want this to be overwritten by user
    @gen_defaults! data begin
        scale = scale
        shape = RECTANGLE
        offset = Vec2f0(0)
    end
    sprites(p, s, data)
end

# There is currently no way to get the two following two signatures
# under one function, which is why we delegate to sprites
_default{Primitive<:Sprites, P<:Point}(p::Tuple{Primitive, VecTypes{P}}, s::Style, data::Dict) =
    sprites(p,s,data)

_default{Primitive<:Sprites, G<:Grid}(p::Tuple{Primitive, G}, s::Style, data::Dict) =
    sprites(p,s,data)

function _default{Pr <: Sprites, G <: Tuple}(
            p::Tuple{Pr, G}, s::Style, data::Dict
        )
        @gen_defaults! data begin
            shape      = primitive_shape(p[1])
            position   = nothing => GLBuffer
            position_x = p[2][1] => GLBuffer
            position_y = p[2][2] => GLBuffer
            position_z = length(p[2]) > 2 ? p[2][3] : 0f0 => GLBuffer
        end
        sprites(p, s, data)
    end

"""
Main assemble functions for sprite particles.
Sprites are anything like distance fields, images and simple geometries
"""
function sprites(p, s, data)
    @gen_defaults! data begin
        shape       = primitive_shape(p[1])
        position    = p[2]    => GLBuffer
        position_x  = nothing => GLBuffer
        position_y  = nothing => GLBuffer
        position_z  = nothing => GLBuffer

        scale       = primitive_scale(p[1])  => GLBuffer
        scale_x     = nothing                => GLBuffer
        scale_y     = nothing                => GLBuffer
        scale_z     = nothing                => GLBuffer

        rotation    = Vec3f0(0,0,1)          => GLBuffer
        offset      = primitive_offset(p[1], scale) => GLBuffer

    end
    inst = _Instances(
        position, position_x, position_y, position_z,
        scale, scale_x, scale_y, scale_z,
        rotation, SimpleRectangle{Float32}(0,0,1,1)
    )
    @gen_defaults! data begin
        intensity = nothing => GLBuffer
        color_map        = nothing => Texture
        color_norm       = nothing
        color            = (color_map == nothing ? default(RGBA, s) : nothing) => GLBuffer

        glow_color       = RGBA{Float32}(0,0,0,0) => GLBuffer
        stroke_color     = RGBA{Float32}(0,0,0,0) => GLBuffer
        stroke_width     = 0f0
        glow_width       = 0f0
        uv_offset_width  = primitive_uv_offset_width(p[1]) => GLBuffer

        image            = nothing => Texture
        distancefield    = primitive_distancefield(p[1]) => Texture
        indices          = const_lift(length, p[2]) => to_indices
        boundingbox      = const_lift(GLBoundingBox, inst)
        # rotation and billboard don't go along
        billboard        = rotation == Vec3f0(0,0,1) => "if `billboard` == true, particles will always face camera"
        preferred_camera = :orthographic_pixel
        shader           = GLVisualizeShader(
            "fragment_output.frag", "util.vert", "sprites.geom",
            "sprites.vert", "distance_shape.frag",
            view=Dict("position_calc"=>position_calc(position, position_x, position_y, position_z, GLBuffer))
        )
        gl_primitive        = GL_POINTS
    end
end

"""
Transforms text into a particle system of sprites, by inferring the
texture coordinates in the texture atlas, widths and positions of the characters.
"""
function _default{S<:AbstractString}(main::TOrSignal{S}, s::Style, data::Dict)

    @gen_defaults! data begin
        relative_scale = Vec2f0(1)
        start_position = Point2f0(0)
        atlas          = get_texture_atlas()
        distancefield  = atlas.images
        stroke_width   = 0f0
        glow_width     = 0f0
        font           = DEFAULT_FONT_FACE
        scale_primitive = true
        position        = const_lift(calc_position, main, start_position, relative_scale, font, atlas)
        offset          = const_lift(calc_offset, main, relative_scale, font, atlas)
        uv_offset_width = const_lift(main) do str
            Vec4f0[glyph_uv_width!(atlas, c, font) for c=str]
        end
        scale           = const_lift(main, relative_scale) do str, s
            Vec2f0[glyph_scale!(atlas, c, font).*s for c=str]
        end
    end

    _default((DISTANCEFIELD, position), s, data)
end


text:
import MobergIO: Moberg, DigitalOut, DigitalIn, AnalogOut, AnalogIn, EncoderIn

Base.@kwdef struct FurutaPendulum <: PhysicalProcess
    m::Moberg                   = Moberg()
    energy_limited::DigitalIn   = DigitalIn(m, Unsigned(40))
    calibrate::DigitalOut       = DigitalOut(m, Unsigned(40))
    system_reset::DigitalOut    = DigitalOut(m, Unsigned(41))
    base_encoder::EncoderIn     = EncoderIn(m, Unsigned(40))
    arm_encoder::EncoderIn      = EncoderIn(m, Unsigned(41))
    base_angle::AnalogIn        = AnalogIn(m, Unsigned(40))
    arm_angle::AnalogIn         = AnalogIn(m, Unsigned(41))
    base_velocity::AnalogIn     = AnalogIn(m, Unsigned(42))
    arm_velocity::AnalogIn      = AnalogIn(m, Unsigned(43))
    total_energy::AnalogIn      = AnalogIn(m, Unsigned(44))
    base_energy::AnalogIn       = AnalogIn(m, Unsigned(45))
    arm_energy::AnalogIn        = AnalogIn(m, Unsigned(46))
    d_current::AnalogIn         = AnalogIn(m, Unsigned(47))
    q_current::AnalogIn         = AnalogIn(m, Unsigned(48))
    torque::AnalogOut           = AnalogOut(m, Unsigned(40))  # Uses internal FOC controller which does not work well
    voltage::AnalogOut          = AnalogOut(m, Unsigned(41))

end


text:
export AbstractStateSpace, ndim, volume, isuniform, encode, decode

abstract type AbstractStateSpace end

@unimpl Base.size(space::AbstractStateSpace)

ndim(space::AbstractStateSpace) = space.ndim

volume(space::AbstractStateSpace) = space.volume
Base.length(space::AbstractStateSpace) = volume(space)

isuniform(space::AbstractStateSpace) = space.isuniform

@unimpl encode(space::AbstractStateSpace, state)

@unimpl decode(space::AbstractStateSpace, code)

@unimpl Base.iterate(space::AbstractStateSpace)
@unimpl Base.iterate(space::AbstractStateSpace, state)


text:
module Plasmons

export dielectric, polarizability, dispersion
export julia_main

using CUDA
using Adapt
using LinearAlgebra
using ArgParse
using HDF5
using Dates
using Logging
using LoggingExtras

include("utilities.jl")

"""
    fermidirac(E; mu, kT) -> f

Return Fermi-Dirac distribution ``f`` at energy `E`, chemical potential `mu`, and temperature
`kT`. Note that `kT` is assumed to be temperature multiplied by the Boltzmann constant, i.e.
physical dimension of `kT` is the same as `E` (e.g. electron-volts).
"""
@inline fermidirac(E; mu, kT) = 1 / (1 + exp((E - mu) / kT))

function _build_g(ħω::Complex{ℝ}, E::AbstractVector{ℝ}; mu::ℝ, kT::ℝ) where {ℝ <: Real}
    f = map(x -> fermidirac(x, mu = mu, kT = kT), E)
    n = length(E)
    G = similar(E, complex(ℝ), n, n)
    @inbounds for j in 1:n
        @. G[:, j] = (f - f[j]) / (E - E[j] - ħω)
    end
    return G
end
function _postprocess_g(
    ::Type{T},
    G::AbstractMatrix,
    blocks::Bool,
) where {T <: Union{Real, Complex}}
    split_if(x) = T <: Real ? (real(x), imag(x)) : x
    compress_if(x) = blocks ? ThreeBlockMatrix(x) : x
    return _map_tuple(compress_if, split_if(G))
end
"""
    _map_tuple(f::Function, x)

Apply `f` to every element of `x`, but only if `x` is a `Tuple`. If not, just apply `f` to
the whole `x`.
"""
_map_tuple(f::Function, x::Tuple) = f.(x)
_map_tuple(f::Function, x::Any) = f(x)

@doc raw"""
    g_matrix(::Type{ℝ}, ħω::ℂ, E; mu::ℝ, kT::ℝ, blocks::Bool) -> (Gᵣ, Gᵢ)
    g_matrix(::Type{ℂ}, ħω::ℂ, E; mu::ℝ, kT::ℝ, blocks::Bool) -> G

Compute matrix ``G``
```math
    G_{ij} = \frac{f(E_i) - f(E_j)}{E_i - E_j - \hbar\omega}
```
where ``f`` is [Fermi-Dirac
distribution](https://en.wikipedia.org/wiki/Fermi%E2%80%93Dirac_statistics) at chemical
potential `mu` (``\mu``) and temperature `kT` (``k_B T``). `E` is a vector of eigenenergies.
`ħω` is a complex frequency including Landau damping (i.e.  ``\hbar\omega + i\eta``). All
arguments are assumed to be in energy units.

Matrix ``G`` is always complex, but the first argument specifies which type should be used
to represent it. For example, if you call `g_matrix(Float32, ...)` then a tuple of two
values will be returned -- representing real and imaginary parts of ``G``. If, however, you
call `g_matrix(ComplexF64, ...)` then ``G`` will not be split into real and imaginary parts
and a single object is returned.

Sometimes one can further exploit the structure of ``G``. For ``E \ll \mu`` or ``E \gg \mu``
Fermi-Dirac distribution is just a constant and ``G`` goes to 0 for all ``\omega``. We can
use this fact to construct a block-sparse version of ``G`` to speed up matrix-matrix
products containing it. To return ``G`` in block-sparse form, set `blocks = true`.

This function works with both CPU and GPU arrays. If vector `E` is stored on the GPU, then
the returned matrix (or all submatrices for block-sparse version) will also reside on GPU.
"""
g_matrix(::Type{T}, ħω::Complex{ℝ}, E::Vector{ℝ}; mu::ℝ, kT::ℝ, blocks::Bool) where {T, ℝ <: Real} =
    _postprocess_g(T, _build_g(ħω, E; mu = mu, kT = kT), blocks)
g_matrix(
    ::Type{T},
    ħω::Complex{ℝ},
    E::CuVector{ℝ};
    mu::ℝ,
    kT::ℝ,
    blocks::Bool,
) where {T, ℝ <: Real} = _map_tuple(
    x -> adapt(CuArray, x),
    _postprocess_g(T, _build_g(ħω, Vector(E); mu = mu, kT = kT), blocks),
)

@doc raw"""
    polarizability(ħω::ℂ, E, ψ; mu::ℝ, kT::ℝ, method = :batched, blocks = true) -> χ

Compute polarizability matrix ``\chi`` using method `method` (either `:simple` or
`:batched`). `ħω` specifies complex frequency including Landau damping (i.e. ``\hbar\omega +
i\eta``). `mu` and `kT` specify chemical potential and temperature respectively. `E` and `ψ`
are eigenvalues and eigenvectors of the system Hamiltonian. They can be either CPU or GPU
arrays.
"""
function polarizability(
    ħω::Complex,
    E::AbstractVector{<:Real},
    ψ::AbstractMatrix{<:Union{Real, Complex}};
    mu::Real,
    kT::Real,
    method::Symbol = :batched,
    blocks::Bool = true,
)
    if size(ψ, 1) != size(ψ, 2)
        throw(DimensionMismatch("'ψ' has wrong shape: $(size(ψ)); expected a square matrix"))
    end
    if size(E, 1) != size(ψ, 1)
        throw(DimensionMismatch(
            "dimensions of 'E' and 'ψ' do not match: $(size(E)) & $(size(ψ)); " *
            "expected 'ψ' to be of the same dimension as 'E'",
        ))
    end
    if method == :simple
        polarizability_simple(ħω, E, ψ; mu = mu, kT = kT)
    elseif method == :batched
        polarizability_batched(ħω, E, ψ; mu = mu, kT = kT, blocks = blocks)
    else
        throw(ArgumentError("invalid 'method': $method; expected either :simple or :batched"))
    end
end

struct Workspace{M <: AbstractMatrix}
    A::M
    temp::M
end

@doc raw"""
    polarizability_batched(ħω::ℂ, E, ψ; mu::ℝ, kT::ℝ, blocks::Bool) -> χ

Compute polarizability matrix ``\chi`` using "batched" method.
"""
function polarizability_batched(
    ħω::Complex,
    E::AbstractVector{ℝ},
    ψ::AbstractMatrix{ℂ};
    mu::ℝ,
    kT::ℝ,
    blocks::Bool,
) where {ℝ <: Real, ℂ <: Union{ℝ, Complex{ℝ}}}
    G = g_matrix(ℂ, ħω, E; mu = mu, kT = kT, blocks = blocks)
    ℂ <: Real ? polarizability_batched(G[1], G[2], ψ) : polarizability_batched(G, ψ)
end
function polarizability_batched(
    Gᵣ::Union{AbstractMatrix{ℝ}, ThreeBlockMatrix{ℝ}},
    Gᵢ::Union{AbstractMatrix{ℝ}, ThreeBlockMatrix{ℝ}},
    ψ::AbstractMatrix{ℝ},
) where {ℝ <: Real}
    N = size(ψ, 1)
    χ = similar(ψ, complex(ℝ))
    fill!(χ, zero(eltype(χ)))
    ws = Workspace(similar(ψ), similar(ψ))
    for b in 1:N
        a = b:N
        _batched_mat_el!(view(χ, a, b), ws, a, b, Gᵣ, Gᵢ, ψ)
    end
    χ .= 2 .* (χ .+ transpose(χ))
    χ[diagind(χ)] ./= 2
    χ
end
function polarizability_batched(
    G::Union{AbstractMatrix{ℂ}, ThreeBlockMatrix{ℂ}},
    ψ::AbstractMatrix{ℂ},
) where {ℂ <: Complex}
    N = size(ψ, 1)
    χ = similar(ψ)
    fill!(χ, zero(eltype(χ)))
    ws = Workspace(similar(ψ), similar(ψ))
    for b in 1:N
        _batched_mat_el!(view(χ, :, b), ws, b, G, ψ)
    end
    χ .*= 2
    χ
end
function _batched_mat_el!(
    out::AbstractVector{Complex{ℝ}},
    ws::Workspace,
    as::UnitRange{Int},
    b::Int,
    Gᵣ::Union{AbstractMatrix{ℝ}, ThreeBlockMatrix{ℝ}},
    Gᵢ::Union{AbstractMatrix{ℝ}, ThreeBlockMatrix{ℝ}},
    ψ::AbstractMatrix{ℝ},
) where {ℝ <: Real}
    batch_size = length(as)
    A = view(ws.A, 1:batch_size, :)
    temp = view(ws.temp, 1:batch_size, :)
    A .= view(ψ, as, :) .* transpose(view(ψ, b, :))
    _mul_with_strides!(temp, A, Gᵣ, one(ℝ), zero(ℝ))
    outᵣ = view(reinterpret(ℝ, out), 1:2:(2 * batch_size))
    dot_batched!(outᵣ, A, temp)
    _mul_with_strides!(temp, A, Gᵢ, one(ℝ), zero(ℝ))
    outᵢ = view(reinterpret(ℝ, out), 2:2:(2 * batch_size))
    dot_batched!(outᵢ, A, temp)
    out
end
function _batched_mat_el!(
    out::AbstractVector{ℂ},
    ws::Workspace,
    b::Int,
    G::Union{AbstractMatrix{ℂ}, ThreeBlockMatrix{ℂ}},
    ψ::AbstractMatrix{ℂ},
) where {ℂ <: Complex}
    A = ws.A
    temp = ws.temp
    A .= ψ .* conj.(transpose(view(ψ, b, :)))
    _mul_with_strides!(temp, A, G, one(ℂ), zero(ℂ))
    dot_batched!(out, A, temp)
    out
end

@doc raw"""
    polarizability_simple(ħω, E, ψ; mu, kT) -> χ

Calculates the polarizability matrix ``\chi`` by direct evaluation of eq. (3). This approach
is single-threaded and doesn't utilize the CPU well. It is provided here for testing
purposes only.
"""
function polarizability_simple(
    ħω::Complex{ℝ},
    E::AbstractVector{ℝ},
    ψ::AbstractMatrix{<:Union{ℝ, Complex{ℝ}}};
    mu::ℝ,
    kT::ℝ,
) where {ℝ <: Real}
    G = _build_g(ħω, E; mu = mu, kT = kT)
    χ = similar(G)
    @inbounds for a in 1:size(χ, 2)
        for b in 1:size(χ, 1)
            acc = zero(eltype(χ))
            for i in 1:length(E)
                for j in 1:length(E)
                    acc += G[i, j] * conj(ψ[a, j]) * ψ[a, i] * conj(ψ[b, i]) * ψ[b, j]
                end
            end
            χ[a, b] = 2 * acc
        end
    end
    return χ
end

@doc raw"""
    dielectric(χ, V) -> ε

Compute dielectric function ``\varepsilon`` given polarizability matrix ``\chi`` and Coulomb
interaction potential ``V``. Dielectric function is simply ``\varepsilon = 1 - \chi V``.
"""
function dielectric(χ::AbstractMatrix{Complex{ℝ}}, V::AbstractMatrix{ℝ}) where {ℝ <: Real}
    ℂ = complex(ℝ)
    if size(χ, 1) != size(χ, 2) || size(χ) != size(V)
        throw(DimensionMismatch(
            "dimensions of χ and V do not match: $(size(χ)) != $(size(V)); " *
            "expected two square matrices of the same size",
        ))
    end
    ε = similar(χ)
    fill!(ε, zero(ℂ))
    @inbounds ε[diagind(ε)] .= one(ℂ)
    mul!(ε, V, χ, -one(ℂ), one(ℂ))
    return ε
end

function _momentum_eigenvectors(
    qs::AbstractVector{NTuple{3, ℝ}},
    rs::AbstractVector{NTuple{3, ℝ}}
) where {ℝ <: Real}
    out = similar(rs, complex(ℝ), length(rs), length(qs))
    @inbounds for j in 1:size(out, 2)
        (qˣ, qʸ, qᶻ) = qs[j]
        @simd for i in 1:size(out, 1)
            (xᵢ, yᵢ, zᵢ) = rs[i]
            out[i, j] = (1 / sqrt(length(rs))) * exp(1im * (qˣ * xᵢ + qʸ * yᵢ + qᶻ * zᵢ))
        end
    end
    out
end
function dispersion_function(
    qs::AbstractVector{NTuple{3, ℝ}},
    rs::AbstractVector{NTuple{3, ℝ}},
) where {ℝ <: Real}
    return let r_q = _momentum_eigenvectors(qs, rs),
        temp = similar(r_q, length(rs), size(r_q, 2))

        A -> begin
            mul!(temp, A, r_q)
            out = similar(r_q, length(qs))
            @inbounds for i in 1:length(out)
                out[i] = dot(view(r_q, :, i), view(temp, :, i))
            end
            out
        end
    end
end
"""
    dispersion(A(ω, r, r'), q, x, y, z; n = 100) -> A(ω, q)

Calculate diagonal elements of the Fourier transform of matrix A.
"""
# function dispersion(As, qs, x, y, z)
#     fn! = _dispersion_function(qs, x, y, z)
#     matrix = Any[]
#     for (i, A) in enumerate(As)
#         out = similar(x, complex(eltype(x)), length(qs))
#         fn!(out, A)
#         push!(matrix, out)
#     end
#     hcat(matrix...)
# end
dispersion(A::AbstractMatrix, qs::AbstractVector, rs::AbstractVector) = dispersion_function(qs, rs)(A)

function main(
    E::AbstractVector{ℝ},
    ψ::AbstractMatrix{ℂ};
    kT::Real,
    μ::Real,
    η::Real,
    ωs::Vector{<:Real},
    out::Union{HDF5.File, HDF5.Group},
    V::Union{AbstractMatrix{ℝ}, Nothing} = nothing,
) where {ℝ <: Real, ℂ <: Union{ℝ, Complex{ℝ}}}
    if kT <= 0
        throw(ArgumentError("invalid 'kT': $kT; expected a positive real number"))
    end
    if η <= 0
        throw(ArgumentError("invalid 'η': $η; expected a positive real number"))
    end
    HDF5.attributes(out)["kT"] = kT
    HDF5.attributes(out)["μ"] = μ
    HDF5.attributes(out)["η"] = η
    group_χ = create_group(out, "χ")
    group_ε::Union{HDF5.Group, Nothing} = isnothing(V) ? nothing : create_group(out, "ε")
    @info "Polarizability matrices χ(ω) will be saved to group 'χ'"
    if !isnothing(V)
        @info "Dielectric functions ε(ω) will be saved to group 'ε'"
    else
        @warn "Coulomb interaction matrix not provided: dielectric function ε(ω) will not be computed"
    end

    for (i, ω) in enumerate(map(x -> x + 1im * η, ωs))
        @info "Calculating χ(ω = $ω) ..."
        name = string(i, pad = 4)
        t₀ = time_ns()
        χ = Array(polarizability(
            convert(complex(ℝ), ω),
            E,
            ψ;
            mu = convert(ℝ, μ),
            kT = convert(ℝ, kT),
        ))
        t₁ = time_ns()
        group_χ[name] = χ
        HDF5.attributes(group_χ[name])["ħω"] = ω
        HDF5.attributes(group_χ[name])["time"] = (t₁ - t₀) / 1e9
        flush(group_χ)
        if !isnothing(V)
            t₀ = time_ns()
            ε = Array(dielectric(χ, V))
            t₁ = time_ns()
            group_ε[name] = ε
            HDF5.attributes(group_ε[name])["ħω"] = ω
            HDF5.attributes(group_ε[name])["time"] = (t₁ - t₀) / 1e9
            flush(group_ε)
        end
    end
end
function main!(H::AbstractMatrix{ℂ}; kwargs...) where {ℂ}
    @warn "Eigenvalues or eigenvectors not provided: diagonalizing the Hamiltonian..."
    E, ψ = _eigen!(Hermitian(H))
    @info "Diagonalization completed"
    main(E, ψ; kwargs...)
end

ArgParse.parse_item(::Type{Vector{Float64}}, x::AbstractString) =
    map(s -> parse(Float64, s), split(x, ","))

function _parse_commandline()
    s = ArgParseSettings()
    #! format: off
    @add_arg_table! s begin
        "--kT"
            help = "Temperature kʙ·T in eletron-volts"
            arg_type = Float64
            required = true
        "--mu"
            help = "Chemical potential μ in eletron-volts"
            arg_type = Float64
            required = true
        "--damping"
            help = "Landau damping η in eletron-volts"
            arg_type = Float64
            required = true
        "--frequency"
            help = "Comma-separated list of ħω values in eletron-volts for which to compute χ and ε"
            arg_type = Vector{Float64}
            required = true
        "--hamiltonian"
            help = "Path to Hamiltonian matrix in the input HDF5 file"
            arg_type = String
        "--coulomb"
            help = "Path to unscreened Coulomb interaction matrix in the input HDF5 file"
            arg_type = String
        "--eigenvalues"
            help = "Path to eigenvalues vector in the input HDF5 file"
            arg_type = String
        "--eigenvectors"
            help = "Path to eigenvectors matrix in the input HDF5 file"
            arg_type = String
        "--cuda"
            help = "Index of CUDA device to use"
            arg_type = Int
        "input_file"
            help = "Input HDF5 file"
            required = true
        "output_file"
            help = "Output HDF5 file"
            required = true
    end
    #! format: on
    return parse_args(s, as_symbols = true)
end

function tryread(io::HDF5.File, path::Union{<:AbstractString, Nothing} = nothing)
    isnothing(path) && return nothing
    if !haskey(io, path)
        @warn "Path '$path' not found in input HDF5 file"
        return nothing
    end
    read(io, path)
end

const date_format = "yyyy-mm-dd HH:MM:SS"

timestamp_logger(logger) =
    TransformerLogger(logger) do log
        merge(log, (; message = "$(Dates.format(now(), date_format))  $(log.message)"))
    end

function julia_main()::Cint
    CUDA.allowscalar(false)
    args = _parse_commandline()
    H, E, ψ, V = h5open(args[:input_file], "r") do io
        tryread(io, get(args, :hamiltonian, nothing)),
        tryread(io, get(args, :eigenvalues, nothing)),
        tryread(io, get(args, :eigenvectors, nothing)),
        tryread(io, get(args, :coulomb, nothing))
    end
    if (isnothing(E) || isnothing(ψ)) && isnothing(H)
        @error "if Hamiltonian 'H' is not specified, both eigenvalues 'E'" *
               " and eigenvectors 'ψ' must be provided"
        return 1
    end
    cuda_device = args[:cuda]
    if !isnothing(cuda_device)
        @info "Setting CUDA device to $cuda_device..."
        device!(cuda_device)
        if !isnothing(H)
            H = adapt(CuArray, H)
        end
        if !isnothing(E)
            E = adapt(CuArray, E)
        end
        if !isnothing(ψ)
            ψ = adapt(CuArray, ψ)
        end
        if !isnothing(V)
            V = adapt(CuArray, V)
        end
    end
    h5open(args[:output_file], "w") do io
        kwargs = (
            kT = args[:kT],
            μ = args[:mu],
            η = args[:damping],
            ωs = args[:frequency],
            out = io,
            V = V,
        )
        with_logger(ConsoleLogger(stdout, Logging.Info) |> timestamp_logger) do
            if (isnothing(E) || isnothing(ψ))
                main!(H; kwargs...)
            else
                main(E, ψ; kwargs...)
            end
        end
    end
    return 0
end

end # module


text:
import PanGraph

cdir = pwd()
home = replace(cdir,"/pangraph/lib/julia"=>"")
root = "$home/data/synthetic"

# usage
PanGraph.main([])                      # main usage
PanGraph.main(["help", "build"])       # build usage
PanGraph.main(["help", "export"])      # export usage
PanGraph.main(["help", "generate"])    # generate usage
PanGraph.main(["help", "marginalize"]) # marginalize usage
PanGraph.main(["help", "polish"])      # polish usage
PanGraph.main(["help", "help"])		   # help usage

# build
PanGraph.main(["build", "-c", "-u", "-b", "0", "-a", "0", "$root/test.fa"])

# export
PanGraph.main(["export", "-o", "$root/export", "$root/test.json"])

# export
PanGraph.main(["generate", "-o", "$root/export.json", "$root/test.fa"])

# marginalize
PanGraph.main(["marginalize", "-o", "$root/pairs", "$root/test.json"])

# polish
PanGraph.main(["polish", "-c", "-l", "10000", "$root/test.json"])

# help
PanGraph.main(["help", "build"])

max_stars_repo_name: neherlab/pangraph

max_issues_repo_name: neherlab/pangraph

max_forks_repo_name: neherlab/pangraph


text:
using NeuralThermalODE

# Test data file.
test_data = read_data("Fig_67.csv_Result.csv")

# Domain.
const N = 10
const tmat = 0.00635
const dx = tmat / N

# Boundary conditions.
const sinkT = test_data[end, 2]
const ambT = test_data[1, 2]
u20 = ambT * ones(N)
const tspan = (Float64(test_data[1, 1]), Float64(test_data[end, 1]))

const pr = BuildThermalODE(dx, 450.0, 7850.0, 42.0, sinkT, ambT, 4.0)

# Solve ODE.
p = read_checkpoint("initial_guess.csv")

du = zeros(N)

# Pass additional arguments through closure.
prob = ODEProblem((du, u20, p, t) -> heat_transfer(pr, du, u20, p, t), u20, tspan, p, saveat = 1.0)

data = Iterators.repeated((), 10000)
opt = ADAM(20.0)

# 180 readings with plus/minus 2.0C each
target_loss = 180.0 * 2.2

# Track loss function.
tot_loss = []

training_cnt = -1
cb = function ()
    early_exit = false
    save_curr = false

    global training_cnt += 1

    if training_cnt % 10 == 0
        save_curr = true
    end

    # Save loss function.
    loss = loss_rd()
    append!(tot_loss, loss)

    print("Loss $loss ")
    println("Iteration $training_cnt")

    if loss <= target_loss
        save_curr = true
        early_exit = true
    end

    sol = solve(remake(prob, p = Flux.data(p)), Tsit5(), saveat = 1.0)

    if save_curr == true
        plot(sol.t[:], sol[end, :], label = "Iteration #$training_cnt",
            color = :black,
            linestyle = :dash)
        plot!(test_data[:, 1], test_data[:, 2], label = "Test data",
            color = :black)
        plot!(fmt = :svg,
            xlims = (0, 200),
            ylims = (-200, 50),
            grid = false,
            legend = :topright)
        savefig("temperature_fit_$training_cnt.svg")

        df = DataFrame(x1 = Float64(test_data[1, 1]):Float64(test_data[end, 1]),
                       x2 = Flux.data(p)[1:end-1])
        write_checkpoint("checkpoint_htc_$training_cnt.csv", df)
        plot(Flux.data(p)[1:end-1], color = :black, label = "Iteration #$training_cnt")
        plot!(fmt = :svg,
            xlims = (0, 200),
            grid = false,
            legend = :topleft)
        savefig("checkpoint_htc_$training_cnt.svg")
    end

    if early_exit == true
        Flux.stop()
    end
end

cb()
Flux.train!(loss_rd, [p], data, opt, cb = cb)

# Save last iteration.
sol = solve(remake(prob, p = Flux.data(p)), Tsit5(), saveat = 1.0)
plot(sol.t[:], sol[end, :], label = "Iteration #$training_cnt",
    color = :black,
    linestyle = :dash)
plot!(test_data[:, 1], test_data[:, 2], label = "Test data",
    color = :black)
plot!(fmt = :svg,
    xlims = (0, 200),
    ylims = (-200, 50),
    grid = false,
    legend = :topright)
savefig("last_temperature_fit_$training_cnt.svg")

df = DataFrame(x1 = Float64(test_data[1, 1]):Float64(test_data[end, 1]),
            x2 = Flux.data(p)[1:end-1])
write_checkpoint("last_checkpoint_htc_$training_cnt.csv", df)
plot(Flux.data(p)[1:end-1], color = :black, label = "Iteration #$training_cnt")
plot!(fmt = :svg,
    xlims = (0, 200),
    grid = false,
    legend = :topleft)
savefig("last_checkpoint_htc_$training_cnt.svg")

# Save loss function.
df_loss = DataFrame(x1 = 1:length(tot_loss),
                    x2 = tot_loss)
write_checkpoint("loss.svg", df_loss)
plot(tot_loss, color = :black, label = "Loss",
    fmt = :svg,
    grid = false)
savefig("loss.svg")

max_stars_repo_name: NikStoyanov/NeuralThermalODE

max_issues_repo_name: NikStoyanov/NeuralThermalODE

max_forks_repo_name: NikStoyanov/NeuralThermalODE


text:
type CharLoader
    filepath::AbstractString
    cur_id::Int
    data::AbstractString
    chars::Set{Char}
    id2char::Dict{Int,Char}
    char2id::Dict{Char, Int}
end

function CharLoader(filepath)
    data = readstring(filepath)
    chars = Set(data)
    id2char = Dict(i=>c for (i,c) in enumerate(chars))
    char2id = Dict(c=>i for (i,c) in enumerate(chars))
    return CharLoader(filepath, 1, data, chars, id2char, char2id)
end

function next_batch(loader::CharLoader, batch_size)
    x = zeros(Float32, batch_size, length(loader.chars))
    y = zeros(Int, batch_size)
    local seq::AbstractString
    if loader.cur_id + batch_size > length(loader.data)
        diff = loader.cur_id + batch_size - length(loader.data)
        padding = repeat(" ", diff)
        seq = loader.data[loader.cur_id:length(loader.data)] * padding
        loader.cur_id = 1
    else
        seq = loader.data[loader.cur_id:loader.cur_id+batch_size]
        loader.cur_id += batch_size
    end
    for j=2:length(seq)
        cx = seq[j-1]
        cy = seq[j]
        x[j-1, loader.char2id[cx]] = 1
        y[j-1] = loader.char2id[cy]
    end
    x, y
end


text:
# NOTE: Gradient of incomplete beta function is not implemented.
#       Using skew normal instead.

import Pkg; Pkg.activate(joinpath(@__DIR__, "../"))
using Turing
using StatsFuns
Turing.setadbackend(:tracker)

include(joinpath(@__DIR__, "Helper.jl"))

collect_finite(y) = y[isfinite.(y)]
count_inf(y) = sum(y .== -Inf)

compute_gammaT_star(gammaC, gammaT, p) = p * gammaT + (1 - p) * gammaC
function compute_etaT_star(etaC, etaT, p, gammaC, gammaT, gammaT_star)
  numer = p * etaT * (1 - gammaT) + (1 - p) * etaC * (1 - gammaC)
  denom = 1 - gammaT_star
  return numer / denom
end
function compute_star(etaC, etaT, p, gammaC, gammaT)
  gammaT_star = compute_gammaT_star(gammaC, gammaT, p)
  etaT_star = compute_etaT_star(etaC, etaT, p, gammaC, gammaT, gammaT_star)
  return (gammaT_star, etaT_star)
end

function mixlpdf(loc, scale, skew, y, w)
  lm = logsumexp(Helper.skewnormal_logpdf.(loc', scale', skew', y) .+
                 log.(w)', dims=2)
  return sum(lm)
end

@model cytofdens(yC_finite, yT_finite, nC, nT, zC, zT, K) = begin
  p ~ Beta(.5, .5)
  gammaC ~ Beta(1, 1)
  gammaT ~ Beta(1, 1)
  etaT ~ Dirichlet(K, 1/K)
  etaC ~ Dirichlet(K, 1/K)

  gammaT_star = compute_gammaT_star(gammaC, gammaT, p)
  etaT_star = compute_etaT_star(etaC, etaT, p, gammaC, gammaT, gammaT_star)

  mu ~ filldist(Normal(0, 3), K)
  sigma ~ filldist(LogNormal(0, 0.5), K)
  # nu ~ filldist(LogNormal(3.5, 0.5), K)
  phi ~ filldist(Normal(0, 3), K)

  zC ~ Binomial(nC, gammaC)
  zT ~ Binomial(nT, gammaT_star)

  Turing.acclogp!(_varinfo, mixlpdf(mu, sigma, phi, yC_finite, etaC))
  Turing.acclogp!(_varinfo, mixlpdf(mu, sigma, phi, yT_finite, etaT_star))
  # yC_finite .~ UnivariateGMM(mu, sigma, Categorical(etaC))
  # yT_finite .~ UnivariateGMM(mu, sigma, Categorical(etaT_star))
end

# NUTS
yC = [randn(100) .+ 1; fill(-Inf, 50)]
yT = [randn(100); fill(-Inf, 100)]
@time nuts_chain = begin
    K = 5
    n_samples = 500
    nadapt = 500
    iterations = n_samples + nadapt
    target_accept_ratio = 0.8

    yC_finite = collect_finite(yC)
    yT_finite = collect_finite(yT)
    nC = length(yC)
    nT = length(yT)
    zC = count_inf(yC)
    zT = count_inf(yT)
    sample(cytofdens(yC_finite, yT_finite, nC, nT, zC, zT, K),
           NUTS(nadapt, target_accept_ratio, max_depth=6),
           iterations);
end;

max_stars_repo_name: luiarthur/CytofDensityEstimation

max_issues_repo_name: luiarthur/CytofDensityEstimation

max_forks_repo_name: luiarthur/CytofDensityEstimation


text:
# coding: utf-8

# In[1]:

#coin change is a straight forward dynamic programming problem
#you are given one note and a set of coins of different values
#assuming you have infinite supply of coins of different values
#we need to compute different ways of making change
#the order of the coins doesnt matter in this case
#more details can be found in the following link
# https://www.geeksforgeeks.org/coin-change-dp-7/
#the script can also be done in recursion
# https://github.com/je-suis-tm/recursion-and-dynamic-programming/blob/master/coin%20change%20recursion.jl

# In[2]:

#to solve this problem via tabulation
#we divide one big problem into two sub problems
#the case where coin of η value is excluded in the solutions
#and the case where at least one coin of η value is included
function coin_change(num,choice)

    #create matrix (num+1)*length(choice)
    #the raison d'être is the computation starts from 0 to num
    #0 is when num is perfectly substituted by coins
    tabulation=[[0 for _ in 1:length(choice)] for _ in 1:num+1]

    #initialize
    #when the remain value happens to be η
    #one solution is found
    for i in 1:length(choice)

        tabulation[1][i]=1

    end

    #since we initialize the null case
    #the outerloop starts from 2
    #i actually refers to i-1
    for i in 2:num+1

        for j in 1:length(choice)

            #annoying part of julia
            #if index starts at zero
            #will be a lot easier
            if i-choice[j]>=1

                #the case where at least one coin of η value is included
                #we just need the computation where η is deducted
                include=tabulation[i-choice[j]][j]

            else

                include=0

            end

            #the case where coin of η value is excluded in the solutions
            if j>=2

                exclude=tabulation[i][j-1]

            else

                exclude=0

            end

            #two sub problems merge into a big one
            tabulation[i][j]=exclude+include

        end

    end

    #get the final answer
    return tabulation[num+1][length(choice)]

end

# In[3]:

coin_change(10,[1,2,5])

max_stars_repo_path: coin change dynamic programming.jl

max_issues_repo_path: coin change dynamic programming.jl

max_forks_repo_path: coin change dynamic programming.jl


text:
### series elastic component

# SEC parameters
struct SECParameters
    k::Float64  # SEC stiffness
end
SECParameters(d::Dict{Symbol,Float64}) = SECParameters(d[:k])

# SEC struct
mutable struct SEC
    τ::Float64  # torque
    θ::Float64  # SEC angle
    ω::Float64  # SEC angular velocity
    sec_p::SECParameters # SEC parameters

    SEC(θ₀, ω₀, p::SECParameters) = new(θ₀ * p.k, θ₀, ω₀, p)
end

# calculate SEC torque
torque(sec::SEC) = sec.sec_p.k * sec.θ
torque!(sec::SEC) = sec.τ =  torque(sec)
torque(θ, p::SECParameters) = p.k * θ