text:
import DFControl.Structures: Orbital, Structure, orbital, size, orbsize

import DFControl.Display.Crayons: @crayon_str
import DFControl.Display: dfprintln, dfprint

setup_ω_grid(ωh, ωv, n_ωh, n_ωv, offset = 0.001) = vcat(range(ωh,             ωh + ωv * 1im,     length = n_ωv)[1:end - 1],
                                                        range(ωh + ωv * 1im,     offset + ωv * 1im, length = n_ωh)[1:end - 1],
                                                        range(offset + ωv * 1im, offset,          length = n_ωv))

struct ExchangeKGrid{T,MT} <: AbstractKGrid{T}
    hamiltonian_kgrid::HamiltonianKGrid{T,MT}
    phases::Vector{Complex{T}}
    D::Matrix{Complex{T}}
end

core_kgrid(x::ExchangeKGrid) = core_kgrid(x.hamiltonian_kgrid)
eigvecs(x::ExchangeKGrid) = eigvecs(x.hamiltonian_kgrid)
eigvals(x::ExchangeKGrid) = eigvals(x.hamiltonian_kgrid)

function ExchangeKGrid(hami::TbHami, kpoints::Vector{Vec3{T}}, R = zero(Vec3{T})) where {T}
    D       = ThreadCache(zeros_block(hami))
    hami_kpoints = HamiltonianKGrid(hami, kpoints, x->D .+= x)
    nk = length(hami_kpoints)
    Ds = gather(D)
    return ExchangeKGrid(hami_kpoints, phases(kpoints, R), Array(Ds[Up()] - Ds[Down()]) / nk)
end

# struct ColinGreensFunction{T <: AbstractFloat}
#     ω::T
#     G::ColinMatrix{Complex{T}} # forward part is always spin up, backward spin down
# end

function calc_greens_functions(ω_grid, kpoints, μ::T) where T
    g_caches = [ThreadCache(fill!(similar(eigvecs(kpoints)[1]), zero(Complex{T}))) for i = 1:3]
    Gs = [fill!(similar(eigvecs(kpoints)[1]), zero(Complex{T})) for i = 1:length(ω_grid)]
    function iGk!(G, ω)
        fill!(G, zero(Complex{T}))
        integrate_Gk!(G, ω, μ, kpoints, cache.(g_caches))
    end
    p = Progress(length(ω_grid), 1, "Calculating contour G(ω)...")
    @threads for j = 1:length(ω_grid)
        iGk!(Gs[j], ω_grid[j])
        next!(p)
    end
    return Gs
end

function integrate_Gk!(G::AbstractMatrix, ω, μ, kpoints, caches)
    dim = blockdim(G)
    cache1, cache2, cache3 = caches

    @inbounds for ik = 1:length(kpoints)
    # Fill here needs to be done because cache1 gets reused for the final result too
        fill!(cache1, zero(eltype(cache1)))
        for x = 1:2dim
            cache1[x, x] = 1.0 / (μ + ω - eigvals(kpoints)[ik][x])
        end
    # Basically Hvecs[ik] * 1/(ω - eigvals[ik]) * Hvecs[ik]'

        mul!(cache2, eigvecs(kpoints)[ik], cache1)
        adjoint!(cache3, eigvecs(kpoints)[ik])
        mul!(cache1, cache2, cache3)
        t = kpoints.phases[ik]
        tp = t'
        for i in 1:dim, j in 1:dim
            G[i, j] += cache1[i, j] * t
            G[i + dim, j + dim] += cache1[i + dim, j + dim] * tp
            G[i+dim, j] = cache1[i + dim, j]
            G[i, j+dim] = cache1[i, j+dim]
        end
    end
    G  ./= length(kpoints)
end

function integrate_Gk!(G::ColinMatrix, ω, μ, kpoints, caches)
    dim = size(G, 1)
    cache1, cache2, cache3 = caches

    @inbounds for ik = 1:length(kpoints)
    # Fill here needs to be done because cache1 gets reused for the final result too
        fill!(cache1, zero(eltype(cache1)))
        for x = 1:dim
            cache1[x, x] = 1.0 / (μ + ω - kpoints.hamiltonian_kgrid.eigvals[ik][x])
            cache1[x, x + dim] = 1.0 / (μ + ω - kpoints.hamiltonian_kgrid.eigvals[ik][x + dim])
        end
    # Basically Hvecs[ik] * 1/(ω - eigvals[ik]) * Hvecs[ik]'

        mul!(cache2, kpoints.hamiltonian_kgrid.eigvecs[ik], cache1)
        adjoint!(cache3, kpoints.hamiltonian_kgrid.eigvecs[ik])
        mul!(cache1, cache2, cache3)
        t = kpoints.phases[ik]
        tp = t'
        for i in 1:dim, j in 1:dim
            G[i, j]     += cache1[i, j] * t
            G[i, j + dim] += cache1[i, j + dim] * tp
        end
    end
    G  ./= length(kpoints)
end

function integrate_Gk!(G_forward::ThreadCache, G_backward::ThreadCache, ω, μ, Hvecs, Hvals, R, kgrid, caches)
    dim = size(G_forward, 1)
    cache1, cache2, cache3 = caches

    @inbounds for ik = 1:length(kgrid)
    # Fill here needs to be done because cache1 gets reused for the final result too
        fill!(cache1, zero(eltype(cache1)))
        for x = 1:dim
            cache1[x, x] = 1.0 / (μ + ω - Hvals[ik][x])
        end
    # Basically Hvecs[ik] * 1/(ω - eigvals[ik]) * Hvecs[ik]'
        mul!(cache2, Hvecs[ik], cache1)
        adjoint!(cache3, Hvecs[ik])
        mul!(cache1, cache2, cache3)
        t = exp(2im * π * dot(R, kgrid[ik]))
        G_forward  .+= cache1 .* t
        G_backward .+= cache1 .* t'
    end
    G_forward.caches  ./= length(kgrid)
    G_backward.caches ./= length(kgrid)
end

abstract type Exchange{T <: AbstractFloat} end
Base.eltype(::Exchange{T}) where T = T
Base.eltype(::Type{Exchange{T}}) where T = T

function (::Type{E})(at1::Atom, at2::Atom; site_diagonal::Bool = false) where {E <: Exchange}
    l1 = length(uprange(at1))
    l2 = length(uprange(at2))
    return site_diagonal ? E(zeros(Float64, l1, l2), at1, at2) : E(zeros(Float64, l1, l1), at1, at2)
end

"""
    Exchange2ndOrder{T <: AbstractFloat}

This holds the exhanges between different orbitals and calculated sites.
Projections and atom datablocks are to be found in the corresponding wannier input file.
It turns out the ordering is first projections, then atom order in the atoms datablock.
"""
mutable struct Exchange2ndOrder{T <: AbstractFloat} <: Exchange{T}
    J::Matrix{T}
    atom1::Atom
    atom2::Atom
end

function Base.show(io::IO, e::Exchange)
    dfprint(io, crayon"red", "atom1:", crayon"reset")
    dfprintln(io, "name: $(e.atom1.name), pos: $(e.atom1.position_cryst)")
    dfprint(io, crayon"red", " atom2:", crayon"reset")
    dfprintln(io, "name: $(e.atom2.name), pos: $(e.atom2.position_cryst)")

    dfprint(io, crayon"red", " J: ", crayon"reset", "$(tr(e.J))")
end

"""
    Exchange4thOrder{T <: AbstractFloat}

This holds the exhanges between different orbitals and calculated sites.
Projections and atom datablocks are to be found in the corresponding wannier input file.
It turns out the ordering is first projections, then atom order in the atoms datablock.
"""
mutable struct Exchange4thOrder{T <: AbstractFloat} <: Exchange{T}
    J::Matrix{T}
    atom1::Atom
    atom2::Atom
end

function calc_exchanges(hami,  atoms, fermi::T, ::Type{E} = Exchange2ndOrder;
                        nk::NTuple{3,Int} = (10, 10, 10),
                        R                  = Vec3(0, 0, 0),
                        ωh::T              = T(-30.), # starting energy
                        ωv::T              = T(0.15), # height of vertical contour
                        n_ωh::Int          = 3000,
                        n_ωv::Int          = 500,
                        temp::T            = T(0.01),
                        site_diagonal      = false) where {T <: AbstractFloat,E <: Exchange}
    R_ = Vec3(R...)
    μ               = fermi
    ω_grid          = setup_ω_grid(ωh, ωv, n_ωh, n_ωv)

    exchanges       = E{T}[]
    for at1 in atoms
        for at2 in atoms
            push!(exchanges, E(at1, at2, site_diagonal = site_diagonal))
        end
    end
    kpoints = ExchangeKGrid(hami, uniform_shifted_kgrid(nk...), R_)

    D_ = site_diagonal ? site_diagonalize(kpoints.D, atoms) : kpoints.D
    calc_exchanges!(exchanges, μ, ω_grid, kpoints, D_)

    return exchanges
end

function site_diagonalize(D::Matrix{Complex{T}}, ats::Vector{DFC.Structures.Atom}) where {T}
    Ts = zeros(D)
    Dvals = zeros(T, size(D, 1))
    for at in ats
        t_vals, t_vecs = eigen(Hermitian(D[at]))
        Ts[at] .= t_vecs
        Dvals[range(at)] .=  real.(t_vals)
    end
    return SiteDiagonalD(Dvals, Ts)
end

function calc_exchanges!(exchanges::Vector{<:Exchange{T}},
                                 μ::T,
                                 ω_grid::AbstractArray{Complex{T}},
                                 kpoints  ,
                                 D::Union{Matrix{Complex{T}},SiteDiagonalD{T}}) where T <: AbstractFloat
    dim     = size(kpoints.hamiltonian_kgrid.eigvecs[1])
    d2      = div(dim[1], 2)
    J_caches = [ThreadCache(zeros(T, size(e.J))) for e in exchanges]
    Gs = calc_greens_functions(ω_grid, kpoints, μ)
    @threads for i = 1:length(Gs) - 1
        for (eid, exch) in enumerate(exchanges)
            J_caches[eid] .+= Jω(exch, D, Gs[i], ω_grid[i + 1] - ω_grid[i])
        end
    end
    for (eid, exch) in enumerate(exchanges)
        exch.J = -1e3 / 2π * gather(J_caches[eid])
    end
end

spin_sign(D) = -sign(real(tr(D))) # up = +1, down = -1. If D_upup > D_dndn, onsite spin will be down and the tr(D) will be positive. Thus explaining the - in front of this.
spin_sign(D::Vector) = sign(real(sum(D))) # up = +1, down = -1

perturbation_bubble(::Exchange2ndOrder, D_site1, G_forward, D_site2, G_backward) =
    D_site1 * G_forward * D_site2 * G_backward

perturbation_bubble(::Exchange4thOrder, D_site1, G_forward, D_site2, G_backward) =
    D_site1 * G_forward * D_site2 * G_backward * D_site1 * G_forward * D_site2 * G_backward

@inline function Jω(exch, D, G, dω)
    if size(D, 1) < size(G, 1)
        ra1 = uprange(exch.atom1)
        ra2 = uprange(exch.atom2)
    else
        ra1 = range(exch.atom1)
        ra2 = range(exch.atom2)
    end
    D_site1    = view(D, ra1, ra1)
    D_site2    = view(D, ra2, ra2)
    G_forward  = view(G, exch.atom1, exch.atom2, Down())
    G_backward = view(G, exch.atom2, exch.atom1, Up())
    return spin_sign(D_site1) .* spin_sign(D_site2) .* imag.(perturbation_bubble(exch,
                                                                                D_site1,
                                                                                G_forward,
                                                                                D_site2,
                                                                                G_backward) * dω)
end

@inline function Jω(exch, D::SiteDiagonalD, G, dω)
    s1 = spin_sign(D.values[exch.atom1])
    s2 = spin_sign(D.values[exch.atom2])
    t  = zeros(exch.J)
    G_forward  = D.T[exch.atom1]' * G[exch.atom1, exch.atom2, Down()] * D.T[exch.atom2]
    G_backward = D.T[exch.atom2]' * G[exch.atom2, exch.atom1, Up()] *  D.T[exch.atom1]
    for j in 1:size(t, 2), i in 1:size(t, 1)
        t[i, j] = s1 * s2 * imag(D.values[exch.atom1][i] * G_forward[i, j] * D.values[exch.atom2][j] * G_backward[j, i] * dω)
    end
    return t
end

mutable struct AnisotropicExchange2ndOrder{T <: AbstractFloat} <: Exchange{T}
    J::Matrix{Matrix{T}}
    atom1::Atom
    atom2::Atom
end

AnisotropicExchange2ndOrder(at1::Atom, at2::Atom) =
    AnisotropicExchange2ndOrder{Float64}([zeros(length(range(at1)), length(range(at1))) for i = 1:3,j = 1:3], at1, at2)

function calc_anisotropic_exchanges(hami,  atoms, fermi::T;
                             nk::NTuple{3,Int} = (10, 10, 10),
                             R                  = Vec3(0, 0, 0),
                             ωh::T              = T(-30.), # starting energy
                             ωv::T              = T(0.1), # height of vertical contour
                             n_ωh::Int          = 3000,
                             n_ωv::Int          = 500,
                             temp::T            = T(0.01)) where T <: AbstractFloat

    μ               = fermi
    k_grid          = uniform_shifted_kgrid(nk...)
    ω_grid          = setup_ω_grid(ωh, ωv, n_ωh, n_ωv)
    exchanges       = setup_anisotropic_exchanges(atoms)

    Hvecs, Hvals, D = DHvecvals(hami, k_grid, atoms)

    calc_anisotropic_exchanges!(exchanges, μ, R, k_grid, ω_grid, Hvecs, Hvals, D)
    return exchanges
end

function calc_anisotropic_exchanges!(exchanges::Vector{AnisotropicExchange2ndOrder{T}},
                                   μ::T,
                                   R::Vec3,
                                   k_grid::AbstractArray{Vec3{T}},
                                   ω_grid::AbstractArray{Complex{T}},
                                   Hvecs::Vector{Matrix{Complex{T}}},
                                   Hvals::Vector{Vector{T}},
                                   D::Vector{Vector{Matrix{Complex{T}}}}) where T <: AbstractFloat
    dim      = size(Hvecs[1])
    J_caches = [ThreadCache([zeros(T, size(e.J[i, j])) for i = 1:3, j = 1:3]) for e in exchanges]
    g_caches = [ThreadCache(zeros(Complex{T}, dim)) for i = 1:3]
    G_forward, G_backward = [ThreadCache(zeros(Complex{T}, dim)) for i = 1:2]

    function iGk!(ω)
        fill!(G_forward, zero(Complex{T}))
        fill!(G_backward, zero(Complex{T}))
        integrate_Gk!(G_forward, G_backward, ω, μ, Hvecs, Hvals, R, k_grid, g_caches)
    end

    for j = 1:length(ω_grid[1:end - 1])
        ω   = ω_grid[j]
        dω  = ω_grid[j + 1] - ω
        iGk!(ω)
    # The two kind of ranges are needed because we calculate D only for the projections we care about
    # whereas G is calculated from the full Hamiltonian, the is needed.
        for (eid, exch) in enumerate(exchanges)
            rm  = range(exch.atom1)
            rn  = range(exch.atom2)
            for i in 1:3, j in 1:3 # x,y,z
                J_caches[eid][i, j] .+=  imag.((view(D[eid][i], 1:length(rm), 1:length(rm)) *
                                  view(G_forward, rm, rn)    *
                                  view(D[eid][j], 1:length(rn), 1:length(rn))  *
                                  view(G_backward, rn, rm)) .* dω)
            end
        end
    end

    for (eid, exch) in enumerate(exchanges)
        exch.J = 1e3 / 2π * gather(J_caches[eid])
    end
