text:
__precompile__(true)

module DEM

using LinearAlgebra, LoopVectorization, Tensors, CSV, DataFrames, StaticArrays

export dem, semiaxes, saveIO, cleantensor!

"""
Axis structure
"""
struct Axis{T}
    a1::T
    a2::T
    a3::T
end
"""
Interaction tensor structure
"""
struct Tcache
    n1::Int16
    n2::Int16
    θ::Vector{Float64}
    ϕ::Vector{Float64}
    T::Array{Float64,4}
    Tv::Array{Float64,4}
    Av::Array{Float64,2}
end

# -- version 1 : one axis
@inline  semiaxes(a1::Float64,a2::Float64,a3::Float64) = [Axis(a1,a2,a3)]

# -- version 1.2 : multiple axes
@inline function semiaxes(a10::Number,a1f::Number, Δ1::Number,
                  a20::Number,a2f::Number, Δ2::Number,
                  a30::Number,a3f::Number, Δ3::Number)

    a1 = range(a10, a1f, step=Δ1); n1 = length(a1)
    a2 = range(a20, a2f, step=Δ2); n2 = length(a2)
    a3 = range(a30, a3f, step=Δ3); n3 = length(a3)
    a1, a2, a3 = promote(a1, a2, a3)
    a  = Axis[]
    for i1 ∈ 1:n1, i2 ∈ 1:n2, i3 ∈ 1:n3
        push!(a,Axis(a1[i1],a2[i2],a3[i3]))
    end
    return a
end

# -- version 2 : multiple axes
@inline function semiaxes(a10,a1f,a20,a2f,a30,a3f,da)
    a1 = range(a10, a1f, step=da); n1 = length(a1)
    a2 = range(a20, a2f, step=da); n2 = length(a2)
    a3 = range(a30, a3f, step=da); n3 = length(a3)

    a  = Axis[]
    for i1 ∈ 1:n1, i2 ∈ 1:n2, i3 ∈ 1:n3
        push!(a,Axis(a1[i1],a2[i2],a3[i3]))
    end
    return a
end

# -- version 3 : given an array
@inline function semiaxes(a1::T,a2::T,a3::T) where T<:AbstractArray
    a  = Axis[]
    for i in eachindex(a1), j in eachindex(a2), k in eachindex(a3)
        push!(a,Axis(a1[i],a2[j],a3[k]))
    end
    return a
end
"""
dem(ηinc, ηmat, ϕf, dϕ, AX::Vector{Semiaxes{T}}; n1 = 100, n2 = 100)

    where ηinc is the viscosity of the inclusion, ηmat is the viscosity of the matrix, ϕf is the maximum volume fraction, dϕ is the increments of volume fraction,
    AX is the inclusion semiaxes, and n1 and n2 are the spatial discretization of θ and ϕ of the interaction tensor T

"""
function dem(ηinc,ηmat,ϕf,dϕ,AX;n1 = 100, n2 = 100)

    # -- Make sure isotropic viscosity is Float64 to ensure type stability
    if eltype(ηinc) ≠ Float64
        ηinc = convert(Float64,ηinc)
    end
    if eltype(ηmat) ≠ Float64
        ηmat = convert(Float64,ηmat)
    end

    # -- Range of volume fractions
    ϕ           = range(dϕ,ϕf,step=dϕ)

    # -- Get some arrays which constant throughtout the whole simulation
    cache       = getTcache(n1,n2)
    β           = βarray()    # constant array of arrays to compute 4th rank tensor inversion

    # -- Deviatoric viscosity tensors
    Js,Jd       = identityFourthOrder()     # Symmetric and Deviatoric 4th rank identity tensors
    Citensor    = 2 * ηinc * Jd             # Inclusion viscosity tensor
    Cmtensor    = 2 * ηmat * Jd             # Matrix viscosity tensor
    c0dem       = copy(Citensor)
    Ctemp       = copy(Cmtensor)

    # -- Voigt notation
    Cinc        = tensor2voigt(Citensor)  # Inclusion viscosity tensor in Voigts reduced notation
    Cmat        = tensor2voigt(Cmtensor)  # Matrix viscosity tensor in Voigts reduced notation

    # -- Output arrays
    CdemVoigt   = Array{Float64,2}[]
    outAx       = Array{Float64,1}[]
    outVol      = Float64[]

    # -- Semi axis iterations
    for i ∈ axes(AX,1)
        # - Solve DEM
        solvedem!(CdemVoigt,outAx,outVol,
            Js,Jd,Cmtensor,Citensor,AX[i],dϕ,ϕ,β,cache)
    end

    return  CdemVoigt, outAx, outVol
end

"""
Parallel version
"""
function parallel_launcher(ηinc,ηmat,ϕf,dϕ,AX)
    nblk = 4*2
    sblk = floor(Int,length(a1)/nblk)
    p    = 1:sblk

    t1   = time()
    @sync for b ∈ 1:nblk
        # C, outAx, outVol =
        Threads.@spawn begin
            main(ηinc,ηmat,ϕf,dϕ,AX[p])
            # main(ηinc,ηmat,ϕf,dϕ,AX[p])
            p = p.+sblk
            if b == nblk-1
                p = p[1]:length(a1)
            end
        end

    end
end

# DEM solver
function solvedem!(CdemVoigt,outAx,outVol,
    Js,Jd,Cmtensor,Citensor,AX,dϕ,ϕ,β,cache)

    c0dem   = copy(Cmtensor)

    # -- store results for ϕ = 0
    push!(CdemVoigt,tensor2voigt(c0dem))
    push!(outAx,[AX.a1, AX.a2, AX.a3 ])
    push!(outVol,0.0)

    ijkl    = @SMatrix Int8[1 6 5; 6 2 4; 5 4 3]    # voigt to tensor mapping
    A       = Array{Float64,2}(undef, 6, 6)

    # -- incremental volume fraction loop
    for ϕi ∈ ϕ
        # -- ~time integration
        c0dem   = RK4(Js,Jd,Citensor,c0dem,β, ϕi, dϕ , AX, cache)
        # -- to Voigt
        tensor2voigt!(A, c0dem, ijkl)
        # -- store
        push!(CdemVoigt, deepcopy(A))
        push!(outAx,[AX.a1, AX.a2, AX.a3])
        push!(outVol,ϕi)
    end

end

# -- Green function
@inline function greenfunction(Js, Jd, Ctemp,Cinc3, β, ϕ, AX,cache)
    T           = interactiontensor(Ctemp, AX, cache)
    Ts          = symmetricinteractiontensor(T)
    dC          = SymmetricTensor{4, 3}(Cinc3 - Ctemp)
    SdC         = Ts ⊡ dC                    # Eshelby tensor := S = Jd*T*Cm = Ts*Cm
    #=
        Add symmetric identity 4th-order tensor
        In reality it should be Ai = Id + SdC, which is however not invertible in the
        cartesian space (Lebenson et al., 1998)
    =#
    Ai          = SymmetricTensor{4, 3}(Js .+ SdC)
    # Ai          = Js + SdC
    invAi       = inv(Ai)                   # overloaded by Tensors.jl
    invAiDev    = Jd ⊡ invAi                # take deviatoric part of Ai
    dc0dem      = dC ⊡ invAiDev             # new composite stiffness tensor
    dc0dem      = map(x -> x * ϕ , dc0dem ) # multiply by volume fraction increment (slightly faster than normal .*)
    # return reinterpret(Float64,dc0dem)
    return SArray{NTuple{4,3}}(dc0dem)
end # end of Green function

"""
Runge-Kutta 4
"""
# -- Runge-Kuta 4
@inline function RK4(Js, Jd, Cinc3, c0dem, β, ϕ, dϕ , AX, cache)
    ϕi      = log(1 + dϕ / (1 - ϕ))
    # c0dem   = c0dem
    Ctemp   = copy(c0dem)
    # Ctemp   = SArray{NTuple{4,3}}(c0dem)
    c0dem   = SArray{NTuple{4,3}}(c0dem)
        # (1) :=
    dC1     = greenfunction(Js, Jd, Ctemp, Cinc3, β, ϕi, AX, cache)
    Ctemp   = c0dem + dC1 / 2
        # (2) :=
    dC2     = greenfunction(Js, Jd, Ctemp, Cinc3, β, ϕi, AX, cache)
    Ctemp   = c0dem + dC2 / 2
        # (3) :=
    dC3     = greenfunction(Js, Jd, Ctemp, Cinc3, β, ϕi, AX, cache)
    Ctemp   = c0dem + dC3
        # (4) :=
    dC4     = greenfunction(Js, Jd, Ctemp, Cinc3, β, ϕi, AX, cache)
        # Final average
    c0dem   = c0dem + (dC1 + 2*dC2 + 2*dC3 + dC4) / 6
    return c0dem
end # end of Runge Kutta integration

"""
Interaction tensor and auxiliar functions
"""
function getTcache(n1,n2)
    step    = π/n1
    dθ,dϕ   = π/n1, 2π/n2
    θ       = rad2deg.(range(dθ, stop = pi,  length = n1))
    ϕ       = rad2deg.(range(dϕ, stop = 2pi, length = n2))
    T       = @MArray zeros(3,3,3,3) #fill(0.0, 3, 3, 3, 3)
    Tv      = @MArray zeros(3,3,3,3) #fill(0.0, 3, 3, 3, 3)
    # Av      = fill(0.0, 4, 4)
    Av      = @MMatrix zeros(4,4)

    return Tcache(n1,n2,θ,ϕ,T,Tv,Av)
end

@inline function interactiontensor(C, AX, cache)

    T       = cache.T # unpack structure
    Av      = @MMatrix zeros(4,4)
    invAv   = similar(Av)
    xi      = Vector{Float64}(undef, 3)
    surface = 0.0
    for p ∈ 1:cache.n1

        # to avoid redundant allocations and operations
        sinθp,cosθp = sincosd(cache.θ[p])

        for q ∈ 1:Int(cache.n2/2)

            # -- to avoid redundant operations
            sinϕq,cosϕq = sincosd(cache.ϕ[q])

            # -- Director cosines
            xi = @SVector [sinθp*cosϕq/AX.a1, sinθp*sinϕq/AX.a2, cosθp/AX.a3]

            # -- Christoffel viscosity tensor
            Christoffel!(Av, C, xi)
            fillAv!(Av, xi)

            # -- Inverse of Christoffel tensor
            # inv4x4!(invAv,Av)
            invAv = inv(SMatrix{4,4}(Av))

            # -- Jiang, 2014, JSG
            tensorT!(T,invAv,xi,sinθp)

            surface += sinθp
        end
    end

    return T ./= surface
end

# -- Christoffel viscosity tensor
@inline function Christoffel!(Av,C,xi)
    @avx for t ∈ 1:3, r ∈ 1:3
        aux = zero(eltype(C))
        for u ∈ 1:3, s ∈ 1:3
            aux += C[r, s, t, u] * xi[s] * xi[u]
        end
        Av[r, t] = aux
    end
end

@inline function tensorT!(T,invAv,xi,sinθp)
    @avx for k ∈ 1:3, i ∈ 1:3
        aux = invAv[i, k]
        for l ∈ 1:3, j ∈ 1:3
            T[i, j, k, l] += aux * xi[j] * xi[l] * sinθp
        end
    end
end

@inline function fillAv!(Av, xi)
    @avx for i ∈ 1:3
        xi0 = xi[i]
        Av[i, 4] = xi0
        Av[4, i] = xi0
    end
end

###################################################
#     Convert C from voigt to tensor notation
###################################################
@inline function voigt2tensor(A::Array{Float64,2})
    ijkl    = Int8[1 6 5; 6 2 4; 5 4 3]    # voigt to tensor mapping
    B       = Array{Float64,4}(undef, 3, 3, 3, 3)
    B       = zero(SymmetricTensor{4, 3, Float32})
    @avx for i = 1:3
        for j = 1:3
            indx = ijkl[i, j]
            for k = 1:3
                for l = 1:3
                    B[i, j, k, l] = A[indx, ijkl[k, l]]
                end
            end
        end
    end
    B = SymmetricTensor{4, 3}(B)
    return B
end

###################################################
#   Convert C from tensor to voigt notation
###################################################
@inline function tensor2voigt(A)
    ijkl    = Int8[1 6 5; 6 2 4; 5 4 3]    # voigt to tensor mapping
    B       = Array{Float64,2}(undef, 6, 6)
    @avx for l = 1:3
        for k = 1:3
            for j = 1:3
                for i = 1:3
                    B[ijkl[i, j], ijkl[k, l]] = A[i, j, k, l]
                end
            end
        end
    end
    return B
end ## End of function

@inline function tensor2voigt!(B,A,ijkl) # zero allocations version
    @avx for l = 1:3
        for k = 1:3
            for j = 1:3
                for i = 1:3
                    B[ijkl[i, j], ijkl[k, l]] = A[i, j, k, l]
                end
            end
        end
    end
end ## End of function

###################################################
# DOUBLE CONTRACTION OF TWO 4TH RANK TENSORS
###################################################
@inline function contraction(A::Array{Float64,4}, B::Array{Float64,4})
    C = fill(0.0, 3, 3, 3, 3)
    @avx for i = 1:3, j = 1:3, k = 1:3, l = 1:3
        aux = zero(eltype(C))
        for m = 1:3, n = 1:3
            aux += A[i, j, m, n] * B[m, n, k, l]
        end
        C[i, j, k, l] = aux
    end
    return C
end ## END OF FUNCTION

