text:
# included from Potentials.jl
# part of the module JuLIP.Potentials

using JuLIP: zerovecs, JVecsF, JVecF, JMatF, neighbourlist

using NeighbourLists

export ZeroPairPotential, PairSitePotential,
         LennardJones, lennardjones,
         Morse, morse, grad, hess

function site_energies(pp::PairPotential, at::AbstractAtoms)
   Es = zeros(length(at))
   for (i,_2,r,_3) in pairs(at, cutoff(pp))
      Es[i] += 0.5 * pp(r)
   end
   return Es
end

import Base.sum

sum(V::PairPotential, r) = sum( V(s) for s in r )

function sum(V::PairPotential, r::Vector{T}) where T <: Real
   E = 0.0
   @simd for n = 1:length(r)
      @inbounds E += V(r[n])
   end
   return E
end

# a simplified way to calculate gradients of pair potentials
@inline grad(V::PairPotential, r::Real, R::JVec) = (evaluate_d(V, r) / r) * R

# function energy(V::PairPotential, at::AbstractAtoms)
#    E = 0.0
#    for (_₁, _₂, r, _₃) in pairs(at, cutoff(V))
#       E += V(r)
#    end
#    return 0.5 * E
# end

function energy(V::PairPotential, at::AbstractAtoms)
   nlist = neighbourlist(at, cutoff(V))::PairList
   return 0.5 * sum(V, nlist.r)
end

# function forces(V::PairPotential, at::AbstractAtoms)
#    dE = zerovecs(length(at))
#    for (i,j,r,R) in pairs(at, cutoff(V))
#       dE[i] += @GRAD V(r, R)
#    end
#    return dE
# end

function forces(V::PairPotential, at::AbstractAtoms)
   nlist = neighbourlist(at, cutoff(V))::PairList
   dE = zerovecs(length(at))
   @simd for n = 1:npairs(nlist)
      @inbounds dE[nlist.i[n]] += grad(V, nlist.r[n], nlist.R[n])
   end
   return dE
end

# TODO: rewrite using generator once bug is fixed (???or maybe decide not to bother???)
function virial(pp::PairPotential, at::AbstractAtoms)
   S = zero(JMatF)
   for (_₁, _₂, r, R) in pairs(at, cutoff(pp))
      S -= 0.5 * grad(pp, r, R) * R'  # (((@D pp(r)) / r) * R) * R'
   end
   return S
end

function hess(V::PairPotential, r, R)
   R̂ = R/r
   P = R̂ * R̂'
   dV = (@D V(r))/r
   return ((@DD V(r)) - dV) * P + dV * eye(JMatF)
end

# an FF preconditioner for pair potentials
function precon(V::PairPotential, r, R)
   dV = @D V(r)
   hV = @DD V(r)
   Id = eye(JMatF)
   S = R/r
   return 0.9 * (abs(hV) * S * S' + abs(dV / r) * (Id - S * S')) +
          0.1 * (abs(hV) + abs(dV / r)) * Id
end

hessian_pos(V::PairPotential, at::AbstractAtoms) =
      _precon_or_hessian_pos(V, at, hess)

#
# this assembles a hessian or preconditioner as a block-matrix
#
function _precon_or_hessian_pos(V::PairPotential, at::AbstractAtoms, hfun)
   nlist = neighbourlist(at, cutoff(V))
   I, J, Z = Int[], Int[], JMatF[]
   for C in (I, J, Z); sizehint!(C, 2*npairs(nlist)); end
   for (i, j, r, R) in pairs(nlist)
      h = 0.5 * hfun(V, r, R)
      append!(I, (i,  i,  j, j))
      append!(J, (i,  j,  i, j))
      append!(Z, (h, -h, -h, h))
   end
   hE = sparse(I, J, Z, length(at), length(at))
   return hE
end

"""
`LennardJones(σ, e0):` constructs the 6-12 Lennard-Jones potential [wiki](https://en.wikipedia.org/wiki/Lennard-Jones_potential)

   e0 * 4 * ( (σ/r)¹² - (σ/r)⁶ )

Constructor with kw arguments: `LennardJones(; kwargs...)`

* `e0, σ` : standard LJ parameters, default `e0 = 1.0, σ = 1.0`
* `r0` : equilibrium distance - if `r0` is specified then `σ` is ignored
* `a0` : FCC lattice parameter, if `a0` is specified then `σ` is ignored

(`r0, a0` cannot both be specified at the same time)
"""
LennardJones(σ, e0) = @analytic r -> e0 * 4.0 * ((σ/r)^(12) - (σ/r)^(6))

function ljparams(; σ=1.0, e0=1.0, r0 = nothing, a0 = nothing)
   if r0 != nothing && a0 != nothing
      error("`LenndardJones`: cannot specify both `r0` and `a0`")
   end
   if a0 != nothing    # r0 = nn-dist = a0/sqrt(2) in FCC
      r0 = a0 / sqrt(2)
   end
   if r0 != nothing    # standard LJ is minimised at r = 2^(1/6)
      σ =  r0 / 2^(1/6)
   end
   return σ, e0
end

LennardJones(; kwargs...) = LennardJones(ljparams(;kwargs...)...)

"""
`lennardjones(; kwargs...)`

simplified constructor for `LennardJones` (note this is type unstable!)

In addition to the `kwargs` of `LennardJones`, this accepts also

* `rcut` : default `:auto` which gives `rcut = (1.9*σ, 2.7*σ)`. Use
`nothing` or `Inf` to specify no cutoff, or specify a tuple or
array with two elements specifying the lower and upper cut-off radii to be
used with `SplineCutoff`.
"""
function lennardjones(; rcut = :auto, kwargs...)
   σ, e0 = ljparams(; kwargs...)
   if (rcut == nothing || rcut == Inf)
      return LennardJones(σ, e0)
   elseif rcut == :auto
      rcut = (1.9*σ, 2.7*σ)
   end
   return SplineCutoff(rcut[1], rcut[2]) * LennardJones(σ, e0)
end

"""
`Morse(A, e0, r0)` or `Morse(;A=4.0, e0=1.0, r0=1.0)`: constructs a
`PairPotential` for
   e0 ( exp( -2 A (r/r0 - 1) ) - 2 exp( - A (r/r0 - 1) ) )
"""
Morse(A, e0, r0) = @analytic(
   r -> e0 * ( exp(-(2.0*A) * (r/r0 - 1.0)) - 2.0 * exp(-A * (r/r0 - 1.0)) ) )
Morse(;A=4.0, e0=1.0, r0=1.0) = Morse(A, e0, r0)

"""
`morse(A=4.0, e0=1.0, r0=1.0, rcut=(1.9*r0, 2.7*r0))`

simplified constructor for `Morse` (type unstable)
"""
morse(;A=4.0, e0=1.0, r0=1.0, rcut=(1.9*r0, 2.7*r0)) = (
   (rcut == nothing || rcut == Inf)
         ?  Morse(A, e0, r0)
         :  SplineCutoff(rcut[1], rcut[2]) * Morse(A, e0, r0) )

@pot struct ZeroPairPotential <: PairPotential end

"""
`ZeroPairPotential()`: creates a potential that just returns zero
""" ZeroPairPotential
evaluate(p::ZeroPairPotential, r::Float64) = 0.0
evaluate_d(p::ZeroPairPotential, r::Float64) = 0.0
evaluate_dd(p::ZeroPairPotential, r::Float64) = 0.0
cutoff(p::ZeroPairPotential) = 0.0

# ========================================================
# wrapping a pair potential in a site potential

SitePotential(pp::PairPotential) = PairSitePotential(pp)

@pot struct PairSitePotential{P} <: SitePotential
   pp::P
end

cutoff(psp::PairSitePotential) = cutoff(psp.pp)

# TODO: get rid of this!
function _sumpair_(pp, r)
   # cant use a generator here since type is not inferred!
   # Watch out for a bugfix
   s = 0.0
   for t in r
      s += pp(t)
   end
   return s
end

# special implementation of site energy and forces for a plain pair potential
evaluate(psp::PairSitePotential, r, R) = 0.5 * _sumpair_(psp.pp, r)

evaluate_d(psp::PairSitePotential, r, R) =
            [ 0.5 * grad(psp.pp, s, S) for (s, S) in zip(r, R) ]

struct MultiPairPotential{TV <: PairPotential}
   V::Matrix{TV}
   Z2V::Dict{Int, Int}
end


text:
using ReinforcementLearningZoo
using Test
using ReinforcementLearningBase
using ReinforcementLearningCore
using ReinforcementLearningEnvironments
using Flux
using Statistics
using Random

@testset "ReinforcementLearningZoo.jl" begin

    @testset "training" begin
        mktempdir() do dir
            for method in (:BasicDQN, :DQN, :PrioritizedDQN, :Rainbow, :IQN)
                res = run(Experiment(
                    Val(:JuliaRL),
                    Val(method),
                    Val(:CartPole),
                    nothing;
                    save_dir = joinpath(dir, string(method)),
                ))
                @info "stats for $method" avg_reward = mean(res.hook[1].rewards) avg_fps =
                    1 / mean(res.hook[2].times)
            end

            for method in (:BasicDQN, :DQN)
                res = run(Experiment(
                    Val(:JuliaRL),
                    Val(method),
                    Val(:MountainCar),
                    nothing;
                    save_dir = joinpath(dir, string(method)),
                ))
                @info "stats for $method" avg_reward = mean(res.hook[1].rewards) avg_fps =
                    1 / mean(res.hook[2].times)
            end

            for method in (:A2C, :A2CGAE, :PPO)
                res = run(Experiment(Val(:JuliaRL), Val(method), Val(:CartPole), nothing))
                @info "stats for $method" avg_reward =
                    mean(Iterators.flatten(res.hook.rewards))
            end

            for method in (:DDPG, :SAC)
                res = run(Experiment(Val(:JuliaRL), Val(method), Val(:Pendulum), nothing))
                @info "stats for $method" avg_reward =
                    mean(Iterators.flatten(res.hook.rewards))
            end
        end
    end

    @testset "run pretrained models" begin
        for x in ("JuliaRL_BasicDQN_CartPole",)
            e = Experiment(x)
            e.agent.policy = load_policy(x)
            Flux.testmode!(e.agent)
            run(e.agent, e.env, StopAfterEpisode(1), e.hook)
            @info "result of evaluating pretrained model: $x for once:" reward =
                e.hook[1].rewards[end]
        end
    end
end

ext: jl


text:
function setup_spacetime(fmax::Number, multiplier::Number, npadding::Int = 0; xlims::NTuple{2,T}, ylims::NTuple{2,T}) where T
    pts_per_wavelength = 16
    factor = c₀*1.2/Sc
    Δ = c₀ /(pts_per_wavelength*fmax)
    Δt = Δ/factor
    x = LinearAxis(xlims[1]-Δ*(npadding+1),xlims[2]+Δ*(npadding+1),Δ) # compensating axis to the new units
    y = LinearAxis(ylims[1]-Δ*(npadding+1),ylims[2]+Δ*(npadding+1),Δ)

    space = Space2D(x,y)

    Nt = round(Int, multiplier*(maximum(domain.(space)))/(c₀*Δt))

    tmin = 0.0 ; tmax = tmin + Δt*Nt
    time = LinearAxis(tmin, tmax, Δt)
    for axis in space
        @assert axis.N > 2*npadding
    end

    return (space, time)
end

function setup_spacetime(fmax::Number, multiplier::Number, npadding::Int = 0; xlims::NTuple{2,T}, ylims::NTuple{2,T}, zlims::NTuple{2,T}) where T
    pts_per_wavelength = 16
    factor = c₀*1.2/Sc
    Δ = c₀ /(pts_per_wavelength*fmax)
    Δt = Δ/factor
    x = LinearAxis(xlims[1]-Δ*(npadding+1),xlims[2]+Δ*(npadding+1),Δ) # compensating axis to the new units
    y = LinearAxis(ylims[1]-Δ*(npadding+1),ylims[2]+Δ*(npadding+1),Δ)
    space = Space2D(x,y,z)

    Nt = round(Int, multiplier*(maximum(domain.(space)))/(c₀*Δt))

    tmin = 0.0 ; tmax = tmin + Δt*Nt
    time = LinearAxis(tmin, tmax, Δt)
    for axis in space
        @assert axis.N > 2*npadding
    end

    return (space, time)
end

    function ⊗(δ::Matrix, F::AbstractArray, dim::Int)
        s = size(F)
        out = similar(F)
        if dim == 1
            for j in 1:s[3]
                for i in 1:s[2]
                    @views out[:,i,j] = δ*F[:,i,j]
                end
            end
        elseif dim == 2
            for j in 1:s[3]
                for i in 1:s[1]
                    @views out[i,:,j] = δ*F[i,:,j]
                end
            end
        elseif dim == 3
            for j in 1:s[2]
                for i in 1:s[1]
                    @views out[i,j,:] = δ*F[i,j,:]
                end
            end
        end
        return out
    end

curl(F, δ, dims) = ⊗(δ[dims[2]], F[dims[1]], dims[2]) .-
                        ⊗(δ[dims[1]], F[dims[2]],dims[1])

@inline function cover3d(A::AbstractVector, dim::Integer)
    extension = ones(Int64, 3)
    if dim == 1
        perm = (1, 2, 3)
    elseif dim == 2
        perm = (2, 1, 3)
    elseif dim == 3
        perm = (2, 3, 1)
    else
        throw("dim must be less than 3")
    end
    return permutedims(repeat(A, 1, 1, 1), perm)
end

@inline function cover2d(A::AbstractVector, dim::Integer)
    extension = ones(Int64, 2)
    if dim == 1
        perm = (1, 2)
    elseif dim == 2
        perm = (2, 1)
    else
        throw("dim must be less than 3")
    end
    return permutedims(repeat(A, 1, 1), perm)
end

perpsize(A::AbstractArray, i::Int) = size(A)[setdiff(1:ndims(A), i)]