end

function setup_anisotropic_exchanges(atoms::Vector{Atom})
    exchanges = AnisotropicExchange2ndOrder{Float64}[]
    for (i, at1) in enumerate(atoms), at2 in atoms[i:end]
        push!(exchanges, AnisotropicExchange2ndOrder(at1, at2))
    end
    return exchanges
end

@doc raw"""
    DHvecvals(hami::TbHami{T, Matrix{T}}, k_grid::Vector{Vec3{T}}, atoms::Atom{T}) where T <: AbstractFloat

Calculates $D(k) = [H(k), J]$, $P(k)$ and $L(k)$ where $H(k) = P(k) L(k) P^{-1}(k)$.
`hami` should be the full Hamiltonian containing both spin-diagonal and off-diagonal blocks.
"""
function DHvecvals(hami, k_grid::AbstractArray{Vec3{T}}, atoms::Vector{Atom}) where T <: AbstractFloat

    nk        = length(k_grid)
    Hvecs     = [zeros_block(hami) for i = 1:nk]
    Hvals     = [Vector{T}(undef, blocksize(hami, 1)) for i = 1:nk]
    δH_onsite = ThreadCache([[zeros(Complex{T}, 2length(range(at)), 2length(range(at))) for i = 1:3] for at in atoms])
    calc_caches = [EigCache(block(hami[1])) for i = 1:nthreads()]
    for i = 1:nk
    # for i=1:nk
        tid = threadid()
        # Hvecs[i] is used as a temporary cache to store H(k) in. Since we
        # don't need H(k) but only Hvecs etc, this is ok.
        Hk!(Hvecs[i], hami, k_grid[i])

        for (δh, at) in zip(δH_onsite, atoms)
            rat = range(at)
            lr  = length(rat)
            δh .+= commutator.((view(Hvecs[i], rat, rat),), at[:operator_block].J) # in reality this should be just range(at)
            # δh .+= commutator.(([Hvecs[i][rat, rat] zeros(Complex{T},lr, lr); zeros(Complex{T}, lr, lr) Hvecs[i][div(blocksize(hami, 1), 2) .+ rat, div(blocksize(hami, 1), 2) .+ rat]],), at[:operator_block].J) #in reality this should be just range(at)
        end
        eigen!(Hvals[i], Hvecs[i], calc_caches[tid])
    end
    return Hvecs, Hvals, gather(δH_onsite) ./ nk
end

commutator(A1, A2) where T = A1 * A2 - A2 * A1

function totocc(Hvals, fermi::T, temp::T) where T
    totocc = zero(Complex{T})
    for i = 1:length(Hvals)
        totocc += gather(1 ./ (exp.((Hvals[i] .- fermi) ./ temp) .+ 1))
    end
    return totocc / length(Hvals)
end

"Generates a Pauli σx matrix with the dimension that is passed through `n`."
σx(::Type{T}, n::Int) where {T} =
    kron(diagm(0 => ones(T, div(n, 2))), SMatrix{2,2}(0, 1, 1, 0)) / 2
σx(n::Int) = σx(Float64, n)

"Generates a Pauli σy matrix with the dimension that is passed through `n`."
σy(::Type{T}, n::Int) where {T} =
    kron(diagm(0 => ones(T, div(n, 2))), SMatrix{2,2}(0, -1im, 1im, 0)) / 2
σy(n::Int) = σy(Float64, n)

"Generates a Pauli σz matrix with the dimension that is passed through `n`."
σz(::Type{T}, n::Int) where {T} =
    kron(diagm(0 => ones(T, div(n, 2))), SMatrix{2,2}(1, 0, 0, -1)) / 2
σz(n::Int) = σz(Float64, n)

for s in (:σx, :σy, :σz)
    @eval @inline $s(m::AbstractArray{T}) where {T} =
    $s(T, size(m, 1))
    @eval $s(m::TbHami) = $s(block(m[1]))
end

σx(m::ColinMatrix{T}) where {T} =
    zeros(m)

σy(m::ColinMatrix{T}) where {T} =
    zeros(m)

σz(m::ColinMatrix{T}) where {T} =
    ColinMatrix(diagm(0 => ones(T, size(m, 1))), diagm(0 => -1.0 .* ones(T, size(m, 1))))

# function calc_onsite_spin(kpoints::AbstractArray{<:KPoint{T}}, atom, fermi=0.0) where {T}
    # 	S = Vec3(σx(kpoints[1].eigvecs[atom]),
    # 	         σy(kpoints[1].eigvecs[atom]),
    # 	         σz(kpoints[1].eigvecs[atom]))

    # 	S_out = zero(Vec3{T})
    # 	for k in kpoints
    # 		for (i, v) in enumerate(k.eigvals)
    # 			if i - fermi <= 0.0
    # 				vec = k.eigvecs[:, i][atom]
    # 			    S_out += real((vec',) .* S .* (vec,))
    # 			end
    # 		end
    # 	end
    # 	return S_out./length(kpoints)
# end

# Old stuff
# function heisenberg_energy(moments::Vector{<:Vec3}, exchanges::Vector{Matrix{T}}, H, anisotropy) where T
#     energy = zero(T)
#     nexch  = length(exchanges)
#     nmom   = length(moments)
#     for i = 1:nmom
#         for j = i+1:nmom
#             for k = 1:nexch
#                 energy += exchanges[k][i,j] * dot(moments[i], moments[j]) + H * moments[i][1] + anisotropy * moments[i][3]^2
#             end
#         end
#     end
#     return energy
# end

# function AFMmap(structure, Rcryst)
#     R = cell(structure)' * Rcryst/2
#     atoms = atoms(structure)
#     map1 = Dict{Atom, Atom}()
#     for at1 in atoms, at2 in atoms
#         if bondlength(at2, at1, R) < 1.0e-7
#             map1[at1] = at2
#         end
#     end
#     return map1
# end

# function rs(R, r1, r2, r3, r4)
#     sig =
#     if sign(R) < 0
#         return r4, r3, r2, r1
#     else
#         return r1, r2, r3, r4
#     end
# end

# function symmetrize!(tb_hamis::NTuple{2, Vector{TbBlock{T, LT}}}, structure::AbstractStructure{T, LT}) where  {T, LT<:Length{T}}
#     forwardmap = AFMmap(structure, Vec3(1,0,0))
#     Hup = getfirst(x -> x.R_cryst == Vec3(0,0,0), tb_hamis[1]).block
#     Hdn = getfirst(x -> x.R_cryst == Vec3(0,0,0), tb_hamis[2]).block

#     for (at1, at2) in forwardmap, (at3, at4) in forwardmap
#         for (r1, r2) in zip(range.(projections(at1)), range.(projections(at2))), (r3, r4) in zip(range.(projections(at3)), range.(projections(at4)))
#             Hup[r1, r3] .= (Hup[r1, r3] .+ Hdn[r2, r4]) ./ 2
#             Hdn[r2, r4] .= Hup[r1, r3]

#             for R=-1:2:1
#                 Hu = getfirst(x -> x.R_cryst == Vec3(R,0,0), tb_hamis[1]).block
#                 Hd = getfirst(x -> x.R_cryst == Vec3(R,0,0), tb_hamis[2]).block
#                 r1_, r2_, r3_, r4_ = rs(R, r1, r2, r3, r4)
#                 Hu[r2_, r3_] .= (Hu[r2_, r3_] .+ Hdn[r1_, r4_]) ./ 2
#                 Hdn[r1_, r4_] .= Hu[r2_, r3_]
#                 Hd[r2_, r3_]  .= (Hd[r2_, r3_] .+ Hup[r1_, r4_]) ./ 2
#                 Hup[r1_, r4_] .= Hd[r2_, r3_]
#             end
#         end
#     end
# end


text:
using AdventOfCode2018
using AdventOfCode2018.Day10

input = get_aoc_input(10)
pos, vel = parse_input(input)
find_min_area(pos, vel)

ext: jl


text:
using DataFrames, LinearAlgebra, MixedModels, RData, StatsBase, Test

if !@isdefined(dat) || !isa(dat, Dict{Symbol, DataFrame})
    dat = Dict(Symbol(k) => v for (k, v) in load(joinpath(dirname(@__FILE__), "dat.rda")))
end

@testset "contra" begin
    contraception = dat[:Contraception]
    contraception[:a2] = abs2.(contraception[:a])
    contraception[:urbdist] = string.(contraception[:urb], contraception[:d])
    gm0 = fit!(GeneralizedLinearMixedModel(@formula(use ~ 1 + a + a2 + urb + l + (1 | urbdist)), contraception,
                    Bernoulli()), fast = true);
    @test gm0.lowerbd == [0.]
    @test isapprox(getθ(gm0)[1], 0.5720734451352923, atol=0.001)
    @test isapprox(deviance(gm0,true), 2361.657188518064, atol=0.001)
    gm1 = fit(GeneralizedLinearMixedModel, @formula(use ~ 1 + a + a2 + urb + l + (1 | urbdist)), contraception,
        Bernoulli());
    @test isapprox(gm1.θ[1], 0.573054, atol=0.005)
    @test lowerbd(gm1) == push!(fill(-Inf, 7), 0.)
    @test isapprox(deviance(gm1,true), 2361.57129, rtol=0.00001)
    @test isapprox(loglikelihood(gm1), -1180.78565, rtol=0.00001)
    @test StatsBase.dof(gm0) == length(gm0.β) + length(gm0.θ)
    @test StatsBase.nobs(gm0) == 1934
    fit!(gm0, fast=true, nAGQ=7)
    @test isapprox(deviance(gm0, 7), 2360.9838, atol=0.001)
    fit!(gm0, nAGQ=7)
    @test isapprox(deviance(gm0, 7), 2360.9002, atol=0.001)
    @test gm0.β == gm0.beta
    @test gm0.θ == gm0.theta
    @test isnan(gm0.σ)
    @test gm0.formula == @formula(use ~ 1 + a + a2 + urb + l + (1 | urbdist))
    @test length(gm0.y) == size(gm0.X, 1)
    @test :θ in propertynames(gm0)
    # the next three values are not well defined in the optimization
    #@test isapprox(logdet(gm1), 75.7275, atol=0.1)
    #@test isapprox(sum(abs2, gm1.u[1]), 48.4747, atol=0.1)
    #@test isapprox(sum(gm1.resp.devresid), 2237.344, atol=0.1)
    show(IOBuffer(), gm1)
end

@testset "cbpp" begin
    cbpp = dat[:cbpp]
    cbpp[:prop] = cbpp[:i] ./ cbpp[:s]
    gm2 = fit!(GeneralizedLinearMixedModel(@formula(prop ~ 1 + p + (1 | h)), cbpp, Binomial(),
                    wt = cbpp[:s]));
    @test isapprox(deviance(gm2,true), 100.09585619324639, atol=0.0001)
    @test isapprox(sum(abs2, gm2.u[1]), 9.723175126731014, atol=0.0001)
    @test isapprox(logdet(gm2), 16.90113, atol=0.0001)
    @test isapprox(sum(gm2.resp.devresid), 73.47179193718736, atol=0.001)
    @test isapprox(loglikelihood(gm2), -92.02628186555876, atol=0.001)
    @test isnan(sdest(gm2))
    @test varest(gm2) == 1
end

@testset "verbagg" begin
    verbagg = dat[:VerbAgg]
    gm3 = fit(GeneralizedLinearMixedModel, @formula(r2 ~ 1 + a + g + b + s + (1 | id) + (1 | item)),
         verbagg, Bernoulli());
    @test isapprox(deviance(gm3), 8151.40, atol=0.01)
    @test lowerbd(gm3) == vcat(fill(-Inf, 6), zeros(2))
    @test fitted(gm3) == predict(gm3)
    # these two values are not well defined at the optimum
    @test isapprox(sum(x -> sum(abs2, x), gm3.u), 273.31563469936697, atol=0.1)
    @test isapprox(sum(gm3.resp.devresid), 7156.558983084621, atol=0.1)
end

@testset "grouseticks" begin
    gm4 = fit!(GeneralizedLinearMixedModel(@formula(t ~ 1 + y + ch + (1|i) + (1|b) + (1|l)),
              dat[:grouseticks], Poisson()), fast=true)  # fails in pirls! with fast=false
    @test isapprox(deviance(gm4), 851.4046, atol=0.001)
    @test lowerbd(gm4) == vcat(fill(-Inf, 4), zeros(3))
    # these two values are not well defined at the optimum
    #@test isapprox(sum(x -> sum(abs2, x), gm4.u), 196.8695297987013, atol=0.1)
    #@test isapprox(sum(gm4.resp.devresid), 220.92685781326136, atol=0.1)
end


text:
using HTTP
using JSON

# Constantes
const ADDR = "0.0.0.0"
const PORT = 3000

const MAX_CATS = 10

const BREEDS = ["Abyssinian",
    "Australian Mist", "Balinese",
    "Bengal", "Birman",
    "Bombay", "British Shorthair",
    "Burmese", "Burmilla",
    "Cornish Rex", "Devon Rex",
    "Egyptian Mau", "Exotic Shorthair",
    "Japanese Bobtail", "Korat",
    "La Perms", "Maine Coon",
    "Manx", "Norwegian Forest",
    "Ocicat", "Oriental Shorthair",
    "Persian Longhair", "Ragdoll",
    "Russian Blue", "Scottish Fold",
    "Siamese", "Siberian Forest",
    "Singapura", "Snowshoe",
    "Somali", "Sphynx",
    "Tiffanie", "Tonkinese",
    "Turkish Van"]
const STATUSES = ["playing", "sleeping", "eating", "drinking", "outside", "unknown"]

"Génère et retourne un chat"
function getcat(agemax=24)
    return Dict("breed" => rand(BREEDS),
                "age" => string(rand(1:agemax)),
                "status" => rand(STATUSES))
end

"Retourne `n` chats sous forme JSON"
function getcats(n, agemax=24)
    return json([getcat() for i in 1:n])
end

"Lance le serveur et écoute sur le port (3000 par défaut)"
function main()
    println("Server listening on port ", PORT)
    HTTP.listen(ADDR, PORT) do http::HTTP.Stream
        HTTP.setstatus(http, 200)
        HTTP.setheader(http, "Content-Type" => "application/json")

        HTTP.startwrite(http)
        write(http, getcats(rand(1:MAX_CATS)))
        HTTP.closewrite(http)
    end
end

main()


text:
for steps in 50:50:100
    results = mean([binaryTokenChallenge(steps, 50, 50, 20, [Benchmarks.accuracyBoost, Benchmarks.accuracyBalanceCorrelation, Benchmarks.registryMean]) for i in 1:40])
    boost = results[1]
    correlation = results[2]
    registry = results[3]
    println("$steps, $boost, $correlation, $registry")
end

max_stars_repo_name: mlodykaz/tcr-simulation

max_issues_repo_name: mlodykaz/tcr-simulation