###################################################
# DOUBLE CONTRACTION OF TWO 4TH RANK TENSORS
###################################################
@inline function identityFourthOrder()
    # deviatoric and symmetric 4th rank identity matrices
    Js = fill(0.0, 3, 3, 3, 3)
    Jd = similar(Js)
    @inbounds for l = 1:3, k = 1:3, j = 1:3, i = 1:3
        Js[i, j, k, l] = 0.5 * (I[i, k] * I[j, l] + I[i, l] * I[j, k])
        Jd[i, j, k, l] = Js[i, j, k, l] - I[i, j] * I[k, l] / 3
    end
    Js = SymmetricTensor{4, 3}(Js)
    Jd = SymmetricTensor{4, 3}(Jd)
    return Js, Jd
end ## END OF FUNCTION

###################################################
# DOUBLE CONTRACTION OF TWO 4TH RANK TENSORS
###################################################
@inline function symmetricinteractiontensor(T)
    Ts = fill(0.0, 3, 3, 3, 3)
    @avx for l = 1:3, k = 1:3, j = 1:3, i = 1:3
        Ts[i, j, k, l] +=
            0.25 * (T[i, j, k, l] + T[j, i, k, l] + T[i, j, l, k] + T[j, i, l, k])
    end
    force_symmetry!(Ts)
    Ts = SymmetricTensor{4, 3}(Ts)
    return Ts
end ## END OF FUNCTION

@inline function force_symmetry!(Ts)
    @inbounds for l = 1:3, k = 1:3, j = 1:3, i = 1:3
        Ts[j, i, k, l] = Ts[i, j, k, l]
        Ts[i, j, l, k] = Ts[i, j, k, l]
    end
end

###################################################
# Inversion of 4th rank tensor
###################################################
@inline function inversefourth(X, β)
    h       = 6
    Cab     = fill(0.0, h, h)
    Cinv    = fill(0.0, 3, 3, 3, 3)
    c1, c2  = 1, 1
    # -- step 1
    @inbounds for A ∈ 1:h, B ∈ 1:h
        βA, βB = β[A], β[B]
        @fastmath for i ∈ 1:3, j ∈ 1:3, k ∈ 1:3, l ∈ 1:3
            Cab[A, B] += X[i, j, k, l] * βA[i, j] * βB[k, l]
        end
    end
    # -- step 2
    invCab = inv(Cab)
    # -- step 3
    @inbounds for A ∈ 1:h, B ∈ 1:h
        βA, βB = β[A], β[B]
        @avx for i ∈ 1:3, j ∈ 1:3, k ∈ 1:3, l ∈ 1:3
            Cinv[i, j, k, l] += invCab[A, B] * βA[i, j] * βB[k, l]
        end
    end
    Cinv = SymmetricTensor{4, 3}(Cinv)
    return Cinv
end ## END OF FUNCTION

@inline function βarray()
    # create β array of arrays
    b1 = [
        -1.0 0.0 0.0
        0.0 -1.0 0.0
        0.0 0.0 2.0
    ] ./ sqrt(6)
    b2 = [
        -1.0 0.0 0.0
        0.0 1.0 0.0
        0.0 0.0 0.0
    ] ./ sqrt(2)
    b3 = [
        0.0 0.0 0.0
        0.0 0.0 1.0
        0.0 1.0 0.0
    ] ./ sqrt(2)
    b4 = [
        0.0 0.0 1.0
        0.0 0.0 0.0
        1.0 0.0 0.0
    ] ./ sqrt(2)
    b5 = [
        0.0 1.0 0.0
        1.0 0.0 0.0
        0.0 0.0 0
    ] / sqrt(3)
    b6 = Diagonal(fill(sqrt(3), 3))
    β  = [b1, b2, b3, b4, b5, b6]
end

"""
Save output in CSV file
"""
function saveIO(C, outAx, outVol,outName)
    n       = length(C)
    # -- Type conversions (F64 -> F32 to lighten file size)
    C       = convert(Array{Array{Float32,2},1},C)
    outAx   = convert(Array{Array{Float32,1},1},outAx)
    outVol  = convert(Array{Float32,1},outVol)
    # -- Allocations
    η11     = Vector{Float32}(undef,n)
    η22     = similar(η11)
    η33     = similar(η11)
    η44     = similar(η11)
    η55     = similar(η11)
    η66     = similar(η11)
    η12     = similar(η11)
    η13     = similar(η11)
    η23     = similar(η11)
    a1      = similar(η11)
    a2      = similar(η11)
    a3      = similar(η11)
    ϕ       = similar(η11)
    # -- Store data in vectors
    Threads.@threads for i ∈ 1:n
        # -- Semi axis
        a1[i]   = outAx[i][1]
        a2[i]   = outAx[i][2]
        a3[i]   = outAx[i][3]
        # -- Volume fraction
        ϕ[i]    = outVol[i]
        # -- Viscous tensor components
        DEM     = C[i]
        η11[i]  = C[i][1,1]
        η22[i]  = C[i][2,2]
        η33[i]  = C[i][3,3]
        η44[i]  = C[i][4,4]
        η55[i]  = C[i][5,5]
        η66[i]  = C[i][6,6]
        η12[i]  = C[i][1,2]
        η13[i]  = C[i][1,3]
        η23[i]  = C[i][2,3]
    end
    # -- Make Data Frame
    df = DataFrame( nu11 = η11,
                    nu22 = η22,
                    nu33 = η33,
                    nu44 = η44,
                    nu55 = η55,
                    nu66 = η66,
                    nu12 = η12,
                    nu13 = η13,
                    nu23 = η23,
                    a1   = a1,
                    a2   = a2,
                    a3   = a3,
                    phi  = ϕ
                    )
    # -- Save data as CSV
    CSV.write(outName,df)
end

## 4 by 4 matrix functions
@inline function det4x4(A)
    detA = @inbounds @fastmath (A[1,1]*A[2,2]*A[3,3]*A[4,4] +
            A[1,1]*A[2,3]*A[3,4]*A[4,2] +
            A[1,1]*A[2,4]*A[3,2]*A[4,3] -
            A[1,1]*A[2,4]*A[3,3]*A[4,2] -
            A[1,1]*A[2,3]*A[3,2]*A[4,4] -
            A[1,1]*A[2,2]*A[3,4]*A[4,3] -
            A[1,2]*A[2,1]*A[3,3]*A[4,4] -
            A[1,3]*A[2,1]*A[3,4]*A[4,2] -
            A[1,3]*A[2,1]*A[3,2]*A[4,3] +
            A[1,4]*A[2,1]*A[3,3]*A[4,2] +
            A[1,3]*A[2,1]*A[3,2]*A[4,4] +
            A[1,2]*A[2,1]*A[3,4]*A[4,3] +
            A[1,2]*A[2,3]*A[3,1]*A[4,4] +
            A[1,3]*A[2,4]*A[3,1]*A[4,2] +
            A[1,4]*A[2,2]*A[3,1]*A[4,3] -
            A[1,4]*A[2,3]*A[3,1]*A[4,2] -
            A[1,3]*A[2,2]*A[3,1]*A[4,4] -
            A[1,2]*A[2,4]*A[3,1]*A[4,3] -
            A[1,2]*A[2,3]*A[3,4]*A[4,1] -
            A[1,3]*A[2,4]*A[3,2]*A[4,1] -
            A[1,4]*A[2,2]*A[3,3]*A[4,1] +
            A[1,4]*A[2,3]*A[3,2]*A[4,1] +
            A[1,3]*A[2,2]*A[3,4]*A[4,1] +
            A[1,2]*A[2,4]*A[3,3]*A[4,1]
            )
end

@fastmath function adjA!(Â,A)
    Â[1,1] =@inbounds (A[2,2]*A[3,3]*A[4,4] + A[2,3]*A[3,4]*A[4,2] + A[2,4]*A[3,2]*A[4,3] -
                     A[2,4]*A[3,3]*A[4,2] - A[2,3]*A[3,2]*A[4,4] - A[2,2]*A[3,4]*A[4,3] )

    Â[1,2] =@inbounds (-A[1,2]*A[3,3]*A[4,4] - A[1,3]*A[3,4]*A[4,2] - A[1,4]*A[3,2]*A[4,3] +
                      A[1,4]*A[3,3]*A[4,2] + A[1,3]*A[3,2]*A[4,4] + A[1,2]*A[3,4]*A[4,3] )

    Â[1,3] =@inbounds (A[1,2]*A[2,3]*A[4,4] + A[1,3]*A[2,4]*A[4,2] + A[1,4]*A[2,2]*A[4,3] -
                     A[1,4]*A[2,3]*A[4,2] - A[1,3]*A[2,2]*A[4,4] - A[1,2]*A[2,4]*A[4,3] )

    Â[1,4] =@inbounds (-A[1,2]*A[2,3]*A[3,4] - A[1,3]*A[2,4]*A[3,2] - A[1,4]*A[2,2]*A[3,3] +
                      A[1,4]*A[2,3]*A[3,2] + A[1,3]*A[2,2]*A[3,4] + A[1,2]*A[2,4]*A[3,3] )

    Â[2,1] =@inbounds (-A[2,1]*A[3,3]*A[4,4] - A[2,3]*A[3,4]*A[4,1] - A[2,4]*A[3,1]*A[4,3] +
                      A[2,4]*A[3,3]*A[4,1] + A[2,3]*A[3,1]*A[4,4] + A[2,1]*A[3,4]*A[4,3] )

    Â[2,2] =@inbounds (A[1,1]*A[3,3]*A[4,4] + A[1,3]*A[3,4]*A[4,2] + A[1,4]*A[3,1]*A[4,3] -
                     A[1,4]*A[3,3]*A[4,1] - A[1,3]*A[3,1]*A[4,4] - A[1,1]*A[3,4]*A[4,3] )

    Â[2,3] =@inbounds (-A[1,1]*A[2,3]*A[4,4] - A[1,3]*A[2,4]*A[4,1] - A[1,4]*A[2,1]*A[4,3] +
                      A[1,4]*A[2,3]*A[4,1] + A[1,3]*A[2,1]*A[4,4] + A[1,1]*A[2,4]*A[4,3] )

    Â[2,4] =@inbounds (A[1,1]*A[2,3]*A[3,4] + A[1,3]*A[2,4]*A[3,1] + A[1,4]*A[2,1]*A[3,3] -
                     A[1,4]*A[2,3]*A[3,1] - A[1,3]*A[2,1]*A[3,4] - A[1,1]*A[2,4]*A[3,3] )

    Â[3,1] =@inbounds (A[2,1]*A[3,2]*A[4,4] + A[2,2]*A[3,4]*A[4,1] + A[2,4]*A[3,1]*A[4,2] -
                     A[2,4]*A[3,2]*A[4,1] - A[2,2]*A[3,1]*A[4,4] - A[2,1]*A[3,4]*A[4,2] )

    Â[3,2] =@inbounds (-A[1,1]*A[3,2]*A[4,4] - A[1,2]*A[3,4]*A[4,1] - A[1,4]*A[3,1]*A[4,2] +
                      A[1,4]*A[3,2]*A[4,1] + A[1,2]*A[3,1]*A[4,4] + A[1,1]*A[3,4]*A[4,2] )

    Â[3,3] =@inbounds (A[1,1]*A[2,2]*A[4,4] + A[1,2]*A[2,4]*A[4,1] + A[1,4]*A[2,1]*A[4,2] -
                     A[1,4]*A[2,2]*A[4,1] - A[1,2]*A[2,1]*A[4,4] - A[1,1]*A[2,4]*A[4,2] )

    Â[3,4] =@inbounds (-A[1,1]*A[2,2]*A[3,4] - A[1,2]*A[2,4]*A[3,1] - A[1,4]*A[2,1]*A[3,2] +
                      A[1,4]*A[2,2]*A[3,1] + A[1,2]*A[2,1]*A[3,4] + A[1,1]*A[2,4]*A[3,2] )

    Â[4,1] =@inbounds (-A[2,1]*A[3,2]*A[4,3] - A[2,2]*A[3,3]*A[4,1] - A[2,3]*A[3,1]*A[4,2] +
                      A[2,3]*A[3,2]*A[4,1] + A[2,2]*A[3,1]*A[4,3] + A[2,1]*A[3,3]*A[4,2] )

    Â[4,2] =@inbounds (A[1,1]*A[3,2]*A[4,3] + A[1,2]*A[3,3]*A[4,1] + A[1,3]*A[3,1]*A[4,2] -
                     A[1,3]*A[3,2]*A[4,1] - A[1,2]*A[3,1]*A[4,3] - A[1,1]*A[3,3]*A[4,2] )

    Â[4,3] =@inbounds (-A[1,1]*A[2,2]*A[4,3] - A[1,2]*A[2,3]*A[4,1] - A[1,3]*A[2,1]*A[4,2] +
                      A[1,3]*A[2,2]*A[4,1] + A[1,2]*A[2,1]*A[4,3] + A[1,1]*A[2,3]*A[4,2] )

    Â[4,4] =@inbounds (A[1,1]*A[2,2]*A[3,3] + A[1,2]*A[2,3]*A[3,1] + A[1,3]*A[2,1]*A[3,2] -
                     A[1,3]*A[2,2]*A[3,1] - A[1,2]*A[2,1]*A[3,3] - A[1,1]*A[2,3]*A[3,2] )
end

function inv4x4!(Â,A)
    adjA!(Â,A)
    d = det4x4(A)
    @avx for j ∈ axes(A,2),i ∈ axes(A,1)
        Â[i,j] /= d
    end
end

function cleantensor!(C)
    mask1 = [1:3,4:6]
    mask2 = [4:6,1:3]
    @inbounds @simd for ll ∈ eachindex(C)
        @views C[ll][mask1[1],mask1[2]] .= 0.0
        @views C[ll][mask2[1],mask2[2]] .= 0.0
        for j ∈ 4:6, i ∈ 4:6
            if i ≠ j
                C[ll][i,j] = 0.0
            end
        end
    end