text:
using LinearAlgebraicRepresentation
using SparseArrays
Lar = LinearAlgebraicRepresentation
L = Lar
using ViewerGL
GL = ViewerGL

function twocubegrids(n,m,p)
    V,(VV,EV,FV,CV) = Lar.cuboidGrid([n,m,p],true)
    mybox = (V,CV,FV,EV)

    twocubs = Lar.Struct([mybox, L.t(.3,.4,.5), L.r(pi/5,0,0), L.r(0,0,pi/12), mybox])

    V,CV,FV,EV = Lar.struct2lar(twocubs)
    GL.VIEW([ GL.GLGrid(V,FV, GL.COLORS[1], 0.5) ]);

    cop_EV = Lar.coboundary_0(EV::Lar.Cells);
    cop_EW = convert(Lar.ChainOp, cop_EV);
    cop_FE = Lar.coboundary_1(V, FV::Lar.Cells, EV::Lar.Cells);
    W = convert(Lar.Points, V');

    V, copEV, copFE, copCF = Lar.space_arrangement( W, cop_EW, cop_FE)

    EV = Lar.cop2lar(copEV)
    FE = Lar.cop2lar(copFE)
    CF = Lar.cop2lar(copCF)
    FV = [union([EV[e] for e in fe]...) for fe in FE]
    CV = [union([FV[f] for f in cf]...) for cf in CF]
    W = convert(Lar.Points, V')
    WW = [[k] for k=1:size(W,2)]

    GL.VIEW(GL.numbering(.5)((W,[WW,EV,FV,CV])));

    V,CVs,FVs,EVs = Lar.pols2tria(W, copEV, copFE, copCF)
    GL.VIEW(GL.GLExplode(V,FVs,1.2,1.2,1.2,99,1));
    GL.VIEW(GL.GLExplode(V,EVs,1.2,1.2,1.2,99,1));
    GL.VIEW(GL.GLExplode(V,CVs[2:end],5,5,5,99,0.5));

    return V, copEV, copFE, copCF
end

twocubegrids(3,3,3);


text:
println("Starting runtests.jl ...")
using ImagePlot
using Base.Test

# using TestImages
# img = plot(testimage("lighthouse"), 30*randn(2,10000).+[80;600])
# img = plot(img, [60 60 90; 580 620 600], radius = 6, fillcolor = "black")

s = (300,200,3)
a = rand(Float32, 300,200)
a3 = rand(Float32, s...)
p = randn(2,10)+50

img = plot(a, p)
@test size(img) == s

img = plot(a3, p)
@test size(img) == s

img = plot(a, p, 1:10)
@test size(img) == s

img = plot(a, p, [1,2,3,1,2,3,1,2,3,4])
@test size(img) == s

println("runtests done.")


text:
using HDF5

"""
Syntactic sugar for h5open
"""
function load_h5(path::String, image_id::String)
    image = h5open(path, "r") do file
        read(file, image_id)
    end
    return image
end

max_stars_repo_name: jjgomezcadenas/JImages

max_issues_repo_name: jjgomezcadenas/JImages

max_forks_repo_name: jjgomezcadenas/JImages


text:
using
    PDFEstimators,
    DataFrames,
    Distributions,
    MLJBase,
    Test

@testset "Overall" begin
    include("core.jl")
    include("actual.jl")
    include("kde.jl")
    include("ash.jl")
    include("mvn.jl")
    include("ffjord.jl")
    include("utils.jl")
end


text:
using CounterfactualRegret
const CFR = CounterfactualRegret
using StaticArrays
using Random
using RecipesBase
using LinearAlgebra
using Test

Random.seed!(1337)

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

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

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

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

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

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

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

max_stars_repo_name: WhiffleFish/HelloCFR

max_issues_repo_name: WhiffleFish/HelloCFR

max_forks_repo_name: WhiffleFish/HelloCFR


text:
# We are going to make the conjugate vectors more conjugate.

mutable struct ConjGradModified{T <: Number}
    A::Function              # Linear opeartor
    b::AbstractArray{T}         # RHS vector
    tensor_size::Tuple       # The size of the tensor the linear operator is acting on.

    x::Vector{T}        # current guess, started with initial guess.
    r::Vector{T}        # previous computed residual
    rnew::Vector{T}     # Current residual
    d::Vector{T}        # Conjugate Direction
    Ad::Vector{T}       # for reducing garbage collector time.

    itr::UInt64              # Iteration count.
    Q::Union{Nothing, Vector{Vector{T}}}       # Re-Orthogonalization basis.

    storage_limit::UInt64    # storage limit for the Q vector, default is n - 1
    reorthogonalize::Bool    # Whether to perform reorthogonalization.
    ortho_period::Int64

    function ConjGradModified(
        A::Function,
        b::AbstractArray,
        x0::Union{AbstractArray,Nothing}=nothing
    )
        x0 = x0===nothing ? b .+ 0.1 : x0
        r = reshape(b - A(x0), :)
        T = typeof(r[1])            # type extraction
        this = new{T}()
        this.r = r
        this.A = A
        this.b = b
        this.tensor_size = size(b)
        this.x = x0
        # this.r = ComputeResidualVec(this, this.x)
        this.rnew = similar(this.r)
        this.d = this.r
        this.Ad = similar(this.r)
        this.x = x0
        this.itr = 0

        this.Q = Vector{Vector{T}}()
        push!(this.Q, this.r/norm(this.r))

        this.storage_limit = length(this.r) - 1       # Maximal limit, if not, we have a problem.
        this.reorthogonalize = true
        this.ortho_period = 1
        return this
    end

    function ConjGradModified(A::AbstractMatrix, b::AbstractVector, x0::Union{AbstractArray,Nothing}=nothing)
        return ConjGradModified((x)-> A*x, b, x0)
    end

end

"""
    Apply the linear transformation on the tensor and return the
    residual vector after transformation.
"""
function ComputeResidualVec(this::ConjGradModified, x::AbstractArray)
    resVec = this.b - this.A(reshape(x, this.tensor_size))
    return reshape(resVec, length(resVec))
end

"""
    Get the norm of the current residual of Modified Conjugate
    Gradient.
"""
function GetResidualNorm(this::ConjGradModified)
    return norm(this.r)
end

"""
    Compute the vector using the linear transformation.
"""
function ComputeVec(this::ConjGradModified, x::AbstractArray)
    vec = this.A(reshape(x, this.tensor_size))
    return reshape(vec, length(vec))
end

"""
    Pre-allocates a vector and get the values to the given
    pre-allocated vector.
"""
function ComputeVec!(this::ConjGradModified, x::AbstractArray, vec::AbstractArray)
    vec .= this.A(reshape(x, this.tensor_size))
    return reshape(vec, :)
end

"""
    Performs One step of conjugate Gradient.
"""
function (this::ConjGradModified)()
    r = this.r
    if norm(r) == 0
        return 0 # The problem is solved already.
    end
    d = this.d
    Ad = this.Ad
    ComputeVec!(this, d, Ad)
    a = dot(r, r)/dot(d, Ad)

    if a < 0
        error("CG got a non-definite matrix")
    end

    this.x += a*d
    this.rnew .= r - a*Ad                    # update rnew
    rnewNorm = norm(this.rnew)

    if this.reorthogonalize
        orErr = abs(dot(this.Q[1], this.rnew/rnewNorm))
        if this.itr%this.ortho_period == 0
            println("Itr = $(this.itr); re-orthog error $orErr")
            for q in this.Q
                this.rnew .-= dot(q, this.rnew)*q
            end
            this.d = this.rnew  # Basically restart
        else
            b = dot(this.rnew, this.rnew)/dot(r, r)
            this.d = this.rnew + b*d
        end

        # manage the vectors.
        if length(this.Q) == this.storage_limit
            newQ = popfirst!(this.Q)
            newQ .= this.rnew/rnewNorm
            push!(this.Q, newQ)
        else
            push!(this.Q, this.rnew/rnewNorm)
        end
    end

    # @assert abs(dot(rnew + β*d, Ad)) < 1e-8 "Not conjugate"

    this.r = copy(this.rnew)                      # Override
    this.itr += 1
    return convert(Float64, rnewNorm)
end

"""
    Performs j iterations of conjugate gradient, collect the
    2norm of the residual and return all of them in a vector.
"""
function (this::ConjGradModified)(j::Int64)
    results = Vector()
    for _ in 1: j
        push!(results, this())
    end
    return results
end

### Using stuff.

using LinearAlgebra


text:
## ------------------------------------------------------------------
# Meta

get_sec(q::RBQuote) = get_parent(q)
set_sec!(q::RBQuote, sec::RBSection) = set_parent!(q, sec)

## ------------------------------------------------------------------
# Data

get_text(q::RBQuote) = getproperty(q, :text, "")
set_text!(q::RBQuote, text::String) = setproperty!(q, :text, strip(text))

## ------------------------------------------------------------------
# show
function Base.show(io::IO, q::RBQuote)
    println(io, "RBQuote(\"", get_label(q), "\")")
    print(io, " - path: '", localrelpath(q), "'")
    text = get_text(q)
    !isempty(text) && print(io, "\n - text: \"", _preview(io, text), "\"")
    tags = join(get_tags(q), ", ")
    !isempty(tags) && print(io, "\n - tags: \"", _preview(io, tags), "\"")
    print(io, "\n", _preview(io, "-"^70))
    return nothing
end


text:
# Function to translate parameter struct vector of parameters

"""
    par2θ(θ, model)
Function two convert parameter struct to vector

* `θ`: Parameter struct
* `model`: Model specification CountModel, INGARCHModel, ...

# Example
```julia-repl
model = CountModel(past_obs = 1) # Specify INARCH(1)
pars = θ2par([10, 0.5], model)
θ = par2θ(pars, model)
```
"""
function par2θ(θ::parameter, model::T where T<:CountModel)
    if model <: INARMA
        nϕ = sum(model.distr .== "NegativeBinomial")
        if nϕ == 2
            if length(model.external) == 0
                nϕ -= 1
            else
                if sum(model.external) == 0
                    nϕ -= 1
                end
            end
        end
    end

    if model <: INGARCH
        nϕ = model.distr == "NegativeBinomial"
    end

    [θ.β0; θ.α; θ.β; θ.η; θ.ϕ[1:nϕ]; ifelse(model.zi, θ.ω, Array{Float64, 1}([]))]
end

function par2θ(θ::parameter, model::T where T<:INGARCH)
    nϕ = model.distr == "NegativeBinomial"

    [θ.β0; θ.α; θ.β; θ.η; θ.ϕ[1:nϕ]; ifelse(model.zi, θ.ω, Array{Float64, 1}([]))]
end

function par2θ(θ::parameter, model::IIDModel)
    nb = model.distr == "NegativeBinomial"
    [θ.β0; θ.η; ifelse(nb, θ.ϕ, Array{Float64, 1}([])); ifelse(model.zi, θ.ω, Array{Float64, 1}([]))]
end

function par2θ(θ::parameter, model::INGARCHModel)
    nb = model.distr == "NegativeBinomial"
    [θ.β0; θ.α; θ.β; θ.η; ifelse(nb, θ.ϕ, Array{Float64, 1}([])); ifelse(model.zi, θ.ω, Array{Float64, 1}([]))]
end

function par2θ(θ::parameter, model::INARCHModel)
    nb = model.distr == "NegativeBinomial"
    [θ.β0; θ.α; θ.η; ifelse(nb, θ.ϕ, Array{Float64, 1}([])); ifelse(model.zi, θ.ω, Array{Float64, 1}([]))]
end

function par2θ(θ::parameter, model::T where T<:INARMA)
    nϕ = sum(model.distr .== "NegativeBinomial")
    nϕ = sum(model.distr .== "NegativeBinomial")
    if nϕ == 2
        if length(model.external) == 0
            nϕ -= 1
        else
            if sum(model.external) == 0
                nϕ -= 1
            end
        end
    end
    [θ.β0; θ.α; θ.β; θ.η; θ.ϕ[1:nϕ]; ifelse(model.zi, θ.ω, Array{Float64, 1}([]))]
end


text:
function bob(stimulus)
    without_spaces = filter(!isspace, stimulus)
    letters = filter(isletter, without_spaces)
    if isempty(without_spaces)
        # silence
        return "Fine. Be that way!"
    elseif last(without_spaces) == '?'
        if !isempty(letters) && all(isuppercase(c) for c in letters)
            return "Calm down, I know what I'm doing!"
        else
            return "Sure."
        end
    elseif !isempty(letters) && all(isuppercase(c) for c in letters)
        return "Whoa, chill out!"
    else
        return "Whatever."
    end
end

max_stars_repo_name: AbstractGeek/exercism-expeditions

max_issues_repo_name: AbstractGeek/exercism-expeditions

max_forks_repo_name: AbstractGeek/exercism-expeditions


text:
"""
    ThreeStateMorse()

Three state morse potential referred to as Model IA here:
[J. Chem. Phys. 150, 244102 (2019)](https://doi.org/10.1063/1.5096276)

Models IB and IC retain the same functional form and need only a change of parameters.
"""
Parameters.@with_kw struct ThreeStateMorse <: DiabaticModel
    d1::Float64 = 0.02
    d2::Float64 = 0.02
    d3::Float64 = 0.003

    α1::Float64 = 0.4
    α2::Float64 = 0.65
    α3::Float64 = 0.65

    r1::Float64 = 4.0
    r2::Float64 = 4.5
    r3::Float64 = 6.0

    c1::Float64 = 0.02
    c2::Float64 = 0.0
    c3::Float64 = 0.02

    a12::Float64 = 0.005
    a13::Float64 = 0.005
    a23::Float64 = 0.0

    α12::Float64 = 32.0
    α13::Float64 = 32.0
    α23::Float64 = 0.0

    r12::Float64 = 3.40
    r13::Float64 = 4.97
    r23::Float64 = 0.0
end

NQCModels.ndofs(::ThreeStateMorse) = 1
NQCModels.nstates(::ThreeStateMorse) = 3

function NQCModels.potential(model::ThreeStateMorse, R::Real)
    V_ii(x, d, α, r, c) = d * (1 - exp(-α*(x-r)))^2 + c
    V_ij(x, a, α, r) = a * exp(-α*(x-r)^2)

    V11 = V_ii(R, model.d1, model.α1, model.r1, model.c1)
    V22 = V_ii(R, model.d2, model.α2, model.r2, model.c2)
    V33 = V_ii(R, model.d3, model.α3, model.r3, model.c3)

    V12 = V_ij(R, model.a12, model.α12, model.r12)
    V13 = V_ij(R, model.a13, model.α13, model.r13)
    V23 = V_ij(R, model.a23, model.α23, model.r23)

    return Hermitian(SMatrix{3,3}(V11, V12, V13, V12, V22, V23, V13, V23, V33))
end

function NQCModels.derivative(model::ThreeStateMorse, R::Real)

    function D_ii(x, d, α, r)
        ex = exp(-α*(x-r))
        return 2 * d * α * (ex - ex^2)
    end

    D_ij(x, a, α, r) = -2 * a * α * (x-r) * exp(-α*(x-r)^2)

    D11 = D_ii(R, model.d1, model.α1, model.r1)
    D22 = D_ii(R, model.d2, model.α2, model.r2)
    D33 = D_ii(R, model.d3, model.α3, model.r3)

    D12 = D_ij(R, model.a12, model.α12, model.r12)
    D13 = D_ij(R, model.a13, model.α13, model.r13)
    D23 = D_ij(R, model.a23, model.α23, model.r23)

    return Hermitian(SMatrix{3,3}(D11, D12, D13, D12, D22, D23, D13, D23, D33))

end


text:
@testset "lz4_compression" begin
    text = b"""
    Lorem ipsum dolor sit amet, consectetur adipiscing elit. Aenean sollicitudin
    mauris non nisi consectetur, a dapibus urna pretium. Vestibulum non posuere
    erat. Donec luctus a turpis eget aliquet. Cras tristique iaculis ex, eu
    malesuada sem interdum sed. Vestibulum ante ipsum primis in faucibus orci luctus
    et ultrices posuere cubilia Curae; Etiam volutpat, risus nec gravida ultricies,
    erat ex bibendum ipsum, sed varius ipsum ipsum vitae dui.
    """

    @testset "Transcoding" begin
        @test LZ4FastCompressorStream <: TranscodingStream
        @test LZ4SafeDecompressorStream <: TranscodingStream

        compressed = transcode(LZ4FastCompressor, text)
        @test sizeof(compressed) < sizeof(text)

        decompressed = transcode(LZ4SafeDecompressor, compressed)
        @test sizeof(decompressed) > sizeof(compressed)
        @test decompressed == Vector{UInt8}(text)

        test_roundtrip_fileio(LZ4FastCompressor, LZ4SafeDecompressor)
        test_roundtrip_transcode(LZ4FastCompressor, LZ4SafeDecompressor)

        file = IOBuffer(text)
        stream = LZ4SafeDecompressorStream(LZ4FastCompressorStream(file))
        flush(stream)

        @test read(stream) == text
        close(stream)
        close(file)

        file = IOBuffer(text)
        stream = LZ4SafeDecompressorStream(LZ4FastCompressorStream(file; acceleration = 5))
        flush(stream)

        @test read(stream) == text
        close(stream)
        close(file)

        file = IOBuffer("")
        stream = LZ4SafeDecompressorStream(LZ4FastCompressorStream(file))
        flush(stream)

        @test read(stream) == b""
        close(stream)
        close(file)

        teststring = rand(UInt8, 10000)
        file = IOBuffer(teststring)
        stream = LZ4SafeDecompressorStream(LZ4FastCompressorStream(file; block_size = 2048); block_size = 2048)
        flush(stream)

        @test read(stream) == teststring
        close(stream)
        close(file)
    end

    @testset "Errors" begin
        @testset "Uninitialized" begin
            input = Memory(Vector{UInt8}(text))
            output = Memory(Vector{UInt8}(undef, 1280))
            not_initialized = LZ4FastCompressor()
            @test TranscodingStreams.startproc(not_initialized, :read, Error()) == :error
            @test TranscodingStreams.process(not_initialized, input, output, Error()) == (0, 0, :error)

            compressed = Memory(transcode(LZ4FastCompressor, text))
            not_initialized = LZ4SafeDecompressor()
            @test TranscodingStreams.startproc(not_initialized, :read, Error()) == :error
            @test TranscodingStreams.process(not_initialized, compressed, output, Error()) == (0, 0, :error)
        end

        @testset "Bad Input" begin
            # Malformed decompression input
            @test_throws CodecLz4.LZ4Exception transcode(LZ4SafeDecompressor, text)
            @test_throws BoundsError transcode(LZ4SafeDecompressor, [0x00])

            # Properly compressed but not formatted as a stream
            compressed = lz4_compress(text)
            @test_throws CodecLz4.LZ4Exception transcode(LZ4SafeDecompressor, text)

            # Block size too large
            @test_throws ArgumentError LZ4FastCompressor(; block_size = CodecLz4.LZ4_MAX_INPUT_SIZE + 1)
        end

        @testset "Bad Buffer Size" begin
            # Decompression with too-small block_size
            output = Memory(Vector{UInt8}(undef, 1024))
            compressed = transcode(LZ4FastCompressor, text)
            decompressor = LZ4SafeDecompressor(; block_size = 200)
            try
                @test_nowarn TranscodingStreams.initialize(decompressor)
                @test TranscodingStreams.startproc(decompressor, :read, Error()) == :ok
                err = Error()
                @test TranscodingStreams.process(decompressor, Memory(compressed), output, err) == (0, 0, :error)

                @test err[] isa CodecLz4.LZ4Exception
                @test err[].msg == "Decompression failed."
            finally
                TranscodingStreams.finalize(decompressor)
            end

            # Compression into too-small buffer
            input = Memory(Vector{UInt8}(text))
            output = Memory(Vector{UInt8}(undef, 1))
            compressor = LZ4FastCompressor()

            try
                @test_nowarn TranscodingStreams.initialize(compressor)
                @test TranscodingStreams.startproc(compressor, :read, Error()) == :ok
                err = Error()
                @test TranscodingStreams.process(compressor, input, output, err) == (0, 0, :error)

                @test err[] isa BoundsError
            finally
                TranscodingStreams.finalize(compressor)
            end

            # Decompression into too-small buffer
            output = Memory(Vector{UInt8}(undef, 1))
            compressed = transcode(LZ4FastCompressor, text)
            decompressor = LZ4SafeDecompressor()
            try
                @test_nowarn TranscodingStreams.initialize(decompressor)
                @test TranscodingStreams.startproc(decompressor, :read, Error()) == :ok
                err = Error()
                @test TranscodingStreams.process(decompressor, Memory(compressed), output, err) == (0, 0, :error)

                @test err[] isa BoundsError
            finally
                TranscodingStreams.finalize(decompressor)
            end
        end
    end

    @testset "dst size fix" begin
        teststring = rand(UInt8, 800000)
        io = IOBuffer(teststring)
        stream = LZ4FastCompressorStream(io)
        result = read(stream)
        @test_nowarn close(stream)
        close(io)

        io = IOBuffer(result)
        stream = LZ4SafeDecompressorStream(io)
        result = read(stream)
        @test result == teststring
        @test_nowarn close(stream)
        close(io)
    end
end


text:
if !isdefined(Main, :CustomStencil)
    include("../module/CustomStencil.jl")
    using Main.CustomStencil
end
include("misc.jl")
coefs = Float64.([1 1 1 1 1])

star_stencil = @def_stencil_expression c[0]D[x,y,z] + @sum(i, 1, 4, c[i]*(
            D[x+i,y,z] + D[x,y+i,z] + D[x,y,z+i] +
            D[x-i,y,z] + D[x,y-i,z] + D[x,y,z-i]))

st_inst = NewStencilInstance(star_stencil, 4, coefs)

radius = 4
nx = 5
ny = 5
nz = 5
dx = 1<<(nx)
dy = 1<<(ny)
dz = 1<<(nz)

t_steps = 5
## Actual time steps
data = CreateData(dx,dy,dz,radius)
dout = ApplyStencil(st_inst, data, t_steps)
#heatmap(dout1[:,:,Int(round(dz/2))])
o = Array(dout)

function cpu_star_stencil(data, radius,coefs)
    out = zeros(Float32,size(data))
    for i = (radius+1):(size(data,1)-radius)
        for j = (radius+1):(size(data,2)-radius)
            for k = (radius+1):(size(data,3)-radius)
                c = coefs[1]*data[i,j,k]
                for r = 1:radius
                    c += coefs[r+1]*(data[i + r,j,k] + data[i - r,j,k] +
                    data[i,j + r,k] + data[i,j - r,k] +
                    data[i,j,k + r] + data[i,j,k - r])
                end
                out[i,j,k] = c
            end
        end
    end
    return out
end

cout = cpu_star_stencil(data, radius, coefs)
for i = 2:t_steps
    global cout
    cout = cpu_star_stencil(cout, radius, coefs)
end

@assert(sum(abs.(cout .- o))./(reduce(*, size(data))) < 0.001)
mean_err = sum(abs.(cout .- o))./(reduce(*, size(data)))
# using Plots
# shift = -13
# heatmap(cout[:,:,Int(round(dz/2))+shift])
# heatmap(o[:,:,Int(round(dz/2))+shift])


text:
import LinearAlgebra: dot, norm
import Base: +, -, inv, *, transpose, sqrt, ones
export dot, norm, add, add!, product, -, +, *, transpose
export hadamard, inv, sqrt, ones, norm

function ones(x::PointE{T, matT}) where {T, matT}
    mats = matT[]
    for dim in x.dims
        mat = matT(zeros(T, dim, dim))
        for i=1:dim
            mat[i, i] = convert(T, 1.0)
        end
        push!(mats, mat)
    end

    return PointE(mats, ones(T, length(x.vec)))
end

function ones(x::PointE{T, matT}, matrixT) where {T, matT}
    mats = matrixT[]
    for dim in x.dims
        mat = matrixT(zeros(T, dim, dim))
        for i=1:dim
            mat[i, i] = convert(T, 1.0)
        end
        push!(mats, mat)
    end

    return PointE(mats, ones(T, length(x.vec)))
end

"""
Scalar product over E, primal space
"""
function dot(pt1::PointE{T}, pt2::PointE{T}) where {T<:Number}
    @assert pt1.dims == pt2.dims

    innerprod = dot(pt1.vec, pt2.vec)

    for n=1:length(pt1.dims)
        innerprod += dot(pt1.mats[n], pt2.mats[n])
    end

    return innerprod
end

function norm(x::PointE{T, matT}, p::Real) where {T<:Number, matT <: AbstractArray}
    if p == 2
        return sqrt(dot(x, x))

    elseif p == Inf
        sup = maximum(x.vec)

        for mat in x.mats
            sup = max(sup, norm(mat, Inf))
        end

        return sup
    end
    @error "norm(): unhandled p"
end

################################################################################
### Additive algebra
################################################################################

function add!(x::PointE{T}, y::PointE{T}) where {T<:Number}
    for (i, mat) in enumerate(y.mats)
        x.mats[i] = add(x.mats[i], mat)
    end

    for (i, yi) in enumerate(y.vec)
        x.vec[i] += yi
    end
    nothing
end

function add!(x::Array{Float64,1}, y::Array{Float64,1})
    x += y
    nothing
end

## Matrix addition
function add(mat1::AbstractArray{T}, mat2::AbstractArray{T}) where {T<:Number}
    return mat1 + mat2
end

function add(mat1::Sparse{T}, mat2::Sparse{T}) where {T<:Number}
    return mat1 + mat2
end

function add(mat1::Symmetric{T}, mat2::Symmetric{T}) where {T<:Number}
    if mat1.uplo == mat2.uplo
        mat1.data += mat2.data
    else
        mat1.data += transpose(mat2.data)
    end
    nothing
end

# function add!(mat1::DenseHerm{T}, mat2::DenseHerm{T}) where {T<:Number}
#     if mat1.uplo == mat2.uplo
#         mat1.data += mat2.data
#     else
#         mat1.data += adjoint(mat2.data)
#     end
#     nothing
# end

function add(x::PointE{T, U}, y::PointE{T, U}) where {T<:Number, U<:AbstractArray}
    point = PointE(x.dims, length(x.vec), T, U)

    for i=1:length(x.dims)
        point.mats[i] = x.mats[i] + y.mats[i]
    end

    point.vec = x.vec .+ y.vec
    return point
end

function add(x::PointE{T, U}, y::PointE{T, V}) where {T<:Number, U<:AbstractArray, V<:AbstractSparseArray}
    point = PointE(x.dims, length(x.vec), T, U)

    for i=1:length(x.dims)
        point.mats[i] = x.mats[i] + y.mats[i]
    end

    point.vec = x.vec .+ y.vec
    return point
end

function add(x::PointE{T, U}, y::PointE{T, V}) where {T<:Number, U<:AbstractSparseArray, V<:AbstractArray}
    return add(y, x)
end

function add(x::PointE{T, U}, y::PointE{T, V}) where {T<:Number, U<:AbstractArray, V<:Symmetric}
    point = PointE(x.dims, length(x.vec), T, U)

    for i=1:length(x.dims)
        point.mats[i] = x.mats[i] + y.mats[i]
    end

    point.vec = x.vec .+ y.vec
    return point
end

function add(x::PointE{T, U}, y::PointE{T, V}) where {T<:Number, U<:Symmetric, V<:AbstractArray}
    return add(y, x)
end

+(x::PointE, y::PointE) = add(x, y)
-(x::PointE, y::PointE) = add(x, -y)

################################################################################
### Multiplicative algebra
################################################################################

function product!(x::PointE{T}, lambda::T) where {T<:Number}
    for (i, mat) in enumerate(x.mats)
        x.mats[i] *= lambda
    end

    x.vec *= lambda
    nothing
end

function product(x::PointE{T}, lambda::T) where {T<:Number}
    xprod = copy(x)
    product!(xprod, lambda)
    return xprod
end

*(x::PointE{T, U}, lambda::T) where {T<:Number, U} = product(x, lambda)
*(lambda::T, x::PointE{T, U}) where {T<:Number, U} = product(x, lambda)

-(x::PointE{T}) where T<:Number = product(x, convert(T, -1))

################################################################################
### Factorization ops: inverse, sqrt
################################################################################
"""
    Compute the inverse of E, assuming `x` lies in strict feasible space
"""
function inv(x::PointE{T}) where T<:Number
    xinv = PointE(x.dims, length(x.vec), T, Dense{T})

    for (i, mat) in enumerate(x.mats)
        xinv.mats[i] = inv(factorize(mat))
    end

    xinv.vec = 1 ./ x.vec
    return xinv
end

"""
    Compute the square root of E, assuming `x` lies in strict feasible space
"""
function sqrt(x::PointE{T, Union{Dense{T}, DenseSym{T}}}) where T<:Number
    xsqrt = PointE(c.dims, length(x.vec), T, Dense{T})

    for (i, mat) in enumerate(x.mats)
        xsqrt.mats[i] = sqrt(factorize(mat))
    end

    xsqrt.vec = sqrt.(x.vec)
    return xsqrt
end

################################################################################
### Hadamard operations
################################################################################
"""
    hadamard(x::PointE{T, U}, y::PointE{T, U}) where {T<:Number, U<:SDCOAbstractMatrix}

Compute the hadamard product of the points x and y, returning a dense matrix object.
"""
function hadamard(x::PointE{T}, y::PointE{T}) where {T<:Number}
    res = PointE(x.dims, length(x.vec), T, Dense{T})

    for (i, mati) in enumerate(x.mats)
        res.mats[i] = x.mats[i] * y.mats[i]
    end

    res.vec = x.vec .* y.vec
    return res
end

function hadamard(x::PointE{T}, y::PointE{T}, z::PointE{T}; transposefirst=false) where T<:Number
    res = PointE(x.dims, length(x.vec), T, Dense{T})

    hadamard!(x, y, z, res, transposefirst=transposefirst)
    return res
end

function hadamard!(x::PointE{T}, y::PointE{T}, z::PointE{T}, out::PointE{T, Dense{T}}; transposefirst=false) where T<:Number
    if transposefirst
        for n=1:length(x.dims)
            out.mats[n] = transpose(x.mats[n])
            out.mats[n] *= y.mats[n]
            out.mats[n] *= z.mats[n]
        end
    else
        for n=1:length(x.dims)
            out.mats[n] = x.mats[n] * y.mats[n] * z.mats[n]
        end
    end

    out.vec = x.vec .* y.vec .* z.vec

    nothing
end

#####################################################################"
## Point Primal Dual
#####################################################################"
function add(z1::PointPrimalDual, z2::PointPrimalDual)
    return PointPrimalDual(z1.x+z2.x, z1.y+z2.y, z1.s+z2.s)
end

function add!(z1::PointPrimalDual, z2::PointPrimalDual)
    add!(z1.x, z2.x)
    z1.y += z2.y
    add!(z1.s, z2.s)
end

+(z1::PointPrimalDual, z2::PointPrimalDual) = add(z1, z2)
-(z1::PointPrimalDual, z2::PointPrimalDual) = add(z1, -z2)

function product(z::PointPrimalDual{T}, lambda::U) where {T<:Number, U<:Number}
    return PointPrimalDual(z.x * convert(T, lambda), z.y * convert(T, lambda), z.s * convert(T, lambda))
end

*(z::PointPrimalDual, lambda::T) where T<:Number = product(z, lambda)
*(lambda::T, z::PointPrimalDual) where T<:Number = product(z, lambda)

function -(z::PointPrimalDual{T}) where T<:Number
    return product(z, convert(T, -1))
end

function norm(z::PointPrimalDual)
    return norm(z.x) + norm(z.y) + norm(z.s)
end

#####################################################################
## transpose
#####################################################################
function transpose(x::PointE{T, Dense{T}}) where T<:Number
    mats = Vector{Dense{T}}([zeros(T, n, n) for n in x.dims])

    for (matind, mat) in enumerate(x.mats)
        mats[matind] = transpose(mat)
    end

    return PointE(mats, x.vec)
end


text:
# included in runtests.jl
using Test
import WebSockets.HTTP
import Base: convert, BufferStream
using WebSockets
import WebSockets:  generate_websocket_key, upgrade

include("logformat.jl")
include("handshaketest_functions.jl")

# Test generate_websocket_key"
@test generate_websocket_key("dGhlIHNhbXBsZSBub25jZQ==") == "s3pPLMBiTxaQ9kYGzzhZRbK+xOo="

# Test reject / switch format"
io = IOBuffer()
const REJECT = "HTTP/1.1 400 Bad Request"
Base.write(io, HTTP.Response(400))
@test takefirstline(io) == REJECT
Base.write(io, HTTP.Response(400))
@test takefirstline(io) == REJECT

const SWITCH = "HTTP/1.1 101 Switching Protocols"
Base.write(io, HTTP.Response(101))
@test takefirstline(io) == SWITCH
Base.write(io, HTTP.Response(101))
@test takefirstline(io) == SWITCH

# "Test simple handshakes that are unacceptable"

for r in templaterequests()
    @test handshakeresponse(r) == REJECT
    sethd(r, "Sec-WebSocket-Version"    => "11")
    @test handshakeresponse(r) == REJECT
    sethd(r, "Sec-WebSocket-Version"    => "13")
    @test handshakeresponse(r) == REJECT
    sethd(r,   "Sec-WebSocket-Key"      => "mumbojumbobo")
    @test handshakeresponse(r) == REJECT
    sethd(r,  "Sec-WebSocket-Key"       => "17 bytes key is not accepted")
    @test handshakeresponse(r) == REJECT
    sethd(r,  "Sec-WebSocket-Key"       => "AQIDBAUGBwgJCgsMDQ4PEA==")
    sethd(r,  "Sec-WebSocket-Protocol"  => "unsupported")
    @test handshakeresponse(r) == REJECT
end

#  Test simple handshakes, acceptable
for r in templaterequests()
    sethd(r, "Sec-WebSocket-Version"    => "13")
    sethd(r,  "Sec-WebSocket-Key"       => "AQIDBAUGBwgJCgsMDQ4PEA==")
    @test handshakeresponse(r) == SWITCH
end

#  Test unacceptable subprotocol handshake subprotocol
for r in templaterequests()
    sethd(r, "Sec-WebSocket-Version"    => "13")
    sethd(r, "Sec-WebSocket-Key"        => "AQIDBAUGBwgJCgsMDQ4PEA==")
    sethd(r, "Sec-WebSocket-Protocol"       => "my.server/json-zmq")
    @test handshakeresponse(r) == REJECT
end

#  add simple subprotocol to acceptable list
@test true == WebSockets.addsubproto("xml")

# add subprotocol with difficult name
@test true == WebSockets.addsubproto("my.server/json-zmq")

# "Test handshake subprotocol now acceptable"
for r in templaterequests()
    sethd(r, "Sec-WebSocket-Version"    => "13")
    sethd(r,  "Sec-WebSocket-Key"        => "AQIDBAUGBwgJCgsMDQ4PEA==")
    sethd(r, "Sec-WebSocket-Protocol"       => "xml")
    @test handshakeresponse(r) == SWITCH
    sethd(r, "Sec-WebSocket-Protocol"       => "my.server/json-zmq")
    @test handshakeresponse(r) == SWITCH
end


text:
using RedBlackGraph
using Test

function test_basic_properties(T::Type)
    zip = zero(T)
    @test typeof(zip) == T
    eins = one(T)
    @test typeof(eins) == T
    r1 = red_one(T)
    @test typeof(r1) == T

    @test r1 != eins
    @test r1 != zip

    @test typemin(T) == zip
    max = typemax(T)
    @test unsigned(max) < unsigned(r1)

    @test eins < zip
    @test r1 < zip
    @test r1 < eins

    @test zip > r1
    @test zip > eins
    @test eins > r1

    second_r1 = red_one(eins)
    @test second_r1 == r1

    @test iseven(r1)
    @test !iseven(eins)
    @test iseven(zip)
end

function test_basic_avus_sum(T1::Type, T2::Type)
    @test T1(0) + T2(1) == T2(1)
    @test T1(1) + T2(1) == T1(1)
    @test T1(1) + T2(2) == T1(1)
end

function test_basic_avus_product(T1::Type, T2::Type)
    @test T1(0) * T2(1) == 0
    @test T1(1) * T2(1) == 1
    @test T1(1) * T2(2) == 2
    @test T1(2) * T2(2) == 4
    @test T1(2) * T2(3) == 5
    @test T1(3) * T2(2) == 6
    @test T1(3) * T2(3) == 7
end

@testset "Basic Avus Tests" begin
    test_basic_properties(AInt8)
    test_basic_properties(AInt16)
    test_basic_properties(AInt32)
    test_basic_properties(AInt64)
    test_basic_properties(AInt128)

    test_basic_avus_sum(AInt8, AInt8)
    test_basic_avus_sum(AInt16, AInt16)
    test_basic_avus_sum(AInt32, AInt32)
    test_basic_avus_sum(AInt64, AInt64)
    test_basic_avus_sum(AInt128, AInt128)

    test_basic_avus_product(AInt8, AInt8)
    test_basic_avus_product(AInt16, AInt16)
    test_basic_avus_product(AInt32, AInt32)
    test_basic_avus_product(AInt64, AInt64)
    test_basic_avus_product(AInt128, AInt128)
end


text:
# ===========================================================================================
# Main API

using GLM
@reexport using GLM
import GLM: glm,
            # Model
            LinPredModel, AbstractGLM, GeneralizedLinearModel, LinearModel,
            LmResp, GlmResp,
            # Pred
            LinPred, DensePred,
            DensePredChol, SparsePredChol, QRCompactWY, SparseMatrixCSC,
            # prediction
            installbeta!, delbeta!, linpred, linpred!,
            updateμ!, cholfactors,
            # other
            FP, BlasReal, Link, dispersion, deviance, dof, dof_residual, nobs

anova(models::Vararg{TableRegressionModel{<: LinearModel, <: AbstractArray}, N};
        test::Type{T} = FTest,
        kwargs...) where {N, T <: GoodnessOfFit} =
    anova(test, models...; kwargs...)

anova(models::Vararg{TableRegressionModel{<: GeneralizedLinearModel, <: AbstractArray}, N};
        test::Type{T} = canonicalgoodnessoffit(models[1].model.rr.d),
        kwargs...) where {N, T <: GoodnessOfFit} =
    anova(test, models...; kwargs...)

# ==================================================================================================================
# ANOVA by F test
# LinearModels

anova(::Type{FTest},
    trm::TableRegressionModel{<: Union{LinearModel, GeneralizedLinearModel{<: GLM.GlmResp{T, <: Normal, IdentityLink}}}};
    type::Int = 1,
    kwargs...) where T = _anova_vcov(trm; type, kwargs...)

function _anova_vcov(trm::TableRegressionModel{<: Union{LinearModel, GeneralizedLinearModel}};
                    type::Int = 1, kwargs...)
    type in [1, 2, 3] || throw(ArgumentError("Invalid type"))

    assign = trm.mm.assign
    df = dof(assign)
    filter!(>(0), df)
    # May exist some floating point error from dof_residual
    push!(df, round(Int, dof_residual(trm)))
    df = tuple(df...)
    if type in [1, 3]
        # vcov methods
        varβ = vcov(trm.model)
        β = trm.model.pp.beta0
        if type == 1
            fs = abs2.(cholesky(Hermitian(inv(varβ))).U * β)
            offset = first(assign) == 1 ? 0 : 1
            fstat = ntuple(last(assign) - offset) do fix
                sum(fs[findall(==(fix + offset), assign)]) / df[fix]
            end
        else
            # calculate block by block
            offset = first(assign) == 1 ? 0 : 1
            fstat = ntuple(last(assign) - offset) do fix
                select = findall(==(fix + offset), assign)
                β[select]' * inv(varβ[select, select]) * β[select] / df[fix]
            end
        end
        σ² = dispersion(trm.model, true)
        devs = (fstat .* σ²..., σ²) .* df
    else
        # refit methods
        devs = deviances(trm; type, kwargs...)
        MSR = devs ./ df
        fstat = MSR[1:end - 1] ./ dispersion(trm.model, true)
    end
    pvalue = (ccdf.(FDist.(df[1:end - 1], last(df)), abs.(fstat))..., NaN)
    AnovaResult{FTest}(trm, type, df, devs, (fstat..., NaN), pvalue, NamedTuple())
end

function anova(::Type{FTest},
                trm::TableRegressionModel{<: Union{LinearModel, GeneralizedLinearModel}};
                type::Int = 1, kwargs...)
    type in [1, 2, 3] || throw(ArgumentError("Invalid type"))

    assign = trm.mm.assign
    devs = deviances(trm; type, kwargs...)
    df = dof(assign)
    filter!(>(0), df)
    # May exist some floating point error from dof_residual
    push!(df, round(Int, dof_residual(trm)))
    length(df) == length(devs) + 1 && popfirst!(df)
    df = tuple(df...)
    msr = devs ./ df
    fstat = msr[1:end - 1] ./ dispersion(trm.model, true)
    pvalue = (ccdf.(FDist.(df[1:end - 1], last(df)), abs.(fstat))..., NaN)
    AnovaResult{FTest}(trm, type, df, devs, (fstat..., NaN), pvalue, NamedTuple())
end

# ----------------------------------------------------------------------------------------
# ANOVA for genaralized linear models
# λ = -2ln(𝓛(̂θ₀)/𝓛(θ)) ~ χ²ₙ , n = difference of predictors

function anova(::Type{LRT},
            trm::TableRegressionModel{<: Union{LinearModel, GeneralizedLinearModel}})
    Δdev = deviances(trm, type = 1)
    df = dof(trm.mm.assign)
    filter!(>(0), df)
    isnullable(trm.model) || popfirst!(df)
    df = tuple(df...)
    # den = last(ss) / (nobs(trm) - dof(trm) + 1)
    # lrstat = ss[1:end - 1] ./ den
    σ² = dispersion(trm.model, true)
    lrstat = Δdev[1:end - 1] ./ σ²
    n = length(lrstat)
    dev = collect(Δdev)
    i = n
    while i > 0
        dev[i] += dev[i + 1]
        i -= 1
    end
    pval = ccdf.(Chisq.(df), abs.(lrstat))
    AnovaResult{LRT}(trm, 1, df, tuple(dev[2:end]...), lrstat, pval, NamedTuple())
end

# =================================================================================================================
# Nested models

function anova(::Type{FTest},
        trms::Vararg{TableRegressionModel{<: Union{LinearModel, GeneralizedLinearModel}}};
        check::Bool = true,
        isnested::Bool = false)
    df = dof.(trms)
    ord = sortperm(collect(df))
    df = df[ord]
    trms = trms[ord]

    # check comparable and nested
    check && @warn "Could not check whether models are nested: results may not be meaningful"

    Δdf = _diff(df)
    # May exist some floating point error from dof_residual
    dfr = round.(Int, dof_residual.(trms))
    dev = deviance.(trms)
    msr = _diffn(dev) ./Δdf
    σ² = dispersion(last(trms).model, true)
    fstat = msr ./ σ²
    pval = map(zip(Δdf, dfr[2:end], fstat)) do (dof, dofr, fs)
        fs > 0 ? ccdf(FDist(dof, dofr), fs) : NaN
    end
    AnovaResult{FTest}(trms, 1, df, dev, (NaN, fstat...), (NaN, pval...), NamedTuple())
end

function anova(::Type{LRT},
        trms::Vararg{<: TableRegressionModel};
        check::Bool = true,
        isnested::Bool = false)
    df = dof.(trms)
    ord = sortperm(collect(df))
    trms = trms[ord]
    # check comparable and nested
    df = df[ord]
    _lrt_nested(trms, df, deviance.(trms), dispersion(last(trms).model, true); nestedwarn = true)
end

# =================================================================================================================================
# Fit new models

"""
    anova_lm(X, y; test::Type{<: GoodnessOfFit} = FTest, <keyword arguments>)

    anova_lm(test::Type{<: GoodnessOfFit}, X, y; <keyword arguments>)

    anova(test::Type{<: GoodnessOfFit}, ::Type{LinearModel}, X, y;
        type::Int = 1,
        <keyword arguments>)

ANOVA for simple linear regression.

The arguments `X` and `y` can be a `Matrix` and a `Vector` or a `Formula` and a `DataFrame`.
* `type` specifies type of anova.
* `dropcollinear` controls whether or not `lm` accepts a model matrix which is less-than-full rank. If true (the default), only the first of each set of linearly-dependent columns is used. The coefficient for redundant linearly dependent columns is 0.0 and all associated statistics are set to NaN.

`anova_lm` generate a `TableRegressionModel` object, which is fitted by `lm`.
"""
anova_lm(X, y;
        test::Type{<: GoodnessOfFit} = FTest,
        kwargs...) =
    anova(test, LinearModel, X, y; kwargs...)

anova_lm(test::Type{<: GoodnessOfFit}, X, y; kwargs...) =
    anova(test, LinearModel, X, y; kwargs...)

function anova(test::Type{<: GoodnessOfFit}, ::Type{LinearModel}, X, y;
        type::Int = 1,
        kwargs...)
    model = lm(X, y; kwargs...)
    anova(test, model; type = type)
end

"""
    anova_glm(X, y, d::UnivariateDistribution, l::Link = canonicallink(d);
            test::Type{<: GoodnessOfFit} = canonicalgoodnessoffit(d), <keyword arguments>)

    anova_glm(test::Type{<: GoodnessOfFit}, X, y, d::UnivariateDistribution, l::Link = canonicallink(d); <keyword arguments>)

    anova(test::Type{<: GoodnessOfFit}, X, y, d::UnivariateDistribution, l::Link = canonicallink(d); <keyword arguments>)

ANOVA for genaralized linear models.

The arguments `X` and `y` can be a `Matrix` and a `Vector` or a `Formula` and a `DataFrame`.
* `d`: a `GLM.UnivariateDistribution`.
* `l`: a `GLM.Link`

For other keyword arguments, see `fit`.
"""
anova_glm(X, y,
        d::UnivariateDistribution, l::Link = canonicallink(d);
        test::Type{<: GoodnessOfFit} = canonicalgoodnessoffit(d),
        kwargs...) =
    anova(test, GeneralizedLinearModel, X, y, d, l; kwargs...)

anova_glm(test::Type{<: GoodnessOfFit}, X, y,
        d::UnivariateDistribution, l::Link = canonicallink(d);
        kwargs...) =
    anova(test, GeneralizedLinearModel, X, y, d, l; kwargs...)

function anova(test::Type{<: GoodnessOfFit}, ::Type{GeneralizedLinearModel}, X, y,
        d::UnivariateDistribution, l::Link = canonicallink(d);
        type::Int = 1,
        kwargs...)
    trm = glm(X, y, d, l; kwargs...)
    anova(test, trm; type, kwargs... )
end

"""
    GLM.glm(f::FormulaTerm, df::DataFrame, d::Binomial, l::GLM.Link, args...; kwargs...)

Automatically transform dependent variable into 0/1 for family `Binomial`.
"""
GLM.glm(f::FormulaTerm, df::DataFrame, d::Binomial, l::Link, args...; kwargs...) =
    fit(GeneralizedLinearModel, f,
        combine(df, :, f.lhs.sym => ByRow(==(last(unique(df[!, f.lhs.sym])))) => f.lhs.sym),
        d, l, args...; kwargs...)


text:
################################################################################
#
#	CONCRETE TYPE
#
#   Site{L,D} <: AbstractSite{L,D}
#   --> L is the label type
#   --> D is the dimension of embedding space
#
#   FILE CONTAINS
#       - concrete struct definition
#       - interface implementation
#
################################################################################

################################################################################
#
#   CONCRETE STRUCT DEFINITION
#
################################################################################
mutable struct Site{L,D} <: AbstractSite{L,D}

	# point
	point	:: Vector{Float64}

	# label
	label	:: L

end

# export the concrete type
export Site

################################################################################
#
#	IMPLEMENTATION OF INTERFACE FOR CONCRETE SITE TYPE
#	(functions that had to be overwritten by concrete type)
#
################################################################################

# default constructor interface
# used for creation of new sites
function newSite(
            :: Type{Site{L,D}},
            point   :: Vector{<:Real},
            label   :: L
        ) :: Site{L,D} where {L,D}

    # check if correct dimension is given
    @assert length(point) == D

    # return a new object
    return Site{L,D}(point, label)
end

# get label
function label(
            s :: Site{L,D}
        ) :: L where {L,D}

    # return the label
    return s.label
end
# set label
function label!(
            s :: Site{L,D},
            l :: L
        ) where {L,D}

    # return the label
    s.label = l
end

# get point
function point(
            s :: Site{L,D}
        ) :: Vector{Float64} where {L,D}

    # return the point
    return s.point
end
# set point
function point!(
            s :: Site{L,D},
            p :: Vector{<:Real}
        ) where {L,D}

    # return the point
    s.point = p
end

function saveSites(
        ss :: Vector{S},
        fn :: AbstractString,
        group :: AbstractString = "sites"
        ;
        append :: Bool = false
    ) where {LS,D,S<:Site{LS,D}}

    # determine the mode based on if one wants to append stuff
    if append
        mode = "r+"
    else
        mode = "w"
    end

    # open the file in mode
    h5open(fn, mode) do ucfile
        # create the group in which the bonds are saved
        group_sites = g_create(ucfile, group)
        # save the parameters
        attrs(group_sites)["D"] = Int64(D)
        attrs(group_sites)["L"] = string(LS)
        # save all Positions (D dimensions)
        if 0 < Int64(D)
            for n in 1:Int64(D)
                group_sites["point_$(n)"] = Float64[point(s)[n] for s in ss]
            end
        end
        # save all labels
        if LS <: Number
            group_sites["label"] = [label(s) for s in ss]
        else
            group_sites["label"] = String[string(label(s)) for s in ss]
        end
    end

    # return nothing
    return nothing
end

function loadSites(
        ::Type{SI},
        fn :: AbstractString,
        group :: AbstractString = "sites"
    ) where {LSI,DI,SI<:Union{Site{LSI,DI},Site}}

    # read attribute data
    attr_data = h5readattr(fn, group)
    # determine D and L based on this
    D  = attr_data["D"]
    LS = Meta.eval(Meta.parse(attr_data["L"]))

    # load all remaining data
    st_label = h5read(fn, group*"/label")

    if D == 0
        st_point = [Float64[] for i in 1:length(st_label)]
    else
        st_point_parts = [
            h5read(fn, group*"/point_"*string(j)) for j in 1:D
        ]
        st_point = Vector{Float64}[Float64[st_point_parts[j][i] for j in 1:D] for i in 1:length(st_label)]
    end

    # create list of sites
    ss = Site{LS,D}[
        newSite(Site{LS,D}, st_point[i], LS(st_label[i]))
        for i in 1:length(st_label)
    ]

    # return nothing
    return ss
end

# convinience function for standard type
function loadSites(
        fn :: AbstractString,
        group :: AbstractString = "sites"
    )

    return loadSites(Site, fn, group)
end


text:
unsafe_mulby_1000(x) = (x<<10) - (x << 5) + (x << 3)

safe_mulby_1000(x::Int64) = x <= 9_223_372_036_854_775 ? mulby_100(x) :
                                                         ArgumentError("$(x) is too large")
safe_mulby_1000(x::Int32) = x <= 2_147_483 ? mulby_100(x) : ArgumentError("$(x) is too large")

# nanoseconds per microsecond
# microseconds per millisecond
# milliseconds per second
# safe for 0 < x < 268435456  (2^28)
unsafe_fld_1000(x::T) where {T} =
    ((x >> 3) * 34_359_739) >> 32

fld_1000(x::T) where {T<:Union{Int64,UInt64}} =
    if (x < 434_934_000) # 2^28 = 268_435_456
    # if x & ~0xffffffff == 0
        unsafe_fld_1000(x)
    else
        fld(x, 1_000)
    end

function fldmod_1000(x::T) where {T<:Union{Int64,UInt64}}
    if abs(x) <= 9_780_955_816
        quotient = if (x < 434_934_000)          # in 2^28 .. 2^29
            unsafe_fld_1000(x)
        else
            fld(x, 1_000)
        end

        remainder = x - unsafe_mulby_1000(quotient)

        if signbit(quotient) && signbit(remainder)
            quotient -= 1
            remainder += 1000
        end
    else
        quotient, remainder = fldmod(x, 1000)
    end

    quotient, remainder
end

# nanoseconds per millisecond
# microseconds per second

function mulby1000000(x)
    y = Int128(x)
    z = (y << 6) - y
    z = (z << 5) - z
    z = (z << 3) + y
    z << 6
end

# safe for x:2^2:Int64 <= 2^27 - 1
unsafe_fld_1_000_000(x::T) where {T<:Union{Int64,UInt64}} =
    ((x >> 6) * 274_878) >> 32

fld_1_000_000(x::T) where {T<:Union{Int64,UInt64}} =
    if (x < 134_217_728)          # 2^27 -1
        unsafe_fld_1_000_000(x)
    else
        fld(x, 1_000_000)
    end

function fldmod_1_000_000(x::T) where {T<:Union{Int64,UInt64}}
    quotient = if (x < 134_217_728)          # 2^27 -1
        unsafe_fld_1_000_000(x)
    else
        fld(x, 1_000_000)
    end

    remainder = x - quotient * 1_000_000
    quotient, remainder
end

# 86_400 (SecondsPerDay)
function mulby86400(x)
    y = x
    z = (y << 4) - y
    z = (z << 4) - z
    z = (z << 2) - z
    z << 7
end

# 86_400_000 (MillisecondsPerDay)
function mulby86400000(x)
    y = x
    z = (y << 5) - y
    z = (z << 2) + y
    z = (z << 4) - z
    z = (z << 4) - z
    z = (z << 2) - z
    z << 10
end

# 60 (SecondsPerMinute, MinutesPerHour)
function mulby60(x)
    y = x
    z = (y << 4) - y
    z << 2
end

# 3600 = 60*60 (SecondsPerHour)
function mulby3600(x)
    y = x
    z = (y << 4) - y
    z = (z << 4) - z
    z << 4
end

# 1440 = 24*60 (MinutesPerDay)
function mulby1440(x)
    y = x
    z = (y << 4) - y
    z = (z << 2) - z
    z << 5
end

# 1_000_000_000 (NanosecondsPerSecond)
function mulby_1000000000(x)
    y = x
    z = (y << 6) - y
    z = (z << 5) - z
    z = (z << 3) + y
    z = (z << 2) + z
    z = (z << 2) + z
    z = (z << 2) + z
    z << 9
end

#=

reference

Euclidean Affine Functions and Applications to Calendar Algorithms
by Cassio Neri, Lorenz Schneider
    arXiv:2102.06959v1 [cs.DS] 13 Feb2021
    https://arxiv.org/pdf/2102.06959.pdf
    pg 11
=#

mod2p32(x) = (x & 0xffff)

unsafe_fld_3600(x) =
    (1193047 * x) >> 32

function fld_3600(x)
    if x < 2_257_200
        unsafe_fld_3600(x)
    else
        fld(x, 3600)
    end
end

unsafe_mod_3600(x) =
    3600 * mod2p32(1193047 * x) >> 32

function mod_3600(x)
    if x < 2_255_762
        unsafe_mod_3600(x)
    else
        mod(x, 3600)
    end
end

function fldmod_3600(x)
    if x < 2_255_762
        (unsafe_fld_3600(x), unsafe_mod_3600(x))
    else
        fldmod(x, 3600)
    end
end

unsafe_fld_60(x) =
    mod2p32(71_582_789 * x) >> 32

function fld_60(x)
    if x < 97_612_920
        unsafe_fld_60(x)
    else
        fld(x, 60)
    end
end

unsafe_mod_60(x) =
    60 * mod2p32(71_582_789 * x) >> 32

function mod_60(x)
    if x < 97_612_920
        unsafe_mod_60(x)
    else
        mod(x, 60)
    end
end

function fldmod_60(x)
    if x < 2_255_762
        (unsafe_fld_60(x), unsafe_mod_60(x))
    else
        fldmod(x, 60)
    end
end

unsafe_fld_10(x) =
    mod2p32(429_496_730 * x) >> 32

function fld_10(x)
    if x < 1_073_741_829
        unsafe_fld_10(x)
    else
        fld(x, 10)
    end
end

unsafe_mod_10(x) =
    10 * mod2p32(429_496_730 * x) >> 32

function mod_10(x)
    if x < 1_073_741_829
        unsafe_mod_10(x)
    else
        mod(x, 10)
    end
end

function fldmod_10(x)
    if x < 1_073_741_829
        (unsafe_fld_10(x), unsafe_mod_10(x))
    else
        fldmod(x, 10)
    end
end

#=

const twopow32 = typemax(UInt32)
mod2p32(x) = (x & 0xffff)

const k3600 =   1_193_047     #  n_div_3600 = mod2p32(k3600 * mod32(n))
const   k60 =  71_582_789     #    n_div_60 = mod2p32(  k60 * mod32(n))
const   k10 = 429_496_730     #    n_div_10 = mod2p32(  k10 * mod32(n))

const k3600 =   1_193_047     #  n_mod_3600 = (k3600 * mod32(n)) ÷ k3600
const   k60 =  71_582_789     #    n_mod_60 = (  k60 * mod32(n)) ÷ k60
const   k10 = 429_496_730     #    n_mod_10 = (  k10 * mod32(n)) ÷ k10

const k3600 =   1_193_047     #  n_div_3600 = mod2p32(k3600 * mod32(n))
const   k60 =  71_582_789     #    n_div_60 = mod2p32(  k60 * mod32(n))
const   k10 = 429_496_730     #    n_div_10 = mod2p32(  k10 * mod32(n))

const k3600 =   1_193_047     #  n_mod_3600 = (k3600 * mod32(n)) ÷ k3600
const   k60 =  71_582_789     #    n_mod_60 = (  k60 * mod32(n)) ÷ k60
const   k10 = 429_496_730     #    n_mod_10 = (  k10 * mod32(n)) ÷ k10

const k3600 =   1_193_047     #  n_div_3600 = mod2p32(k3600 * mod32(n))
const   k60 =  71_582_789     #    n_div_60 = mod2p32(  k60 * mod32(n))
const   k10 = 429_496_730     #    n_div_10 = mod2p32(  k10 * mod32(n))

const k3600 =   1_193_047     #  n_mod_3600 = (k3600 * mod32(n)) ÷ k3600
const   k60 =  71_582_789     #    n_mod_60 = (  k60 * mod32(n)) ÷ k60
const   k10 = 429_496_730     #    n_mod_10 = (  k10 * mod32(n)) ÷ k10

=#

#=
    divmod10(k)   -> (k ÷   10, k %   10)     ∀n ∈ [0,     2_257_200]
    divmod60(k)   -> (k ÷   60, k %   60)     ∀n ∈ [0,    97_612_920]
    divmod3600(k) -> (k ÷ 3600, k % 3600)     ∀n ∈ [0, 1_073_741_830]

    div10(k)    -> k ÷   10   ↤  ( 429_496_730 * n ) & twopow32
    div60(k)    -> k ÷   60   ↤  (  71_582_789 * n ) & twopow32
    div3600(k)  -> k ÷ 3600   ↤  (   1_193_047 * n ) & twopow32

    mod10(k)    -> k %   10   ↤  ( 429_496_730 * n ) % twopow32
    mod60(k)    -> k %   60
    mod3600(k)  -> k % 3600
=#


text:
# Autogenerated wrapper script for CompilerSupportLibraries_jll for i686-w64-mingw32-libgfortran4
export libatomic, libgcc_s, libgfortran, libgomp, libstdcxx

JLLWrappers.@generate_wrapper_header("CompilerSupportLibraries")
JLLWrappers.@declare_library_product(libatomic, "libatomic-1.dll")
JLLWrappers.@declare_library_product(libgcc_s, "libgcc_s_sjlj-1.dll")
JLLWrappers.@declare_library_product(libgfortran, "libgfortran-4.dll")
JLLWrappers.@declare_library_product(libgomp, "libgomp-1.dll")
JLLWrappers.@declare_library_product(libstdcxx, "libstdc++-6.dll")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libatomic,
        "bin\\libatomic-1.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgcc_s,
        "bin\\libgcc_s_sjlj-1.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgfortran,
        "bin\\libgfortran-4.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgomp,
        "bin\\libgomp-1.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libstdcxx,
        "bin\\libstdc++-6.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# This file contains code that was formerly part of Julia. License is MIT: https://julialang.org/license

using AbstractFFTs
using AbstractFFTs: Plan
using LinearAlgebra
using Test

mutable struct TestPlan{T} <: Plan{T}
    region
    pinv::Plan{T}
    TestPlan{T}(region) where {T} = new{T}(region)
end

mutable struct InverseTestPlan{T} <: Plan{T}
    region
    pinv::Plan{T}
    InverseTestPlan{T}(region) where {T} = new{T}(region)
end

AbstractFFTs.plan_fft(x::Vector{T}, region; kwargs...) where {T} = TestPlan{T}(region)
AbstractFFTs.plan_bfft(x::Vector{T}, region; kwargs...) where {T} = InverseTestPlan{T}(region)
AbstractFFTs.plan_inv(p::TestPlan{T}) where {T} = InverseTestPlan{T}

# Just a helper function since forward and backward are nearly identical
function dft!(y::Vector, x::Vector, sign::Int)
    n = length(x)
    length(y) == n || throw(DimensionMismatch())
    fill!(y, zero(complex(float(eltype(x)))))
    c = sign * 2π / n
    @inbounds for j = 0:n-1, k = 0:n-1
        y[k+1] += x[j+1] * cis(c*j*k)
    end
    return y
end

mul!(y::Vector, p::TestPlan, x::Vector) = dft!(y, x, -1)
mul!(y::Vector, p::InverseTestPlan, x::Vector) = dft!(y, x, 1)

Base.:*(p::TestPlan, x::Vector) = mul!(copy(x), p, x)
Base.:*(p::InverseTestPlan, x::Vector) = mul!(copy(x), p, x)

@testset "Custom Plan" begin
    x = AbstractFFTs.fft(collect(1:8))
    # Result computed using FFTW
    fftw_fft = [36.0 + 0.0im,
                -4.0 + 9.65685424949238im,
                -4.0 + 4.0im,
                -4.0 + 1.6568542494923806im,
                -4.0 + 0.0im,
                -4.0 - 1.6568542494923806im,
                -4.0 - 4.0im,
                -4.0 - 9.65685424949238im]
    @test x ≈ fftw_fft

    fftw_bfft = [Complex{Float64}(8i, 0) for i in 1:8]
    @test AbstractFFTs.bfft(x) ≈ fftw_bfft

    fftw_ifft = [Complex{Float64}(i, 0) for i in 1:8]
    @test AbstractFFTs.ifft(x) ≈ fftw_ifft

    @test eltype(plan_fft(collect(1:8))) == Int
end

@testset "Shift functions" begin
    @test AbstractFFTs.fftshift([1 2 3]) == [3 1 2]
    @test AbstractFFTs.fftshift([1, 2, 3]) == [3, 1, 2]
    @test AbstractFFTs.fftshift([1 2 3; 4 5 6]) == [6 4 5; 3 1 2]

    @test AbstractFFTs.fftshift([1 2 3; 4 5 6], 1) == [4 5 6; 1 2 3]
    @test AbstractFFTs.fftshift([1 2 3; 4 5 6], ()) == [1 2 3; 4 5 6]
    @test AbstractFFTs.fftshift([1 2 3; 4 5 6], (1,2)) == [6 4 5; 3 1 2]
    @test AbstractFFTs.fftshift([1 2 3; 4 5 6], 1:2) == [6 4 5; 3 1 2]

    @test AbstractFFTs.ifftshift([1 2 3]) == [2 3 1]
    @test AbstractFFTs.ifftshift([1, 2, 3]) == [2, 3, 1]
    @test AbstractFFTs.ifftshift([1 2 3; 4 5 6]) == [5 6 4; 2 3 1]

    @test AbstractFFTs.ifftshift([1 2 3; 4 5 6], 1) == [4 5 6; 1 2 3]
    @test AbstractFFTs.ifftshift([1 2 3; 4 5 6], ()) == [1 2 3; 4 5 6]
    @test AbstractFFTs.ifftshift([1 2 3; 4 5 6], (1,2)) == [5 6 4; 2 3 1]
    @test AbstractFFTs.ifftshift([1 2 3; 4 5 6], 1:2) == [5 6 4; 2 3 1]
end

@testset "FFT Frequencies" begin
    @test fftfreq(8) isa Frequencies
    @test copy(fftfreq(8)) isa Frequencies

    # N even
    @test fftfreq(8) == [0.0, 0.125, 0.25, 0.375, -0.5, -0.375, -0.25, -0.125]
    @test rfftfreq(8) == [0.0, 0.125, 0.25, 0.375, 0.5]
    @test fftshift(fftfreq(8)) == -0.5:0.125:0.375

    # N odd
    @test fftfreq(5) == [0.0, 0.2, 0.4, -0.4, -0.2]
    @test rfftfreq(5) == [0.0, 0.2, 0.4]
    @test fftshift(fftfreq(5)) == -0.4:0.2:0.4

    # Sampling Frequency
    @test fftfreq(5, 2) == [0.0, 0.4, 0.8, -0.8, -0.4]
    # <:Number type compatibility
    @test eltype(fftfreq(5, ComplexF64(2))) == ComplexF64

    @test_throws ArgumentError Frequencies(12, 10, 1)
end

@testset "normalization" begin
    # normalization should be inferable even if region is only inferred as ::Any,
    # need to wrap in another function to test this (note that p.region::Any for
    # p::TestPlan)
    f9(p::Plan{T}, sz) where {T} = AbstractFFTs.normalization(real(T), sz, p.region)
    @test @inferred(f9(plan_fft(zeros(10), 1), 10)) == 1/10
end


text:
import Distributed

function add_processes(count)
    count -= Distributed.nprocs()
    if count > 0
        @info "Starting $count more processes ..."
        Distributed.addprocs(count)
    end
end

function verify(::Type{SpawnProcesses}, sig::Model.Signature, pkhashpairs)
    # kick off workers in parallel
    future_rhside = Distributed.@spawn Curve.curve_miller(Config.G2, sig.sig)
    lhside = Distributed.@distributed (*) for pair in pkhashpairs
        Curve.curve_miller(pair.first.pk, pair.second.hash)
    end
    res = Curve.field_final_exp(fetch(lhside) // fetch(future_rhside))
    return res == one(typeof(res))
end


text:
module TruncatedMVN

using Distributions
using Compat
using Convex
using SCS

import Base.rand

export TruncatedNormalSampler,
       GibbsTMVN

include("univariate.jl")
include("multivariate.jl")

end # module


text:
module SimpleTrees

using ..Traits: getcontent, parametertype, efficientoperations, rawtype

import ..Traits: contentnames, dissolve

export simpletreedepth, simpletreewidth
export SimpleTreeCore, AbstractSimpleTree, SimpleTree
export root, children, ancestor, descendants, siblings, leaves, subtree
export isleaf, level, addnode!, deletenode!, move!

"""
    SimpleTreeCore()

The core of a simple tree.
"""
mutable struct SimpleTreeCore{N, D}
    root::Union{N, Nothing}
    contents::Dict{N, D}
    parent::Dict{N, N}
    children::Dict{N, Vector{N}}
    SimpleTreeCore{N, D}() where {N, D} = new{N, D}(nothing, Dict{N, D}(), Dict{N, N}(), Dict{N, Vector{N}}())
end
@inline Base.:(==)(tc1::TC, tc2::TC) where {TC<:SimpleTreeCore} = ==(efficientoperations, tc1, tc2)
@inline Base.isequal(tc1::TC, tc2::TC) where {TC<:SimpleTreeCore} = isequal(efficientoperations, tc1, tc2)

abstract type SimpleTreeIteration end
struct SimpleTreeDepth <: SimpleTreeIteration end
struct SimpleTreeWidth <: SimpleTreeIteration end

"""
    simpletreedepth

Indicate that the iteration over a tree is depth-first.
"""
const simpletreedepth = SimpleTreeDepth()

"""
    simpletreewidth

Indicate that the iteration over a tree is width-first.
"""
const simpletreewidth = SimpleTreeWidth()

"""
    AbstractSimpleTree{Node, Data}

Abstract type for all concrete trees.
"""
abstract type AbstractSimpleTree{N, D} end
@inline contentnames(::Type{<:AbstractSimpleTree}) = (:TREECORE,)
@inline Base.:(==)(t₁::T, t₂::T) where {T<:AbstractSimpleTree} = ==(efficientoperations, t₁, t₂)
@inline Base.isequal(t₁::T, t₂::T) where {T<:AbstractSimpleTree} = isequal(efficientoperations, t₁, t₂)

"""
    SimpleTreeCore(tree::AbstractSimpleTree) -> SimpleTreeCore

Get the core of a simple tree.
"""
@inline SimpleTreeCore(tree::AbstractSimpleTree) = getcontent(tree, :TREECORE)

"""
    keytype(tree::AbstractSimpleTree)
    keytype(::Type{T}) where {T<:AbstractSimpleTree}

Get a tree's node type.
"""
@inline Base.keytype(tree::AbstractSimpleTree) = keytype(typeof(tree))
@inline @generated Base.keytype(::Type{T}) where {T<:AbstractSimpleTree} = parametertype(supertype(T, :AbstractSimpleTree), 1)

"""
    valtype(tree::AbstractSimpleTree)
    valtype(::Type{T}) where {T<:AbstractSimpleTree}

Get a tree's data type.
"""
@inline Base.valtype(tree::AbstractSimpleTree) = valtype(typeof(tree))
@inline @generated Base.valtype(::Type{T}) where {T<:AbstractSimpleTree} = parametertype(supertype(T, :AbstractSimpleTree), 2)

"""
    eltype(tree::AbstractSimpleTree)
    eltype(::Type{T}) where {T<:AbstractSimpleTree}

Get the eltype of a tree.
"""
@inline Base.eltype(tree::AbstractSimpleTree) = eltype(typeof(tree))
@inline Base.eltype(::Type{T}) where {T<:AbstractSimpleTree} = Pair{keytype(T), valtype(T)}

"""
    root(tree::AbstractSimpleTree) -> Union{keytype(tree), Nothing}

Get a tree's root node.
"""
@inline root(tree::AbstractSimpleTree) = SimpleTreeCore(tree).root

"""
    haskey(tree::AbstractSimpleTree{N}, node::N) where N -> Bool

Check whether a node is in a tree.
"""
@inline Base.haskey(tree::AbstractSimpleTree{N}, node::N) where N = haskey(SimpleTreeCore(tree).contents, node)

"""
    length(tree::AbstractSimpleTree) -> Int

Get the number of a tree's nodes.
"""
@inline Base.length(tree::AbstractSimpleTree) = length(SimpleTreeCore(tree).contents)

"""
    parent(tree::AbstractSimpleTree{N}, node::N, superparent::Union{N, Nothing}=nothing) where N -> Union{N, Nothing}

Get the parent of a tree's node. When `node` is the tree's root, return `superparent`.
"""
@inline Base.parent(tree::AbstractSimpleTree{N}, node::N, superparent::Union{N, Nothing}=nothing) where {N} = (node == root(tree) ? superparent : SimpleTreeCore(tree).parent[node])

"""
    children(tree::AbstractSimpleTree) -> Vector{keytype(tree)}
    children(tree::AbstractSimpleTree, ::Nothing) -> Vector{keytype(tree)}
    children(tree::AbstractSimpleTree{N}, node::N) where N -> Vector{N}

Get the children of a tree's node.
"""
@inline children(tree::AbstractSimpleTree) = children(tree, nothing)
@inline children(tree::AbstractSimpleTree, ::Nothing) = (root(tree) === nothing) ? error("children error: empty tree!") : [root(tree)]
@inline children(tree::AbstractSimpleTree{N}, node::N) where N = SimpleTreeCore(tree).children[node]

"""
    addnode!(tree::AbstractSimpleTree{N}, node::N) where N} -> typeof(tree)
    addnode!(tree::AbstractSimpleTree{N}, ::Nothing, node::N) where N -> typeof(tree)
    addnode!(tree::AbstractSimpleTree{N}, parent::N, node::N) where N -> typeof(tree)

Update the structure of a tree by adding a node. When the parent is `nothing`, the input tree must be empty and the input node becomes the tree's root.
"""
@inline addnode!(tree::AbstractSimpleTree{N}, node::N) where {N} = addnode!(tree, nothing, node)
function addnode!(tree::AbstractSimpleTree{N}, ::Nothing, node::N) where N
    @assert root(tree) === nothing "addnode! error: not empty tree."
    SimpleTreeCore(tree).root = node
    SimpleTreeCore(tree).children[node] = N[]
    return tree
end
function addnode!(tree::AbstractSimpleTree{N}, parent::N, node::N) where N
    @assert haskey(tree, parent) "addnode! error: parent($parent) not in tree."
    @assert !haskey(tree, node) "addnode! error: node($node) already in tree."
    push!(SimpleTreeCore(tree).children[parent], node)
    SimpleTreeCore(tree).parent[node] = parent
    SimpleTreeCore(tree).children[node] = N[]
    return tree
end

"""
    deletenode!(tree::AbstractSimpleTree{N}, node::N) where N -> typeof(tree)

Update the structure of a tree by deleting a node.
"""
function deletenode!(tree::AbstractSimpleTree{N}, node::N) where N
    @assert haskey(tree, node) "deletenode! error: node($node) not in tree."
    if node == root(tree)
        SimpleTreeCore(tree).root = nothing
    else
        pnode = parent(tree, node)
        haskey(SimpleTreeCore(tree).children, pnode) && filter!(!=(node), SimpleTreeCore(tree).children[pnode])
    end
    pop!(SimpleTreeCore(tree).contents, node)
    pop!(SimpleTreeCore(tree).parent, node, nothing)
    pop!(SimpleTreeCore(tree).children, node)
    return tree
end

"""
    getindex(tree::AbstractSimpleTree{N}, node::N) where N -> N

Get the data of a tree's node.
"""
@inline Base.getindex(tree::AbstractSimpleTree{N}, node::N) where {N} = SimpleTreeCore(tree).contents[node]

"""
    setindex!(tree::AbstractSimpleTree{N, D}, data::D, node::N) where {N, D}

Set the data of a tree's node.
"""
@inline Base.setindex!(tree::AbstractSimpleTree{N, D}, data::D, node::N) where {N, D} = (SimpleTreeCore(tree).contents[node] = data)

"""
    empty(tree::AbstractSimpleTree)

Construct an empty tree of the same type with the input one.
"""
@inline Base.empty(tree::AbstractSimpleTree) = rawtype(typeof(tree))(dissolve(tree, empty)...)
@inline dissolve(tree::AbstractSimpleTree{N, D}, ::Val{:TREECORE}, ::typeof(empty), args::Tuple, kwargs::NamedTuple) where {N, D} = SimpleTreeCore{N, D}()

"""
    keys(tree::AbstractSimpleTree{N}, ::SimpleTreeDepth, node::Union{N, Nothing}=root(tree)) where N
    keys(tree::AbstractSimpleTree{N}, ::SimpleTreeWidth, node::Union{N, Nothing}=root(tree)) where N

Iterate over a tree's nodes starting from a certain `node` by depth first search or width first search.
"""
function Base.keys(tree::AbstractSimpleTree{N}, ti::SimpleTreeIteration, node::Union{N, Nothing}=root(tree)) where N
    TreeKeys{typeof(ti), N, valtype(tree), typeof(tree)}(tree, node)
end
struct TreeKeys{P<:SimpleTreeIteration, N, D, T<:AbstractSimpleTree{N, D}}
    tree::T
    node::Union{N, Nothing}
end
Base.eltype(::Type{<:TreeKeys{<:SimpleTreeIteration, N, D, <:AbstractSimpleTree{N, D}} where D}) where {N} = N
Base.IteratorSize(::Type{<:TreeKeys}) = Base.SizeUnknown()
Base.iterate(tk::TreeKeys) = (root(tk.tree) === nothing) ? nothing : (tk.node, copy(children(tk.tree, tk.node)))
function Base.iterate(tk::TreeKeys{SimpleTreeDepth, N}, state::Vector{N}) where N
    (length(state) == 0) ? nothing : (node = popfirst!(state); prepend!(state, children(tk.tree, node)); (node, state))
end
function Base.iterate(tk::TreeKeys{SimpleTreeWidth, N}, state::Vector{N}) where N
    (length(state) == 0) ? nothing : (node = popfirst!(state); append!(state, children(tk.tree, node)); (node, state))
end

"""
    values(tree::AbstractSimpleTree{N}, ::SimpleTreeDepth, node::Union{N, Nothing}=root(tree)) where N
    values(tree::AbstractSimpleTree{N}, ::SimpleTreeWidth, node::Union{N, Nothing}=root(tree)) where N

Iterate over a tree's data starting from a certain `node` by depth first search or width first search.
"""
@inline Base.values(tree::AbstractSimpleTree{N}, ti::SimpleTreeIteration, node::Union{N, Nothing}=root(tree)) where {N} = (tree[key] for key in keys(tree, ti, node))

"""
    pairs(tree::AbstractSimpleTree{N}, ::SimpleTreeDepth, node::Union{N, Nothing}=root(tree)) where N
    pairs(tree::AbstractSimpleTree{N}, ::SimpleTreeWidth, node::Union{N, Nothing}=root(tree)) where N

Iterate over a tree's (node, data) pairs starting from a certain `node` by depth first search or width first search.
"""
@inline Base.pairs(tree::AbstractSimpleTree{N}, ti::SimpleTreeIteration, node::Union{N, Nothing}=root(tree)) where {N} = ((key, tree[key]) for key in keys(tree, ti, node))

"""
    isleaf(tree::AbstractSimpleTree{N}, node::N) where N -> Bool

Judge whether a tree's node is a leaf (a node without children) or not.
"""
@inline isleaf(tree::AbstractSimpleTree{N}, node::N) where {N} = length(children(tree, node)) == 0

"""
    level(tree::AbstractSimpleTree{N}, node::N) where N -> Int

Get the level of tree's node.
"""
@inline function level(tree::AbstractSimpleTree{N}, node::N) where N
    result = 1
    while node != root(tree)
        result += 1
        node = parent(tree, node)
    end
    return result
end

"""
    ancestor(tree::AbstractSimpleTree{N}, node::N, generation::Int=1) where N -> N

Get the ancestor of a tree's node of the n-th generation.
"""
@inline function ancestor(tree::AbstractSimpleTree{N}, node::N, generation::Int=1) where N
    @assert generation >= 0 "ancestor error: generation($generation) must be non-negative."
    result = node
    for i = 1:generation
        result = parent(tree, result)
    end
    result
end

"""
    descendants(tree::AbstractSimpleTree{N}, node::N, generation::Int=1) where N -> Vector{N}

Get the descendants of a tree's node of the nth generation.
"""
function descendants(tree::AbstractSimpleTree{N}, node::N, generation::Int=1) where N
    @assert generation >= 0 "descendants error: generation($generation) must be non-negative."
    result = N[node]
    for i = 1:generation
        result = vcat((children(tree, node) for node in result)...)
    end
    result
end

"""
    siblings(tree::AbstractSimpleTree{N}, node::N) where N -> Vector{N}

Get the siblings (other nodes sharing the same parent) of a tree's node.
"""
@inline siblings(tree::AbstractSimpleTree{N}, node::N) where N = filter(!=(node), children(tree, parent(tree, node)))

"""
    leaves(tree::AbstractSimpleTree) -> Vector{keytype(tree)}

Get a tree's leaves.
"""
@inline leaves(tree::AbstractSimpleTree) = keytype(tree)[node for node in keys(tree, simpletreedepth) if isleaf(tree, node)]

"""
    push!(tree::AbstractSimpleTree{N, D}, node::N, data::D) where {N, D} -> typeof(tree)
    push!(tree::AbstractSimpleTree{N, D}, parent::Union{N, Nothing}, node::N, data::D) where {N, D} -> typeof(tree)

Push a new node to a tree. When `parent` is `nothing`, this function set the root node of an empty tree.
"""
Base.push!(tree::AbstractSimpleTree{N, D}, node::N, data::D) where {N, D} = push!(tree, nothing, node, data)
function Base.push!(tree::AbstractSimpleTree{N, D}, parent::Union{N, Nothing}, node::N, data::D) where {N, D}
    addnode!(tree, parent, node)
    tree[node] = data
    return tree
end

"""
    append!(tree::AbstractSimpleTree{N, D}, subtree::AbstractSimpleTree{N, D}) where {N, D} -> typeof(tree)
    append!(tree::AbstractSimpleTree{N, D}, node::Union{N, Nothing}, subtree::AbstractSimpleTree{N, D}) where {N, D} -> typeof(tree)

Append a subtree to a tree.
"""
Base.append!(tree::AbstractSimpleTree{N, D}, subtree::AbstractSimpleTree{N, D}) where {N, D} = append!(tree, nothing, subtree)
function Base.append!(tree::AbstractSimpleTree{N, D}, node::Union{N, Nothing}, subtree::AbstractSimpleTree{N, D}) where {N, D}
    for (key, value) in pairs(subtree, simpletreewidth)
        push!(tree, parent(subtree, key, node), key, value)
    end
    return tree
end

"""
    delete!(tree::AbstractSimpleTree{N}, node::N) where N -> typeof(tree)

Delete a node and all its descendants from a tree.
"""
function Base.delete!(tree::AbstractSimpleTree{N}, node::N) where N
    for key in collect(N, keys(tree, simpletreedepth, node))
        deletenode!(tree, key)
    end
    return tree
end

"""
    empty!(tree::AbstractSimpleTree) -> typeof(tree)

Empty a tree.
"""
@inline Base.empty!(tree::AbstractSimpleTree) = delete!(tree, root(tree))

"""
    subtree(tree::AbstractSimpleTree{N}, node::N) where N -> typeof(tree)

Get a subtree whose root is `node`.
"""
function subtree(tree::AbstractSimpleTree{N}, node::N) where N
    result = empty(tree)
    for (i, (key, value)) in enumerate(pairs(tree, simpletreedepth, node))
        push!(result, (i == 1) ? nothing : parent(tree, key), key, value)
    end
    return result
end

"""
    move!(tree::AbstractSimpleTree{N}, node::N, parent::N) where N -> typeof(tree)

Move a subtree to a new position.
"""
function move!(tree::AbstractSimpleTree{N}, node::N, parent::N) where N
    sub = subtree(tree, node)
    delete!(tree, node)
    append!(tree, parent, sub)
    return tree
end

"""
    SimpleTree{N, D} <: AbstractSimpleTree{N, D}

The minimum tree structure that implements all the default tree methods.
"""
struct SimpleTree{N, D} <: AbstractSimpleTree{N, D}
    TREECORE::SimpleTreeCore{N, D}
end

"""
    SimpleTree{N, D}() where {N, D}

Construct an empty simple tree.
"""
SimpleTree{N, D}() where {N, D} = SimpleTree(SimpleTreeCore{N, D}())

end #module


text:
using Test
using MPI
using CLIMA
using CLIMA.Mesh.Topologies
using CLIMA.Mesh.Grids
using CLIMA.Mesh.Grids: VerticalDirection, HorizontalDirection, EveryDirection
using CLIMA.VTK
using Logging
using Printf
using LinearAlgebra

const ArrayType = CLIMA.array_type()

let
  # boiler plate MPI stuff
  CLIMA.init()
  mpicomm = MPI.COMM_WORLD
  ll = uppercase(get(ENV, "JULIA_LOG_LEVEL", "INFO"))
  loglevel = ll == "DEBUG" ? Logging.Debug :
    ll == "WARN"  ? Logging.Warn  :
    ll == "ERROR" ? Logging.Error : Logging.Info
  logger_stream = MPI.Comm_rank(mpicomm) == 0 ? stderr : devnull
  global_logger(ConsoleLogger(logger_stream, loglevel))

  # Mesh generation parameters
  N = 4
  Nq = N+1
  Neh = 10
  Nev = 4

  @testset "$(@__FILE__) DGModel matrix" begin
    for FT in (Float64, Float32)
      for dim = (2, 3)
        if dim == 2
          brickrange = (range(FT(0); length=Neh+1, stop=1),
                        range(FT(1); length=Nev+1, stop=2))
        elseif dim == 3
          brickrange = (range(FT(0); length=Neh+1, stop=1),
                        range(FT(0); length=Neh+1, stop=1),
                        range(FT(1); length=Nev+1, stop=2))
        end

        topl = StackedBrickTopology(mpicomm, brickrange)

        function warpfun(ξ1, ξ2, ξ3)
          FT = eltype(ξ1)

          ξ1 ≥ FT(1//2) && (ξ1 = FT(1//2) + 2(ξ1 - FT(1//2)))
          if dim == 2
            ξ2 ≥ FT(3//2) && (ξ2 = FT(3//2) + 2(ξ2 - FT(3//2)))
          elseif dim == 3
            ξ2 ≥ FT(1//2) && (ξ2 = FT(1//2) + 2(ξ2 - FT(1//2)))
            ξ3 ≥ FT(3//2) && (ξ3 = FT(3//2) + 2(ξ3 - FT(3//2)))
          end
          (ξ1, ξ2, ξ3)
        end

        grid = DiscontinuousSpectralElementGrid(topl,
                                                FloatType = FT,
                                                DeviceArray = ArrayType,
                                                polynomialorder = N,
                                                meshwarp = warpfun)

        # testname = "grid_poly$(N)_dim$(dim)_$(ArrayType)_$(FT)"
        # filename(rank) = @sprintf("%s_mpirank%04d", testname, rank)
        # writevtk(filename(MPI.Comm_rank(mpicomm)), grid)
        # if MPI.Comm_rank(mpicomm) == 0
        #   writepvtu(testname, filename.(0:MPI.Comm_size(mpicomm)-1), ())
        # end

        ξ = referencepoints(grid)
        hmnd = (ξ[2] - ξ[1]) / (2Neh)
        vmnd = (ξ[2] - ξ[1]) / (2Nev)

        @test hmnd ≈ min_node_distance(grid, EveryDirection())
        @test vmnd ≈ min_node_distance(grid, VerticalDirection())
        @test hmnd ≈ min_node_distance(grid, HorizontalDirection())

      end
    end
end
end

nothing

max_stars_repo_name: emassoud/CLIMA

max_issues_repo_name: emassoud/CLIMA

max_forks_repo_name: emassoud/CLIMA


text:
include("lasso.jl")
include("utils.jl")

using Test
using UnicodePlots
using Plots

Test1Results = 0

function Test1()
    """
        Test 1, some basic test on a random matrix, copied columns of matrix.
    """
    deg = 3
    N = 2000
    B = rand(N, deg)
    A = fill(0.0, (N, 2deg))
    A[:, 1:deg] = B
    A[:, deg + 1:end] = 2B.^2 + B + 0.1.*randn(size(B))

    b = A*ones(2deg, 1)
    instance = LassoSCOP(A, b)
    results, λs = LassoPath(instance)
    display(results)
    display(λs)

    # Let's plot this
    global Test1Results = results

    ylim = (minimum(Test1Results), maximum(Test1Results))
    ThePlot = lineplot(log2.(λs), Test1Results[1, :], ylim=ylim, title="Lasso Path")

    for II ∈ 2:size(Test1Results, 1)
        lineplot!(ThePlot, log2.(λs), Test1Results[II, :])
    end
    display(ThePlot)
    ThePlot = UnicodePlots.heatmap(results[end:-1:begin, :])
    display(ThePlot)

    # Visualize the path.
    VisualizeLassoPath(instance)

    # Finding the important weights for it.
    PrintTitle("These are the important weights: ")
    WeightsIndices = CaptureImportantWeights(instance)
    println(WeightsIndices)
    return true

end

@test Test1()

function Test2()
    N = 10
    deg = 10
    xgrid = collect(LinRange(-2π, 2π, N))
    f(x) = sin.(x)
    y = f(xgrid) + 0.1.*rand(N)
    V = VanderMonde(reshape(xgrid, (N, 1)), deg)
    instance = LassoSCOP(V, reshape(y, (N, 1)))
    results, λs = LassoPath(instance)
    λs = log2.(λs)
    display(results)

    ylim = (minimum(results), maximum(results))
    ThePlot = lineplot(λs, results[1, :], ylim=ylim)
    for II ∈ 2:size(results, 1)
        lineplot!(ThePlot, λs ,results[II, :])
    end
    display(ThePlot)

    YHat = V*results[:, end]
    ThePlot = lineplot(xgrid, y)
    scatterplot!(ThePlot, xgrid, YHat)
    display(ThePlot)
    display(UnicodePlots.heatmap(results[end:-1:begin, :]))

    return true
end

@test Test2()


text:
import Base: /
using ..YAXTools: PickAxisArray
struct YAXSlice{T,N,P}
    c::P
    sliceaxes::Any
    otheraxes::Any
end
function YAXSlice(p, dims)
    N = ndims(p)
    ax = caxes(p)
    isliceax = (findAxis.(dims, Ref(ax))...,)
    any(isnothing, isliceax) && error("Axis not in cube")
    iotherax = (filter(i -> !in(i, isliceax), 1:N)...,)
    YAXSlice{Array{eltype(p),length(dims)},N - length(dims),typeof(p)}(
        p,
        (isliceax, getindex.(Ref(ax), isliceax)),
        (iotherax, getindex.(Ref(ax), iotherax)),
    )
end
/(c::YAXArray, s::Union{String,Symbol}) = YAXSlice(c, (s,))
/(c::YAXArray, s) = YAXSlice(c, s)

dimvals(x::YAXSlice, i) = x.c.axes[x.otheraxes[1][i]].values

function dimname(x::YAXSlice, i)
    axsym(x.otheraxes[2][i])
end

getattributes(x::YAXSlice) = x.c.properties

iscontdim(x::YAXSlice, i) = isa(x.c.axes[x.otheraxes[1][i]], RangeAxis)

function getdata(x::YAXSlice)
    m = fill!(Array{Union{Colon,Bool}}(undef, ndims(x.c)), true)
    m[collect(x.sliceaxes[1])] .= Colon()
    PickAxisArray(x.c.data, m)
end

getCubeDes(s::YAXSlice) =
    string(join(axname.(s.sliceaxes[2]), " x "), " slices over an ", getCubeDes(s.c))
cubesize(s::YAXSlice) = cubesize(s.c)
Base.ndims(::YAXSlice{<:Any,N}) where {N} = N

Base.show(io::IO, s::YAXSlice) = YAXArrays.Cubes.show_yax(io, s)

max_forks_repo_name: meggart/YAXArrays


text:
# PML coordinate streching functions
function s_PML(x,σs,Hd,k,LHp,LHn,dpml)
    σ = (x[2]<0)&&(x[2]+Hd>0) ? σs[2] : σs[1]
    xf = [x[1] x[2]]
    u = @. ifelse(xf>0 , xf-LHp , -xf-LHn)
    return @. ifelse(u > 0,  1-(1im*σ/k)*(u/dpml)^2, $(1.0+0im))
end

function Λ(x,σs,Hd,k,LHp,LHn,dpml)
    s_x,s_y = s_PML(x,σs,Hd,k,LHp,LHn,dpml)
    return TensorValue(1/s_x,0,0,1/s_y)
end

max_stars_repo_name: Tiamat-Tech/hpinn

max_issues_repo_name: louhz/hpinn

max_forks_repo_name: louhz/hpinn


text:
export bennett, bennett!

function direct_emulate(step, x0::T, args...; N::Int, kwargs...) where T
    xpre = copy(x0)
    local x
    for i=1:N
        x = _zero(xpre)
        res = step(x, xpre, args...; kwargs...)
        xpre = res[1]
        args = res[3:end]
    end
    return xpre
end

struct BennettLog
    fcalls::Vector{NTuple{3,Any}}  # depth, function index f_i := s_{i-1} -> s_{i}, length should be `(2k-1)^n` and function
    peak_mem::Base.RefValue{Int}  # should be `n*(k-1)+2`
    depth::Base.RefValue{Int}
end
BennettLog() = BennettLog(NTuple{3,Any}[], Ref(0), Ref(0))

# hacking the reversible program
function logfcall(l::BennettLog, i, f)
    push!(l.fcalls, (l.depth[], i, f))
    l, i, f
end
function ilogfcall(l::BennettLog, i, f)
    push!(l.fcalls, (l.depth[], i, ~f))
    l, i, f
end

@dual logfcall ilogfcall

Base.show(io::IO, ::MIME"text/plain", logger::BennettLog) = Base.show(io, logger)
function Base.show(io::IO, logger::BennettLog)
    nreverse = count(x->x[3] isa Inv, logger.fcalls)
    print(io, """Bennett log
| peak memory usage = $(logger.peak_mem[])
| number of function forward/backward calls = $(length(logger.fcalls)-nreverse)/$nreverse""")
end

"""
    bennett(step, y, x, args...; k, N, logger=BennettLog(), kwargs...)

* `step` is a reversible step function,
* `y` is the output state,
* `x` is the input state,
* `k` is the number of steps in each Bennett's recursion,
* `N` is the total number of steps,
* `logger=BennettLog()` is the logging of Bennett's algorithm,
* `args...` and `kwargs...` are additional arguments for steps.
"""
@i function bennett(step, y::T, x::T, args...; k::Int, N::Int, logger=BennettLog(), kwargs...) where T
    state ← Dict{Int, T}()
    state[1] ← _zero(x)
    state[1] +=  x
    bennett!((@skip! step), state, k, 1, N, args...; do_uncomputing=true, logger=logger, kwargs...)
    SWAP(y, state[N+1])
    state[1] -= x
    state[1] → _zero(x)
    state[N+1] → _zero(x)
    state → Dict{Int, T}()
end

"""
    bennett!(step, state::Dict, args...; k, N, logger=BennettLog(), do_uncomputing=false, kwargs...)

* `step` is a reversible step function,
* `state` is the dictionary state, with `state[1]` the input state, the return value is stored in `state[N+1]`,
* `k` is the number of steps in each Bennett's recursion,
* `N` is the total number of steps,
* `logger=BennettLog()` is the logging of Bennett's algorithm,
* `args...` and `kwargs...` are additional arguments for steps.
"""
@i function bennett!(step, state::Dict{Int,T}, args...; k::Int, N::Int, logger=BennettLog(), do_uncomputing=false, kwargs...) where T
    bennett!(step, state, k, 1, N, args...; logger=logger, do_uncomputing=do_uncomputing, kwargs...)
end

@i function bennett!(step, state::Dict{Int,T}, k::Int, base, len, args...; logger, do_uncomputing, kwargs...) where T
    @safe logger !== nothing && (logger.depth[] += 1)
    @invcheckoff if len == 1
        state[base+1] ← _zero(state[base])
        @safe logger !== nothing && (logger.peak_mem[] = max(logger.peak_mem[], length(state)))
        getf(step, base)(state[base+1], state[base], args...; kwargs...)
        if logger !== nothing
            logfcall(logger, (@const base+1), (@const getf(step, base)))
        end
    else
        @routine begin
            @zeros Int nstep n
            n += ceil((@skip! Int), (@const len / k))
            nstep += ceil((@skip! Int), (@const len / n))
        end
        for j=1:nstep
            bennett!(step, state, k, (@const base+n*(j-1)), (@const min(n,len-n*(j-1))), args...; logger=logger, do_uncomputing=true, kwargs...)
        end
        if do_uncomputing
            for j=nstep-1:-1:1
                ~bennett!(step, state, k, (@const base+n*(j-1)), n, args...; logger=logger, do_uncomputing=true, kwargs...)
            end
        end
        ~@routine
    end
end

getf(f, i::Int) = f
getf(f::AbstractArray, i::Int) = f[i]