max_forks_repo_name: mlodykaz/tcr-simulation


text:
struct Size{T}
  height::T
  width::T
end
Size(size::T) where T<:Int = Size(size, size)

width(s::Size) = s.width
height(s::Size) = s.height

"""
genpatch!(patch::AbstractArray, set::Vector)

Fill an array with random values selected from a given set
"""
function genpatch!(patch::AbstractArray, set::Vector)
  for i in 1:length(patch)
    @inbounds patch[i] = rand(set)
  end
  return patch
end
genpatch!(patch::AbstractArray) = genpatch!(patch, UInt8[0,255])

"Return an array of size n x n filled with random values
selected from a given set"
function genpatch(set::Vector, n::Integer)
  res = Matrix{eltype(set)}(undef, n,n)
  genpatch!(res, set)
end
"Return an array of size n x n filled with random values of 0/1"
genpatch(n::Integer) = genpatch(UInt8[0,255], n)

"""
Return an Array of given size, filled with a specified number of
randomly filled patches of a given size of which a certain amount are
identical
"""
function gensample(T::Type, size::Size, patchsize::Size,
                   nbpatches::Int, nbsame::Int)
  @assert nbpatches >= nbsame "Number of identical patches has to be <=
            number of total patches!"
  # Generate image
  # res = zeros(T, width(size), height(size))
  res = fill(convert(T, 128), (width(size), height(size)))
  # Gather non-overlapping positions for patches
  patchpos = nonoverlappingrandrect(width(size), height(size),
                width(patchsize), height(patchsize), nbpatches)
  # Fill in the identical patches
  patch = genpatch(width(patchsize))
  for i in 1:nbsame
    pos = pop!(patchpos)
    paste!(res, patch, left(pos), top(pos))
  end
  # Fill in non-identical patches
  while !isempty(patchpos)
    genpatch!(patch)
    pos = pop!(patchpos)
    paste!(res, patch, left(pos), top(pos))
  end
  return res
end

gensample(size::Size, patchsize::Size,
  nbpatches::Int, nbsame::Int) =
    gensample(UInt8, size, patchsize, nbpatches, nbsame)


text:
use HTTP

function make_API_call(url)
    try
        response = HTTP.request("GET", url)
        return String(response.body)
    catch e
        return "Error occured: $e"
    end
end


text:
# The purpose of this script is to compare the filtering capability of the single snapshots
# estimates and the gas filter on a given

#function that Generates the dynamical fitnesses and estimates single snapshot
using DynNets,JLD
 using PyCall; pygui(:qt); using PyPlot
 Nsample_2 = 3
 save_fold = "./data/estimatesTest/asympTest/"
 @load(save_fold*"FilterTestDynamicIncreasingSize_extraLarge_denseAndSparse_$(Nsample_2)_scaling.jld", dynFitDgp,
        dynDegsSam,estFitSS,indsTVnodes,storeGasPar,rmseSSandGas,filFitGas,convGasFlag)

 dynDegsSam_2 = dynDegsSam
 rmseSSandGas_2 = rmseSSandGas
 indsTVnodes_2 = indsTVnodes
 Nvals_2 = [round(Int,length(dynDegsSam[i,1,1][:,1])/2) for i=1:size(dynDegsSam)[1]]

Nsample= 10
 @load(save_fold*"FilterTestDynamicIncreasingSize_large_denseAndSparse_$(Nsample)_scaling.jld", dynFitDgp,
       dynDegsSam,estFitSS,indsTVnodes,storeGasPar,rmseSSandGas,filFitGas,convGasFlag)

 Nvals = [round(Int,length(dynDegsSam[i,1,1][:,1])/2) for i=1:size(dynDegsSam)[1]]
 Nsample = size(dynDegsSam)[2]
 T = length(dynDegsSam[1,1,1][1,:])
 Nvals_tot = [Nvals;Nvals_2]

# TWO PLOTS FOR THE PAPER

close("all")
 figure(figsize=(9,7))
 lineW = 5
 ax = axes()
 yPlot= [mean( [sum(dynDegsSam[n,s,1]/(T*Nvals[n]*(Nvals[n]-1))) for s=1:Nsample]) for n=1:length(Nvals)]
 append!(yPlot, [mean( [sum(dynDegsSam_2[n,s,1]/(T*Nvals_2[n]*(Nvals_2[n]-1))) for s=1:Nsample_2]) for n=1:length(Nvals_2)])
 plot(Nvals_tot,yPlot,".-k" ,linewidth=lineW)
 yPlot =  [mean( [sum(dynDegsSam[n,s,2]/(T*Nvals[n]*(Nvals[n]-1))) for s=1:Nsample]) for n=1:length(Nvals)]
 append!(yPlot,[mean( [sum(dynDegsSam_2[n,s,2]/(T*Nvals_2[n]*(Nvals_2[n]-1))) for s=1:Nsample_2]) for n=1:length(Nvals_2)])
 plot(Nvals_tot,yPlot, ".--k",linewidth=lineW)
 xlabel("Number of Nodes ",fontsize=38)
 ylabel("Density",fontsize=38)
 grid()
 ax[:tick_params]("both", labelsize=18)
 legend(["Sparse"; "Dense"],fontsize = 28)
 tight_layout()

close("all")
 figure(figsize=(9,7))
 d=1
 ax = axes()
 plot(Nvals_tot,append!([mean(rmseSSandGas[i,d][:,1,:]) for i=1:length(Nvals)],[mean(rmseSSandGas_2[i,d][:,1,:]) for i=1:length(Nvals_2)]),".-b",linewidth=lineW)
 plot(Nvals_tot,append!([mean(rmseSSandGas[i,d][:,2,:]) for i=1:length(Nvals)],[mean(rmseSSandGas_2[i,d][:,2,:]) for i=1:length(Nvals_2)]),".-r",linewidth=lineW)
 d=2
 ax = axes()
 plot(Nvals_tot,append!([mean(rmseSSandGas[i,d][:,1,:]) for i=1:length(Nvals)],[mean(rmseSSandGas_2[i,d][:,1,:]) for i=1:length(Nvals_2)]),".--b",linewidth=lineW)
 plot(Nvals_tot,append!([mean(rmseSSandGas[i,d][:,2,:]) for i=1:length(Nvals)],[mean(rmseSSandGas_2[i,d][:,2,:]) for i=1:length(Nvals_2)]),".--r",linewidth=lineW)

 #the commented lines plot the rmse only on the TV nodes
 #plot(Nvals_tot,append!([mean(rmseSSandGas[i,d][!indsTVnodes[i,1],1,:]) for i=1:length(Nvals)],[mean(rmseSSandGas_2[i,d][!indsTVnodes_2[i,1],1,:]) for i=1:length(Nvals_2)]),".--b")
 #plot(Nvals_tot,append!([mean(rmseSSandGas[i,d][!indsTVnodes[i,1],2,:]) for i=1:length(Nvals)],[mean(rmseSSandGas_2[i,d][!indsTVnodes_2[i,1],2,:]) for i=1:length(Nvals_2)]),".--r")
 xlabel("Number of Nodes",fontsize=38)
 ylabel("Average RMSE",fontsize=38)
 grid()
 ax[:tick_params]("both", labelsize=18)
 legend(["Beta Model Sparse";"SD-Beta Model Sparse";"Beta Model Dense";"SD-Beta Model Dense"],fontsize = 28)
 tight_layout()

#THE PLOTS FOR THE PAER FINISH HERE

for i=1:3
 println(storeGasPar[1,i,:][1][3])
end

indN,s,d = 1,1,1
 N = Nvals[indN]

 indN = 1
 indsPlot = [1,Nvals[indN]-1]