end

end

max_stars_repo_name: albert-de-montserrat/DEM

max_issues_repo_name: albert-de-montserrat/DEM

max_forks_repo_name: albert-de-montserrat/DEM


text:
using LHEF
using Test

const SAMPLES_DIR = joinpath(@__DIR__, "")

@testset "laziness" begin
    lazyevents = parse_lhe(joinpath(SAMPLES_DIR, "ft.lhe"))
    @test lazyevents isa Base.Generator
    events = collect(lazyevents)
    @test events isa Vector
    @test length(events) == 45
end

@testset "parsing" begin
    events = parse_lhe(joinpath(SAMPLES_DIR, "ft.lhe")) |> collect
    event = events[1]
    particles = event.particles
    header = event.header
    @test values(header) == (6, 0, 0.1182900E-04, 0.2556536E+03, 0.7546771E-02, 0.1112889E+00)
    @test length(particles) == 6
    @test values(particles[1]) == (0, 21, -1, 0, 0, 502, 503, 0.00000000000E+00 , 0.00000000000E+00 , 0.10709531583E+04 , 0.10709531583E+04, 0.00000000000E+00, 0., -1.)
    @test values(particles[2]) == (1, 21, -1, 0, 0, 501, 504, 0.00000000000E+00 , 0.00000000000E+00 , -0.77476002582E+03, 0.77476002582E+03, 0.00000000000E+00, 0., 1.)
    @test values(particles[3]) == (2, 6 , 1 , 1, 2, 501, 0  , 0.11337785248E+03 , 0.11416185862E+03 , -0.41887649846E+02, 0.23993966451E+03, 0.17300000000E+03, 0., 1.)
    @test values(particles[4]) == (3, 6 , 1 , 1, 2, 502, 0  , 0.34597641987E+02 , -0.27246642769E+03,-0.24576811815E+03 , 0.40714360973E+03, 0.17300000000E+03, 0., 1.)
    @test values(particles[5]) == (4, -6, 1 , 1, 2, 0  , 503, 0.15534573574E+02 , 0.18289123966E+03 , 0.82271340950E+03 , 0.86050966450E+03, 0.17300000000E+03, 0., -1.)
    @test values(particles[6]) == (5, -6, 1 , 1, 2, 0  , 504, -0.16351006804E+03,-0.24586670591E+02 ,-0.23886450899E+03 , 0.33812024543E+03, 0.17300000000E+03, 0., -1.)
end

@testset "gzip" begin
    e1 = parse_lhe(joinpath(SAMPLES_DIR, "ft.lhe")) |> collect
    e2 = parse_lhe(joinpath(SAMPLES_DIR, "ft.lhe.gz")) |> collect
    @test all([e1[i].particles == e2[i].particles for i in 1:length(e1)])
end

@testset "flattened particles" begin
    particles = flatparticles(joinpath(SAMPLES_DIR, "ft.lhe.gz"))
    @test length(particles) == 270
    @test maximum(p.eventnum for p in particles) == 45
end

ext: jl


text:
module DifferentialAlgebra

end # End of module


text:
function GMT_Create_Session(tag::String="GMT", pad=2, mode=0, print_func::Ptr{Cvoid}=C_NULL)
	API = ccall((:GMT_Create_Session, libgmt), Ptr{Cvoid}, (Ptr{UInt8}, UInt32, UInt32, Ptr{Cvoid}), tag, pad, mode, print_func)
	(API == C_NULL) && error("Failure to create a GMT Session")
	API
end

function GMT_Create_Data(API::Ptr{Cvoid}, family::Integer, geometry, mode, dim=NULL, wesn=NULL,
                         inc=NULL, registration=0, pad=2, data::Ptr{Cvoid}=NULL)
	ptr = ccall((:GMT_Create_Data, libgmt), Ptr{Cvoid}, (Cstring, UInt32, UInt32, UInt32, Ptr{UInt64},
		Ptr{Cdouble}, Ptr{Cdouble}, UInt32, Cint, Ptr{Cvoid}), API, family, geometry, mode, dim, wesn, inc,
		registration, pad, data)

	(ptr == C_NULL) && error("Failure to allocate GMT resource")
	ptr
end

GMT_Destroy_Session(API::Ptr{Cvoid}) = ccall((:GMT_Destroy_Session, libgmt), Cint, (Cstring,), API)

#=		Not used yet, so comment
function GMT_Read_Data(API::Ptr{Cvoid}, family, method, geometry, mode, wesn, input=C_NULL, data=C_NULL)

	if (family == GMT_IS_DATASET)
		ret_type = Ptr{GMT_DATASET}
	elseif (family == GMT_IS_GRID)
		ret_type = Ptr{GMT_GRID}
	elseif (family == GMT_IS_PALETTE)
		ret_type = Ptr{GMT_PALETTE}
	elseif (family == GMT_IS_IMAGE)
		ret_type = Ptr{GMT_IMAGE}
	elseif (family == GMT_IS_MATRIX)
		ret_type = Ptr{GMT_MATRIX}
	elseif (family == GMT_IS_VECTOR)
		ret_type = Ptr{GMT_VECTOR}
	else
		ret_type = Ptr{Cvoid}			# Should be error instead
	end

	ptr = ccall( (:GMT_Read_Data, libgmt), Ptr{Cvoid}, (Ptr{Cvoid}, UInt32, UInt32, UInt32, UInt32, Ptr{Cdouble},
		Ptr{UInt8}, Ptr{Cvoid}), API, family, method, geometry, mode, wesn, input, data)

	convert(ret_type, ptr)
end

function GMT_Insert_Data(API::Ptr{Cvoid}, object_ID::Integer, data)
	ccall((:GMT_Insert_Data, libgmt), Cint, (Ptr{Cvoid}, Cint, Ptr{Cvoid}), API, object_ID, data)
end

function GMT_Duplicate_Data(API::Ptr{Cvoid}, family::Integer, mode::Integer, data::Ptr{Cvoid})
	ccall( (:GMT_Duplicate_Data, libgmt), Ptr{Cvoid}, (Ptr{Cvoid}, UInt32, UInt32, Ptr{Cvoid}), API, family, mode, data)
end

function GMT_Get_Record(API::Ptr{Cvoid}, mode::Integer, retval::Ptr{Cint})
	ccall( (:GMT_Get_Record, libgmt), Ptr{Cvoid}, (Ptr{Cvoid}, UInt32, Ptr{Cint}), API, mode, retval)
end

function GMT_Manage_Session(API::Ptr{Cvoid}, mode::Integer, args::Ptr{Cvoid})
	ccall( (:GMT_Manage_Session, libgmt), Cint, (Ptr{Cvoid}, UInt32, Ptr{Cvoid}), API, mode, args)
end

function GMT_Register_IO(API, family::Integer, method::Integer, geometry::Integer, direction::Integer=0,
		wesn=C_NULL, resource=C_NULL)
	err = ccall((:GMT_Register_IO, libgmt), Cint, (Ptr{Cvoid}, UInt32, UInt32, UInt32, UInt32, Ptr{Cdouble},
		Ptr{Cvoid}), API, family, method, geometry, direction, wesn, resource)
end

function GMT_Init_IO(API::Ptr{Cvoid}, family::UInt32, geometry::UInt32, direction::UInt32, mode::UInt32, n_args::UInt32, args::Ptr{Cvoid})
	ccall( (:GMT_Init_IO, libgmt), Cint, (Ptr{Cvoid}, UInt32, UInt32, UInt32, UInt32, UInt32, Ptr{Cvoid}), API, family, geometry, direction, mode, n_args, args)
end
function GMT_Begin_IO(API::Ptr{Cvoid}, family::UInt32, direction::UInt32, header::UInt32)
	ccall( (:GMT_Begin_IO, libgmt), Cint, (Ptr{Cvoid}, UInt32, UInt32, UInt32), API, family, direction, header)
end
function GMT_Status_IO(API::Ptr{Cvoid}, mode::UInt32)
	ccall( (:GMT_Status_IO, libgmt), Cint, (Ptr{Cvoid}, UInt32), API, mode)
end
function GMT_End_IO(API::Ptr{Cvoid}, direction::UInt32, mode::UInt32)
	ccall( (:GMT_End_IO, libgmt), Cint, (Ptr{Cvoid}, UInt32, UInt32), API, direction, mode)
end

function GMT_Write_Data(API::Ptr{Cvoid}, family::Integer, method::Integer, geometry::Integer, mode::Integer,
	wesn, output::String, data)
	err = ccall((:GMT_Write_Data, libgmt), Cint, (Ptr{Cvoid}, UInt32, UInt32, UInt32, UInt32, Ptr{Cdouble},
		Ptr{UInt8}, Ptr{Cvoid}), API, family, method, geometry, mode, wesn, output, data)
end =#

function GMT_Destroy_Data(API::Ptr{Cvoid}, object)
	ccall((:GMT_Destroy_Data, libgmt), Cint, (Cstring, Ptr{Cvoid}), API, object)
end

function GMT_Set_Comment(API::Ptr{Cvoid}, family::Integer, mode, arg::Ptr{Cvoid}, data::Ptr{Cvoid})
	ccall((:GMT_Set_Comment, libgmt), Cint, (Cstring, UInt32, UInt32, Ptr{Cvoid}, Ptr{Cvoid}), API, family, mode, arg, data)
end

#=
function GMT_Put_Record(API::Ptr{Cvoid}, mode::UInt32, record::Ptr{Cvoid})
	ccall( (:GMT_Put_Record, libgmt), Cint, (Ptr{Cvoid}, UInt32, Ptr{Cvoid}), API, mode, record)
end

function GMT_Encode_ID(API::Ptr{Cvoid}, fname::String, object_ID::Integer)
	err = ccall((:GMT_Encode_ID, libgmt), Cint, (Ptr{Cvoid}, Ptr{UInt8}, Cint), API, fname, object_ID)
end

function GMT_Get_Row(API::Ptr{Cvoid}, rec_no::Cint, G::Ptr{GMT_GRID}, row::Ptr{Cfloat})
	ccall( (:GMT_Get_Row, libgmt), Cint, (Ptr{Cvoid}, Cint, Ptr{GMT_GRID}, Ptr{Cfloat}), API, rec_no, G, row)
end
function GMT_Put_Row(API::Ptr{Cvoid}, rec_no::Cint, G::Ptr{GMT_GRID}, row::Ptr{Cfloat})
	ccall( (:GMT_Put_Row, libgmt), Cint, (Ptr{Cvoid}, Cint, Ptr{GMT_GRID}, Ptr{Cfloat}), API, rec_no, G, row)
end

function GMT_Get_ID(API::Ptr{Cvoid}, family::Integer, dir::Integer, resource=C_NULL)
	ccall((:GMT_Get_ID, libgmt), Cint, (Ptr{Cvoid}, UInt32, UInt32, Ptr{Cvoid}), API, family, dir, resource)
end

function GMT_Get_Family(API::Ptr{Cvoid}, dir::Integer, head::Ptr{GMT_OPTION})
	ccall((:GMT_Get_Family, libgmt), Cint, (Ptr{Cvoid}, UInt32, Ptr{GMT_OPTION}), API, dir, head)
end

function GMT_Get_Index(API::Ptr{Cvoid}, header::Ptr{GMT_GRID_HEADER}, row::Cint, col::Cint)
	ccall((:GMT_Get_Index, libgmt), Clonglong, (Ptr{Cvoid}, Ptr{GMT_GRID_HEADER}, Cint, Cint), API, header, row, col)
end
function GMT_Get_Coord(API::Ptr{Cvoid}, family::Integer, dim::Integer, container::Ptr{Cvoid})
	ccall((:GMT_Get_Coord, libgmt), Ptr{Cdouble}, (Ptr{Cvoid}, UInt32, UInt32, Ptr{Cvoid}), API, family, dim, container)
end

function GMT_Option(API::Ptr{Cvoid}, options)
	ccall((:GMT_Option, libgmt), Cint, (Ptr{Cvoid}, Ptr{UInt8}), API, options)
end

function GMT_Get_Value(API::Ptr{Cvoid}, arg::String, par::Ptr{Cdouble})
	ccall((:GMT_Get_Value, libgmt), Cint, (Ptr{Cvoid}, Ptr{UInt8}, Ptr{Cdouble}), API, arg, par)
end
function GMT_Get_Values(API::Ptr{Cvoid}, arg::String, par::Ptr{Cdouble}, maxpar::Integer)
	ccall((:GMT_Get_Values, libgmt), Cint, (Ptr{Cvoid}, Ptr{UInt8}, Ptr{Cdouble}, Cint), API, arg, par, maxpar)
end =#

function GMT_Get_Common(API::Ptr{Cvoid}, option::Char)
	buffer = Vector{Float64}(undef, 4)
	n_par = ccall((:GMT_Get_Common, libgmt), Cint, (Ptr{Cvoid}, UInt32, Ptr{Cdouble}), API, option, buffer)
	return buffer, n_par
end

function GMT_Get_Default(API::Ptr{Cvoid}, keyword::String, value)
    ccall((:GMT_Get_Default, libgmt), Cint, (Cstring, Ptr{UInt8}, Ptr{UInt8}), API, keyword, value)
end

function GMT_Call_Module(API::Ptr{Cvoid}, _module::AbstractString, mode=0, args=C_NULL)
	#(isa(args, String)) && (args = pointer(args))
	ccall((:GMT_Call_Module, libgmt), Cint, (Cstring, Ptr{UInt8}, Cint, Ptr{Cvoid}), API, _module, mode, args)
end

function GMT_Create_Options(API::Ptr{Cvoid}, argc::Int, args)
	# VERSATILIZAR PARA O CASO DE ARGS SER STRING OU ARRAY DE STRINGS
	# ccall((:GMT_Create_Options, libgmt), Ptr{GMT_OPTION}, (Ptr{Cvoid}, Cint, Ptr{Cvoid}), API, argc, args)
	ccall((:GMT_Create_Options, libgmt), Ptr{GMT_OPTION}, (Cstring, Cint, Cstring), API, argc, args)
end
# GMT_Create_Options(API::Ptr{Cvoid}, argc::Integer, args::String) =
#                   GMT_Create_Options(API, argc, convert(Ptr{Cvoid},pointer(args)))

function GMT_Destroy_Options(API::Ptr{Cvoid}, head::Ref{Ptr{GMT_OPTION}})
	ccall((:GMT_Destroy_Options, libgmt), Cint, (Cstring, Ref{Ptr{GMT_OPTION}}), API, head)
end

GMT_Create_Cmd(API::Ptr{Cvoid}, head::Ptr{GMT_OPTION}) =
	ccall( (:GMT_Create_Cmd, libgmt), Ptr{UInt8}, (Ptr{Cvoid}, Ptr{GMT_OPTION}), API, head)

#=
function GMT_Make_Option(API::Ptr{Cvoid}, option::UInt8, arg::Ptr{UInt8})
	ccall((:GMT_Make_Option, libgmt), Ptr{GMT_OPTION}, (Ptr{Cvoid}, UInt8, Ptr{UInt8}), API, option, arg)
end

function GMT_Find_Option(API::Ptr{Cvoid}, option::UInt8, head::Ptr{GMT_OPTION})
	ccall((:GMT_Find_Option, libgmt), Ptr{GMT_OPTION}, (Ptr{Cvoid}, UInt8, Ptr{GMT_OPTION}), API, option, head)
end

function GMT_Append_Option(API::Ptr{Cvoid}, current::Ptr{GMT_OPTION}, head::Ptr{GMT_OPTION})
	ccall((:GMT_Append_Option, libgmt), Ptr{GMT_OPTION}, (Ptr{Cvoid}, Ptr{GMT_OPTION}, Ptr{GMT_OPTION}), API, current, head)
end
function GMT_Create_Args(API::Ptr{Cvoid}, argc::Ptr{Int}, head::Ptr{GMT_OPTION})
	ccall( (:GMT_Create_Args, libgmt), Ptr{Ptr{UInt8}}, (Ptr{Cvoid}, Ptr{Cint}, Ptr{GMT_OPTION}), API, argc, head)
end

function GMT_Destroy_Args(API::Ptr{Cvoid}, argc::Cint, argv::Ptr{Ptr{Ptr{UInt8}}})
	ccall( (:GMT_Destroy_Args, libgmt), Cint, (Ptr{Cvoid}, Cint, Ptr{Ptr{Ptr{UInt8}}}), API, argc, argv)
end
function GMT_Destroy_Cmd(API::Ptr{Cvoid}, cmd::Ptr{Ptr{UInt8}})
	ccall( (:GMT_Destroy_Cmd, libgmt), Cint, (Ptr{Cvoid}, Ptr{Ptr{UInt8}}), API, cmd)
end
function GMT_Update_Option(API::Ptr{Cvoid}, current::Ptr{GMT_OPTION}, arg::Ptr{UInt8})
	ccall( (:GMT_Update_Option, libgmt), Cint, (Ptr{Cvoid}, Ptr{GMT_OPTION}, Ptr{UInt8}), API, current, arg)
end
function GMT_Delete_Option(API::Ptr{Cvoid}, current::Ptr{GMT_OPTION}, head::Ptr{Ptr{GMT_OPTION}})
	ccall( (:GMT_Delete_Option, libgmt), Cint, (Ptr{Cvoid}, Ptr{GMT_OPTION}, Ptr{Ptr{GMT_OPTION}}), API, current, head)
end
function GMT_Parse_Common(API::Ptr{Cvoid}, given_options::Ptr{UInt8}, options::Ptr{GMT_OPTION})
	ccall( (:GMT_Parse_Common, libgmt), Cint, (Ptr{Cvoid}, Ptr{UInt8}, Ptr{GMT_OPTION}), API, given_options, options)
end

function GMT_FFT_Option(API::Ptr{Cvoid}, option::UInt8, dim::UInt32, string::Ptr{UInt8})
	ccall( (:GMT_FFT_Option, libgmt), UInt32, (Ptr{Cvoid}, UInt8, UInt32, Ptr{UInt8}), API, option, dim, string)
end
function GMT_FFT_Parse(API::Ptr{Cvoid}, option::UInt8, dim::UInt32, args::Ptr{UInt8})
	ccall( (:GMT_FFT_Parse, libgmt), Ptr{Cvoid}, (Ptr{Cvoid}, UInt8, UInt32, Ptr{UInt8}), API, option, dim, args)
end
function GMT_FFT_Create(API::Ptr{Cvoid}, X::Ptr{Cvoid}, dim::UInt32, mode::UInt32, F::Ptr{Cvoid})
	ccall( (:GMT_FFT_Create, libgmt), Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{Cvoid}, UInt32, UInt32, Ptr{Cvoid}), API, X, dim, mode, F)
end
function GMT_FFT_Wavenumber(API::Ptr{Cvoid}, k::UInt16, mode::UInt32, K::Ptr{Cvoid})
	ccall( (:GMT_FFT_Wavenumber, libgmt), Cdouble, (Ptr{Cvoid}, UInt16, UInt32, Ptr{Cvoid}), API, k, mode, K)
end
function GMT_FFT(API::Ptr{Cvoid}, X::Ptr{Cvoid}, direction::Cint, mode::UInt32, K::Ptr{Cvoid})
	ccall( (:GMT_FFT, libgmt), Cint, (Ptr{Cvoid}, Ptr{Cvoid}, Cint, UInt32, Ptr{Cvoid}), API, X, direction, mode, K)
end
function GMT_FFT_Destroy(API::Ptr{Cvoid}, K::Ptr{Cvoid})
	ccall( (:GMT_FFT_Destroy, libgmt), Cint, (Ptr{Cvoid}, Ptr{Cvoid}), API, K)
end
function GMT_FFT_1D(API::Ptr{Cvoid}, data::Ptr{Cfloat}, n::UInt16, direction::Cint, mode::UInt32)
	ccall( (:GMT_FFT_1D, libgmt), Cint, (Ptr{Cvoid}, Ptr{Cfloat}, UInt16, Cint, UInt32), API, data, n, direction, mode)
end
function GMT_FFT_2D(API::Ptr{Cvoid}, data::Ptr{Cfloat}, nx::UInt32, ny::UInt32, direction::Cint, mode::UInt32)
	ccall( (:GMT_FFT_2D, libgmt), Cint, (Ptr{Cvoid}, Ptr{Cfloat}, UInt32, UInt32, Cint, UInt32), API, data, nx, ny, direction, mode)
end

function GMT_Report(API, vlevel::Integer, txt)
	ccall((:GMT_Report, libgmt), Cvoid, (Cstring, Cint, Ptr{UInt8}), API, vlevel, txt)
end =#

function GMT_Encode_Options(V_API::Ptr{Cvoid}, _module, n_argin::Int, head::Ref{Ptr{GMT_OPTION}}, n::Ptr{Int})
	ccall((:GMT_Encode_Options, libgmt), Ptr{GMT_RESOURCE}, (Cstring, Ptr{UInt8}, Int32, Ref{Ptr{GMT_OPTION}},
					Ptr{UInt32}), V_API, _module, n_argin, head, n)
end
function GMT_Encode_Options(V_API::Ptr{Cvoid}, _module, n_argin::Int, head, n::Ptr{Int})
	ccall((:GMT_Encode_Options, libgmt), Ptr{GMT_RESOURCE}, (Cstring, Ptr{UInt8}, Int32, Ptr{Ptr{Cvoid}}, Ptr{UInt32}),
					V_API, _module, n_argin, head, n)
end

function GMT_Expand_Option(V_API::Ptr{Cvoid}, opt::Ptr{GMT_OPTION}, arg)
	ccall((:GMT_Expand_Option, libgmt), Cint, (Cstring, Ptr{GMT_OPTION}, Ptr{UInt8}), V_API, opt, arg)
end

#=
function gmt_core_module_info(API, candidate)
	ccall((:gmt_core_module_info, libgmt), Ptr{UInt8}, (Ptr{Cvoid}, Ptr{UInt8}), API, candidate)
end

function gmtlib_grd_flip_vertical(gridp, n_cols::Integer, n_rows::Integer, n_stride::Integer=0, cell_size::Integer=1)
	ccall((:gmtlib_grd_flip_vertical, libgmt), Cvoid, (Ptr{Cvoid}, UInt32, UInt32, UInt32, Csize_t),
				 gridp, n_cols, n_rows, n_stride, cell_size)
end

function GMT_set_mem_layout(API, mem_layout)
	ccall((:GMT_set_mem_layout, libgmt), Cvoid, (Ptr{Cvoid}, Ptr{UInt8}), API, mem_layout)
end =#

GMT_Set_Default(API::Ptr{Cvoid}, keyword, value) =
	ccall((:GMT_Set_Default, libgmt), Cvoid, (Cstring, Ptr{UInt8}, Ptr{UInt8}), API, keyword, value)

function GMT_blind_change_struct(API::Ptr{Cvoid}, X, what, keyword::String, off::Integer)
	(GMTver > v"6.0") ?		# Use this construct to cheat Coverage
		ccall((:gmtlib_blind_change_struct, libgmt), Cint, (Cstring, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{UInt8}, Csize_t), API, X, what, keyword, off) : ccall((:GMT_blind_change_struct, libgmt), Cint, (Cstring, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{UInt8}, Csize_t), API, X, what, keyword, off)
end

#=
function GMT_Convert_Data(API::Ptr{Cvoid}, In::Ptr{Cvoid}, family_in::Integer, out::Ptr{Cvoid}, family_out::Integer, flag)
	ccall((:GMT_Convert_Data, libgmt), Ptr{Cvoid}, (Cstring, Ptr{Cvoid}, UInt32, Ptr{Cvoid}, UInt32, Ptr{UInt32}), API, In,
				 family_in, out, family_out, flag)
end

function GMT_blind_change_struct_(API::Ptr{Cvoid}, X, what, keyword::String)
	ccall((:GMT_blind_change_struct_, libgmt), Cint, (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{UInt8}),
				 API, X, what, keyword)
end

function GMT_Read_Group(API::Ptr{Cvoid}, family::Integer, method::Integer, geometry::Integer, mode::Integer,
												wesn, sources::Ptr{Cvoid}, n_items::Ptr{UInt32}, data::Ptr{Cvoid})
	ccall((:GMT_Read_Group, libgmt), Ptr{Cvoid}, (Ptr{Cvoid}, UInt32, UInt32, UInt32, UInt32, Ptr{Cdouble},
				 Ptr{Cvoid}, Ptr{UInt32}, Ptr{Cvoid}), API, family, method, geometry, mode, wesn, sources, n_items, data)
end

function GMT_Get_Index(API::Ptr{Cvoid}, header::Ptr{GMT_GRID_HEADER}, row::Integer, col::Integer)
	ccall((:GMT_Get_Index, libgmt), UInt64, (Ptr{Cvoid}, Ptr{GMT_GRID_HEADER}, Cint, Cint), API, header, row, col)
end

function GMT_Get_Pixel(API::Ptr{Cvoid}, header::Ptr{GMT_GRID_HEADER}, row::Integer, col::Integer, layer::Integer)
	ccall((:GMT_Get_Pixel, libgmt), UInt64, (Ptr{Cvoid}, Ptr{GMT_GRID_HEADER}, Cint, Cint, Cint), API, header, row, col, layer)
end

function GMT_Set_Index(API::Ptr{Cvoid}, header::Ptr{GMT_GRID_HEADER}, code)
	ccall((:GMT_Set_Index, libgmt), Cint, (Ptr{Cvoid}, Ptr{GMT_GRID_HEADER}, Ptr{UInt8}), API, header, code)
end =#

function GMT_Alloc_Segment(API::Ptr{Cvoid}, family::Integer, n_rows::Integer, n_columns::Integer, header, S::Ptr{Cvoid})
	(family != GMT_IS_DATASET && family != GMT_WITH_STRINGS) && error("Bad family type")

	ptr = ccall((:GMT_Alloc_Segment, libgmt), Ptr{Cvoid}, (Cstring, UInt32, UInt64, UInt64, Ptr{UInt8}, Ptr{Cvoid}),
							 API, family, n_rows, n_columns, header, S)

	convert(Ptr{GMT_DATASEGMENT}, ptr)
end

#=
function GMT_Set_Columns(API::Ptr{Cvoid}, n_columns::Integer, mode)
	ccall((:GMT_Set_Columns, libgmt), Cint, (Ptr{Cvoid}, UInt32, UInt32), API, n_columns, mode)
end

function GMT_Destroy_Group(API::Ptr{Cvoid}, obj::Ptr{Cvoid}, n_items::Integer)
	ccall((:GMT_Destroy_Group, libgmt), Cint, (Ptr{Cvoid}, Ptr{Cvoid}, Cuint), API, obj, n_items)
end =#

function GMT_Change_Layout(API::Ptr{Cvoid}, family::Integer, code::String, mode::Integer, obj, data=C_NULL, alpha=C_NULL)
	ccall((:GMT_Change_Layout, libgmt), Cint, (Ptr{Cvoid}, Cuint, Ptr{UInt8}, Cuint, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}),
	      API, family, code, mode, obj, data, alpha)