T =250
close("all")
 ax = axes()
 plot(dynFitDgp[indN,s,d][indsTVnodes[indN,d],:][indsPlot,:]')
 plot(estFitSS[indN,s,d][indsTVnodes[indN,d],:][indsPlot,:]',".")
 ylabel(L"\overleftarrow{\theta}",fontsize = 38)
 xlabel("Time",fontsize=38)
 ylim([-4.5,1.5])
 ylimvals = ylim()
 xlim([0,T])
 grid()
 ax[:tick_params]("both", labelsize=18)
 title("Single Snapshots N = $(Nvals[indN])",fontsize = 38)

close("all")
   ax = axes()
    plot(dynFitDgp[indN,s,d][indsTVnodes[indN,d],:][indsPlot,:]')
   plot(filFitGas[indN,s,d][indsTVnodes[indN,d],:][indsPlot,:]',".")
   ylabel(L"\overleftarrow{\theta}",fontsize = 38)
   xlabel("Time",fontsize=38)
   xlim([0,T])
   ylim(ylimvals)
   grid()
   ylimvals = ylim()
   ax[:tick_params]("both", labelsize=18)
   title("Gas N = $(Nvals[indN])",fontsize = 38)

close("all")
 ax = axes()
 plot(dynFitDgp[indN,s,d][1:N,:][indsTVnodes[indN,d][1:N],:]')
 ylabel(L"\overleftarrow{\theta}",fontsize = 38)
 xlabel("Time",fontsize=38)
 grid()
 xlim([0,T])
  ylim(ylimvals)
 ax[:tick_params]("both", labelsize=18)

close("all")
 ax = axes()
 plot(dynFitDgp[indN,s,d][indsTVnodes[indN,d],:][[1,17],:]')
 ylabel(L"\overleftarrow{\theta}",fontsize = 38)
 xlabel("Time",fontsize=38)
 xlim([0,T])
 ylim(ylimvals)
 grid()
 ylimvals = ylim()
 ax[:tick_params]("both", labelsize=18)

fig, ax = subplots(2,2)
 for f=1:2
     for d=1:2
         ax[1,d][:plot](Nvals,[mean([sum(dynDegsSam[n,s,d]/(T*Nvals[n]*(Nvals[n]-1))) for s=1:Nsample]) for n=1:length(Nvals)] )
         ax[1,d][:set_ylim]([0,1])
         ax[2,d][:plot](Nvals,[mean(rmseSSandGas[i,d][indsTVnodes[i,d],f,:]) for i=1:length(Nvals)])
         ax[2,d][:plot](Nvals,[mean(rmseSSandGas[i,d][.!indsTVnodes[i,d],f,:]) for i=1:length(Nvals)])
          ax[2,d][:set_ylim]([0.1,0.5])
     end
     for i=1:length(ax[:,1]),j=1:length(ax[1,:]) ax[i,j][:grid]() end
 end
ax[2,1][:plot](estFitSS[indsTVnodes,:]',".")

storeGasPar[5,:,:]
rmseSSandGas[i,d][indsTVnodes[i,d],1,:]
rmseSSandGas[i,d][indsTVnodes[i,d],2,:]

indN,s,d = 2,2,2
 parGas,~ =  DynNets.estimateTarg(modGasDirBin1;SSest = estFitSS[indN,s,d])
 degsIO_T = dynDegsSam[indN,s,d]
 modGasDirBin1 = DynNets.SdErgmDirBin1(degsIO_T)
 parGas,~ =  DynNets.estimateTarg(modGasDirBin1;SSest = estFitSS[indN,s,d])
 gasFiltFit,~ = DynNets.score_driven_filter_or_dgp( DynNets.SdErgmDirBin1(degsIO_T),[parGas[1];parGas[2];0.05])
 sqrt.(meanSq((dynFitDgp[indN,s,d] - gasFiltFit).^2,2))

plot(gasFiltFit')
plot(dynFitDgp[indN,s,d]')
plot(dynDegsSam[indN,s,d]')
close("all")

for i=[1,20]
    println(i)
end
#
degsDgp = zeros(dynFitDgp)
 for t = 1:T degsDgp[:,t] =  expValStats(fooErgmDirBin1,dynFitDgp[:,t] ) end
 indPlot = indsTVnodes# [1,round(Int,N/2),N]
 fig, ax = subplots(2,1)
 ax[1,1][:plot](dynFitDgp[indPlot,:]')
 ax[2,1][:plot](degsDgp[indPlot,:]')
sampledNets = sampl(fooSnapSeqNetDirBin1,dynFitDgp;Nsample= Nsample)
sampledDegsIO = zeros(Nsample,2N,T)
estSSFit = zeros(Nsample,2N,T)
# compute the degrees and estimate single snapshot
for n=1:Nsample
    sampledDegsIO[n,:,:] = [sumSq(sampledNets[n,:,:,:],2);sumSq(sampledNets[n,:,:,:],1)]
    estSSFit[n,:,:] = estimate(fooSnapSeqNetDirBin1;degsIO_T= sampledDegsIO[n,:,:] )
end

indPlot = [1,round(Int,N/2),N]
 fig, ax = subplots(2,1)
 ax[1,1][:plot](estSSFit[1,indPlot,:]')
 ax[2,1][:plot](sampledDegsIO[1,indPlot,:]')

##

 ##

max_stars_repo_name: domenicodigangi/SDERGM

max_issues_repo_name: domenicodigangi/SDERGM

max_forks_repo_name: domenicodigangi/SDERGM


text:
function render!(
    rendermodel::RenderModel,
    veh::Vehicle,
    color::Colorant=RGB(rand(), rand(), rand())
    )

    p = veh.state.posG
    add_instruction!(rendermodel, render_vehicle, (p.x, p.y, p.θ, veh.def.length, veh.def.width, color))
    rm
end


text:
"""
See also: [`ContinuousImpulse`](@ref), [`frequency_to_time`](@ref), [`DiscreteGaussianImpulse`](@ref)

    DiscreteImpulse{T<:AbstractFloat}

A struct used to represent a numerical impulse. Only the fields: `in_freq` which is the frequency response vector, and the frequency vector `ω` are required to use this struct to use in the function `frequency_to_time`.
"""
struct DiscreteImpulse{T<:AbstractFloat}
    t::Vector{T}
    in_time::Vector{Complex{T}}
    ω::Vector{T}
    in_freq::Vector{Complex{T}}
    function DiscreteImpulse{T}(t_vec::AbstractArray{T}, in_time::AbstractVector{Complex{T}}, ω_vec::AbstractArray{T}, in_freq::AbstractVector{Complex{T}}, do_self_test=true) where {T}
        impulse = new{T}(t_vec,collect(in_time),ω_vec,collect(in_freq))
        if do_self_test self_test(impulse) end
        return impulse
    end
end

"""
Check that the discrete impulse vectors are the right sizes
"""
function self_test(impulse::DiscreteImpulse{T}) where {T}

    right_size = length(impulse.t) == length(impulse.in_time)
    right_size = right_size && length(impulse.ω) == length(impulse.in_freq)

    return right_size
end

function DiscreteImpulse(t_vec::AbstractArray{T}, in_time::Union{AbstractArray{Complex{T}},AbstractArray{T}},
        ω_vec::AbstractArray{T} = t_to_ω(t_vec),
        in_freq::Union{AbstractArray{Complex{T}},AbstractArray{T}} = Complex{T}[];
    kws... ) where {T}

    if isempty(in_freq)
        in_freq = collect(time_to_frequency(Vector{T}(in_time), Vector{T}(t_vec), ω_vec; kws...))[:]
    end

    return DiscreteImpulse{T}(
        Vector{T}(t_vec), Vector{Complex{T}}(in_time), Vector{T}(ω_vec), Vector{Complex{T}}(in_freq)
    )
end

"""
See also: [`DiscreteImpulse`](@ref), [`frequency_to_time`](@ref)

    ContinuousImpulse{T<:AbstractFloat}

A struct used to represent an analytic impulse function. Has two fields: `in_time` a function of time `t`, and `in_freq` a function of the angular frequency `ω`. `in_freq` should be the Fourier transform of `in_time`, though this is not enforced.

We use the Fourier transform convention:
F(ω) =  ∫ f(t)*exp(im*ω*t) dt,
f(t) = (2π)^(-1) * ∫ F(ω)*exp(-im*ω*t) dω.

An impluse f(t) is convoluted in time with the field u(t), however we avoid the convolution by working with the fourier transform F(ω) of the impulse f(t), which results in

frequency to time: (2π)^(-1) * ∫ F(ω)*U(ω)*exp(-im*ω*t) dω
"""
struct ContinuousImpulse{T<:AbstractFloat}
    in_time::Function
    in_freq::Function
    # Enforce that the Types are the same
    function ContinuousImpulse{T}(in_time::Function, in_freq::Function, do_self_test=true) where {T}
        impulse = new{T}(in_time,in_freq)
        if do_self_test self_test(impulse) end
        return impulse
    end
end

"""
    continuous_to_discrete_impulse(impulse::ContinuousImpulse, t_vec, ω_vec = t_to_ω(t_vec); t_shift = 0.0, ω_shift = 0.0)

Returns a [`DiscreteImpulse`](@ref) by sampling `impulse`. The signal can be shifted in time and frequency by choosing `t_shit` and `ω_shift`.
"""
function continuous_to_discrete_impulse(impulse::ContinuousImpulse{T}, t_vec::AbstractArray{T}, ω_vec::AbstractArray{T} = t_to_ω(t_vec); t_shift::T = zero(T), ω_shift::T = zero(T)) where {T}

    # shift in frequency then in time
    in_time = impulse.in_time.(t_vec .- t_shift) .* cos.(ω_shift .* (t_vec .- t_shift))

    # shift in frequency then in time
    in_freq = (impulse.in_freq.(ω_vec .- ω_shift) .+ impulse.in_freq.(ω_vec .+ ω_shift)) .* exp.(im .* ω_vec .* t_shift) ./ T(2)

    return DiscreteImpulse(t_vec, in_time, ω_vec, in_freq)
end

"""
Check that the continuous impulse functions return the correct types
"""
function self_test(impulse::ContinuousImpulse{T}) where {T}

    t = one(T)
    impulse.in_time(t)::Union{T,Complex{T}}

    ω = one(T)
    impulse.in_freq(ω)::Union{T,Complex{T}}

    return true
end

import Base.(+)
function +(i1::ContinuousImpulse{T},i2::ContinuousImpulse{T})::ContinuousImpulse{T} where {T}
    in_time(t) = i1.in_time(t) + i2.in_time(t)
    in_freq(ω) = i1.in_freq(ω) + i2.in_freq(ω)

    ContinuousImpulse{T}(in_time,in_freq)
end

import Base.(*)
function *(a,impulse::ContinuousImpulse{T})::ContinuousImpulse{T} where {T}

    if typeof(one(Complex{T})*a) != Complex{T}
        throw(DomainError(a, "Multiplying impulse by $a would cause impulse type to change, please explicitly cast $a to same type as RegularSource ($T)"))
    end

    in_time(t) = a*impulse.in_time(t)
    in_freq(ω) = a*impulse.in_freq(ω)

    ContinuousImpulse{T}(in_time,in_freq)
end

*(i::ContinuousImpulse,a) = *(a,i)

"""
    TimeDiracImpulse(t0::T)

Dirac Delta function of unit area in the time domain centred at t=t0.

Warning: in the time domain this is a singuarity and so may lead to unexpected behaviour.
"""
function TimeDiracImpulse(t0::T) where {T<:AbstractFloat}
    in_time(t::T) = (t==t0) ? T(Inf) : zero(T)
    in_freq(ω::T) = exp(im*ω*t0)
    ContinuousImpulse{T}(in_time, in_freq, false)
end

function DiscreteTimeDiracImpulse(t0::T, t_vec::AbstractArray{T}, ω_vec::AbstractArray{T} = t_to_ω(t_vec)) where {T}

    return continuous_to_discrete_impulse(TimeDiracImpulse(t0), t_vec, ω_vec)
end

"""
Dirac Delta function of unit area in the frequency domain centred at ω=ω0.

Warning: in frequency space this is a singuarity and so may lead to unexpected behaviour.
"""
function FreqDiracImpulse(ω0::T, dω::T = one(T)) where {T<:AbstractFloat}
    in_time(t::T) = exp(-im*ω0*t) / (T(2)*π)
    in_freq(ω::T) = (ω==ω0) ? T(Inf) : zero(ω)
    ContinuousImpulse{T}(in_time, in_freq, false)
end

"""
See also: [`ContinuousImpulse`](@ref), [`TimeDiracImpulse`](@ref)

    GaussianImpulse(maxω[; σ = 3.0/maxω^2])

Returns a gaussian impulse function, which in the frequency domain is `exp(-σ*ω^2)*(2sqrt(σ*pi))`.
"""
function GaussianImpulse(maxω::T; σ::T = T(3.0)/maxω^2, t_shift = zero(T), ω_shift = zero(T)) where T<:AbstractFloat
    # shift in frequency then in time
    in_time(t::T) = exp(-(t - t_shift)^2 / (4σ)) * cos(ω_shift * (t - t_shift))

    # shift in frequency then in time
    gauss_freq(ω::T) = exp(-σ*ω^2) * (2sqrt(σ*pi))
    in_freq(ω::T) = (gauss_freq(ω - ω_shift) + gauss_freq(ω + ω_shift)) * exp.(im * ω * t_shift) / 2

    ContinuousImpulse{T}(in_time, in_freq)
end

"""
See also: [`ContinuousImpulse`](@ref), [`TimeDiracImpulse`](@ref)

    DiscreteGaussianImpulse(t_vec[, ω_vec])

Returns a discretised gaussian impulse.
"""
function DiscreteGaussianImpulse(t_vec::AbstractArray{T}, ω_vec::AbstractArray{T} = t_to_ω(t_vec);
        kws...) where {T}

    return continuous_to_discrete_impulse(GaussianImpulse(one(T); kws...), t_vec, ω_vec)
end


text:
using QuantumESPRESSOCommands: pw, set_path, get_path
using Test

@testset "Test `pw`" begin
    @test pw("scf.in", "scf.out"; np = 8, nimage = 4, npool = 2).cmd.cmd.exec ==
          ["mpiexec", "-np", "8", "pw.x", "-nimage", "4", "-npool", "2"]
    @test pw("/home/test/ph.in", "/home/test/ph.out").cmd.cmd.dir == "/home/test"
end

@testset "Test `set_path`" begin
    get_path("pw") == "pw.x"
    touch("pw2.x")
    set_path("pw", "pw2.x")
    @test pw("scf.in", "scf.out"; np = 8).cmd.cmd.exec == ["mpiexec", "-np", "8", "pw2.x"]
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 7697c424-7277-11eb-105d-150285a535c4
begin
	using Pkg
	Pkg.activate(mktempdir())
	Pkg.add("Yao")
	Pkg.add("YaoPlots")
	Pkg.add("YaoExtensions")
end

# ╔═╡ a49017b7-64b9-4519-9a91-457f2a303e1c
using Yao, YaoPlots

# ╔═╡ 9b314a8a-d521-4b97-a11b-c5226196c6d4
using YaoExtensions

# ╔═╡ 8cfe674b-a624-465f-9caa-f839ee83380e
md"""
# Yao Basics

First let's import the module by "using" it:
"""

# ╔═╡ 62c75c87-eaa5-4320-9b49-3c2f8a33397b
md"""
Yao represents quantum circuits/gates using **Yao Blocks**, they are a collection of Julia objects.

For example, you can implement a quantum Fourier transformation circuit as following:
"""

# ╔═╡ 763c233b-95cf-4a58-b3fb-ae3485a09f52
A(i, j) = control(i, j=>shift(2π/(1<<(i-j+1))))

# ╔═╡ 22075649-4798-4a89-8a37-02d6be489f3f
B(n, k) = chain(n, j==k ? put(k=>H) : A(j, k) for j in k:n)

# ╔═╡ 250eac58-7c75-482d-8ce7-a423e281fee5
qft(n) = chain(B(n, k) for k in 1:n)

# ╔═╡ 35f0bb9a-6676-4afa-864e-14e15f56ff0d
plot(qft(3))

# ╔═╡ 20d1d437-ceb3-4fb8-ab81-51a9aebe1bb3
md"""
here we use `plot` function to plot the generated quantum circuit, you can also use it to check what are the block `A` and block `B`.

The `chain` function is used to chain two blocks of same size together:
"""

# ╔═╡ cd837663-02bf-4540-9d75-1c0a456940c7
plot(chain(X, Y, H))

# ╔═╡ 76d40c22-3307-4ff8-b0cd-19d8cb815c72
md"the `put` function is used to put a gate on a specific location, it thus creates a larger block"

# ╔═╡ 81a4bb41-5790-42de-bf8c-ee909cc816df
plot(put(5, 2=>H))

# ╔═╡ 960629fa-1409-441d-a3c0-f5276d291c6c
md"""the control gates are defined using `control` block with another block as its input.

- the 1st argument is the number of qubits
- the 2nd argument is the controlled gate and its location
"""

# ╔═╡ a6aa360a-632a-46d7-b061-f48841042b13
plot(control(5, 3, 2=>H))

# ╔═╡ 575c593a-1957-4d97-84b6-622b9a99fa42
md"the quantum blocks defined for a quantum circuit eventually form a tree-like structure, they are also printed in this way:"

# ╔═╡ da1f61c6-f17a-4944-8c23-05ce1af83ca2
qft(3)

# ╔═╡ 3a7e4835-555d-4938-bba7-a0c0339bda5e
md"now we can define registers, Yao provides a builtin register type `ArrayReg` for full amplitude simulation:"

# ╔═╡ 3123fc39-a783-4ac9-b8ca-9c60c72a8a73
@doc ArrayReg

# ╔═╡ 01afe151-a6aa-4121-8db7-66560d3eec94
md"you can also create some common states using convenient functions, such as `rand_state`, `zero_state`:"

# ╔═╡ 43a17b19-16c3-4637-bcad-47bf0c562251
md"""
There are some implemented common components for building a quantum algorithm in the `YaoExtensions` package
"""

# ╔═╡ fc647d13-cfce-4093-a5af-6e049a8f66f5
md"""
for example you can build a variational quantum circuit quickly using `variational_circuit`
"""

# ╔═╡ de04aa00-70a7-4d12-8ff9-622c08f5eadf
variational_circuit(5, 3, [1=>2, 2=>3, 3=>4, 4=>5, 5=>1]) |> plot

# ╔═╡ 829d1f6d-7ee3-45d3-ba10-62529de07c54
md"now we can execute the circuit by on a register"

# ╔═╡ a41c3329-a644-48f4-a10d-c3c4887310ed
r = zero_state(5)

# ╔═╡ 3834a873-c821-4275-843a-11ee24b237b5
r |> variational_circuit(5, 3, [1=>2, 2=>3, 3=>4, 4=>5, 5=>1])

# ╔═╡ 260176ac-b108-4881-97c3-4c6aa6c52c6a
md"and calculate the expectation of a Hamiltonian"

# ╔═╡ 54ae27c2-8d35-4c8c-9fda-5715e1da1b1f
expect(heisenberg(5), r)

# ╔═╡ 4c9198b0-98d2-44ab-b4e6-a086342ab4c4
md"or we can also write it in a more compact form"

# ╔═╡ a3b75abb-03c0-4290-b42a-9b5b991ab85a
expect(heisenberg(5), r=>variational_circuit(5, 3, [1=>2, 2=>3, 3=>4, 4=>5, 5=>1]))

# ╔═╡ 170b8dd4-dbab-4a03-bd27-9e7fa0a9770e
md"now you can add the adjoint `'` mark to get the gradient together with the register using reverse mode differentiation"

# ╔═╡ 1d4839d5-8794-493f-a7aa-fcf20b6b3078
reg, ∇θ = expect'(heisenberg(5), zero_state(5)=>variational_circuit(5, 3, [1=>2, 2=>3, 3=>4, 4=>5, 5=>1]))

# ╔═╡ 2d02eb13-736b-49f3-988a-b37f5465d55a
md"""
GPU is supported by the `CuYao` package, you can install it using Pkg.add("CuYao"),
and switch your CPU register to a CUDA register using `cu` function
"""

# ╔═╡ Cell order:
# ╠═7697c424-7277-11eb-105d-150285a535c4
# ╟─8cfe674b-a624-465f-9caa-f839ee83380e
# ╠═a49017b7-64b9-4519-9a91-457f2a303e1c
# ╟─62c75c87-eaa5-4320-9b49-3c2f8a33397b
# ╠═763c233b-95cf-4a58-b3fb-ae3485a09f52
# ╠═22075649-4798-4a89-8a37-02d6be489f3f
# ╠═250eac58-7c75-482d-8ce7-a423e281fee5
# ╠═35f0bb9a-6676-4afa-864e-14e15f56ff0d
# ╟─20d1d437-ceb3-4fb8-ab81-51a9aebe1bb3
# ╠═cd837663-02bf-4540-9d75-1c0a456940c7
# ╟─76d40c22-3307-4ff8-b0cd-19d8cb815c72
# ╠═81a4bb41-5790-42de-bf8c-ee909cc816df
# ╟─960629fa-1409-441d-a3c0-f5276d291c6c
# ╠═a6aa360a-632a-46d7-b061-f48841042b13
# ╟─575c593a-1957-4d97-84b6-622b9a99fa42
# ╠═da1f61c6-f17a-4944-8c23-05ce1af83ca2
# ╟─3a7e4835-555d-4938-bba7-a0c0339bda5e
# ╠═3123fc39-a783-4ac9-b8ca-9c60c72a8a73
# ╟─01afe151-a6aa-4121-8db7-66560d3eec94
# ╟─43a17b19-16c3-4637-bcad-47bf0c562251
# ╠═9b314a8a-d521-4b97-a11b-c5226196c6d4
# ╟─fc647d13-cfce-4093-a5af-6e049a8f66f5
# ╠═de04aa00-70a7-4d12-8ff9-622c08f5eadf
# ╟─829d1f6d-7ee3-45d3-ba10-62529de07c54
# ╠═a41c3329-a644-48f4-a10d-c3c4887310ed
# ╠═3834a873-c821-4275-843a-11ee24b237b5
# ╟─260176ac-b108-4881-97c3-4c6aa6c52c6a
# ╠═54ae27c2-8d35-4c8c-9fda-5715e1da1b1f
# ╟─4c9198b0-98d2-44ab-b4e6-a086342ab4c4
# ╠═a3b75abb-03c0-4290-b42a-9b5b991ab85a
# ╟─170b8dd4-dbab-4a03-bd27-9e7fa0a9770e
# ╠═1d4839d5-8794-493f-a7aa-fcf20b6b3078
# ╟─2d02eb13-736b-49f3-988a-b37f5465d55a


text:
#---

using ACE, JuLIP, Test, ForwardDiff, LinearAlgebra
using ACEexperimental
using ACEexperimental.Combinations: get_params, set_params!, set_params
using JuLIP.Potentials: evaluate, evaluate_d
using ForwardDiff
using ForwardDiff: Dual

#---

basis = ACE.Utils.ace_basis(; species = :W, N = 4, maxdeg = 6)
Ffun = ρ -> ρ[1] + exp(-ρ[2]^2)
V = FitCombiPotential(basis, Ffun, 2)
set_params!(V, rand(length(get_params(V))))

#---

at = rattle!( bulk(:W, cubic=true, pbc = false) * 2, 0.1 )
energy(V, at)
forces(V, at)

#small change
p0 = get_params(V)
# write the energy  as a vector containing a single value
Efun = p -> [ energy(set_params(V, p), at) ]
Efun(p0)
# we want to use the jacobian
ForwardDiff.jacobian(Efun, p0)

# write forces (vector of vectors) as a single long vector
Ffun = p -> mat(forces(set_params(V, p), at))[:]
Ffun(p0)
ForwardDiff.jacobian(Ffun, p0)


text:
"""
    runno(iter)

runno(1:3) runs
   julia "test_1.jl"   and waits for finish
   julia "test_2.jl"  and does not wait but starts:
   julia "test_3.jl"

NOTE: All process output goes to the same file, and can be overwritten by
    following processes. Check the output of the function for failed processes,
    and run those separately.
"""
function runno(iter)
    procs = [run(pipeline(`julia -E include\(\"test_$(first(iter)).jl\"\)`; stderr = "stderr_$(first(iter)).txt", stdout="stderr_$(first(iter)).txt"); wait = false)]
    for i in iter[2:end]
        push!(procs, run(pipeline(
                                  `julia -E include\(\"test_$i.jl\"\)`;
                                  stderr = "stderr_$i.txt", stdout="stderr_$i.txt")
               ; wait = false))
    end
    procs
end
runno


text:
import DataFrames
import StatsModels

"""
"""
function transform(
        transformer::ImmutableFeatureArrayTransposerTransformer,
        featuresarray::AbstractMatrix,
        labelsarray::AbstractArray;
        kwargs...
        )
    featuresarraytransposed = transpose(featuresarray)
    return featuresarraytransposed, labelsarray
end

"""
"""
function transform(
        transformer::ImmutableFeatureArrayTransposerTransformer,
        featuresarray::AbstractMatrix;
        kwargs...
        )
    featuresarraytransposed = transpose(featuresarray)
    return featuresarraytransposed
end

"""
"""
function fit!(
        transformer::ImmutableFeatureArrayTransposerTransformer,
        featuresarray::AbstractMatrix,
        labelsarray::AbstractArray;
        kwargs...
        )
    return transform(transformer, featuresarray, labelsarray)
end

"""
"""
function predict(
        transformer::ImmutableFeatureArrayTransposerTransformer,
        featuresarray::AbstractMatrix,
        varargs...;
        kwargs...
        )
    return (transform(transformer, featuresarray), varargs...)
end

"""
"""
function predict_proba(
        transformer::ImmutableFeatureArrayTransposerTransformer,
        featuresarray::AbstractMatrix;
        kwargs...
        )
    return transform(transformer, featuresarray)
end

"""
"""
function DataFrame2LIBSVMTransformer(
        feature_names::AbstractVector,
        single_label_name::Symbol;
        levels::AbstractVector = [],
        )
    df2decisiontreetransformer = MutableDataFrame2DecisionTreeTransformer(
        feature_names,
        single_label_name;
        levels = levels,
        )
    featuretransposetransformer =
        ImmutableFeatureArrayTransposerTransformer()
    result = df2decisiontreetransformer |> featuretransposetransformer
    return result
end


text:
#!/usr/bin/env julia

function collatz_len(n::Int64)
    count = 1
    while n > 1
        if n % 2 == 0
            n = div(n,2)
        else
            n = 3 * n + 1
        end
        count += 1
    end
    return count
end

function main()
    max_len = 0
    num = 0
    for i in 1:999999
        l = collatz_len(i)
        if l > max_len
            max_len = l
            num = i
        end
    end
    @printf("The number with the longest Collatz len is %d\n", num)
end

main()


text:
"""
    iris_data()

Returns a dictionary `D` with the Iris flower data set, discussed in
sections 14.2.1 and 14.3.2.  The items in the dictionary are:
* `D["setosa"]`:  a 50 x 4 matrix with 50 examples of Iris Setosa
* `D["versicolor"]`: a 50 x 4 matrix with 50 examples of Iris Versicolor
* `D["virginica"]`: a 50 x 4 matrix with 50 examples of Iris Virginica.

The columns give values for four features:
1. sepal length in cm
2. sepal width in cm
3. petal length in cm
4. petal width in cm.
"""
function iris_data()

    # 50 examples of Iris Setosa
    setosa = [
    5.1 3.5 1.4 0.2;
    4.9 3.0 1.4 0.2;
    4.7 3.2 1.3 0.2;
    4.6 3.1 1.5 0.2;
    5.0 3.6 1.4 0.2;
    5.4 3.9 1.7 0.4;
    4.6 3.4 1.4 0.3;
    5.0 3.4 1.5 0.2;
    4.4 2.9 1.4 0.2;
    4.9 3.1 1.5 0.1;
    5.4 3.7 1.5 0.2;
    4.8 3.4 1.6 0.2;
    4.8 3.0 1.4 0.1;
    4.3 3.0 1.1 0.1;
    5.8 4.0 1.2 0.2;
    5.7 4.4 1.5 0.4;
    5.4 3.9 1.3 0.4;
    5.1 3.5 1.4 0.3;
    5.7 3.8 1.7 0.3;
    5.1 3.8 1.5 0.3;
    5.4 3.4 1.7 0.2;
    5.1 3.7 1.5 0.4;
    4.6 3.6 1.0 0.2;
    5.1 3.3 1.7 0.5;
    4.8 3.4 1.9 0.2;
    5.0 3.0 1.6 0.2;
    5.0 3.4 1.6 0.4;
    5.2 3.5 1.5 0.2;
    5.2 3.4 1.4 0.2;
    4.7 3.2 1.6 0.2;
    4.8 3.1 1.6 0.2;
    5.4 3.4 1.5 0.4;
    5.2 4.1 1.5 0.1;
    5.5 4.2 1.4 0.2;
    4.9 3.1 1.5 0.2;
    5.0 3.2 1.2 0.2;
    5.5 3.5 1.3 0.2;
    4.9 3.6 1.4 0.1;
    4.4 3.0 1.3 0.2;
    5.1 3.4 1.5 0.2;
    5.0 3.5 1.3 0.3;
    4.5 2.3 1.3 0.3;
    4.4 3.2 1.3 0.2;
    5.0 3.5 1.6 0.6;
    5.1 3.8 1.9 0.4;
    4.8 3.0 1.4 0.3;
    5.1 3.8 1.6 0.2;
    4.6 3.2 1.4 0.2;
    5.3 3.7 1.5 0.2;
    5.0 3.3 1.4 0.2
    ];
    # 50 examples of Iris Versicolor
    versicolor = [
    7.0 3.2 4.7 1.4;
    6.4 3.2 4.5 1.5;
    6.9 3.1 4.9 1.5;
    5.5 2.3 4.0 1.3;
    6.5 2.8 4.6 1.5;
    5.7 2.8 4.5 1.3;
    6.3 3.3 4.7 1.6;
    4.9 2.4 3.3 1.0;
    6.6 2.9 4.6 1.3;
    5.2 2.7 3.9 1.4;
    5.0 2.0 3.5 1.0;
    5.9 3.0 4.2 1.5;
    6.0 2.2 4.0 1.0;
    6.1 2.9 4.7 1.4;
    5.6 2.9 3.6 1.3;
    6.7 3.1 4.4 1.4;
    5.6 3.0 4.5 1.5;
    5.8 2.7 4.1 1.0;
    6.2 2.2 4.5 1.5;
    5.6 2.5 3.9 1.1;
    5.9 3.2 4.8 1.8;
    6.1 2.8 4.0 1.3;
    6.3 2.5 4.9 1.5;
    6.1 2.8 4.7 1.2;
    6.4 2.9 4.3 1.3;
    6.6 3.0 4.4 1.4;
    6.8 2.8 4.8 1.4;
    6.7 3.0 5.0 1.7;
    6.0 2.9 4.5 1.5;
    5.7 2.6 3.5 1.0;
    5.5 2.4 3.8 1.1;
    5.5 2.4 3.7 1.0;
    5.8 2.7 3.9 1.2;
    6.0 2.7 5.1 1.6;
    5.4 3.0 4.5 1.5;
    6.0 3.4 4.5 1.6;
    6.7 3.1 4.7 1.5;
    6.3 2.3 4.4 1.3;
    5.6 3.0 4.1 1.3;
    5.5 2.5 4.0 1.3;
    5.5 2.6 4.4 1.2;
    6.1 3.0 4.6 1.4;
    5.8 2.6 4.0 1.2;
    5.0 2.3 3.3 1.0;
    5.6 2.7 4.2 1.3;
    5.7 3.0 4.2 1.2;
    5.7 2.9 4.2 1.3;
    6.2 2.9 4.3 1.3;
    5.1 2.5 3.0 1.1;
    5.7 2.8 4.1 1.3
    ];
    # 50 examples of Iris Virginica
    virginica = [
    6.3 3.3 6.0 2.5;
    5.8 2.7 5.1 1.9;
    7.1 3.0 5.9 2.1;
    6.3 2.9 5.6 1.8;
    6.5 3.0 5.8 2.2;
    7.6 3.0 6.6 2.1;
    4.9 2.5 4.5 1.7;
    7.3 2.9 6.3 1.8;
    6.7 2.5 5.8 1.8;
    7.2 3.6 6.1 2.5;
    6.5 3.2 5.1 2.0;
    6.4 2.7 5.3 1.9;
    6.8 3.0 5.5 2.1;
    5.7 2.5 5.0 2.0;
    5.8 2.8 5.1 2.4;
    6.4 3.2 5.3 2.3;
    6.5 3.0 5.5 1.8;
    7.7 3.8 6.7 2.2;
    7.7 2.6 6.9 2.3;
    6.0 2.2 5.0 1.5;
    6.9 3.2 5.7 2.3;
    5.6 2.8 4.9 2.0;
    7.7 2.8 6.7 2.0;
    6.3 2.7 4.9 1.8;
    6.7 3.3 5.7 2.1;
    7.2 3.2 6.0 1.8;
    6.2 2.8 4.8 1.8;
    6.1 3.0 4.9 1.8;
    6.4 2.8 5.6 2.1;
    7.2 3.0 5.8 1.6;
    7.4 2.8 6.1 1.9;
    7.9 3.8 6.4 2.0;
    6.4 2.8 5.6 2.2;
    6.3 2.8 5.1 1.5;
    6.1 2.6 5.6 1.4;
    7.7 3.0 6.1 2.3;
    6.3 3.4 5.6 2.4;
    6.4 3.1 5.5 1.8;
    6.0 3.0 4.8 1.8;
    6.9 3.1 5.4 2.1;
    6.7 3.1 5.6 2.4;
    6.9 3.1 5.1 2.3;
    5.8 2.7 5.1 1.9;
    6.8 3.2 5.9 2.3;
    6.7 3.3 5.7 2.5;
    6.7 3.0 5.2 2.3;
    6.3 2.5 5.0 1.9;
    6.5 3.0 5.2 2.0;
    6.2 3.4 5.4 2.3;
    5.9 3.0 5.1 1.8
    ]
    return Dict([ ("setosa", setosa), ("versicolor", versicolor),
	("virginica", virginica) ])
end;

max_issues_repo_name: lynnlangit/julia-linear-algebra


text:
module GaPLACTests

using ReTest
using GaPLAC

@testset "CLI formula parseing" begin
    @test true
end

end # module

max_stars_repo_name: biobakery/GaPLAC

max_issues_repo_name: biobakery/GaPLAC

max_forks_repo_name: biobakery/GaPLAC


text:
using YYZRecipes
using Test

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


text:
import RowEchelon
import LinearAlgebra

#=
Homework 7
Section 3.2
Problems 26
=#

# Create a matrix with the vectors as columns
A = [3 2 -2 0 ;
     5 -6 -1 0 ;
     -6 0 3 0 ;
     4 7 0 -3]

# Find the detA
detA = LinearAlgebra.det(A)

println("Original Matrix:")
println(A)
println("Determinate:")
println(detA)

# See if the Determine was 0
# Dont forget floating point error
if (detA > 0.000001 || detA < -0.000001)
    println("Hence the original matrix is invertable, and its columns form a linearly independent set.")
else
    println("Hence the original matrix is NOT invertable, and its columns form a linearly dependent set.")
end


text:
function t1t2model(lx::Int,
                   t1::T, t2::T, mu::Float64;
                   boundary::Symbol=:OBC) where {T<:Number}

    hop1 = FermionQModelInteraction{1, 2, T}(
        t1,
        (1, 1),
        ((0,), (1,)), [FermionHop])
    hop2 = FermionQModelInteraction{1, 2, T}(
        t2,
        (1, 1),
        ((0,), (2,)), [FermionHop])

    chemical = FermionQModelInteraction{1, 1, Float64}(
        mu, (1,), (0,), [FermionNum])

    terms = FermionQModelInteraction[]
    t1 != 0 && push!(terms, hop1)
    t2 != 0 && push!(terms, hop2)

    UnitCellQModel(fermion,
                   QLattice(triangularunitcell, lx, boundary),
                   Z2Charge,
                   terms)
end

function triangularhopping(ls::Tuple{Int, Int},
                           t1::T,
                           t2::T,
                           t3::T,
                           mu::Float64;
                           boundary::Tuple{Symbol, Symbol}=(:PBC, :OBC)) where{T}

    ly, lx = ls

    fermion = Fermion{:f}()
    ft, f = fermionoperators(fermion)

    fhop = [(ft, f), (f, ft)]
    hop1 = FermionQModelInteraction{2, 2, T}(
        t1,
        (1, 1),
        ((0, 0), (0, 1)), fhop)

    hop2 = FermionQModelInteraction{2, 2, T}(
        t2,
        (1, 1),
        ((0, 0), (1, 0)), fhop)

    hop3 = FermionQModelInteraction{2, 2, T}(
        t3,
        (1, 1),
        ((0, 0), (1, -1)), fhop)

    terms = FermionQModelInteraction[]
    t1 != 0 && push!(terms, hop1)
    t2 != 0 && push!(terms, hop2)
    t3 != 0 && push!(terms, hop3)

    UnitCellQModel(fermion,
                   QLattice(triangularunitcell, (ly, lx), boundary),
                   Z2Charge,
                   terms)
end

function honeycombhopping(ls::Tuple{Int, Int},
                          t1::T,
                          t2::T,
                          t3::T,
                          mu1::Float64=0.0,
                          mu2::Float64=0.0;
                          boundary::Tuple{Symbol, Symbol}=(:PBC, :OBC)) where{T}

    ly, lx = ls

    fermion = Fermion{:f}()
    ft, f = fermionoperators(fermion)

    fhop = [(ft, f), (f, ft)]
    hop1 = FermionQModelInteraction{2, 2, T}(
        t1,
        (1, 2),
        ((0, 0), (0, 0)), fhop)

    hop2 = FermionQModelInteraction{2, 2, T}(
        t2,
        (2, 1),
        ((0, 0), (1, 0)), fhop)

    hop3 = FermionQModelInteraction{2, 2, T}(
        t3,
        (2, 1),
        ((0, 0), (1, -1)), fhop)

    terms = FermionQModelInteraction[]
    t1 != 0 && push!(terms, hop1)
    t2 != 0 && push!(terms, hop2)
    t3 != 0 && push!(terms, hop3)

    UnitCellQModel(fermion,
                   QLattice(honeycombunitcell, (ly, lx), boundary),
                   Z2Charge,
                   terms)
end

function kagomestriphopping(lx::Int,
                            tl::T,
                            tc::T,
                            mu::Float64;
                            boundary::Symbol=:OBC) where{T}

    fermion = Fermion{:f}()
    ft, f = fermionoperators(fermion)

    fhop = [(ft, f), (f, ft)]
    hop1 = FermionQModelInteraction{1, 2, T}(
        tl,
        (1, 1),
        ((0,), (1,)), fhop)

    hop2 = FermionQModelInteraction{1, 2, T}(
        tl,
        (3, 3),
        ((0,), (1,)), fhop)

    hop3 = FermionQModelInteraction{1, 2, T}(
        tc,
        (1, 2),
        ((0,), (0,)), fhop)

    hop4 = FermionQModelInteraction{1, 2, T}(
        tc,
        (3, 2),
        ((0,), (0,)), fhop)

    hop5 = FermionQModelInteraction{1, 2, T}(
        tc,
        (1, 2),
        ((0,), (1,)), fhop)

    hop6 = FermionQModelInteraction{1, 2, T}(
        tc,
        (3, 2),
        ((0,), (1,)), fhop)

    terms = FermionQModelInteraction[]
    if tl != zero(T)
        push!(terms, hop1)
        push!(terms, hop2)
    end
    if tc != zero(T)
        push!(terms, hop3)
        push!(terms, hop4)
        push!(terms, hop5)
        push!(terms, hop6)
    end
    UnitCellQModel(fermion,
                   QLattice(kagomestripunitcell, lx, boundary),
                   Z2Charge,
                   terms)
end

max_stars_repo_name: amiragha/mps

max_issues_repo_name: amiragha/mps

max_forks_repo_name: amiragha/mps


text:
module RPS
    include("../../examples/rock_paper_scissors.jl")
end

@testset "rock paper scissors" begin
    g = RPS.RockPaperScissors()
    reset!(g)
    while !terminated(g)
        @test sum(act!(g, rand(actions(g)))) == 0
    end
end


text:
module MyLinInterp

  export LinearInterpolation

  immutable LinearInterpolation # Inplementation as a type class
    grid::Array
    vals::Array
  end

  function call(a::LinearInterpolation, xs)
    map(xs) do x
      #a.grid[1] <= x <= a.grid[end] || throw(DomainError()) # throw an error if input x is out of domain
        if x < a.grid[1]
            return a.vals[1]
        elseif x > a.grid[end]
            return a.vals[end]
        end

      i = searchsortedlast(a.grid, x) # find the first index i in grid where grid[i] >= x
      if i == length(a.grid)
        return a.vals[end]
      else
        return a.vals[i] + (a.vals[i+1] - a.vals[i]) / (a.grid[i+1] - a.grid[i]) * (x - a.grid[i])
      end
    end
  end

end

max_stars_repo_path: Julia/Untitled Folder 3/MyLinInterp.jl

max_issues_repo_path: Julia/Untitled Folder 3/MyLinInterp.jl

max_forks_repo_path: Julia/Untitled Folder 3/MyLinInterp.jl


text:
const registry = Dict{String, NamedTuple{(:org, :url), Tuple{String, String}}}(
    # in alphabetic order
    "BFSU" => (;
        org="北京外国语大学开源软件镜像站",
        url="https://mirrors.bfsu.edu.cn/julia"
    ),
    "NJU" => (;
        org="南京大学开源镜像站",
        url="https://mirrors.nju.edu.cn/julia")
    ,
    "OpenTUNA" => (;
        org="OpenTUNA开源镜像站 -- TUNA 协会",
        url="https://opentuna.cn/julia"
    ),
    "SJTUG" => (;
        org="上海交通大学Linux用户组 (SJTUG) 软件源镜像服务",
        url="https://mirrors.sjtug.sjtu.edu.cn/julia"
    ),
    "SUSTech" => (;
        org="南方科技大学开源镜像站",
        url="http://mirrors.sustech.edu.cn/julia"
    ),
    "TUNA" => (;
        org="清华大学开源软件镜像站 -- TUNA 协会",
        url="https://mirrors.tuna.tsinghua.edu.cn/julia"
    ),
    "USTC" => (;
        org="中国科学技术大学开源软件镜像",
        url="https://mirrors.ustc.edu.cn/julia"
    ),
)

# always put "JuliaLang" to the first
registry_str = join(["- `\"$k\"`: $(registry[k].org)\n" for k in sort(collect(keys(registry)))])
registry["JuliaLang"] = (; org="The official Julia pkg server", url="https://pkg.julialang.org")
registry_str = "- `\"JuliaLang\"`: The official Julia pkg server\n" * registry_str


text:
using GLAbstraction, GLWindow, GLVisualize
using FileIO, GeometryTypes, Reactive, Images
using Colors, Plots, StaticArrays
using Plots; glvisualize(size=(800, 300))

function benchscatter(timings, names, images; n = 30)
    yposs, xposs = Vector{Float64}[], Vector{Float64}[]
    binss, scales = Vector{Float64}[], Vector{Float64}[]
    x = 1; dist = -0.5:0.0001:0.5
    min_res = 100

    for ts in timings
        bins = zeros(n); mini, maxi = minimum(ts), maximum(ts)
        w = max(maxi - mini, 1)
        min_res = min(mini, min_res)
        for t in ts
            t == 0 && continue # ignore the zeros we added
            i = round(Int, ((t-mini)/w)*(n-1) + 1)
            bins[i] += 1
        end
        bins .= (/).(bins, maximum(bins))
        xpos = similar(ts)
        scale = similar(ts)
        for (j,t) in enumerate(ts)
            if t != 0 # ignore the zeros we added
                i = round(Int, ((t-mini)/w)*(n-1) + 1)
                xpos[j] = x + rand(dist)*bins[i]*3
                scale[j] = bins[i]
            else
                xpos[j] = x
                scale[j] = 0.0
            end
        end
        push!(scales, scale)
        push!(xposs, xpos)

        gap = w/n
        push!(yposs, [t+rand(dist)*gap for t in ts])
        push!(binss, bins)
        x += 5
    end

    mscale = map(hcat(scales...)) do s
        s==0 && return s # preserve 0.0
        clamp(s*3.0, 1.5, 3.5) # else, don't make things too small
    end
    p = scatter(
        xposs, yposs, marker_z=hcat(scales...),
        shape=:circle, leg=false, ms=mscale,
        markerstrokewidth=0f0, markeralpha=0.15f0,
        markercolor=cgrad(:inferno, scale=:exp),
        title="Frame rendering times",
        ylabel="time in ms",
    )
    x = [mean(x) for x in xposs] .- 2.5
    scatter!(
        x, fill(min_res-1, length(x)),
        shape=images, hover=names, ms=5
    )
    p
end

function summarize(window, config)
    # These are allowed to fail, since they depend on not installed packages
    allowed_to_fail = ("mesh_edit.jl", "billiard.jl")
    # open(homedir()*"/results.jls", "w") do io
    #     serialize(io, config.attributes)
    # end
    allowed_failures = filter(config.attributes) do k, dict
        !dict[:success] && (basename(k) in allowed_to_fail)
    end
    failures = filter(config.attributes) do k, dict
        !dict[:success] && !(basename(k) in allowed_to_fail)
    end
    successfull = filter(config.attributes) do k, dict
        dict[:success]
    end

    resize!(window, 800, 700)

    empty!(window)

    area = Signal(SimpleRectangle(0, 0, 800, 300))
    area2 = Signal(SimpleRectangle(0, 300, 800, 400))

    plot_screen = Screen(window, name=:plots, area=area)
    glvis_screen = Screen(window, name=:glvis, area=area2)
    GLVisualize.add_screen(plot_screen) # make main screen for Plots

    # ystat = [length(failures), length(allowed_failures), length(successfull)]
    # failur_plt = bar(
    #     ["failures", "allowed failures", "passed"],
    #     ystat,
    #     markerstrokewidth=0f0, leg=false,
    #     title="Test Statistic",
    #     color=[RGBA(0.8, 0.1, 0.2, 0.6), RGBA(0.8, 0.6, 0.1, 0.6),  RGBA(0.1, 0.5, 0.4, 0.6)],
    #     ylabel="number of tests",
    #     hover=map(string, ystat)
    # )
    benchmark_names = String[]; times = Vector{Float64}[]
    success_thumbs = Matrix{RGB{N0f8}}[]
    for (k,v) in successfull
        if haskey(v, :timings) && haskey(v, :thumbnail)
            push!(benchmark_names, basename(k))
            ts = v[:timings] .* 1000.0
            ts = if length(ts) < 300
                append!(ts, fill(0.0, 300 - length(ts)))
            elseif length(ts) > 300
                resize!(ts, 300) # resampling would be more accurate, I suppose
            end
            # take 300 samples, not including the first (JIT, yo)
            push!(times, ts) # for m
            img = rotl90(v[:thumbnail])
            for x in 1:size(img, 2)
                img[:, x] = reverse(view(img, 1:size(img, 1), x))
            end
            push!(success_thumbs, img)
        end
    end
    benchplot = benchscatter(times, benchmark_names, success_thumbs)
    gui()
    rows = 11
    len = length(success_thumbs)-1
    w = 64
    positions = Point2f0[((i%rows)*w*1.05, div(i, rows)*w*1.05) for i=0:len]
    positions = positions .+ Scalar(Point2f0(w/2, 0))

    imgs = visualize(
        (success_thumbs, positions),
        scale=Vec2f0(w), stroke_width=2f0
    )
    _view(imgs, glvis_screen)

    Plots.hover(imgs, benchmark_names, glvis_screen)
    GLWindow.waiting_renderloop(window)

    if !isempty(failures)
        error("Tests not passed with $(length(failures)) failures")
    end
end

# config gets inserted by ExampleRunner
summarize(GLWindow.rootscreen(config.window), config)


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

Base.ACTIVE_PROJECT[] = joinpath(@__DIR__, "..")
using Documenter
using Pkg

makedocs(
    modules = [Pkg],
    format = :html,
    sitename = "Pkg.jl",
    pages = Any[
        "Pkg" => "index.md",
    ]
)

deploydocs(
    repo = "github.com/JuliaLang/Pkg.jl",
    target = "build",
    julia = "nightly",
    deps = nothing,
    make = nothing,
)

max_stars_repo_name: greimel/julia

max_issues_repo_name: greimel/julia

max_forks_repo_name: greimel/julia


text:
"""
    Implements all needed parts for constructing DAGMM model
        (https://sites.cs.ucsb.edu/~bzong/doc/iclr18-dagmm.pdf) for detecting anomalies.
    Code is inspired by multiple reimplementations
    - https://github.com/danieltan07/dagmm
    - https://github.com/tnakae/DAGMM/
    - https://github.com/Newcomer520/tf-dagmm
    - https://github.com/mperezcarrasco/PyTorch-DAGMM
"""

using Flux: mse, softmax, unsqueeze, stack
using LinearAlgebra
using StatsBase
using Statistics
using Distributions
using DistributionsAD
using MLDataPattern: RandomBatches

struct DAGMM{E,D,S}
    encoder::E
    decoder::D
    estimator::S
end

Flux.@functor DAGMM

function DAGMM(;idim::Int=1, activation="tanh", hdim=60, zdim=1, ncomp=2,
        nlayers::Int=3, dropout=0.5, init_seed=nothing, kwargs...)
    # if seed is given, set it
    (init_seed !== nothing) ? Random.seed!(init_seed) : nothing

    encoder = build_mlp(idim, hdim, zdim, nlayers, activation=activation, lastlayer="linear")
    decoder = build_mlp(zdim, hdim, idim, nlayers, activation=activation, lastlayer="linear")
    estimator = build_mlp(zdim+2, hdim, ncomp, nlayers, activation=activation, lastlayer="linear")

    if dropout > 0.0
        estimator = Chain(estimator[1:end-1]..., Dropout(dropout), estimator[end])
    end

    # reset seed
    (init_seed !== nothing) ? Random.seed!() : nothing

    DAGMM(encoder, decoder, estimator)
end

function Base.show(io::IO, model::DAGMM)
    # to avoid the show explosion
    print(io, "DAGMM(...)")
end

norm2(x; dims=1) = sqrt.(sum(abs2, x; dims=dims))

function cosine_similarity(x₁, x₂; dims=1, ε=1e-8)
    T = eltype(x₁)
    nx₁ = norm2(x₁; dims=dims)
    nx₂ = norm2(x₂; dims=dims)
    sum(x₁ .* x₂; dims=dims) ./ max.(nx₁ .* nx₂, T(ε))
end

function compute_reconstruction(x, x̂)
    relrec = norm2(x .- x̂; dims=1) ./ norm2(x; dims=1)
    cosim = cosine_similarity(x, x̂; dims=1)
    relrec, cosim
end

function (model::DAGMM)(x)
    z_c = model.encoder(x)
    x̂ = model.decoder(z_c)
    rec_1, rec_2 = compute_reconstruction(x, x̂)
    z = cat([z_c, rec_1, rec_2]..., dims=1)
    gamma = softmax(model.estimator(z))
    z_c, x̂, z, gamma
end

"""
Computing the parameters phi, mu and gamma for sample energy function

K: number of Gaussian mixture components
N: Number of samples
D: Latent dimension

z = DxN
gamma = KxN

this will be a little bit different in Julia
"""
function compute_params(z, gamma)
    #phi = D
    gamma_sum = sum(gamma; dims=2)
    phi = gamma_sum ./ size(gamma, 2)

    #mu = D x K x 1
    mu = sum(unsqueeze(gamma,1) .* unsqueeze(z,2); dims=3) ./ gamma_sum'

    # z_mu = D x K x N
    z_mu = unsqueeze(z,2) .- mu

    # z_mu = D x D x K x N
    z_mu_z_mu_t = unsqueeze(z_mu, 2) .* unsqueeze(z_mu, 1)

    #cov = K x D x D
    cov = sum(unsqueeze(unsqueeze(gamma, 1), 1) .* z_mu_z_mu_t; dims=4)
    cov = dropdims(cov; dims=4) ./ unsqueeze(gamma_sum', 1)

    phi, mu, cov
end

function compute_energy(z, phi, mu, cov; ε=1e-7)
    T = eltype(z)
    z_mu = unsqueeze(z,2) .- mu
    D, K, _ = size(mu)

    cov = cov .+ Diagonal(ones(Float32, D) .* T(ε))

    cov_inverse = stack([inv(cov[:,:,k]) for k in 1:K], 3)
    det_cov = unsqueeze(unsqueeze([det(cov[:,:,k] .* 2 * T.(π)) for k in 1:K],1),1)
    cov_diag = sum([sum(1 ./ diag(cov[:,:,k])) for k in 1:K])

    E_z = -T(0.5) .* sum(sum(unsqueeze(z_mu, 1) .* cov_inverse, dims=2) .* unsqueeze(z_mu, 2), dims=1)
    E_z = exp.(E_z)
    E_z = -log.(sum(unsqueeze(phi', 1).* E_z ./ sqrt.(det_cov), dims=3) .+ T(ε))

    E_z, cov_diag
end

function loss(model::DAGMM, x, λ₁, λ₂)
    _, x̂, z, gamma = model(x)
    reconst_loss = mse(x, x̂)
    phi, mu, cov = compute_params(z, gamma)
    sample_energy, cov_diag = compute_energy(z, phi, mu, cov)
    reconst_loss + λ₁ * mean(sample_energy) + λ₂ * cov_diag
end

loss(model::DAGMM, x, λ₁, λ₂, batchsize::Int) =
    mean(map(y -> loss(model, y, λ₁, λ₂), Flux.Data.DataLoader(x, batchsize=batchsize)))

# during testing we use params fitted from training
# testmode is not really necessary because dropout is not part if the autoencoder
function StatsBase.predict(model::DAGMM, x, phi, mu, cov)
    testmode!(model, true)
    _, _, z, _ = model(x)
    testmode!(model, false)
    compute_energy(z, phi, mu, cov)[1][:]
end

function StatsBase.fit!(model::DAGMM, data::Tuple; max_train_time=82800,
                        batchsize=64, patience=200, check_interval::Int=1,
                        wreg=1f-6, lambda_e=0.1f0, lambda_d=0.005f0, lr=1f-4, kwargs...)
    # add regularization through weight decay in optimizer
    opt = (wreg > 0) ? ADAMW(lr, (0.9, 0.999), wreg) : Flux.ADAM(lr)

    trn_model = deepcopy(model)
    ps = Flux.params(trn_model);

    X = data[1][1]
    # filter only normal data from validation set
    X_val = data[2][1][:, data[2][2] .== 0.0f0]
    val_N = size(X_val, 2)

    history = MVHistory()
    _patience = patience

    best_val_loss = Inf
    i = 1
    start_time = time()
    frmt(v) = round(v, digits=4)
    for batch in RandomBatches(X, batchsize)
        batch_loss = 0f0

        grad_time = @elapsed begin
            gs = gradient(() -> begin
                batch_loss = loss(trn_model, batch, lambda_e, lambda_d)
            end, ps)
            Flux.update!(opt, ps, gs)
        end

        push!(history, :batch_loss, i, batch_loss)
        push!(history, :grad_time, i, grad_time)

        if (i%check_interval == 0)
            testmode!(trn_model, true)
            val_loss_time = @elapsed val_loss = loss(trn_model, X_val, lambda_e, lambda_d)
            # val_loss_time = @elapsed val_loss = val_N > 5000 ? loss(trn_model, X_val, lambda_e, lambda_d, 512) : loss(trn_model, X_val, lambda_e, lambda_d)
            testmode!(trn_model, false)

            @info "$i - loss: $(frmt(batch_loss)) (batch) | $(frmt(val_loss)) (validation) || $(frmt(grad_time)) (t_grad) | $(frmt(val_loss_time)) (t_val)"

            if isnan(val_loss) || isinf(val_loss) || isnan(batch_loss) || isinf(batch_loss)
                @info "Stopped training after $(i) iterations, $((time() - start_time)/3600) hours due to invalid values."
                error("Encountered invalid values in loss function.")
            end

            push!(history, :validation_loss, i, val_loss)
            push!(history, :val_loss_time, i, val_loss_time)

            if val_loss < best_val_loss
                best_val_loss = val_loss
                _patience = patience

                model = deepcopy(trn_model)
            else # else stop if the model has not improved for `patience` iterations
                _patience -= 1
                if _patience == 0
                    @info "Stopped training after $(i) iterations, $((time() - start_time)/3600) hours."
                    break
                end
            end
        end

        # time limit for training
        if time() - start_time > max_train_time
            @info "Stopped training after $(i) iterations, $((time() - start_time)/3600) hours due to time constraints."
            model = deepcopy(trn_model)
            break
        end

        i += 1
    end

    (history=history, niter=i, model=model, npars=sum(map(p->length(p), ps)))
end


text:
#
# methods.jl -
#
# Implementation of methods for the Julia interface to the FLI SDK.
#
#------------------------------------------------------------------------------
handle(obj::Camera) = getfield(obj, :handle)

isnull(obj::Camera) = isnull(handle(obj))
isnull(ptr::Ptr{T}) where {T} = (ptr == null_pointer(T))

null_pointer(::Type{T}) where {T} = Ptr{T}(0)
null_pointer(x) = null_pointer(typeof(x))

_clear_handle!(obj::Camera) =
    setfield!(obj, :handle, fieldtype(typeof(obj), :handle)(0))

_convert(p::Int64) = convert(Float64,p)

function _finalize(obj::Camera)
    ptr = handle(obj)
    if !isnull(ptr)
        _clear_handle!(obj)

    end
    return nothing
end

#------------------------------------------------------------------------------
# Specialize `getproperty` and `setproperty!` in the name of the member (for
# type-stability and faster code).

getproperty(obj::Camera, sym::Symbol) = getproperty(obj, Val(sym))
setproperty!(obj::Camera, sym::Symbol, val) =
    setproperty!(obj, Val(sym), val)

# The following methods are to deal with errors.
getproperty(obj::Camera, ::Val{M}) where M = throw_unknown_field(Camera, M)

setproperty!(obj::Camera, ::Val{M}, val) where M =
    if M in propertynames(obj)
        throw_read_only_field(Camera, M)
    else
        throw_unknown_field(Camera, M)
    end

@noinline throw_unknown_field(T::Type, sym::Union{Symbol,AbstractString}) =
    throw(ErrorException("objects of type $T have no field `$sym`"))

@noinline throw_read_only_field(T::Type, sym::Union{Symbol,AbstractString}) =
    throw(ErrorException("field `$sym` of objects of type $T is read-only"))


text:
using Requests
using DataFrames
using JLD

import Requests: get, post, put, delete, options

# File that defines `api_key` as my FantasyData api key.
# Excluded from git repo for obvious reasons.
include("api_key.jl")

key_header = Dict("Ocp-Apim-Subscription-Key" => api_key)

######
# Get all game dates
######

date_output = get("https://api.fantasydata.net/v3/nfl/stats/JSON/Schedules/2017REG";
            headers=key_header)
parsed_dates = Requests.json(date_output)
game_dates = Set()
key_format = DateFormat("YYYYMMdd")
for d in parsed_dates
   if typeof(d["Date"]) != Void
      date = Date(d["Date"], DateFormat("y-m-dTH:M:S"))
      date = uppercase(Dates.format(date, "YYYY-u-dd"))
      push!(game_dates, (d["Week"], date))
   end
end

###
# Seed entry
###
week, game_date = pop!(game_dates)

raw = get("https://api.fantasydata.net/v3/nfl/stats/JSON/DailyFantasyPoints/$game_date";
          headers=key_header)
processed = Requests.json(raw)

filtered = map(z -> filter((x,v) -> x in ["Name", "PlayerID", "Position", "Team", "FantasyPointsPPR"], z), processed)
map(z -> merge!(z, Dict("Week" => week)), filtered)
df = DataFrame(filtered[1])
for j in 2:length(filtered)
   push!(df, filtered[j])
end

###
# Full loops
###
for date_tuple in game_dates
       game_date = date_tuple[2]
       week = date_tuple[1]
       raw = get("https://api.fantasydata.net/v3/nfl/stats/JSON/DailyFantasyPoints/$game_date";
                 headers=key_header)
       processed = Requests.json(raw)
       filtered = map(z -> filter((x,v) -> x in ["Name", "PlayerID", "Position", "Team", "FantasyPointsPPR"], z), processed)
       map(z -> merge!(z, Dict("Week" => week)), filtered)

       for j in 1:length(filtered)
          push!(df, filtered[j])
       end
end
unique!(df)

JLD.save("player_stats_2017.jld", "players", df)


text:
module EliminateGraphs
import LightGraphs
include("utils.jl")
include("Core.jl")
include("Graph.jl")
include("iterset.jl")
include("generateset.jl")
include("degrees.jl")
include("eliminate.jl")
include("graphlib.jl")
include("algorithms/algorithms.jl")
end # module


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    Quantile(distribution=Normal())

The quantile transform to a given `distribution`.

# Examples

```julia
Quantile()
Quantile(Normal())
```
"""
struct Quantile{D} <: Colwise
  dist::D
end

Quantile() = Quantile(Normal())

assertions(::Type{Quantile}) = [assert_continuous]

isrevertible(::Type{<:Quantile}) = true

colcache(::Quantile, x) = EmpiricalDistribution(x)

function colapply(transform::Quantile, x, c)
  origin, target = c, transform.dist
  qqtransform(x, origin, target)
end

function colrevert(transform::Quantile, y, c)
  origin, target = transform.dist, c
  qqtransform(y, origin, target)
end

# transform samples from original to target distribution
function qqtransform(samples, origin, target)
  # avoid evaluating the quantile at 0 or 1
  T = eltype(samples)
  low = zero(T) + T(0.001)
  high = one(T) - T(0.001)
  map(samples) do sample
    prob = cdf(origin, sample)
    quantile(target, clamp(prob, low, high))
  end
end


text:
function solve()
    a, b = [parse(Int, x) for x in split(readline())]
    2a + 100 - b
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
import BFloat16s.BFloat16

"""The BFloat16 + stochastic rounding type."""
primitive type BFloat16sr <: AbstractFloat 16 end

Base.sign_mask(::Type{BFloat16sr}) = 0x8000
Base.exponent_mask(::Type{BFloat16sr}) = 0x7f80
Base.significand_mask(::Type{BFloat16sr}) = 0x007f

Base.iszero(x::BFloat16sr) = iszero(Float32(x))
Base.isfinite(x::BFloat16sr) = isfinite(Float32(x))
Base.isnan(x::BFloat16sr) = isnan(Float32(x))

Base.precision(::Type{BFloat16sr}) = 8
Base.one(::Type{BFloat16sr}) = reinterpret(BFloat16sr,0x3f80)
Base.zero(::Type{BFloat16sr}) = reinterpret(BFloat16sr,0x0000)

const InfB16sr = reinterpret(BFloat16sr, 0x7f80)
const NaNB16sr = reinterpret(BFloat16sr, 0x7fc0)

Base.typemin(::Type{BFloat16sr}) = -InfB16sr
Base.typemax(::Type{BFloat16sr}) = InfB16sr
Base.floatmin(::Type{BFloat16sr}) = reinterpret(BFloat16sr,0x0080)
Base.floatmax(::Type{BFloat16sr}) = reinterpret(BFloat16sr,0x7f7f)
minpos(::Type{BFloat16sr}) = reinterpret(BFloat16sr,0x0001)

Base.typemin(::BFloat16sr) = typemin(BFloat16sr)
Base.typemax(::BFloat16sr) = typemax(BFloat16sr)
Base.floatmin(::BFloat16sr) = floatmin(BFloat16sr)
Base.floatmax(::BFloat16sr) = floatmax(BFloat16sr)
minpos(::BFloat16sr) = minpos(BFloat16sr)

# Truncation from Float32
Base.uinttype(::Type{BFloat16sr}) = UInt16
Base.trunc(::Type{BFloat16sr}, x::Float32) = reinterpret(BFloat16sr,
        (reinterpret(UInt32, x) >> 16) % UInt16)

# same for BFloat16sr, but do not apply stochastic rounding to avoid InexactError
Base.round(x::BFloat16sr, r::RoundingMode{:Up}) = BFloat16sr(ceil(Float32(x)))
Base.round(x::BFloat16sr, r::RoundingMode{:Down}) = BFloat16sr(floor(Float32(x)))
Base.round(x::BFloat16sr, r::RoundingMode{:Nearest}) = BFloat16sr(round(Float32(x)))

# Conversions
"""Convert BFloat16sr to Float32 by padding trailing zeros."""
Base.Float32(x::BFloat16sr) = reinterpret(Float32, UInt32(reinterpret(UInt16, x)) << 16)

BFloat16sr(x::Float64) = BFloat16sr(Float32(x))
BFloat16sr(x::Float16) = BFloat16sr(Float32(x))
BFloat16sr(x::Integer) = BFloat16sr(Float32(x))

Base.Float64(x::BFloat16sr) = Float64(Float32(x))
Base.Float16(x::BFloat16sr) = Float16(Float32(x))

# conversion between BFloat16 and BFloat16sr
BFloat16(x::BFloat16sr) = reinterpret(BFloat16,x)
BFloat16sr(x::BFloat16) = reinterpret(BFloat16sr,x)

# conversion to integer
(::Type{T})(x::BFloat16sr) where {T<:Integer} = T(Float32(x))

"""Convert to BFloat16sr from Float32 via round-to-nearest
and tie to even. Identical to BFloat16(::Float32)."""
function BFloat16sr(x::Float32)
    isnan(x) && return NaNB16sr
    h = reinterpret(UInt32, x)
    h += 0x7fff + ((h >> 16) & 1)
    return reinterpret(BFloat16sr, (h >> 16) % UInt16)
end

"""Stochastically round x::Float32 to BFloat16 with distance-proportional probabilities."""
function BFloat16_stochastic_round(x::Float32)
	ui = reinterpret(UInt32,x)
	ui += rand(Xor128[],UInt16)							# add stochastic perturbation in [0,u)
	return reinterpret(BFloat16sr,(ui >> 16) % UInt16)	# round to zero and convert to BFloat16
end

"""Chance that x::Float32 is round up when converted to BFloat16sr."""
function BFloat16_chance_roundup(x::Float32)
	xround = Float32(BFloat16(x))
	xround == x && return zero(Float32)
	xround_down, xround_up = xround < x ? (xround,Float32(nextfloat(BFloat16sr(xround)))) :
		(Float32(prevfloat(BFloat16sr(xround))),xround)

	return (x-xround_down)/(xround_up-xround_down)
end

# Basic arithmetic
for f in (:+, :-, :*, :/, :^)
	@eval Base.$f(x::BFloat16sr, y::BFloat16sr) = BFloat16_stochastic_round($(f)(Float32(x), Float32(y)))
end

# negation via signbit flip
Base.:(-)(x::BFloat16sr) = reinterpret(BFloat16sr, reinterpret(UInt16, x) ⊻ 0x8000)

# absolute value by setting the signbit to zero
Base.abs(x::BFloat16sr) = reinterpret(BFloat16sr, reinterpret(UInt16, x) & 0x7fff)

for func in (:sin,:cos,:tan,:asin,:acos,:atan,:sinh,:cosh,:tanh,:asinh,:acosh,
             :atanh,:exp,:exp2,:exp10,:expm1,:log,:log2,:log10,:sqrt,:cbrt,:log1p)
    @eval begin
        Base.$func(a::BFloat16sr) = BFloat16_stochastic_round($func(Float32(a)))
    end
end

for func in (:atan,:hypot)
    @eval begin
        Base.$func(a::BFloat16sr,b::BFloat16sr) = BFloat16_stochastic_round($func(Float32(a),Float32(b)))
    end
end

# Floating point comparison
function Base.:(==)(x::BFloat16sr, y::BFloat16sr)
    return BFloat16(x) == BFloat16(y)
end

for op in (:<, :<=, :isless)
    @eval Base.$op(a::BFloat16sr, b::BFloat16sr) = ($op)(Float32(a), Float32(b))
end

Base.promote_rule(::Type{Float32}, ::Type{BFloat16sr}) = Float32
Base.promote_rule(::Type{Float64}, ::Type{BFloat16sr}) = Float64

for t in (Int8, Int16, Int32, Int64, Int128, UInt8, UInt16, UInt32, UInt64, UInt128)
	@eval Base.promote_rule(::Type{BFloat16sr}, ::Type{$t}) = BFloat16sr
end

# Showing
function Base.show(io::IO, x::BFloat16sr)
    if isinf(x)
        print(io, x < 0 ? "-InfB16" : "InfB16")
    elseif isnan(x)
        print(io, "NaNB16")
    else
		io2 = IOBuffer()
        print(io2,Float32(x))
        f = String(take!(io2))
        print(io,"BFloat16sr("*f*")")
    end
end

Base.bitstring(x::BFloat16sr) = bitstring(reinterpret(UInt16,x))

function Base.bitstring(x::BFloat16sr,mode::Symbol)
    if mode == :split	# split into sign, exponent, signficand
        s = bitstring(x)
		return "$(s[1]) $(s[2:9]) $(s[10:end])"
    else
        return bitstring(x)
    end
end

function Base.nextfloat(x::BFloat16sr)
    if isfinite(x)
		ui = reinterpret(UInt16,x)
		if ui < 0x8000	# positive numbers
			return reinterpret(BFloat16sr,ui+0x0001)
		elseif ui == 0x8000		# =-zero(T)
			return reinterpret(BFloat16sr,0x0001)
		else				# negative numbers
			return reinterpret(BFloat16sr,ui-0x0001)
		end
	else	# NaN / Inf case
		return x
	end
end

function Base.prevfloat(x::BFloat16sr)
    if isfinite(x)
		ui = reinterpret(UInt16,x)
		if ui == 0x0000		# =zero(T)
			return reinterpret(BFloat16sr,0x8001)
		elseif ui < 0x8000	# positive numbers
			return reinterpret(BFloat16sr,ui-0x0001)
		else				# negative numbers
			return reinterpret(BFloat16sr,ui+0x0001)
		end
	else	# NaN / Inf case
		return x
	end
end


text:
module TestAutoscale
using AWS
using AWS.Autoscale
using Base.Test

function check_resp(resp, resp_obj_type=Any)
    info("response: ", resp.obj)
    @test 200 <= resp.http_code <= 206
    @test isa(resp, AWS.Autoscale.AutoScalingResponse)
    @test isa(resp.obj, resp_obj_type)
    @test !isa(resp.obj, AutoScalingError)
end

function runtests(env, config)
    println("getting autoscale limits")
    resp = DescribeAccountLimits(env)
    check_resp(resp, DescribeAccountLimitsResult)
    @test resp.obj.numberOfLaunchConfigurations >= 0
    @test resp.obj.maxNumberOfLaunchConfigurations >= 0
    @test resp.obj.numberOfAutoScalingGroups >= 0
    @test resp.obj.maxNumberOfAutoScalingGroups >= 0

    println("listing autoscale groups")
    resp = DescribeAutoScalingGroups(env)
    check_resp(resp, DescribeAutoScalingGroupsResult)
    valid_autoscale_groups = Vector{String}()
    for grp in resp.obj.autoScalingGroups
        println(grp.autoScalingGroupName)
        @test isa(grp.autoScalingGroupName, String)
        push!(valid_autoscale_groups, grp.autoScalingGroupName)
    end

    println("listing launch configurations")
    resp = DescribeLaunchConfigurations(env)
    check_resp(resp, DescribeLaunchConfigurationsResult)
    for lc in resp.obj.launchConfigurations
        println(lc.launchConfigurationName)
        @test isa(lc.launchConfigurationName, String)
    end

    println("listing load balancers")
    for grp in valid_autoscale_groups
        req = DescribeLoadBalancersRequest(; autoScalingGroupName=grp)
        resp = DescribeLoadBalancers(env, req)
        check_resp(resp, DescribeLoadBalancersResult)
        for lb in resp.obj.loadBalancers
            println(lc.loadBalancerName)
            @test isa(lc.loadBalancerName, String)
        end
    end

    ami = get(config, "ami-generic", "")
    insttype = get(config, "insttype", "")
    keyname = get(config, "keyname", "")
    ownertag = get(config, "ownertag", "")
    nametag = get(config, "nametag", "")
    az = get(config, "availability-zones", [])

    for v in (ami, insttype, keyname, ownertag, nametag, az)
        if isempty(v)
            warn("skipping further autoscaling tests, not all parameters configured")
            return
        end
    end

    println("create a launch configuration")
    req = CreateLaunchConfigurationRequest(; launchConfigurationName="awstest", imageId=ami, instanceType=insttype, instanceMonitoring=InstanceMonitoring(;enabled=false), keyName=keyname)
    resp = CreateLaunchConfiguration(env, req)
    check_resp(resp)

    println("create a autoscaling group")
    tags = [Tag(;key="Name", value=nametag), Tag(;key="Owner", value=ownertag)]
    req = CreateAutoScalingGroupRequest(; autoScalingGroupName="awstest", availabilityZones=az, desiredCapacity=0, launchConfigurationName="awstest", maxSize=2, minSize=0, tags=tags)
    resp = CreateAutoScalingGroup(env, req)
    check_resp(resp)

    println("scale the group up")
    req = SetDesiredCapacityRequest(; autoScalingGroupName="awstest", desiredCapacity=1)
    resp = SetDesiredCapacity(env, req)
    check_resp(resp)

    println("wait for group to scale up")
    while true
        req = DescribeAutoScalingGroupsRequest(; autoScalingGroupNames=["awstest"])
        resp = DescribeAutoScalingGroups(env, req)
        check_resp(resp, DescribeAutoScalingGroupsResult)
        grps = resp.obj.autoScalingGroups
        @test isa(grps, Vector{AutoScalingGroup})
        grp = grps[1]
        @test grp.autoScalingGroupName == "awstest"
        @test grp.desiredCapacity == 1
        @test grp.maxSize == 2
        @test grp.minSize == 0
        instances = grp.instances
        @test isa(instances, Vector{Instance})
        healthy = 0
        for instance in instances
            println("instance ", instance.instanceId, " - ", instance.healthStatus)
            (instance.healthStatus == "Healthy") && (healthy += 1)
        end
        (healthy >= 1) && break
        println("only ", healthy, "/1 instances up")
    end

    println("scale the group down")
    req = SetDesiredCapacityRequest(; autoScalingGroupName="awstest", desiredCapacity=0)
    resp = SetDesiredCapacity(env, req)
    check_resp(resp)

    println("wait for group to scale down")
    while true
        req = DescribeAutoScalingGroupsRequest(; autoScalingGroupNames=["awstest"])
        resp = DescribeAutoScalingGroups(env, req)
        check_resp(resp, DescribeAutoScalingGroupsResult)
        grps = resp.obj.autoScalingGroups
        @test isa(grps, Vector{AutoScalingGroup})
        grp = grps[1]
        @test grp.autoScalingGroupName == "awstest"
        @test grp.desiredCapacity == 0
        @test grp.maxSize == 2
        @test grp.minSize == 0
        instances = grp.instances
        @test isa(instances, Vector{Instance})
        healthy = 0
        for instance in instances
            println("instance ", instance.instanceId, " - ", instance.healthStatus)
            (instance.healthStatus == "Healthy") && (healthy += 1)
        end
        (healthy == 0) && break
        println("still ", healthy, "/1 instances up")
    end

    println("delete autoscaling group")
    req = DeleteAutoScalingGroupRequest(; autoScalingGroupName="awstest", forceDelete=true)
    resp = DeleteAutoScalingGroup(env, req)
    check_resp(resp)

    println("delete launch configuration")
    req = DeleteLaunchConfigurationRequest(; launchConfigurationName="awstest")
    resp = DeleteLaunchConfiguration(env, req)
    check_resp(resp)
end

end # module TestAutoscale


text:
# Doing epsilon greedy
# https://youtu.be/PnHCvfgC_ZA?t=3013

using Game2048
using Flux, Statistics, StatsBase
using Flux: crossentropy
using BSON: @save, @load
using CUDA:cu
using CUDA
CUDA.allowscalar(false)

const EPSILON = 0.1
const LAMBDA = 0.95 # the discount fcator

function mirrors_rotations(board)
    tboard = transpose(board) |> collect

    return board,  rotl90(board),  rot180(board),  rotr90(board),
           tboard, rotl90(tboard), rot180(tboard), rotr90(tboard)
end

actor_model = Chain(
    normalised_board -> reshape(normalised_board, 4, 4, 1, 1),
    Conv((2,2), 1=>256, relu),
    Conv((2,2), 256=>512, relu),
    Conv((2,2), 512=>1024, relu),
    x->reshape(x, 1024),
    Dense(1024, 4),
    softmax
) |> gpu

function loss_actor(board_float32, prob)
    crossentropy(actor_model(board_float32), prob)
end

critic_model1 = Chain(
    board -> reshape(board, 4, 4, 1, 1),
    Conv((2,2), 1=>256, relu),
    Conv((2,2), 256=>512, relu),
    Conv((2,2), 512=>1024, relu),
    x->reshape(x, 1024)
) |> gpu

critic_model2 = Chain(
    Dense(1024+4, 512, relu),
    Dense(512, 256, relu),
    Dense(256, 1, exp)
) |> gpu

function critic_score(board, dir)
    v1 = critic_model1(cu(Float32.(board)))
    v2 = cu(action2vec(dir))
    res = critic_model2(vcat(v1, v2)) |> collect
    res[1]
end

function action2vec(dir)
    a = zeros(Float32, 4)
    a[Int(dir)+1] = 1
    a
end

function loss_qa((board, v2), expected_value)
    v1 = critic_model1(board)
    v3 = critic_model2(vcat(v1, v2))
    sum((v3 .- expected_value).^2)
end

# figure out the possible actions
# and sample a move from the actor
function choose_action(board)
    res = move.(Ref(board), DIRS)
    possible = map(r->r[2][2], res)

    if !any(possible)
        return false, left, 0, board
    end

    if rand() <= EPSILON
        i = wsample(1:4, possible)
    else
        prob_dir = collect(actor_model(cu(Float32.(board))))
        if all(isnan, prob_dir)
            error("wtf")
        end
        prob_dir .= prob_dir .* possible
        i = wsample(1:4, prob_dir)
    end

    reward = res[i][2][1]
    dir = DIRS[i]
    new_board = res[i][1]

    true, dir, reward, new_board
end

function train_one_episode() #p, opt)
    # keep track of actions taken
    historiesX = []
    board = initboard()

    cstmp = critic_score.(Ref(board), DIRS)
    println("*****************my abilities now***********************")
    println(cstmp |> mean)
    display(cstmp)
    display(softmax(cstmp./100))
    println("****************************************")

    # choose an action
    possible, dir, reward, new_board = choose_action(board)

    # save the history
    push!(historiesX, (cu(Float32.(board)), cu(action2vec(dir))))

    max_reward = 2^8
    tot_reward = 0
    rewards = Int[]
    max_critic_scores = 0
    while possible
        # move according to choosen action
        board = new_board
        gen_new_tile!(board)

        # update the critic network with latest info
        possible, dir, reward, new_board = choose_action(board)

        tot_reward += reward
        push!(rewards, reward)

        if possible
            cs = critic_score(board, dir)
            y = zeros(Float32, length(historiesX))
            yscore = cs*LAMBDA
            for i in length(y):-1:1
                yscore += rewards[i]
                y[i] = yscore
                yscore *= LAMBDA
            end

            if reward >= max_reward
                max_reward = reward
                display(board)
                println(dir)
                println(y)
            end

            training_data = [(X, cu([y1])) for (X, y1) in zip(historiesX, y)]

            # update the critic network
            Flux.train!(loss_qa, p, training_data, opt)

            # use the updated critic network to update actor network
            critic_scores = critic_score.(Ref(board), DIRS)
            critic_scores .= softmax(critic_scores ./ 100)
            if maximum(critic_scores) > max_critic_scores
                _, pos = findmax(critic_scores)
                max_critic_scores = maximum(critic_scores)
                println("*************best score******************")
                println("$(DIRS[pos]) $max_critic_scores")
                display(board)
            end
            nboard = Float32.(board)
            nboard = cu(nboard ./ maximum(nboard))

            Flux.train!(loss_actor, pam, [(nboard, cu(critic_scores))], opt_pam)

            # display(board)
            # println(dir)
            # println(cs)
            # to save on computional power just reuse the previous chosen action
            push!(historiesX, (cu(Float32.(board)), action2vec(dir) |> cu))
        else
            return board, tot_reward, rewards
        end

    end
end

if 1 == 2
    p = Flux.params(critic_model1, critic_model2)
    opt = ADAM()

    pam = Flux.params(actor_model)
    opt_pam = ADAM()
elseif 2 == 3
    using Zygote
    @load "actor"
    @load "critic1"
    @load "critic2"
    @load "actor_p"
    @load "critic1_p"
    @load "critic2_p"
    @load "opt"
    @load "opt_pam"
    Flux.loadparams!(actor_model, pam)
    Flux.loadparams!(critic_model1, c1)
    Flux.loadparams!(critic_model2, c2)
end

@time res = train_one_episode()

for i in 1:600
    @time println(train_one_episode())
    if i % 10 == 0
        @save "actor" actor_model
        @save "critic1" critic_model1
        @save "critic2" critic_model2

        @save "opt" opt
        @save "opt_pam" opt_pam

        @save "actor_p" pam
        c1 = Flux.params(critic_model1)
        @save "critic1_p" c1

        c2 = Flux.params(critic_model2)
        @save "critic2_p" c2
    end
end


text:
function createmediatable!(data::Dict{Symbol, Any})
  push!(data, :media=>Dict{Symbol, Float64}())
  push!(data[:media], :mud=>11.35)
  push!(data[:media], :water=>8.339)
  push!(data[:media], :air=>0.0)
  push!(data[:media], :heavymud=>12.0)
  push!(data[:media], :lightmud=>10.0)
end


text:
#
#--------------------------------------#
export semmesh
#--------------------------------------#
#
# creates 1D sem mesh of poly. order n-1
# and e elements [-1,1].
#
function semmesh(E,n)

z0,w0 = FastGaussQuadrature.gausslobatto(n); # [-1,1]

z0 = 0.5 .* (z0.+1); # [0,1]
w0 = 0.5 .*  w0;

ze = range(start=-1,stop=1,length=E+1) |> Array # element mesh
dz = diff(ze);

z = kron(dz,z0) + kron(ze[1:end-1],ones(n,1));

w = kron(dz,w0);

z = reshape(z,E*n);
w = reshape(w,E*n);

return z,w
end
#--------------------------------------#
export semreshape
#--------------------------------------#
function semreshape(u,nr,ns,Ex,Ey)

    v = zeros(nr,ns,Ex*Ey)

    Ix = 1:nr
    Iy = 1:ns

    for i=1:Ex
        for j=1:Ey

            ie = j + (i-1)*Ey

            ix = Ix .+ (i-1)*nr
            iy = Iy .+ (j-1)*ns

            v[:,:,ie] = @view u[ix,iy]
        end
    end

    return v
end
#--------------------------------------#