end

function GMT_Duplicate_String(API::Ptr{Cvoid}, str)
	ccall((:GMT_Duplicate_String, libgmt), Ptr{UInt8}, (Cstring, Ptr{UInt8}), API, str)
end

function GMT_Open_VirtualFile(API::Ptr{Cvoid}, family::Integer, geometry::Integer, dir::Integer, data, name)
	if (GMTver >= v"6.1")  dir |= GMT_Get_Enum(API, "GMT_IS_REFERENCE")  end
	ccall((:GMT_Open_VirtualFile, libgmt), Cint, (Cstring, UInt32, UInt32, UInt32, Ptr{Cvoid}, Ptr{UInt8}), API, family, geometry, dir, data, name)
end
function GMT_Close_VirtualFile(API::Ptr{Cvoid}, str)
	ccall((:GMT_Close_VirtualFile, libgmt), Cint, (Cstring, Ptr{UInt8}), API, str)
end

function GMT_Read_VirtualFile(API::Ptr{Cvoid}, str)
	ccall((:GMT_Read_VirtualFile, libgmt), Ptr{Cvoid}, (Cstring, Ptr{UInt8}), API, str)
end

#=
function GMT_Init_VirtualFile(API::Ptr{Cvoid}, mode::Integer, name)
	ccall((:GMT_Init_VirtualFile, libgmt), Ptr{Cvoid}, (Ptr{Cvoid}, UInt32, Ptr{UInt8}), API, mode, name)
end

function GMT_Put_Vector(API::Ptr{Cvoid}, V::GMT_VECTOR, col::Integer, tipo::Integer, vector)
	ccall((:GMT_Put_Vector, libgmt), Cint, (Ptr{Cvoid}, GMT_VECTOR, UInt32, UInt32, Ptr{Cvoid}), API, V, col, tipo, vector)
end
function GMT_Get_Vector(API::Ptr{Cvoid}, V::Ptr{GMT_VECTOR}, col::Integer)
	ccall((:GMT_Get_Vector, libgmt), Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{GMT_VECTOR}, UInt32), API, V, col)
end

function GMT_Put_Matrix(API::Ptr{Cvoid}, M::GMT_MATRIX, tipo::Integer, matrix)
	ccall((:GMT_Put_Matrix, libgmt), Cint, (Ptr{Cvoid}, GMT_MATRIX, UInt32, Ptr{Cvoid}), API, M, tipo, matrix)
end
function GMT_Get_Matrix(API::Ptr{Cvoid}, M::Ptr{GMT_MATRIX})
	ccall((:GMT_Get_Matrix, libgmt), Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{GMT_MATRIX}), API, M)
end

function GMT_Error_Message(API::Ptr{Cvoid})
	ccall((:GMT_Error_Message, libgmt), Ptr{UInt8}, (Ptr{Cvoid},), API)
end =#

# -------------------------------------------------------------------------------------------------------------
function GMT_Set_AllocMode(API::Ptr{Cvoid}, family::Integer, object)
	ccall((:GMT_Set_AllocMode, libgmt), Cint, (Cstring, UInt32, Ptr{Cvoid}), API, family, object)
end

function gmt_manage_workflow(API::Ptr{Cvoid}, mode::Int, texto)
	ccall((:gmt_manage_workflow, libgmt), Cint, (Cstring, UInt32, Cstring), API, mode, texto)
end

function GMT_Get_Enum(API::Ptr{Cvoid}, enum_name::String)
	ccall((:GMT_Get_Enum, libgmt), Cint, (Cstring, Ptr{UInt8}), API, enum_name)
end

function GMT_Get_Version()
	ver = ccall((:GMT_Get_Version, libgmt), Cfloat, (Ptr{Cvoid}, Ptr{Cuint}, Ptr{Cuint}, Ptr{Cuint}), C_NULL, C_NULL, C_NULL, C_NULL)
end
function GMT_Get_Version(major, minor, patch)
	ver = ccall((:GMT_Get_Version, libgmt), Cfloat, (Ptr{Cvoid}, Ptr{Cuint}, Ptr{Cuint}, Ptr{Cuint}), C_NULL, major, minor, patch)
end

function GMT_Get_Ctrl(API::Ptr{Cvoid})
	#if (GMTver > v"6.0")  ccall((:gmtlib_get_ctrl, libgmt), Ptr{Cvoid}, (Cstring,), API)
	#else                  ccall((:GMT_Get_Ctrl, libgmt), Ptr{Cvoid}, (Cstring,), API)
	#end
	(GMTver > v"6.0") ? ccall((:gmtlib_get_ctrl, libgmt), Ptr{Cvoid}, (Cstring,), API) :
	                    ccall((:GMT_Get_Ctrl, libgmt), Ptr{Cvoid}, (Cstring,), API)
end

#=
function gmt_getpen(API::Ptr{Cvoid}, buffer, P)
	GMT_ = GMT_Get_Ctrl(API)
	ccall((:gmt_getpen, libgmt), Cint, (Cstring, Ptr{Cuint}, Ref{GMT_PEN}), GMT_, buffer, P)
end =#

function gmtlib_setparameter(API, keyword::String, value::String)
	(!isa(API, Ptr{Nothing}) || API == C_NULL) && return UInt32(1)
	ccall((:gmtlib_setparameter, libgmt), Cuint, (Cstring, Ptr{UInt8}, Ptr{UInt8}, Bool), GMT_Get_Ctrl(API), keyword, value, true)
end

function reset_defaults(API::Ptr{Cvoid})
	(GMTver > v"6.1.1") ? ccall((:gmt_conf_SI, libgmt), Cvoid, (Cstring,), GMT_Get_Ctrl(API)) :
	                      ccall((:gmt_conf, libgmt), Cvoid, (Cstring,), GMT_Get_Ctrl(API))
end

function gmt_ogrread(API::Ptr{Cvoid}, fname::String, region=C_NULL)
	GMT_ = GMT_Get_Ctrl(API)
	if (GMTver == v"6.0")  ccall((:gmt_ogrread, libgmt), Ptr{OGR_FEATURES}, (Cstring, Ptr{UInt8}), GMT_, fname)
	else                   ccall((:gmt_ogrread, libgmt), Ptr{OGR_FEATURES}, (Cstring, Ptr{UInt8}, Ptr{Cdouble}), GMT_, fname, region)
	end
end

function gmt_ogrread(API::Ptr{Cvoid}, X)
	ccall((:gmt_ogrread2, libgmt), Ptr{OGR_FEATURES}, (Cstring, Ptr{Cvoid}), GMT_Get_Ctrl(API), X)
end

#=
function gmt_put_history(API::Ptr{Cvoid})
	ccall((:gmt_put_history, libgmt), Cint, (Cstring,), GMT_Get_Ctrl(API))
end =#

function GMT_Put_Strings(API::Ptr{Cvoid}, family::Integer, object::Ptr{Cvoid}, txt::Vector{String})
	ccall((:GMT_Put_Strings, libgmt), Cint, (Cstring, UInt32, Cstring, Ptr{Ptr{UInt8}}), API, family, object, txt)
end

function gmt_get_rgb_from_z(API::Ptr{Cvoid}, P::Ptr{GMT.GMT_PALETTE}, value::Cdouble, rgb::Vector{Float64})
	GMT_ = GMT_Get_Ctrl(API)
	ccall((:gmt_get_rgb_from_z, libgmt), Cint, (Cstring, Ptr{Cvoid}, Cdouble, Ptr{Cdouble}), GMT_, P, value, rgb)
end

function gmt_free_mem(API::Ptr{Cvoid}, mem)
	GMT_ = GMT_Get_Ctrl(API)
	ccall((:gmt_free_func, libgmt), Cvoid, (Cstring, Ptr{Cvoid}, Bool, Cstring), GMT_, mem, true, "Julia")
end

function sprintf(format::String, x...)
	strp = Ref{Ptr{Cchar}}(0)
	if (length(x) == 1)
		len = ccall(:asprintf, Cint, (Ptr{Ptr{Cchar}}, Cstring, Cdouble...), strp, format, x[1])
	elseif (length(x) == 2)
		len = ccall(:asprintf, Cint, (Ptr{Ptr{Cchar}}, Cstring, Cdouble, Cdouble...), strp, format, x[1], x[2])
	elseif (length(x) == 4)
		len = ccall(:asprintf, Cint, (Ptr{Ptr{Cchar}}, Cstring, Cdouble, Cdouble, Cdouble, Cdouble...), strp, format, x[1], x[2], x[3], x[4])
	end
	str = unsafe_string(strp[],len)
	Libc.free(strp[])
	return str
end

#=
function get_common_R(API::Ptr{Cvoid})
	R = COMMON_R((false,false,false,false), false, 0, 0, 0, (0., 0., 0., 0., 0., 0.), (0., 0., 0., 0.), (0., 0.), map(UInt8, (string(repeat(" ",256))...,)))
	Rp = pointer([R])
	ccall((:gmtlib_get_common_R, libgmt), Cint, (Cstring, Ptr{COMMON_R}), API, Rp)
	return unsafe_load(Rp)
end
=#

terrain_filter(data, detail, nrows, ncols, xinc, yinc, coord_type, center_lat=0.0, progress=C_NULL) =
	ccall((:terrain_filter, libgmt), Cint, (Ptr{Cfloat}, Cdouble, Cint, Cint, Cdouble, Cdouble, Cint, Cdouble, Ptr{Cvoid}), data, detail, nrows, ncols, xinc, yinc, coord_type, center_lat, progress)

terrain_image_data(data, contrast, nrows, ncols, image_min=0., image_max=65535.0) =
	ccall((:terrain_image_data, libgmt), Cint, (Ptr{Cfloat}, Cint, Cint, Cdouble, Cdouble, Cdouble), data, nrows, ncols, contrast, image_min, image_max)

fix_mercator(data, detail, nrows, ncols, lat1, lat2) =
	ccall((:fix_mercator, libgmt), Cvoid, (Ptr{Cfloat}, Cdouble, Cint, Cint, Cdouble, Cdouble), data, nrows, ncols, detail, lat1, lat2)

fix_polar_stereographic(data, detail, nrows, ncols, center_res) =
	ccall((:fix_polar_stereographic, libgmt), Cvoid, (Ptr{Cfloat}, Cdouble, Cint, Cint, Cdouble), data, nrows, ncols, detail, center_res)


text:
module perftools
  const _ProtoBuf_Top_ = @static isdefined(parentmodule(@__MODULE__), :_ProtoBuf_Top_) ? (parentmodule(@__MODULE__))._ProtoBuf_Top_ : parentmodule(@__MODULE__)
  module profiles
    const _ProtoBuf_Top_ = @static isdefined(parentmodule(@__MODULE__), :_ProtoBuf_Top_) ? (parentmodule(@__MODULE__))._ProtoBuf_Top_ : parentmodule(@__MODULE__)
    include("profile_pb.jl")
  end
end


text:
#__precompile__()

module ExperimentalDesign

using Random
using Logging
using Primes
using DataFrames
using Distributions
using StatsModels
using DocStringExtensions
using LinearAlgebra

import Random: rand

# Types

export    AbstractDesign,   AbstractScreeningDesign,    AbstractFactorialDesign,
    AbstractOptimalDesign,  PlackettBurman, FullFactorial,  FractionalFactorial,
    OptimalDesign, DesignDistribution, CategoricalFactor

# Methods

export fullfactorial,  explicit_fullfactorial, plackettburman, fold!, isplackettburman,
    paley, rand, random_design!, next_offset_divisible_prime, kl_exchange, d_criterion

# Pretty printing

export show

include("design.jl")
include("plackett_burman.jl")
include("factorial.jl")
include("random.jl")
include("kl_exchange.jl")
include("custom_show.jl")
include("categorical.jl")

end


text:
@testset "nystrom" begin
    dims = [10,5]
    X = rand(dims...)
    k = SqExponentialKernel()
    for obsdim in [1, 2]
        @test kernelmatrix(k, X; obsdim=obsdim) ≈ kernelmatrix(nystrom(k, X, 1.0; obsdim=obsdim))
        @test kernelmatrix(k, X; obsdim=obsdim) ≈ kernelmatrix(nystrom(k, X, collect(1:dims[obsdim]); obsdim=obsdim))
    end
end


text:
# Parts of this code were taken / derived from Graphs.jl. See LICENSE for
# licensing details.

# Maximum adjacency visit / traversal

#################################################
#
#  Maximum adjacency visit
#
#################################################

struct MaximumAdjacency <: SimpleGraphVisitAlgorithm end

abstract type AbstractMASVisitor <: SimpleGraphVisitor end

function maximum_adjacency_visit_impl!(
    g::AGraphOrDiGraph,	                      # the graph
    pq::PriorityQueue{Int, T},               # priority queue
    visitor::AbstractMASVisitor,                      # the visitor
    colormap::Vector{Int}) where T                            # traversal status

    while !isempty(pq)
        u = dequeue!(pq)
        discover_vertex!(visitor, u)
        for v in out_neighbors(g, u)
            examine_neighbor!(visitor, u, v, 0, 0, 0)

            if haskey(pq,v)
                ed = visitor.distmx[u, v]
                pq[v] += ed
            end
        end
        close_vertex!(visitor, u)
    end

end

function traverse_graph!(
    g::AGraphOrDiGraph,
    T::DataType,
    alg::MaximumAdjacency,
    s,
    visitor::AbstractMASVisitor,
    colormap::Vector{Int})

    pq = PriorityQueue{Int,T,typeof(Base.Order.Reverse)}(Base.Order.Reverse)

    # Set number of visited neighbors for all vertices to 0
    for v in vertices(g)
        pq[v] = zero(T)
    end

    @assert haskey(pq,s)
    @assert nv(g) >= 2

    #Give the starting vertex high priority
    pq[s] = one(T)

    #start traversing the graph
    maximum_adjacency_visit_impl!(g, pq, visitor, colormap)
end

#################################################
#
#  Visitors
#
#################################################

#################################################
#
#  Minimum Cut Visitor
#
#################################################

mutable struct MinCutVisitor{G,T,EM} <: AbstractMASVisitor
    g::G
    parities::BitVector
    colormap::Vector{Int}
    bestweight::T
    cutweight::T
    visited::Int
    distmx::EM
    vertices::Vector{Int}
end

function MinCutVisitor(g::AGraphOrDiGraph, distmx::AEdgeMap)
    n = nv(g)
    parities = falses(n)
    T = valtype(distmx)
    return MinCutVisitor(
        g,
        falses(n),
        zeros(Int,n),
        typemax(T),
        zero(T),
        0,
        distmx,
        Vector{Int}()
    )
end

function discover_vertex!(vis::MinCutVisitor, v)
    vis.parities[v] = false
    vis.colormap[v] = 1
    push!(vis.vertices,v)
    return true
end

function examine_neighbor!(vis::MinCutVisitor, u, v, ucolor, vcolor, ecolor)
    ew = vis.distmx[u, v]

    # if the target of e is already marked then decrease cutweight
    # otherwise, increase it

    if vis.colormap[v] != vcolor # here vcolor is 0
        vis.cutweight -= ew
    else
        vis.cutweight += ew
    end
    return true
end

function close_vertex!(vis::MinCutVisitor, v)
    vis.colormap[v] = 2
    vis.visited += 1

    if vis.cutweight < vis.bestweight && vis.visited < nv(vis.g)
        vis.bestweight = vis.cutweight
        for u in vertices(vis.g)
            vis.parities[u] = ( vis.colormap[u] == 2)
        end
    end
    return true
end

#################################################
#
#  MAS Visitor
#
#################################################

mutable struct MASVisitor{EM,I<:IO} <: AbstractMASVisitor
    io::I
    vertices::Vector{Int}
    distmx::EM
    log::Bool
end

function discover_vertex!(visitor::MASVisitor{T}, v) where T
    push!(visitor.vertices,v)
    visitor.log ? println(visitor.io, "discover vertex: $v") : nothing
    return true
end

function examine_neighbor!(visitor::MASVisitor, u, v, ucolor, vcolor, ecolor)
    visitor.log ? println(visitor.io, " -- examine neighbor from $u to $v") : nothing
    return true
end

function close_vertex!(visitor::MASVisitor, v)
    visitor.log ? println(visitor.io, "close vertex: $v") : nothing
    return true
end

#################################################
#
#  Minimum Cut
#
#################################################

"""
    minimum_cut(g, dist_edge map=weights(g))

Finds the `cut` of minimum total weight.

Returns a tuple `(f, cut, labels)`, where `f` is the weight of the cut,
`cut` is a vector of the edges in the cut, and `labels` gives a partitioning
of the vertices in two sets, according to the cut.
An optional `dist_matrix` edge map maybe specified; if omitted, edge distances are assumed to be 1.
"""
function minimum_cut(
        g::AGraphOrDiGraph,
        distmx::AEdgeMap
    )
    T = valtype(distmx)
    visitor = MinCutVisitor(g, distmx)
    colormap = zeros(Int, nv(g))
    traverse_graph!(g, T, MaximumAdjacency(), 1, visitor, colormap)
    labels = visitor.parities + 1
    E = edgetype(g)
    cut = Vector{E}()
    for e in edges(g)
        if labels[src(e)] != labels[dst(e)]
            push!(cut, e)
        end
    end
    return visitor.bestweight, cut, labels
end

minimum_cut(g::AGraphOrDiGraph) = minimum_cut(g, weights(g))

"""
    maximum_adjacency_visit(
        g,
        distmx::AEdgeMap=weights(g);
        log::Bool=false,
        io::IO=STDOUT
    )

Returns the vertices in `g` traversed by maximum adjacency search. An optional
`distmx` edge map may be specified; if omitted, edge distances are assumed to
be 1. If `log` (default `false`) is `true`, visitor events will be printed to
`io`, which defaults to `STDOUT`; otherwise, no event information will be
displayed.
"""
function maximum_adjacency_visit(
        g::AGraphOrDiGraph,
        distmx::AEdgeMap=weights(g);
        log::Bool=false,
        io::IO=stdout
    )
    visitor = MASVisitor(io, Vector{Int}(), distmx, log)
    traverse_graph!(g, valtype(distmx), MaximumAdjacency(), 1, visitor, zeros(Int, nv(g)))
    return visitor.vertices
end


text:
# # Hovercraft Path Planning
# In this case study, we seek to determine an optimal control policy for the
# trajectory of a hovercraft that travels to a set of dynamic waypoints while
# trying to minimize the thrust input.

# ## Formulation

# The corresponding dynamic optimization problem is expressed:
# ```math
# \begin{aligned}
# 	&&\underset{x(t), v(t), u(t)}{\text{min}} &&& \int_{t \in T} |u(t)|_2^2 dt  \\
# 	&&\text{s.t.} &&& v(0) = v0\\
# 	&&&&& \frac{dx}{dt} = v(t), && t \in T\\
#     &&&&& \frac{dv}{dt} = u(t), && t \in T\\
#     &&&&& x(t_i) = xw_i, && i \in I
# \end{aligned}
# ```
# where ``x(t)`` is the Cartesian position, ``v(t)`` is the velocity, ``u(t)`` is
# the thrust input, ``xw_i, \ i \in I,`` are the waypoints, and ``T`` is the time
# horizon.

# ## Model Definition

# Let's implement this in `InfiniteOpt` and first import the packages we need:
using InfiniteOpt, Ipopt

# Next we'll specify our waypoint data:
xw = [1 4 6 1; 1 3 0 1] # positions
tw = [0, 25, 50, 60];    # times

# We initialize the infinite model and opt to use the Ipopt solver:
m = InfiniteModel(optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0));

# Let's specify our infinite parameter which is time ``t \in [0, 60]``:
@infinite_parameter(m, t in [0, 60], num_supports = 61)

# Now let's specify the decision variables:
@variables(m, begin
    ## state variables
    x[1:2], Infinite(t)
    v[1:2], Infinite(t)
    ## control variables
    u[1:2], Infinite(t), (start = 0)
end)

# Specify the objective:
@objective(m, Min, ∫(u[1]^2 + u[2]^2, t))

# Set the initial conditions with respect to the velocity:
@constraint(m, [i = 1:2], v[i] == 0, DomainRestrictions(t => 0))

# Define the point physics ODEs which serve as our system model:
@constraint(m, [i = 1:2], ∂(x[i], t) == v[i])
@constraint(m, [i = 1:2], ∂(v[i], t) == u[i])

# Ensure we hit all the waypoints:
@constraint(m, [i = 1:2, j = eachindex(tw)], x[i] == xw[i, j], DomainRestrictions(t => tw[j]))

# ## Problem Solution

# Optimize the model:
optimize!(m)

# Extract the results:
x_opt = value.(x);

# Plot the results:
using Plots
scatter(xw[1,:], xw[2,:], label = "Waypoints", background_color = :transparent)
plot!(x_opt[1], x_opt[2], label = "Trajectory")
xlabel!("x_1")
ylabel!("x_2")

# That's it, now we have our optimal trajectory!

# ### Maintenance Tests
# These are here to ensure this example stays up to date.
using Test
@test termination_status(m) == MOI.LOCALLY_SOLVED
@test has_values(m)
@test x_opt isa Vector{<:Vector{<:Real}}


text:
abstract type AbstractVectorQuantity{T,N} <: AbstractArray{T,N} end

struct VectorQuantity end

struct VectorField{N,M,T,D<:AbstractArray{T,M},G<:AbstractAxisGrid} <: AbstractPICDataStructure{T,N,G}
    data::D
    grid::G
    name::String
end

struct VectorVariable{T,D<:AbstractVector{T},G<:ParticlePositions} <: AbstractPICDataStructure{T,1,G}
    data::D
    grid::G
    name::String
end

@inline scalarness(::Type{<:VectorField}) = VectorQuantity()
@inline scalarness(::Type{<:VectorVariable}) = VectorQuantity()

function VectorField{N}(data::D, grid::G, name) where {N, M, T, D <: AbstractArray{T,M}, G}
    VectorField{N, M, T, D, G}(data, grid, name)
end

function VectorField(data::AbstractArray{T,M}, grid::G, name) where {N, M, P, G, S<:NTuple{N}, T<:NamedTuple{P,S}}
    @debug "Array of NamedTuple"
    VectorField{N}(data, grid, name)
end

VectorField(data, grid; name="") = VectorField(data, grid, name)

function VectorField(row_data::AbstractArray{T,M}, grid::G; name="") where {N,M,T<:SVector{N},G}
    component_names = Val(staticnames(G))
    VectorField(row_data, grid, component_names, name)
end

function VectorField(row_data::AbstractArray{T,M}, grid::G, ::Val{names}, name="") where {N, M, G, names, T<:SVector{N}}
    # @debug "Building $(N)D VectorField row-wise from $T"
    ElType = NamedTuple{names, NTuple{N,recursive_bottom_eltype(row_data)}}
    # @debug "Got ElType $ElType"
    data = StructArray{ElType}(undef, size(row_data))
    fillvec!(data, row_data)
    VectorField{N}(data, grid, name)
end

@inbounds function fillvec!(target, data)
    for i in eachindex(target, data)
        target[i] = vector2nt(target, data[i])
    end
end

function VectorVariable(data::D, grid::G, name) where {M, T, D <: AbstractArray{T,M}, G}
    VectorVariable{M, T, D, G}(data, grid, name)
end

VectorVariable(data, grid; name="") = VectorVariable(data, grid, name)

@inline function vector2nt(f, v::SArray{Tuple{N},T}) where {N,T}
    # @debug "Data type $(typeof(f.data)) and $(typeof(v))"
    names = propertynames(f)
    NamedTuple{names, NTuple{N,T}}(v)
end

# Broadcasting

Base.BroadcastStyle(::VectorQuantity, ::Type{<:AbstractPICDataStructure}) = Broadcast.ArrayStyle{AbstractVectorQuantity}()

function similar_data(data::StructArray{T}, ElType, dims) where T
    @debug "Building similar StructArray with type $(typeof(data)) end eltype $ElType"

    props = propertynames(data)
    N = length(props)
    eltypes = NTuple{N, eltype(ElType)}
    @debug "Element eltypes: $eltypes"
    output_eltype = NamedTuple{props, eltypes}
    similar(similar_structarray(data, output_eltype), dims)
end

function Base.similar(bc::Broadcasted{ArrayStyle{AbstractVectorQuantity}}, ::Type{ElType}) where ElType
    # Scan the inputs for the AbstractPICDataStructure:
    f = find_field(bc)
    @debug "Building datastructure similar to $(typeof(f)) with eltype $ElType"
    grid = getdomain(f)
    # Keep the same grid for the output
    data = similar_data(unwrapdata(f), ElType, axes(bc))
    @debug "Output data type: $(typeof(data))"
    newstruct(f, data, grid)
end

function Base.similar(bc::Broadcasted{ArrayStyle{AbstractVectorQuantity}}, ::Type{ElType}) where ElType <: Number
    # Scan the inputs for the AbstractPICDataStructure:
    f = find_field(bc)
    @debug "Building datastructure similar to $(typeof(f)) with eltype $ElType"
    grid = getdomain(f)
    # Keep the same grid for the output
    data = similar(unwrapdata(f), ElType, axes(bc))
    @debug "Output data type: $(typeof(data))"
    scalar_from(typeof(f))(data, grid)
end

# Indexing
@propagate_inbounds function Base.getindex(::VectorQuantity, v::T, i) where T
    N = dimensionality(eltype(T))
    SVector{N}(unwrapdata(v)[i]...)
end
@propagate_inbounds function Base.setindex!(::VectorQuantity, f, v::SVector, i::Int)
    # @debug typeof(f.data)
    unwrapdata(f)[i] = vector2nt(f, v)
end

function Base.eltype(::VectorQuantity, v::Type{T}) where T
    SVector{dimensionality(T),recursive_bottom_eltype(T)}
end

function Base.similar(::VectorQuantity, f, ::Type{S}, dims::Dims) where S
    # @debug "Building similar VectorQuantity of type $S"
    data = StructArray(similar(unwrapdata(f), S, dims))
    grid = getdomain(f)
    newstruct(f, data, grid)
end

# Acessing the internal data storage by column names
function get_componenet(f::T, key) where T
    data = getproperty(unwrapdata(f), key)
    grid = getdomain(f)

    scalar_from(T)(data, grid)
end

Base.getproperty(::VectorQuantity, f::AbstractPICDataStructure, key::Symbol) = get_componenet(f, key)
Base.propertynames(::VectorQuantity, f::AbstractPICDataStructure) = propertynames(unwrapdata(f))

vector_from(::Type{<:ScalarField}) = VectorField
vector_from(::Type{<:ScalarVariable}) = VectorVariable
scalar_from(::Type{<:VectorField}) = ScalarField
scalar_from(::Type{<:VectorVariable}) = ScalarVariable

# build_vector must be called with at least one component. NTuple{N,T} can have N==0
function build_vector(components::Tuple{T, Vararg{T,N}}, names::Tuple{Symbol, Vararg{Symbol,N}}; name="") where {N, T}
    # We cannot have StructArrays with ScalarField components because we cannot
    # correctly build them through similar because StructArrays doens't have a
    # similar(::StructArray, ElType)
    # We fake this by creating the field in getproperty
    data_components = NamedTuple(name=>unwrapdata(c) for (name,c) in zip(names, components))
    data = StructArray(data_components)
    x = first(components)

    for c in components
        if c.grid ≠ x.grid
            @warn "Grids for vector variable may not be compatible"
        end
    end

    vectortype = vector_from(T)
    vectortype(data, x.grid; name)
end


text:
"""
    AsynchronousExecution

Functor object for using asynchronous execution in a progressive-hedging algorithm (assuming multiple Julia cores are available). Create by supplying an [`Asynchronous`](@ref) object through `execution` in the `ProgressiveHedgingSolver` factory function and then pass to a `StochasticPrograms.jl` model.

"""
struct AsynchronousExecution{T <: AbstractFloat,
                             A <: AbstractVector,
                             PT <: AbstractPenaltyterm} <: AbstractProgressiveHedgingExecution
    subworkers::Vector{SubWorker{T,A,PT}}
    work::Vector{Work}
    finalize::Vector{Work}
    progressqueue::ProgressQueue{T}
    x̄::Vector{RemoteRunningAverage{A}}
    δ::Vector{RemoteRunningAverage{T}}
    iterates::RemoteIterates{A}
    r::IteratedValue{T}
    active_workers::Vector{Future}
    # Bookkeeping
    subobjectives::Vector{A}
    finished::Vector{Int}
    # Parameters
    κ::T

    function AsynchronousExecution(κ::T,
                                   ::Type{T}, ::Type{A}, ::Type{PT}) where {T <: AbstractFloat,
                                                                            A <: AbstractVector,
                                                                            PT <: AbstractPenaltyterm}
        return new{T,A,PT}(Vector{SubWorker{T,A,PT}}(undef, nworkers()),
                           Vector{Work}(undef, nworkers()),
                           Vector{Work}(undef, nworkers()),
                           RemoteChannel(() -> Channel{Progress{T}}(4 * nworkers())),
                           Vector{RemoteRunningAverage{A}}(undef, nworkers()),
                           Vector{RemoteRunningAverage{T}}(undef, nworkers()),
                           RemoteChannel(() -> IterationChannel(Dict{Int,A}())),
                           RemoteChannel(() -> IterationChannel(Dict{Int,T}())),
                           Vector{Future}(undef, nworkers()),
                           Vector{A}(),
                           Vector{Int}(),
                           κ)
    end
end

function initialize_subproblems!(ph::AbstractProgressiveHedging,
                                 execution::AsynchronousExecution{T,A},
                                 scenarioproblems::DistributedScenarioProblems,
                                 penaltyterm::AbstractPenaltyterm) where {T <: AbstractFloat, A <: AbstractVector}
    # Create subproblems on worker processes
    initialize_subproblems!(ph,
                            execution.subworkers,
                            scenarioproblems,
                            penaltyterm)
    # Continue preparation
    @sync begin
        for w in workers()
            execution.work[w-1] = RemoteChannel(() -> Channel{Int}(round(Int,10/execution.κ)), w)
            execution.finalize[w-1] = RemoteChannel(() -> Channel{Int}(1), w)
            execution.x̄[w-1] = RemoteChannel(() -> Channel{RunningAverage{A}}(1), w)
            execution.δ[w-1] = RemoteChannel(() -> Channel{RunningAverage{T}}(1), w)
            @async remotecall_fetch(
                w,
                execution.subworkers[w-1],
                execution.x̄[w-1],
                length(ph.ξ)) do sw, average, xdim
                    subproblems = fetch(sw)
                    x̄ = mapreduce(+, subproblems, init = zeros(T, xdim)) do subproblem
                        x = subproblem.x
                        π = subproblem.probability
                        return π * x
                    end
                    running_average = RunningAverage(x̄, [s.x for s in subproblems])
                    put!(average, running_average)
                    return x̄
                end
            @async remotecall_fetch(
                w,
                execution.subworkers[w-1],
                execution.δ[w-1],
                length(ph.ξ)) do sw, average, xdim
                    running_average =
                        RunningAverage(zero(T),
                                       fill(zero(T), length(fetch(sw))))
                    put!(average, running_average)
                end
            put!(execution.work[w-1], 1)
        end
        # Prepare memory
        push!(execution.subobjectives, zeros(num_subproblems(ph)))
        push!(execution.finished, 0)
        log_val = ph.parameters.log
        ph.parameters.log = false
        log!(ph)
        ph.parameters.log = log_val
    end
    # Initial reductions
    update_iterate!(ph)
    # Init δ₂
    @sync begin
        for w in workers()
            @async remotecall_fetch(
                w,
                execution.subworkers[w-1],
                ph.ξ,
                execution.δ[w-1]) do sw, ξ, δ
                    for (i, subproblem) in enumerate(fetch(sw))
                        subtract!(fetch(δ), i)
                        x = subproblem.x
                        π = subproblem.probability
                        add!(fetch(δ), i, norm(x - ξ, 2) ^ 2, π)
                    end
                end
        end
    end
    update_dual_gap!(ph)
    return nothing
end

function iterate!(ph::AbstractProgressiveHedging, execution::AsynchronousExecution{T}) where T <: AbstractFloat
    wait(execution.progressqueue)
    while isready(execution.progressqueue)
        # Add new cuts from subworkers
        t::Int, i::Int, Q::SubproblemSolution{T} = take!(execution.progressqueue)
        if !(Q.status ∈ AcceptableTermination)
            # Early termination log
            log!(ph; status = Q.status)
            return Q.status
        end
        execution.subobjectives[t][i] = Q.value
        execution.finished[t] += 1
        if execution.finished[t] == num_subproblems(ph)
            # Update objective
            ph.Q_history[t] = current_objective_value(ph, execution.subobjectives[t])
            ph.data.Q = ph.Q_history[t]
        end
    end
    # Project and generate new iterate
    t = ph.data.iterations
    if execution.finished[t] >= execution.κ * num_subproblems(ph)
        # Get dual gap
        update_dual_gap!(ph)
        # Update progress
        @unpack δ₁, δ₂ = ph.data
        ph.primal_gaps[t] = δ₁
        ph.dual_gaps[t] = δ₂
        # Check if optimal
        if check_optimality(ph)
            # Optimal, final log
            log!(ph, optimal = true)
            return MOI.OPTIMAL
        end
        # Update penalty (if applicable)
        update_penalty!(ph)
        # Update iterate
        update_iterate!(ph)
        # Send new work to workers
        put!(execution.iterates, t + 1, ph.ξ)
        put!(execution.r, t + 1, penalty(ph))
        map((w, aw) -> !isready(aw) && put!(w, t + 1), execution.work, execution.active_workers)
        # Prepare memory for next iteration
        push!(execution.subobjectives, zeros(num_subproblems(ph)))
        push!(execution.finished, 0)
        # Log progress
        log!(ph)
    end
    # Dont return a status as procedure should continue
    return nothing
end

function start_workers!(ph::AbstractProgressiveHedging, execution::AsynchronousExecution)
    # Load initial decision
    put!(execution.iterates, 1, ph.ξ)
    put!(execution.r, 1, penalty(ph))
    for w in workers()
        execution.active_workers[w-1] = remotecall(work_on_subproblems!,
                                                   w,
                                                   execution.subworkers[w-1],
                                                   execution.work[w-1],
                                                   execution.finalize[w-1],
                                                   execution.progressqueue,
                                                   execution.x̄[w-1],
                                                   execution.δ[w-1],
                                                   execution.iterates,
                                                   execution.r)
    end
    return nothing
end

function close_workers!(ph::AbstractProgressiveHedging, execution::AsynchronousExecution)
    t = ph.data.iterations-1
    map((w, aw)->!isready(aw) && put!(w, t), execution.finalize, execution.active_workers)
    map((w, aw)->!isready(aw) && put!(w, -1), execution.work, execution.active_workers)
    map(wait, execution.active_workers)
end

function restore_subproblems!(::AbstractProgressiveHedging, execution::AsynchronousExecution)
    restore_subproblems!(execution.subworkers)
    return nothing
end

function resolve_subproblems!(ph::AbstractProgressiveHedging, execution::AsynchronousExecution)
    return nothing
end

function update_iterate!(ph::AbstractProgressiveHedging, execution::AsynchronousExecution)
    ξ_prev = copy(ph.ξ)
    ph.ξ .= sum(average.(fetch.(execution.x̄)))
    # Update δ₁
    ph.data.δ₁ = norm(ph.ξ - ξ_prev, 2) ^ 2
    return nothing
end

function update_subproblems!(ph::AbstractProgressiveHedging, execution::AsynchronousExecution)
    return nothing
end

function update_dual_gap!(ph::AbstractProgressiveHedging, execution::AsynchronousExecution)
    ph.data.δ₂ = sum(average.(fetch.(execution.δ)))
    return nothing
end

function calculate_objective_value(ph::AbstractProgressiveHedging, execution::AsynchronousExecution)
    return calculate_objective_value(execution.subworkers)
end

# API
# ------------------------------------------------------------
function (execution::Asynchronous)(::Type{T}, ::Type{A}, ::Type{PT}) where {T <: AbstractFloat,
                                                                            A <: AbstractVector,
                                                                            PT <: AbstractPenaltyterm}
    return AsynchronousExecution(execution.κ, T, A, PT)
end

function str(::Asynchronous)
    return "Asynchronous "
end


text:
using Documenter, GoogleLanguageService

makedocs(;
    modules=[GoogleLanguageService],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/sivakon/grpc-julia-gcp-example/blob/{commit}{path}#L{line}",
    sitename="Google Language Service",
    authors="Sivaram Konanki",
    assets=String[],
)

deploydocs(;
    repo="github.com/sivakon/grpc-julia-gcp-example",
)

max_stars_repo_name: sivakon/grpc-julia-gcp-example

max_issues_repo_name: sivakon/grpc-julia-gcp-example

max_forks_repo_name: sivakon/grpc-julia-gcp-example


text:
@doc raw"""
conductance based channel, generic format

$$i = g\cdot \text{krule}(\text{var}) \cdot (V - E)$$

"""
mutable struct GenericIonChannel <: AbstractIonChannel #TODO: optimize the structure for better and easier use.
    name::String
    ion::Symbol
    g::Real

    krule::Function
    kvars::Vector{Kinetics}
end

function dof(ch::GenericIonChannel)
    _ans = zeros(Int, length(ch.kvars))
    for (idx, item) in enumerate(ch.kvars)
        if item._type == :evolving
            _ans[idx] = 1
        end
    end
    _ans
end

function update!(ch::GenericIonChannel)
    for item in ch.kvars
        if ~isnothing(item.Vhalf)
            if item._state == :activation
                item.infty = (V) -> (1 - item.zeta) / (1 + exp((item.Vhalf - V) / item.k)) + item.zeta
            elseif item._state == :inactivation
                item.infty = (V) -> (1 - item.zeta) / (1 + exp(-(item.Vhalf - V) / item.k)) + item.zeta
            else
                nothing
            end
        end

        if ~isnothing(item.Cbase)
            item.tau = (V) -> item.Cbase + item.Camp * exp(- ((item.Vmax - V) / item.sigma)^2)
        end
    end
    ch
end

function step(ch::GenericIonChannel; V::Real, var::Vector{T}, E::Real) where {T<:Real}
    _var_idx = 1
    derivative = zeros(size(var))

    current = ch.g * ch.krule(ch.kvars, var) * (V - E)
    for item in ch.kvars
        if item._type == :evolving
            derivative[_var_idx] = (item.infty(V) - var[_var_idx]) / item.tau(V)
            _var_idx += 1
        else
            nothing
        end
    end

    return (current, derivative)
end

function itr_kinetics(ch::GenericIonChannel)
    ch.kvars
end

#NOTE: not tested
function current(ch::GenericIonChannel; V::Vector{T}, var::Array{T, 2}, E::T) where {T<:Real}
    _var_idx = 1
    _krule = (var_item) -> ch.krule(ch.kvars, var_item)
    _var_reform = [Tuple(var[:, idx]) for idx =1:size(var, 2)]
    current = ch.g .* _krule.(_var_reform) .* (V .- E)
    return current
end


text:
# utilities

half(x::Float64) = x * 0.5
half(x::Float32) = x *0.5f0

nonneg(x::Float64) = x < 0.0 ? 0.0 : x  # this also properly takes care of NaN
nonneg(x::Float32) = x < 0.0f0 ? 0.0f0 : x

macro _checkdims(cond)
    quote
        ($cond) || throw(DimensionMismatch())
    end
end

function axpby!{T<:BlasReal}(a::T, x::StridedVector{T}, b::T, y::StridedVector{T})
    if b == zero(T)
        scale!(y, x, a)
    elseif b == one(T)
        axpy!(a, x, y)
    else
        n = length(x)
        @_checkdims length(y) == n
        @inbounds for i = 1:n
            y[i] = a * x[i] + b * y[i]
        end
    end
    y
end

gets(x::StridedVector, i::Int) = x[i]
gets(x::StridedMatrix, i::Int) = view(x,:,i)
gets{T}(x::StridedArray{T,3}, i::Int) = view(x,:,:,i)

shrink{T<:AbstractFloat}(x::T, t::T) = (x > t ? x - t : x < -t ? x + t : zero(T))
shrink{T<:AbstractFloat}(x::StridedVector{T}, t::T) = T[shrink(v, t) for v in x]

no_op(args...) = nothing


text:
# # Selecting boundary faces

# In this tutorial, we will learn
#
#    -  How to use some low-level functionality of the MeshCore library.
#    -  How to find the boundary of the mesh.
#    -  How to select pieces of the boundary surfaces.

# The tutorial will produce files for mesh visualization in the
# [Paraview](https://www.paraview.org/) format (VTK). One can display this
# information by loading the file with `paraview.exe`. When the tutorial is
# executed in `mybinder.org`, the graphics file needs to be downloaded to your
# desktop, and then visualized locally.

# As an example we generate a regular tetrahedral mesh from hardwired arrays for
# the locations of the nodes and the connectivities of the elements.
include("samplet4.jl")
using Main.samplet4: samplet4mesh
xyz, cc = samplet4mesh()
# The two arrays `xyz`, `cc`, define the locations of the vertices (`xyz`, one
# node per row), and  the numbers of the vertices connected by the tetrahedral
# elements(again, one element per row of the array `cc`).

# Here we show some low-level operations using the core mesh library. The
# locations of the vertices are stored in static arrays in order to achieve fast
# access and processing.
using StaticArrays

# First we find out how many space coordinate dimensions there are (there should be 3).
N, T = size(xyz, 2), eltype(xyz)
# Now we can create an attribute for the vertex  shape collection to represent
# the locations of the vertices.
using MeshCore: VecAttrib
locs =  VecAttrib([SVector{N, T}(xyz[i, :]) for i in 1:size(xyz, 1)])

# We can create two shape collections, one for the vertices, and one for the
# tetrahedral elements. We also attach the geometry attribute to the vertex
# shape collection.
using MeshCore: P1, T4, ShapeColl
vrts = ShapeColl(P1, length(locs))
vrts.attributes["geom"] = locs
tets = ShapeColl(T4, size(cc, 1))

# Finally we connect the two shaped collections in an incidence relation.
using MeshCore: IncRel
ir = IncRel(tets, vrts, cc)

# The incidence relation (which really represents what is meant by "mesh"), can
# be written out to a file for visualization.
using MeshSteward: vtkwrite
vtkwrite("samplet4-elements", ir)

# The boundary of the tetrahedral mesh may be derived by a topological operation
# that constructs the skeleton of the mesh and identifies faces that are on the
# boundary.
using MeshCore: ir_boundary
bir = ir_boundary(ir)

# Now we select the boundary faces which are inside a box. The box is defined by
# the six numbers, for each of the  three coordinates we give the lower and
# upper limit, and then the box is increased (inflated) in all directions by a
# small amount in order to be able to capture vertices whose location is not
# precisely in the original box. This is important because the original box is
# in fact of zero volume,  being of zero thickness in the Y direction (both the
# lower and upper bound in Y are 8.0).
using MeshSteward: eselect
@show el = eselect(bir; box = [0.0, 3.0, 8.0, 8.0, 0.0, 5.0], inflate = 0.01)

# A subset of the boundary within the box (the faces given by the list `el`) may
# now be extracted and written out for visualization into a file. Both the file
# exported above and this one can now be loaded into  visualization software
# (we assume `paraview` here, at any visualization software capable of reading
# VTK files will do), and the correctness of the selection may be verified
# visually.
using MeshCore: ir_subset
using MeshSteward: vtkwrite
vtkwrite("samplet4-boundary-y=8", ir_subset(bir, el))


text:
"""
	ControlAffineSystem(system::String)

Construct a control affine system from a string.
"""

function ControlAffineSystem(system::String)
    if system == "single_integrator2D"
        return single_integrator2D()
    elseif system == "single_integrator3D"
        return single_integrator3D()
    elseif system == "double_integrator1D"
        return double_integrator1D()
    elseif system == "double_integrator2D"
        return double_integrator2D()
    elseif system == "double_integrator3D"
        return double_integrator3D()
    elseif system == "unicycle_3D"
        return unicycle_3D()
    elseif system == "unicycle_4D"
        return unicycle_4D()
    end
end

"""
	single_integrator2D()

Construct a ControlAffineSystem with 2D single integrator dynamics.
"""
function single_integrator2D()
    n = 2
    m = 2
    f(x) = zeros(n)
    g(x) = Matrix(1.0I, m, m)

    return ControlAffineSystem(n, m, f, g)
end

"""
	single_integrator3D()

Construct a ControlAffineSystem with 3D single integrator dynamics.
"""
function single_integrator3D()
    n = 3
    m = 3
    f(x) = zeros(n)
    g(x) = Matrix(1.0I, m, m)

    return ControlAffineSystem(n, m, f, g)
end

"""
	double_integrator1D()

Construct a ControlAffineSystem with 1D double integrator dynamics.
"""
function double_integrator1D()
    n = 2
    m = 1
    f(x) = [x[2], 0.0]
    g(x) = [0.0, 1.0]

    return ControlAffineSystem(n, m, f, g)
end

"""
	double_integrator2D()

Construct a ControlAffineSystem with 2D double integrator dynamics.
"""
function double_integrator2D()
    n = 4
    m = 2
    f(x) = vcat(x[(m + 1):end], zeros(m))
    g(x) = vcat(zeros(m, m), Matrix(1.0I, m, m))

    return ControlAffineSystem(n, m, f, g)
end

"""
	double_integrator3D()

Construct a ControlAffineSystem with 3D double integrator dynamics.
"""
function double_integrator3D()
    n = 6
    m = 3
    f(x) = vcat(x[(m + 1):end], zeros(m))
    g(x) = vcat(zeros(m, m), Matrix(1.0I, m, m))

    return ControlAffineSystem(n, m, f, g)
end

"""
    unicycle_3D()

Construct a ControlAffineSystem with 3D unicycle dynamics.
"""
function unicycle_3D()
    n = 3
    m = 2
    f(x) = zeros(n)
    g(x) = [cos(x[3]) 0.0; sin(x[3]) 0.0; 0.0 1.0]

    return ControlAffineSystem(n, m, f, g)
end

"""
    unicycle_4D()

Construct a ControlAffineSystem with 4D unicycle dynamics.
"""
function unicycle_4D()
    n = 4
    m = 2
    f(x) = [x[4] * cos(x[3]), x[4] * sin(x[3]), 0.0, 0.0]
    g(x) = [0.0 0.0; 0.0 0.0; 1.0 0.0; 0.0 1.0]

    return ControlAffineSystem(n, m, f, g)
end


text:
module AtomIO
using AtomsBase
using PyCall
using Unitful

export load_system, save_system

function load_system(file::AbstractString)
    ase = pyimport_e("ase")
    ispynull(ase) && error("Install ASE to load data from exteral files")
    ase_atoms = pyimport("ase.io").read(file)

    T = Float64
    cell_julia = convert(Array, ase_atoms.cell)  # Array of arrays
    box = [[T(cell_julia[j][i]) * u"Å" for i = 1:3] for j = 1:3]

    atoms = map(ase_atoms) do at
        atnum = convert(Int, at.number)
        Atom(atnum, at.position * u"Å"; magnetic_moment=at.magmom)
    end
    bcs = [p ? Periodic() : DirichletZero() for p in ase_atoms.pbc]
    atomic_system(atoms, box, bcs)
end

function save_system(filename::AbstractString, system::AbstractSystem)
    error("Not yet implemented.")
end

end


text:
a = 1
c = -5
limit = 10

mf = SigmoidMF(a, c, limit)

# Evaluation tests
@assert Fuzzy.eval(mf, c) == 0.5
@assert Fuzzy.eval(mf, c - 5) < Fuzzy.eval(mf, c + 5)

# Mean finding tests
@assert Fuzzy.mean_at(mf, 1) > Fuzzy.mean_at(mf, 0.6) > Fuzzy.mean_at(mf, 0.3) > Fuzzy.mean_at(mf, 0)


text:
using ..CUDASupport
using ..TopOpt: @init_cuda
@init_cuda()
using ..GPUUtils
import ..TopOpt: whichdevice

get_f(problem, vars::CuArray) = f = zeros(typeof(vars), ndofs(problem.ch.dh))

function update_f!(
    f::CuVector{T}, fes, fixedload, dof_cells, black, white, penalty, vars, varind, xmin
) where {T}
    args = (
        f,
        fes,
        fixedload,
        dof_cells.offsets,
        dof_cells.values,
        black,
        white,
        penalty,
        vars,
        varind,
        xmin,
        length(f),
    )
    callkernel(dev, assemble_kernel1, args)
    return CUDAdrv.synchronize(ctx)
end

function assemble_kernel1(
    f,
    fes,
    fixedload,
    dof_cells_offsets,
    dof_cells_values,
    black,
    white,
    penalty,
    vars,
    varind,
    xmin,
    ndofs,
)
    dofidx = @thread_global_index()
    offset = @total_threads()

    while dofidx <= ndofs
        f[dofidx] = fixedload[dofidx]
        r = dof_cells_offsets[dofidx]:(dof_cells_offsets[dofidx + 1] - 1)
        for i in r
            cellidx, localidx = dof_cells_values[i]
            if black[cellidx]
                f[dofidx] += fes[cellidx][localidx]
            elseif white[cellidx]
                px = xmin
                f[dofidx] += px * fes[cellidx][localidx]
            else
                if PENALTY_BEFORE_INTERPOLATION
                    px = density(penalty(vars[varind[cellidx]]), xmin)
                else
                    px = penalty(density(vars[varind[cellidx]], xmin))
                end
                f[dofidx] += px * fes[cellidx][localidx]
            end
        end
        dofidx += offset
    end

    return nothing
end

whichdevice(p::StiffnessTopOptProblem) = whichdevice(p.ch)
whichdevice(ch::ConstraintHandler) = whichdevice(ch.dh)
whichdevice(dh::DofHandler) = whichdevice(dh.grid)
whichdevice(g::Ferrite.Grid) = whichdevice(g.cells)

@define_cu(
    ElementFEAInfo,
    :Kes,
    :fes,
    :fixedload,
    :cellvolumes,
    :metadata,
    :black,
    :white,
    :varind,
    :cells
)
@define_cu(TopOptProblems.Metadata, :cell_dofs, :dof_cells, :node_cells, :node_dofs)
@define_cu(Ferrite.ConstraintHandler, :values, :prescribed_dofs, :dh)
@define_cu(Ferrite.DofHandler, :grid)
@define_cu(Ferrite.Grid, :cells)
for T in (PointLoadCantilever, HalfMBB, LBeam, TieBeam, InpStiffness)
    @eval @define_cu($T, :ch, :black, :white, :varind)
end


text:
using JLD2, LinearAlgebra, Oceananigans, Printf
record_interaction = false
include(pwd() * "/scripts/vizinanigans.jl")
include(pwd() * "/scripts/states.jl")
include(pwd() * "/scripts/compare.jl")
include(pwd() * "/scripts/zonalstates.jl")
include(pwd() * "/analysis_scripts/" * "post_analysis.jl")

file = pwd() * "/Abernathy_16_zonal_averages.jld2"
file = pwd() * "/Abernathy_16_τ0.05_Q10_zonal_averages.jld2"
zonalstatistics = jldopen(file)
tkeys = keys(zonalstatistics["timeseries"]["t"])

i = 10

label = "vb"
vb = zonalstatistics["timeseries"][label][tkeys[i]][1,:,:]

##
filename = pwd() * "/Abernathy_16_checkpoint_iteration1158859.jld2"
filename = pwd() * "/Abernathy_16_τ0.05_Q10_zonal_averages.jld2"
states, statenames, units = grabstates(filename)
v = states[3]
b = states[5]
vb = mean(v .* b, dims = 1)[1,:,:]
vpbp = vb - mean(v, dims=1)[1,:,:] .* mean(b, dims=1)[1,:,:]
zstates = [vb, vpbp]
zstatenames = ["vb", "v'b'"]

heatmap(mean(b,dims=1)[1,:,:], colormap = :balance, colorrange = (0,2e-3))
scene = visualize(zstates, statenames = zstatenames)

max_stars_repo_name: sandreza/Mesoscale

max_issues_repo_name: sandreza/Mesoscale

max_forks_repo_name: sandreza/Mesoscale


text:
using DensityEstimationDatasets
using Test
using DataFrames: DataFrame

@testset "API Tests" begin

    file_train = twenty_datasets_file("plants", "train")
    file_valid =  twenty_datasets_file("plants", "valid")
    file_test =  twenty_datasets_file("plants", "test")

    @test file_train isa String && isfile(file_train)
    @test file_valid isa String && isfile(file_valid)
    @test file_test isa String && isfile(file_test)

    train, valid, test = twenty_datasets("plants");

    @test train isa DataFrame
    @test valid isa DataFrame
    @test test isa DataFrame

    @test size(train) == (17412, 69)
    @test size(valid) == (2321, 69)
    @test size(test) == (3482, 69)

end


text:
requiredTypes = [ "AbstractLevel" ]

for reqType in requiredTypes
    if !isdefined( Symbol( uppercase( string( reqType[ 1 ] ) ) * reqType[ 2:end ] ) )
        include(  reqType * ".jl" )
    end  # if !isdefined( Symbol( ...
end

export Affiliation
type Affiliation <: AbstractLevel

    # Name of the Affiliation ex: AC, AS
    Name::String

    # Maximum and minimum stay in this Affiliation
    MinStay::Int
    MaxStay::Int

    # A list of allowed Jobs
    NextLevel::Array{String}

    Affiliation() = new("", 0, 0, Array{String}())
    Affiliation(nm::String;
                min::Int = 0,
                max::Int = 0,
                next::Array{String} = Array{String}()) = new(nm, min, max, next)

end

max_stars_repo_name: Omazaria/CareerPathBeDef

max_issues_repo_name: Omazaria/CareerPathBeDef

max_forks_repo_name: Omazaria/CareerPathBeDef